2013高教社杯全国大学生数学建模竞赛
承 诺 书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）：
参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)：
（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容
请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。）
日期： 2013 年 09 月 13 日 赛区评阅编号（由赛区组委会评阅前进行编号）：
2013高教社杯全国大学生数学建模竞赛
编 号 专 用 页
评 阅 人 评 分 备 注 赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
古塔的变形
摘要
中国古塔---是中国五千年文明史的载体之一。古塔为祖国城市山林增光添彩，被佛教界人士尊为佛塔。矗立在大江南北的古塔，被誉为中国古代杰出的高层建筑。但是，由于长时间承受自重、气温、风力等各种作用，偶然还要受到地震等自然灾害的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。
现在已知某测绘公司先后四次对某古塔进行测量的数据，要求通过这些数据来分析古塔变形的情况以及变形的趋势。
针对问题一，根据古塔每层已知的八个测量点，本文利用最小二乘法，通过Mathematica数学软件中的拟合曲线fit[]函数将这八个点的x、y坐标拟合成一个圆，则圆心就是各层的平面中心。最后对各测量点的z坐标求平均值，即可得到各层的空间中心。
针对问题二,倾斜程度可以通过古塔塔尖坐标相对于底层坐标在竖直方向偏移的角度来判断；弯曲度通过古塔各相邻两层之间的中心轴线的夹角大小反映；扭曲度通过古塔顶层各测量点相对底层相对应的测量点在水平方向偏移的角度大小反映。
针对问题三，通过问题二得出的数据分析古塔各方面变形程度的强弱。进而判断该古塔变形的趋势。最后制定出相应的保护措施，防止古塔继续形变、遭到破坏。
关键词: 最小二乘法 Mathematica Matlab 数据拟合
1
一、问题重述
因为长期承受自重、气温、风力等各种作用，偶然还要受到地震、飓风的影响。所以古塔会产生各种形变。出于对古塔保护的目的，文物部门委托测绘公司对某塔适时进行测量。然后通过对测量数据的分析，总结出古塔变形的情况以及变形的趋势。最后制定出相应的保护措施。防止古塔继续形变，遭到破坏。
二、问题分析
2.1 问题一分析
利用Matlab软件将附件1中的坐标画在坐标轴上，然后从俯视图中我们可以清楚得看出，该塔为八边形（如图1）。又因为在测量数据中，1986年和1996年该塔第13层的第5个数据空白，而2009年和2011年该数据有测量值。则分析可知，该塔在1986年到1996年之间，第13层的第5个角遭到破坏，1996年到2009年之间对该塔进行过维修，并将被破坏的角修好。这也很好的解释了附件1中1986年和1996年第13层第5个数据空白，而2009年和2011年该数据有测量值。同时可判断，对该塔是以每层的八个角为测量点进行测量。于是我们利用Mathmatic软件拟合曲线中的fit[]函数将每层的八个点拟合为圆。然后再联系圆的圆心和八个角的纵坐标的平均值，得出每层的中心坐标。
2
图1 1986古塔坐标俯视图
2.2 问题二分析
1古塔长时间受到各种因素的影响，使其主轴线发生偏移。从上往下○
投影，塔顶与塔底的中心不在一个点上。于是可以利用坐标求出塔顶与塔底中心在水平面上的距离L，然后结合塔顶的高度H构建直角三角形。最后在三角形中解得古塔中轴线与垂线的夹角?，即该塔的倾斜度。如图2所示：
图2 古塔中轴线倾斜度示意图
2根据第一问得出的中心坐标，将相邻三个坐标相连得到一个三角形。○
如?kgf，然后利用matlab求出三角形中间角度的补角，即?。则可知道上一层塔相对于下一层塔偏移的角度。即这两层塔弯曲的程度。如图3所示：
3
图3 古塔各层弯曲示意图
3利用matlab将古塔的坐标画在坐标轴上，从俯视图中可以看见，○
同一方位各层的坐标相连是一条曲线。说明古塔已有扭曲形变。那么将同一方位古塔最底层坐标、最高层坐标、底层中轴线坐标相连得到一个三角形，如?cba。然后利用matlab求出角度?，即可知道顶层相对于底层扭曲了多少。如图4所示：
图4 古塔扭曲示意图
4
三、模型假设
1.假设古塔每层以八个角为测量点； 2.假设每次测量选择的坐标原点为定点； 3.假设古塔不受到人为破坏。
四、符号说明
H：塔高；
D：塔顶投影到塔底距离塔底中心的距离；
?：古塔倾斜角度；
xi：第i次测量古塔顶层x坐标； yi：第i次测量古塔顶层y坐标； zi：第i次测量古塔顶层z坐标； xi，：第i次测量古塔底层x坐标； yi，：第i次测量古塔底层y坐标； zi，：第i次测量古塔底层z坐标。
五、模型的建立和求解
1由于古塔是八边形，每层的八个坐标都已知。所以可以直接利用○
Mathematica的fit[]函数将古塔每层八个点的x、y坐标拟合成圆。则圆心就是各层水平面的中心位置。然后再用平均法求出各层竖直方向的中心位置。由此即可确定古塔各层的中心位置。以计算1986年第1层中心位置为例。程序如下:
t=List[{565.454,528.012},{562.058,525.544},{561.395,521.447},{563.782,518.108},{567.941,517.407},{571.255,519.857},{571.938,523.
5
953},{569.5,527.356}];
tz=Table[Flatten[{t[[i]],t[[i,1]]^2 +t[[i,2]]^2}],{i,1,8}]; Fit[tz,{1,x,y},{x,y}]
通过以上程序拟合结果为：z=-594306+1133.33x+1045.42y
其中：2a=1133.33 2b=1045.42。圆心为（a，b），即：（566.665，522.71）。 然后用平均法求出1986年第1层中心点的纵坐标。算法如下： c=（1.792+1.818+1.783+1.769+1.772+1.77+1.794+1.801+）/8 结算结果为c=1.787
所以古塔1986年第1层的中心坐标为（566.665，522.71，1.787） 用相同的方法求出各次测量的古塔各层中心坐标如下表：
表1 1986年各层中心位置 坐标 层数 第一层 第二层 第三层 第四层 第五层 第六层 第七层 第八层 第九层 第十层 第十一层 第十二层 第十三层 塔尖 X （566.665） （566.825） （566.780） （566.870） （566.870） （566.920） （566.950） （566.985） （567.020） （567.050） （567.105） （567.145） （567.200） （567.245） Y （522.710） （522.615） （522.635） （522.675） （522.575） （522.545） （522.525） （522.510） （522.495） （522.480） （522.445） （522.380） （522.380） （522.240） Z （1.787） （7.320） （12.755） （17.078） （21.721） （26.235） （29.837） （33.351） （36.855） （40.172） （44.441） （48.712） （52.834） （54.550）
表2 1996年各层中心位置 坐标 层数 第一层 第二层
X （566.665） （566.725） 6
Y （522.710） （522.670） Z （1.783） （7.320） 第三层 第四层 第五层 第六层 第七层 第八层 第九层 第十层 第十一层 第十二层 第十三层 塔尖 （566.780） （566.825） （566.875） （566.920） （566.955） （566.990） （566.525） （567.055） （567.125） （567.160） （567.210） （567.250） （522.635） （522.605） （522.570） （522.540） （522.520） （522.505） （522.990） （522.475） （522.425） （522.390） （522.375） （522.235） （12.751） （17.075） （21.716） （26.230） （29.832） （33.345） （36.848） （40.168） （44.435） （48.707） （52.834） （54.550）
表3 2009年各层中心位置 坐标 层数 第一层 第二层 第三层 第四层 第五层 第六层 第七层 第八层 第九层 第十层 第十一层 第十二层 第十三层 塔尖 坐标 层数 第一层 第二层 第三层 第四层 第五层 第六层 第七层 第八层 第九层
X （566.745） （566.780） （566.810） （566.840） （566.865） （566.860） （566.990） （567.030） （567.095） （567.150） （567.190） （567.235） （567.275） （567.336） Y （522.700） （522.670） （522.645） （522.625） （522.600） （522.910） （522.530） （522.470） （522.465） （522.410） （522.360） （522.330） （522.290） （522.215） Z （1.765） （7.291） （12.732） （17.070） （21.709） （26.211） （29.625） （33.340） （36.844） （40.161） （44.432） （48.684） （52.818） （55.091） 表4 2011年各层中心位置 X （566.745） （566.780） （566.810） （566.840） （566.865） （566.960） （566.990） （567.045） （567.095） 7
Y （522.700） （522.670） （522.645） （522.625） （522.605） （522.540） （522.530） （522.495） （522.465） Z （1.763） （7.291） （12.727） （17.052） （21.704） （26.205） （29.817） （33.337） （36.823） 第十层 第十一层 第十二层 第十三层 塔尖 （567.150） （567.190） （567.235） （567.285） （567.338） （522.420） （522.370） （522.330） （522.285） （522.214） （40.144） （44.425） （48.684） （52.813） （55.087）
2古塔的倾斜度即塔顶在底面的投影偏离底面中心的程度。可以用偏离的○
角度?反映其偏离程度大小。
设1986年、1996年、2009年、2011年塔尖中心位置坐标分别为:（x1,y1,z1）;(x2,y2,z2);(x3,y3,z3);(x4,y4,z4)。
1986年、1996年、2009年、2011年塔底中心位置坐标分别为: （x1,,y1,,z1,）;(x2,,y2,,z2,);(x3,,y3,,z3,);(x4,,y4,,z4,)。 塔高H=zi-zi,；投影L=?xi?xi,???yi?yi,?；?=arctan(L/H)。
22对各年古塔倾斜程度计算结果如下表：
表5 各年份古塔倾斜程度表 项目 年份 塔高H 52.763 53.337 53.326 53.324 投影长L 0.7465 0.7536 0.7645 0.7667 倾斜角度? 0.831 0.808 0.819 0.825 1986 1996 2009 2011
将古塔各层的中心位置坐标依次连接起来既可得到该古塔各层的中心轴线。相邻中心轴线的夹角的补角即为这两层古塔弯曲的度数。度数越大，说明弯曲越严重。具体求解过程可借助Matlab解三角形来实现。以1986年第1层和第2层为例，程序如下：
8
j=566.665;k=522.710;l=1.787;d=566.825;e=522.615;f=7.320;g=566.780;h=522.635;i=12.755;
a=sqrt((j-g)^2+(k-h)^2+(l-i)^2); c=sqrt((j-d)^2+(k-e)^2+(l-f)^2); b=sqrt((d-g)^2+(e-h)^2+(f-i)^2); A=acos((b^2+c^2+a^2)/(2*b*c)) 求得结果为A=0+1.7625i
说明1986年该古塔第1层和第2层弯曲的度数为（0+1.7625i）rad。 用相同的方法求得各年各相邻层古塔的弯曲度，结果如下：
表6 各年份古塔相邻层数弯曲程度表 相邻层数 1-2 2-3 3-4 4-5 5-6 6-7 7-8 8-9 9-10 10-11 年份 1986 1986 1986 1986 1986 1986 1986 1986 1986 1986 弯曲度 0+1.7625i 0+1.7810i 0+1.7643i 0+1.7630i 0+1.7807i 0+1.7630i 0+1.7628i 0+1.7638i 0+1.7851i 0+1.7691i 年份 1996 1996 1996 1996 1996 1996 1996 1996 1996 1996 弯曲度 0+1.7629i 0+3.0410i 0+1.7644i 0+1.7630i 0+1.7807i 0+1.7630i 0+1.7558i 0+1.7341i 0+1.7740i 0+1.7646i 9
11-12 12-13 1-2 2-3 3-4 4-5 5-6 6-7 7-8 8-9 9-10 10-11 11-12 12-13
1986 1986 2009 2009 2009 2009 2009 2009 2009 2009 2009 2009 2009 2009 0+1.7632i 0+2.0125i 0+1.88323i 0+1.7632i 0+1.7627i 0+1.7827i 0+1.7567i 0+1.9287i 0+1.7630i 0+1.7799i 0+1.7631i 0+1.7647i 0+1.7810i 0+1.7629i 1996 1996 2011 2011 2011 2011 2011 2011 2011 2011 2011 2011 2011 2011 0+1.7632i 0+1.8796i 0+1.7628i 0+1.7811i 0+1.7646i 0+1.7631i 0+1.7798i 0+1.7630i 0+1.7628i 0+1.7636i 0+1.7854i 0+1.7628i 0+1.7631i 0+1.8839i 古塔顶层某一点在底层的投影相对于底层对应的点在平面上偏移的量与古塔的扭曲有关。扭曲越大，偏移越大。扭曲越小，偏移越小。于是我们可以通过顶层对底层偏移的角度来描述古塔的扭曲程度。为了减小误差，应在顶层选多组数据计算，最后求平均值即可。我们先将顶层的点投影到底面。然后与底面相对应的点及中心位置的点相连接得到一个三角形，最后通过Matlab解出古塔扭曲的角度。如图5所示：
10
图5 1986年古塔扭曲示意图
现以求1986年古塔扭曲程度为例。程序如下所示：
j=562.508;k=528.012;d=566.933;e=522.545;g=564.716;h=523.616; a=sqrt((j-g)^2+(k-h)^2); c=sqrt((j-d)^2+(k-e)^2); b=sqrt((d-g)^2+(e-h)^2); A=acos((b^2+c^2+a^2)/(2*b*c))
解得A = 0 + 1.4760i，即1986年该塔顶层相对于底层扭曲（0 + 1.4760i）rad。
采用相同的方法求得古塔各年份顶层的扭曲度。计算结果如下表：
表7 各年份古塔顶层扭曲程度表 点 2 4 6 8 平均值 1986年 0+1.1046i 0+0.9330i 0+0.6767i 0+0.9309i 0+0.9113i 1996年 0+1.1683i 0+0.9846i 0+0.6741i 0+0.8743i 0+0.9253i 2009年 0+0.9171i 0+0.6883i 0+0.9355i 0+1.1736i 0+0.9286i 2011年 0+0.9067i 0+0.6739i 0+0.9430i 0+1.1778i 0+0.9254i 11
3通过分析第二问得出的数据发现，该古塔的倾斜程度从1986年到1996○
年在减小。但从1996年到2011年在一直增大。并且倾斜的度数明显大于古塔的弯曲度和扭曲度。而古塔的弯曲度和扭曲度始终在一个范围内波动，并没有呈现一直增大或一直减小的趋势。所以从以上分析可知，该古塔的的主要变形是倾斜。并且倾斜程度会逐年增大。因此文物保护部门应采取相关措施防止古塔的进一步倾斜。
六、模型检验
为了得到古塔每层精确的中心位置，我们采用Mathematica拟合的方法。除此之外还可以运用平均值的方法求中心位置。如1986年第1层中心位置的坐标为：
x=（565.454+562.058+561.39+563.782+567.941+571.255+571.938+569.5）/8=566.6645；
y=(528.012+525.544+521.447+518.108+517.407+519.857+523.953+527.356)/8=522.7105；
z=(1.792+1.818+1.783+1.769+1.772+1.770+1.794+1.801)/8=1.7874。 由此可见，其值与拟合法非常相近，所以模型可行。
七、模型评价
优点
1. 利用Mathematica将古塔各层八个点拟合成圆，再根据圆心是中心确定每层的中心位置，使得精确度较高。
2. 根据中心轴线来计算古塔的偏移量，使得计算结果更接近于实际。
12
3. 整个计算方法简便易行，实际操作意义强。 缺点
1. 对形变趋势分析不够全面。
2. 古塔坐标图像太少，不能全面的反应形变情况。
八、参考文献
[1]赵静，但琦，数学建模与实验[M].3版，北京：高等教育出版社，2008年。
[2]宋兆基，徐流美，MATLAB在科学计算机中的应用[M]，北京：清华出版社，2005年1月。
[3]张志勇，精通MATLAB6.5版[M]，北京：北京航空航天大学出版社，2004. [4]韩中庚，数学建模方法及其应用[M].2版，北京：解放军信息工程大学，2009年。
[5]姜启源，谢金星，叶俊，数学模型[M].3版，北京：高等教育出版社，2003.
[6]王文波，数学建模及其基础知识详解[M]，武汉：武汉大学出版社。
13
