2013数学建模C题古塔问题探究
古塔的变形监测与预测
摘要：本文针对古塔长时间承受自重和外力作用引起的组合变形问题，采用数
据处理、几何分析及曲线拟合的方法，利用matlab 、lingo 等数学软件编程、计算，给出了确定古塔各层中心位置的通用方法，根据所给的各次测量数据计算出古塔各层的中心坐标（见表9），分析了该塔的倾斜、弯曲、扭曲等变形情况并对变形趋势做了预测。得出如下结论：
1、随着时间的增长，古塔的倾斜弧度与时间呈二次抛物线
(y =0. 0004x (y =4. 6277x
2
-1. 7468x +1753. 2分布，总的变形是先减小后增大趋势。
)
2、随着时间的增长，古塔的弯曲曲率与时间呈三次函数
2
-27713x 2+5. 5E +7x -3. 7E +10分布，曲率随时间先增大，后减
)
小，再增大的趋势。
3、随着时间的增长，古塔的扭曲角弧度与时间呈二次抛物线
(y =6,. 5E -5x
2
-0. 2566x +255. 95分布，总的扭曲变形是先减小后增大趋势。
)
关键词：中心坐标 投影 组合变形 倾斜 弯曲 扭曲
一、问题重述
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，防止其变形而导致破坏，文物部门需适时对古塔进行观测，通过测量所得数据可以了解各种变形量，以制定必要的保护措施。
某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。请根据提供的4次观测数据，讨论以下问题：1. 列表给出各次测量的古塔各层中心坐标，给出确定古塔各层中心位置的通用方法。2. 分析该塔倾斜、弯曲、扭曲等变形情况。3. 分析该塔的变形趋势。
二、模型假设
1、假设塔的每层是一个平面，不考虑楼层厚度；
2、假设古塔在2011年3月前没有进行过大的维修工程，如地基处理；古塔加固；3、假设观测过程中测的数据精确度高，误差很小，可以忽略不计； 4、假设把各种外界因素对塔的影响看做一个整体来研究； 5、假设古塔在受到外界因素干扰时，未受到严重破坏；
三、符号说明
x i ：第i 层观测点的横坐标； y i ：第i 层观测点的纵坐标； z i ：第i 层观测点的竖坐标；
G (x 0, y 0, z 0) ：表示的各层中心坐标；
θ：第1层塔中心塔尖中心连线的倾斜角；
∆h ：塔尖到第1层竖坐标差；
k ：曲线弯曲曲率；
ϕ：扭转角；
四、模型准备
观察所给的数据，发现 1986年和1996年的测量数据中缺失了第13层观测点5的数据，对此我们观察当年已有的13层观测点5数据，计算各观测点的平均坐标，发现观测值与平均值之差在允许范围内，因此我们选取平均值作为缺失数据的补充，1986年观测点5补充数据（566.195,520.132,52.795），1996年观测点5补充数据（566.2314,520.1482.52.731）。
2009年和2011年观测点只有1个，我们只能将此观测数据作为塔尖的中心坐标。
五、模型建立与求解
（一）问题一
古塔有13个塔层和1个塔尖，通过对数据的分析以及matlab 作图（如图1）
可以知道古塔1到13层的结构相似，而塔尖和塔层观测点数目不同，而且结构各异，由此我们分别讨论塔层和塔尖，找到它们的中心。
图1
1 第1层到13层 1.1 问题分析
每层塔的8个观测点的z 不全相同，它们构成一个不规则的曲面。首先，将8个观测点
(x i , y i , z i )
投影到xoy 面，投影点分别为
(x i , y i )
。利用matlab 拟合一
'
'
'
(x , y , z )
个斜平面Z =Ax +By +C ，满足：8个观测点与斜平面的交点i i i 与对应
观测点的距离之和最小。依次连接这8个交点，所得图形的中心就是要找的塔层中心。再将xoy 面的投影点依次连接，利用lingo 软件拟合出一个到8个顶点的
22
(x -x 0)+(y -y 0)=r 2(x , y )
距离之和最小的圆，圆心G 00可以近似看做该八边形
的中心。最后，过圆心作平行于z 轴的直线，它与斜平面的交点
(x 0, y 0, C +Ax 0+By 0)就是我们所要找的塔层的中心点。
1.2模型建立
设拟合的斜平面方程为：
z =Ax +By +C （1）
建立模型： min d =∑(z i -z )
i =1
'
i
8
2
z i ' =Ax i +By i +C （i =1, ,8）
设xoy 面内拟合的圆的方程为：
(x -x 0)2+(y -y 0)2=r 2 （2） 建立模型：
min d =∑r ) 2
i =18
1.3模型求解
根据已知数据，利用matlab 拟合可以算出A ，B ，C ，得到斜平面的方程，做出图像（如图2）；利用lingo 软件拟合，可以得到圆的方程，如图3。将圆心坐标(x 0, y 0)，代入斜平面的方程可解出
z 0
，即得到中心坐标的通式
(x 0, y 0, C +
Ax 0+By 0)
图2
图3
将四年测量的数据代入之后得到1到13层的中心坐标（如表1）
1.4模型分析与评价
该模型的建立，采用投影、拟合的方法将不规则转化为规则，将空间转化到平面，近似找到塔层的中心坐标，由于第1到13层的结构相似，因此方法具有通用性。即使1986年和1996年的测量数据中缺失了第13层观测点5的数据，但是我们根据其他各层的数据猜测得到，不影响方法的使用。但是也因为猜测和 两次拟合，塔层中心坐标的计算产生误差。
2 塔尖 2.1 问题分析
考虑到塔尖的横截面面积相对塔层横截面较小，且已知的测点数据表明，塔尖各测点变形相对较小，因此我们采用求平均值的方法来计算塔尖的中心坐标
G (x 0y 0, z 0)
。
2.2 模型建立
将塔尖的4个观测点坐标设为中心坐标的计算式为：
(x i , y i , z i )，
4
⎧
x i ∑⎪
⎪x 0=i =1
4⎪
4⎪
y i ⎪∑⎪
⎨y 0=i =1 (i =1. 2. 3. 4)
4⎪
4
⎪
z i ∑⎪
⎪z 0=i =1
4⎪
⎪⎩
2.3 模型求解
将四次观测数据分别代入平均值计算的中心坐标公式：
44
⎛4⎫ ∑x i ∑y i ∑z i ⎪i =1
, i =1, i =1⎪ G 444⎪ ⎪⎝⎭
算出塔尖的中心坐标（如表2）
2.4 模型分析与评价
该模型可以简单、方便地计算出塔尖的中心坐标，但是2009年和2011年观测点只有1个，我们只能将此观测数据作为塔尖的中心坐标。而且平均值计算的方法有它的缺点，观测点的数目越少，结果的误差就会越大。这是该模型有待改进的地方。
（二）问题二 1 问题分析
在古塔刚修建成功时，不考虑修建过程中的误差，各层中心的连线（中性轴）应该是垂直于底平面的一条垂线。但是由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生组合变形，包括倾斜、弯曲、扭曲等，此时其中性轴也会发生相应的变形（见图4）。因此，我们可以通过对古塔中性轴的变形来讨论古塔的变形，分别单独讨论倾斜、弯曲和扭曲。
图4（中心点连线的图）
1.1 倾斜变形 1.1.1问题分析
塔的倾斜变形是一种塑性变形，是塔基固定，塔身和塔尖产生倾斜的一种变形。这种变形反映了整座塔的变形，因此我们考虑塔尖中心相对于第1层中心的倾斜。连接第1层的中心G 1和塔尖直线l 上，
G 0
G 0
的中心，将其投影到过G 1且平行于z 轴的
的投影点到G 1的距离记作∆h ，G 0到l 的距离记作∆d ，倾斜角
θ=arctan(
∆d
（如图5） ) 。∆h
图5
1.1.2模型建立
∆h =z 塔尖
-z 1
∆d =θ=arctan(
∆d ) ∆h
1.1.3 模型求解
将第一题中计算出来的数据带到以上公式，分别计算出四年古塔的倾斜角度（表3）
1.1.4模型分析与评价
根据上述表格所得的数据可知，古塔的倾斜角度是很微小的。在古塔开始发生倾斜时，由于地基基础下的地质分布不同及其他外界不利因素，导致古塔刚开始是向一边倾斜。但是随着时间的变化，以及各种作用下，古塔开始向其另一边发生倾斜，从而抑制塔向一边倾斜的趋势。由二力平衡可知，两边对塔产生力的作用，这样使塔处于相对平衡的状态。由此我们可以从上述计算出的结果看出，在没有大的事故作用下，塔在一定的范围时间内不会被外界作用所毁坏。
2.1 弯曲变形
2.1.1问题分析
曲线的弯曲程度可以用曲率来度量，依次连接塔层中心点，得到的是空间折
⎧x =x (t ) 线段，将它拟合成一条曲线，参数方程为：⎪。代入曲率公式
⎨y =y (t )
⎪z =z (t ) ⎩
k =
可以计算出各点的曲率，由于曲线各点的弯曲程度不同，我们用曲线的最大曲率来度量古塔的弯曲变形程度。
2.1.2模型建立
⎧x =k 1t 2+k 2t +k 3⎪
曲线方程：⎨y =m 1t 2+m 2t +m 3
⎪z =t ⎩
max
k 2.1.3模型求解
利用第一题计算的中心坐标，运用matlab 软件拟合出曲线（见图6）的方程，
计算出参数方程的一阶导数、二阶导数代入曲率公式，得到
k =
-m 2
再利用matlab 对k 求导数，计算出当t =时，k 取得最大值。
2m 1
图6
将第一题求得的中心坐标带入，按照上述方法计算，可以分别计算不同年古塔的最大曲率（见表4）。
3.1扭曲变形 3.1.1问题分析
扭曲变形是指各塔层横截面绕轴线的相对转动，转动后与其原来的位置形成一定的角度，称这个角为相对扭转角。考虑古塔的扭曲是绕着中性轴的转动，而每层扭转的程度不一样，因此分层讨论。假设第1层为基准面，不发生扭曲，考虑其他各层相对于它绕着中性轴的扭曲程度。而塔尖的数据不全，且近似看作在中性轴上，因此只考虑第2层到第13层的扭曲。
每一个塔层有8个观测点，任选一个作为扭曲变形的研究对象，不妨假设观测点1。将第1层的观测点1投影到第i 层中心点G (x Gi , y Gi ) 所在的水平面，投影点记作b (x i -1, y i -1) ，分别连接G (x Gi , y Gi ) 与b (x i -1, y i -1) 以及第i 层观测点1（b /(x i , y i ) ）与G (x Gi , y Gi ) ，它们之间的夹角ϕ就是要找的相对扭转角（如图7）。
图7
3.1.2模型建立
已知坐标G (x Gi , y Gi ) ，b (x i -1, y i -1) , b /(x i , y i ) ，利用余弦定理列出计算式(i =2、3.......13) 。
Gb =(x i -1-x Gi , y i -1-y Gi ) ，Gb /=(x i -x Gi , y i -y Gi ) ，
bb ' =x i -x i -12+y i -y i -12 φ=arc cos(
3.1.3模型求解
将题中所给的数据以及第一题计算的中心坐标代入上式，利用matlab 可以计算出每一层相对于第一层的扭转角（见表5）。
Gb +Gb /-bb /
2Gb Gb /
2
22
)
（三）问题三 3.1问题分析
变形趋势就是一种惯性，没有外力不会轻易改变状态。由于古塔的变形
是多种影响的组合变形，如倾斜、弯曲、扭曲等等。这些影响来自各个方面，不便于综合确定各种影响的程度，这样就造成了无法整体去研究它的变形趋势。因此我们还是分开考虑，分别考虑倾斜、弯曲、扭曲对塔的变形趋势。这些变形趋势都是建立在问题二的基础之上，模型相同，因此在以下变形趋势中同样适用。
3.1.1倾斜变形趋势分析
在问题二倾斜变形模型求解中，我们计算出了塔的整体倾斜角（见表3），我们再用时间（单位是年）与倾斜角（单位是弧度）建立平面直角坐标系，以O 为坐标原点，时间（单位是年）为X 轴, 倾斜角（单位是弧度）为Y 轴，坐标系表示倾斜角与时间变化的关系。所以我们可以得出四个坐标点
M (x i , y i )(i =1, 2, 3, 4)，再通过Excel 处理拟合曲线作图（见表6、图8）。
图8
根据上述图表预测，随着时间的增长，古塔的倾斜弧度与时间呈二次抛物线
2
(y =0. 0004x
-1. 7468x +1753. 2分布，当x =2184（年）时塔的倾斜角最小，即
)
倾斜程度最小，塔的变形最小。总的变形是先减小后增大趋势。
3.1.2 弯曲变形趋势分析
在问题二弯曲变形模型中，我们求解出了曲线的最大曲率，曲线的最大曲率能够帮助我们分析曲线的弯曲程度。基于此情况，我们在问题二的弯曲变形模型基础上去分析弯曲变形趋势。首先，建立平面直角坐标系，O 为坐标原点，OX 轴代表时间（年），OY 轴代表曲率(1⨯10-10)，此坐标系是表示曲率与时间的变
化关系。
所以我们可以得出四个坐标点N (x i , y i )(i =1, 2, 3, 4)，再通过Excel 处理拟合曲线
作图（见表7、图9）。
图9
根据上述图表预测，随着时间的增长，古塔的弯曲曲率与时间呈三次函数
(y =4. 6277x
2
-27713x 2+5. 5E +7x -3. 7E +10分布，曲率随时间先增大，后减
)
小，再增大的趋势。
3.1.3扭曲变形趋势分析
在问题二扭曲变形模型中，已经求解出了各层相对于第1层的扭曲角ϕ（见表
5），观察表的数据，可以找出一规律：随着塔层的增高，扭曲角增大，即扭曲程度增大。最大扭曲角来最能够代表整座塔的变形，因此我们还是选取它来研究扭曲变形趋势。同样建立平面直角坐标系，O 为坐标原点，OX 轴代表时间（年），OY 轴代表最大扭曲角（弧度）。此坐标系是表示扭曲角随时间的变化关系。
所以我们可以得出四个坐标点W (x i , y i )(i =1, 2, 3, 4)，再通过Excel 处理拟合曲线作图（见表8、）。
表8
图10
根据上述图表预测，随着时间的增长，古塔的扭曲角弧度与时间呈二次抛物线(y =6,. 5E -5x
2
-0. 2566x +255. 95分布，当x =1970（年）时塔的扭曲角最小，
)
即是扭曲程度最小。总的扭曲变形是先减小后增大趋势。
六、模型推广与评价
一、优点：
本文通过对数据的处理、分析，利用MATLAB 和LINGO 软件对古塔每年每层的测量点数据进行拟合，有效地将不规则立体转化为规则的平面图形，利用规则平面图形的中心来找到塔层的中心。这样的处理很方便、明了，便于实现。此模型可用于大多数文物建筑，通用性较强。
二、缺点：
在数据的处理中，由于对数据采用了“四舍五入”方法，并且预测补充缺失数据，多次拟合，从而使数据误差增大，降低了其精确性。
七、参考文献
[1] 赵静 但琦，《数学建模与数学实验》第三版，北京：高等教育出版社，2008年1月。 [2] 王兵团，《数学建模基础》，北京:清华大学出版社，2004年11月。
[3] 姜启源, 谢金星, 叶俊, 《数学模型》第三版，北京：高等教育出版社,2003年。 [4] 于润伟 朱晓慧，《MTLAB 基础及应用》第二版，北京：机械工业出版社，2011年1月。
[5] 同济大学航空航天与力学学院基础力学教学研究部，《材料力学》，同济大学出版社，2008年8月。 [6] 同济大学数学系，《高等数学》第六版上•下册，高等教育出版社，2007年4月。
[7] 廖武鹏 邓俊晔 王丹，《管道切片的三维重建》，《工程数学学报》，第19卷：22-34，2002年2月。 [8] 梁海奎，《古塔变形测量方法探讨》，《城市勘测》，第3期：114-118页,2011年6月。
[9] 胡志晓，《古塔倾斜观测和数据分析》，《江苏建筑》，第6期：34-44页，2011年6月。
[10] 卢俊龙 张萌 姚谦峰，《砖石古塔纠偏工程的有限元分析》
八、附录
1、用matlab 实现空间斜平面拟合：
x=[565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5]'; y=[528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356]'; z=[1.792 1.818 1.783 1.769 1.772 1.77 1.794 1.801]'; scatter3(x,y,z,'filled') hold on
X = [ones(8,1) x y]; b = regress(z,X) ；
xfit = min(x):0.1:max(x); yfit = min(y):0.1:max(y);
[XFIT,YFIT]= meshgrid (xfit,yfit);
ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT; mesh (XFIT,YFIT,ZFIT) x=566.6650; y=522.7090;
Z=b(1)+b(2)*x+b(3)*y
2、运用lingo 软件对数据处理的程序：
model: sets:
dianji/1..8/:x,y,d; endsets data:
x=565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5; y=528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356; enddata
min =@sum(dianji(i):(@sqrt((x(i)-x0)^2+(y(i)-y0)^2)-r)^2); end
通过此程序可以到：
1986、1996、2009、2011年每层的中心坐标，下面是1986年所得的数据：
1、 Variable Value Reduced Cost X0 566.6650 0.000000 Y0 522.7090 0.000000 R 5.427619 0.000000
2、 Variable Value Reduced Cost X0 566.7225 0.000000 Y0 522.6714 0.000000 R 5.226233 0.000000
3、 Variable Value Reduced Cost X0 566.7787 0.000000 Y0 522.6348 0.000000 R 5.028751 0.000000
4、 Variable Value Reduced Cost X0 566.8227 0.000000 Y0 522.6058 0.000000 R 4.872247 0.000000
5、 Variable Value Reduced Cost X0 566.8700 0.000000 Y0 522.5746 0.000000 R 4.704324 0.000000
6、 Variable Value Reduced Cost X0 566.9188 0.000000 Y0 522.5443 0.000000 R 4.541165 0.000000 7、 Variable Value Reduced Cost X0 566.9522 0.000000 Y0 522.5271 0.000000 R 4.272584 0.000000 8、 Variable Value Reduced Cost X0 566.9847 0.000000 Y0 522.5103 0.000000 R 4.011286 0.000000 9、 Variable Value Reduced Cost X0 567.0176 0.000000 Y0 522.4934 0.000000 R 3.750459 0.000000 10、 X0 567.0476 0.000000
Variable Value Reduced Cost
Y0 522.4792 0.000000 R 3.503796 0.000000 11、 Variable Value Reduced Cost X0 567.1014 0.000000 Y0 522.4388 0.000000 R 3.276693 0.000000 12、 Variable Value Reduced Cost X0 567.1554 0.000000 Y0 522.3984 0.000000 R 3.049347 0.000000 13、 Variable Value Reduced Cost X0 567.2020 0.000000 Y0 522.3800 0.000000 R 2.819595 0.000000
按照同样的方法我们可以的到1996，2009,2011年古塔的每层中心坐标（见表9）：
表9
⎧x =x (t )
3、根据空间参数方程：⎪ 通过matlab 软件对上述数据进行拟合： ⎨y =y (t )
⎪z =z (t ) ⎩
x=data(:,1); >> z=data(:,3); >> polyfit(z,x,2);
>> y=data(:,2); >> z=data(:,3); >> polyfit(z,y,2);
得到具体的空间曲线方程（以1986年的为例）：
⎧x(t)=1.[**************]e -05t 2+0.0095t +566.6518
⎪2
⎨y(t)=-3.[**************]e -06t -0.0061t + 522.7157 ⎪z(t)=t ⎩
通过下面这个程序，画出1986年的空间曲线：
>> t=0:100;
>> x=1.[**************]e-05*t.^2+0.0095*t+566.6518; >> y=-3.[**************]e-06*t.^2 -0.0061*t+ 522.7157; >> z=t;
>> plot3(x,y,z,'r-');
然后用此方法分别对每年的数据进行拟合，这样我们就可以分别得到每年空间曲线方程
和空间曲线。
用diff 命令对空间曲线方程求导：
syms t
diff(1.[**************]e-05*t^2+0.0095*t+566.6518)
diff(1.[**************]e-05*t^2+0.0095*t+566.6518,2)
代入曲率公式，这样便可求出此曲线的曲率。
4、运用matlab 软件求最优半径程序：
function [R,A,B]=circ(x,y,N)
x1 = 0;
x2 = 0; x3 = 0; y1 = 0; y2 = 0; y3 = 0; x1y1 = 0; x1y2 = 0; x2y1 = 0;
for i = 1 : N
x1 = x1 + x(i);
x2 = x2 + x(i)*x(i);
x3 = x3 + x(i)*x(i)*x(i); y1 = y1 + y(i);
y2 = y2 + y(i)*y(i);
y3 = y3 + y(i)*y(i)*y(i);
x1y1 = x1y1 + x(i)*y(i);
x1y2 = x1y2 + x(i)*y(i)*y(i); x2y1 = x2y1 + x(i)*x(i)*y(i); end
C = N * x2 - x1 * x1;
D = N * x1y1 - x1 * y1;
E = N * x3 + N * x1y2 - (x2 + y2) * x1; G = N * y2 - y1 * y1; H = N * x2y1 + N * y3 - (x2 + y2) * y1;
a = (H * D - E * G)/(C * G - D * D); b = (H * C - E * D)/(D * D - G * C); c = -(a * x1 + b * y1 + x2 + y2)/N;
A = a/(-2); %x 坐标
B = b/(-2); %y 坐标
R = sqrt(a * a + b * b - 4 * c)/2;
5、用matlab 求解曲线方程曲率：
clear
syms k1 k2 m1 m2 t
>>
diff(((2*m1)^2+(2*k1)^2+((2*m1)*(2*k1*t+k2)-(2*k1)*(2*m1*t+m2))^2)/((2*k1+k2)^2+(2*m1*t+m2)^2+1)^3)
ans =
-(12*m1*(m2 + 2*m1*t)*((2*m1*(k2 + 2*k1*t) - 2*k1*(m2 + 2*m1*t))^2 + 4*k1^2 + 4*m1^2))/((m2 + 2*m1*t)^2 + (2*k1 + k2)^2 + 1)^4
>> s=-(12*m1*(m2 + 2*m1*t)*((2*m1*(k2 + 2*k1*t) - 2*k1*(m2 + 2*m1*t))^2 + 4*k1^2 + 4*m1^2))/((m2 + 2*m1*t)^2 + (2*k1 + k2)^2 + 1)^4;
>> solve(s)
ans =
-m2/(2*m1)
>> t=-m2/(2*m1);
>>k=((2*m1)^2+(2*k1)^2+((2*m1)*(2*k1*t+k2)-(2*k1)*(2*m1*t+m2))^2)/((2*k1+k2)^2+(2*m1*t+m2)^2+1)^3
k =
(4*m1^2*(k2 - (k1*m2)/m1)^2 + 4*k1^2 + 4*m1^2)/((2*k1 + k2)^2 + 1)^3
k1=1.56*10^(-5);k2=0.0095;m1=-3.78*10^(-6);m2=(-0.0061);
>> k=(4*m1^2*(k2 - (k1*m2)/m1)^2 + 4*k1^2 + 4*m1^2)/((2*k1 + k2)^2 + 1)^3
k =
1.0303e-09