2013数学建模-古塔变形

承诺书我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。

如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）：我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：参赛队员(打印并签名) ：1.2.3.指导教师或指导教师组负责人(打印并签名)：（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。

以上内容请仔细核对，提交后将不再允许做任何修改。

如填写错误，论文可能被取消评奖资格。

）日期： 2013 年 09 月 13 日赛区评阅编号（由赛区组委会评阅前进行编号）：编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：古塔的变形摘要中国古塔---是中国五千年文明史的载体之一。

古塔为祖国城市山林增光添彩，被佛教界人士尊为佛塔。

矗立在大江南北的古塔，被誉为中国古代杰出的高层建筑。

但是，由于长时间承受自重、气温、风力等各种作用，偶然还要受到地震等自然灾害的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。

现在已知某测绘公司先后四次对某古塔进行测量的数据，要求通过这些数据来分析古塔变形的情况以及变形的趋势。

针对问题一，根据古塔每层已知的八个测量点，本文利用最小二乘法，通过Mathematica数学软件中的拟合曲线fit[]函数将这八个点的x、y坐标拟合成一个圆，则圆心就是各层的平面中心。

最后对各测量点的z坐标求平均值，即可得到各层的空间中心。

针对问题二,倾斜程度可以通过古塔塔尖坐标相对于底层坐标在竖直方向偏移的角度来判断；弯曲度通过古塔各相邻两层之间的中心轴线的夹角大小反映；扭曲度通过古塔顶层各测量点相对底层相对应的测量点在水平方向偏移的角度大小反映。

针对问题三，通过问题二得出的数据分析古塔各方面变形程度的强弱。

进而判断该古塔变形的趋势。

最后制定出相应的保护措施，防止古塔继续形变、遭到破坏。

关键词: 最小二乘法 Mathematica Matlab 数据拟合一、问题重述因为长期承受自重、气温、风力等各种作用，偶然还要受到地震、飓风的影响。

所以古塔会产生各种形变。

出于对古塔保护的目的，文物部门委托测绘公司对某塔适时进行测量。

然后通过对测量数据的分析，总结出古塔变形的情况以及变形的趋势。

最后制定出相应的保护措施。

防止古塔继续形变，遭到破坏。

二、问题分析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[]函数将每层的八个点拟合为圆。

然后再联系圆的圆心和八个角的纵坐标的平均值，得出每层的中心坐标。

图1 1986古塔坐标俯视图2.2 问题二分析○1古塔长时间受到各种因素的影响，使其主轴线发生偏移。

从上往下投影，塔顶与塔底的中心不在一个点上。

于是可以利用坐标求出塔顶与塔底中心在水平面上的距离L，然后结合塔顶的高度H构建直角三角形。

最后在三角形中解得古塔中轴线与垂线的夹角θ，即该塔的倾斜度。

如图2所示：图2 古塔中轴线倾斜度示意图○2根据第一问得出的中心坐标，将相邻三个坐标相连得到一个三角形。

如∆kgf，然后利用matlab求出三角形中间角度的补角，即θ。

则可知道上一层塔相对于下一层塔偏移的角度。

即这两层塔弯曲的程度。

如图3所示：图3 古塔各层弯曲示意图○3利用matlab将古塔的坐标画在坐标轴上，从俯视图中可以看见，同一方位各层的坐标相连是一条曲线。

说明古塔已有扭曲形变。

那么将同一方位古塔最底层坐标、最高层坐标、底层中轴线坐标相连得到一个三角形，如∆cba。

然后利用matlab求出角度θ，即可知道顶层相对于底层扭曲了多少。

如图4所示：图4 古塔扭曲示意图三、模型假设1.假设古塔每层以八个角为测量点；2.假设每次测量选择的坐标原点为定点；3.假设古塔不受到人为破坏。

四、符号说明H：塔高；D：塔顶投影到塔底距离塔底中心的距离；：古塔倾斜角度；x i：第i次测量古塔顶层x坐标；y i：第i次测量古塔顶层y坐标；z i：第i次测量古塔顶层z坐标；x i，：第i次测量古塔底层x坐标；y i，：第i次测量古塔底层y坐标；z i，：第i次测量古塔底层z坐标。

五、模型的建立和求解○1由于古塔是八边形，每层的八个坐标都已知。

所以可以直接利用Mathematica的fit[]函数将古塔每层八个点的x、y坐标拟合成圆。

则圆心就是各层水平面的中心位置。

然后再用平均法求出各层竖直方向的中心位置。

由此即可确定古塔各层的中心位置。

以计算1986年第1层中心位置为例。

程序如下:t=List[{565.454,528.012},{562.058,525.544},{561.395,521.447},{56 3.782,518.108},{567.941,517.407},{571.255,519.857},{571.938,523.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）用相同的方法求出各次测量的古塔各层中心坐标如下表：○2古塔的倾斜度即塔顶在底面的投影偏离底面中心的程度。

可以用偏离的角度θ反映其偏离程度大小。

设1986年、1996年、2009年、2011年塔尖中心位置坐标分别为:（x 1,y 1,z 1）;(x 2,y 2,z 2);(x 3,y 3,z 3);(x 4,y 4,z 4)。

1986年、1996年、2009年、2011年塔底中心位置坐标分别为: （x 1,,y 1,,z 1,）;(x 2,,y 2,,z 2,);(x 3,,y 3,,z 3,);(x 4,,y 4,,z 4,)。

塔高H=z i -z i ,；投影L=()()2,2,i i i i y y x x -+-；θ=arctan(L/H)。

对各年古塔倾斜程度计算结果如下表：将古塔各层的中心位置坐标依次连接起来既可得到该古塔各层的中心轴线。

相邻中心轴线的夹角的补角即为这两层古塔弯曲的度数。

度数越大，说明弯曲越严重。

具体求解过程可借助Matlab 解三角形来实现。

以1986年第1层和第2层为例，程序如下：j=566.665;k=522.710;l=1.787;d=566.825;e=522.615;f=7.320;g=566.78 0;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。

用相同的方法求得各年各相邻层古塔的弯曲度，结果如下：古塔顶层某一点在底层的投影相对于底层对应的点在平面上偏移的量与古塔的扭曲有关。

扭曲越大，偏移越大。

扭曲越小，偏移越小。

于是我们可以通过顶层对底层偏移的角度来描述古塔的扭曲程度。

为了减小误差，应在顶层选多组数据计算，最后求平均值即可。

我们先将顶层的点投影到底面。

然后与底面相对应的点及中心位置的点相连接得到一个三角形，最后通过Matlab解出古塔扭曲的角度。

如图5所示：图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。

采用相同的方法求得古塔各年份顶层的扭曲度。

计算结果如下表：点1986年1996年2009年2011年2 0+1.1046i 0+1.1683i 0+0.9171i 0+0.9067i4 0+0.9330i 0+0.9846i 0+0.6883i 0+0.6739i6 0+0.6767i 0+0.6741i 0+0.9355i 0+0.9430i8 0+0.9309i 0+0.8743i 0+1.1736i 0+1.1778i平均值0+0.9113i 0+0.9253i 0+0.9286i 0+0.9254i○3通过分析第二问得出的数据发现，该古塔的倾斜程度从1986年到1996年在减小。