古塔的形变(共21页)

1、精选优质文档-倾情为你奉上承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）

2、： 43 所属学校（请填写完整的全名）：山西财经大学华商学院 参赛队员 (打印并签名) ：1. 胡珺怡 2. 赵萌 3. 陈茜 指导教师或指导教师组负责人 (打印并签名)： 日期： 2013 年 9 月 16 日赛区评阅编号（由赛区组委会评阅前进行编号）：2013高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：古塔的形变摘 要 本文针对古塔的变形这一问题，用的matlab软件来建立圆形塔模型、斜率倾斜模型和曲率弯曲

3、模型求古塔各层中心位置的问题可转化为求圆心问题。当3个观测点位于一个圆周上时，可以用图解或解析的方法确定圆心的位置；当观测点数大于3时，可以用最小二乘进行曲线拟合，计算出圆心的坐标。下面讨论用最小二乘法进行曲线拟合的问题。，综合运用最小二乘法，机器学习算法，线性回归法，曲线拟合法，梯度下降法来模型求解。使古塔的形变问题得到很好的解决。 通过对问题一所给数据的处理分析，用matlab对其中已知数据的作图发现所给的数据呈圆形，从而建立圆形塔模型，画出相应图形。发现题设所给的数据所画的图形拟合程度极高。针对问题二，我们分别建立斜率倾斜模型和曲率弯曲模型，并用matlab进行了有效的验证。 针对问题三

4、，我们发现圆形塔倾斜程度大扭曲程度高，需要有关部门进行维修保养关键词 最小二乘法 线性回归法 圆形塔模型 斜率倾斜模型 曲率弯曲模型一、问题的提出需要进行安全性检测的古建筑一般是文物，因为文物的保护，旅游开发等需要将其建设成文物保护区或风景区。这类区域的特点是植被密度大，各种附属构筑物的遮挡导致测量时通视条件不好。由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。在这里以一组古塔的变形四次监测数据为例进行探讨。现在就此现象，建立数学模型，回答以下问题。

5、1. 找到圆塔的中心便于受力分析，对古塔的维修和改造有很大的帮助2. 针对古塔的倾斜、弯曲、扭曲，可以探测到古塔的耐用程度3. 根据古塔的变形趋势来判定古塔的变形二、问题的分析关于古塔的现状，我们在数据处理研究时应该怎么做。古塔各层中心点和古中心坐标应该这样计量确定，塔的倾斜角，曲率，扭距应该如何度量。对于问题一，我们利用圆形塔模型，顺利找到塔的中心点。根据观测的数据可以看到每一层有八个点，可以把它看成一个八边形，八个点在八边形的边上。我们利用八边形的外接圆寻找其中心点。因为八边形外接圆的圆心和八边形的中心重合，所以我们提出的圆形塔模型是合理的对于问题二，由于长时间承受风力、温度、自然灾害的影

6、响，面临的问题是塔结构的变形、扭曲、倾斜等。这些问题可由问题一中的圆形塔模型以及拟合出来的曲线图观测出来。再利用斜率倾斜模型和曲率弯曲模型可分析得出结果。对于问题三，通过分析塔中心的变化以及塔的倾斜、弯曲、扭曲的情况图表可得出古塔变形趋势，从而得出古塔有受损情况。三、基本假设观测到的数据都具有真实性。不考虑风力、湿度、温度对古塔变形的影响按Z轴的平均值计算题目所给数据真实可靠四、符号说明及名词定义A古塔各层中心x坐标值B古塔各层中心y坐标值Z古塔各层中心z坐标值R古塔各层半径a,b,c问题1中圆拟合参数di样本集中的点到塔层中心的距离Xi样本集中第i条记录的x坐标Yi样本集中第i条记录的y坐标

7、h(x)问题2中目标拟合直线/曲线方程J()损耗函数0,1,2问题2中的拟合参数古塔倾斜角k拟合曲线的曲率 五、模型的分析、建立及求解5.1 问题1模型的提出及分析使用Matlab对数据进行可视化分析，根据观测数据的二维平面俯瞰分布如图1，提出如下假设：古塔的水平横截面为圆形：图1 观测数据的二维平面俯瞰分布图5.1.1 问题1模型建立与求解求古塔各层中心位置的问题可转化为求圆心问题。当3个观测点位于一个圆周上时，可以用图解或解析的方法确定圆心的位置；当观测点数大于3时，可以用最小二乘进行曲线拟合，计算出圆心的坐标。下面讨论用最小二乘法进行曲线拟合的问题。（1）模型建立： 从解析几何中可知，中

8、心在A,B、半径为R的圆方程是：R2=(x-A)2+(y-B)2 (1)将方程（1）中的括号展开并移项，得x2+y2-2Ax-2By+A2+B2-r2 (2)如果令a=-2A;b=-2B;c=A2+B2-R2则方程（2）变成x2+y2+ax+by+c=0 (3)只要求出参数a, b, c就可以求得圆心半径的参数A=a-2 (4)B=b-2 (5)R=12a2+b2-4c (6)圆心的Z坐标值C可通过求观测坐标Z分量的均值得到：C=ZiN样本集Xi, Yi i(1,2,3N)中点到圆心的距离为di:di2=(Xi-A)2+(Yi-B)2 (7)点Xi, Yi到圆边缘的距离的平方和半径的平方的差为

9、：i=di2-R2=Xi2+Yi2+aXi+bYi+c (8)令Qa,b,c为i的平方和：Qa,b,c=i2=Xi2+Yi2+aXi+bYi+c2 (9)能够使得Q(a,b,c)的值最小的a,b,c即为所要求的解。图2为模型示意图： (Xi,Yi)(A,B) diR图2 模型示意图（2）模型求解：平方差Q(a,b,c)大于0，因此函数存在大于或等于0的极小值。Q(a,b,c)对a,b,c求偏导，令偏导等于0，可求得极值点，比较所有极值点的函数值即可得到最小值。偏导方程组：Q(a,b,c)a=2(Xi2+Yi2+aXi+bYi+c)Xi=0 (10)Q(a,b,c)b=2(Xi2+Yi2+aXi

10、+bYi+c)Yi=0 (11)Q(a,b,c)c=2(Xi2+Yi2+aXi+bYi+c)=0 (12)求解方程组：10*N-12*Xi得：NXi2-XiXia+NXiYi-XiYib+NXi3+NXiYi2-Xi2+Yi2Xi=0 (13)11*N-12*Yi得：NXiYi-XiYia+NYi2-YiYib+NYi3+NXi2Yi-Xi2+Yi2Yi=0 (14)令C=NXi2-XiXiD=NXiYi-XiYiE=NXi3+NXiYi2-Xi2+Yi2XiG=NYi2-YiYiH=NYi3+NXi2Yi-Xi2+Yi2Yi可得方程组：Ca+Db+E=0 (15)Da+Gb+H=0 (16)

11、可解得：a=(HD-EG)(CG-D2) (17)b=(HC-ED)(D2-CG) (18)c=-Xi2+Yi2+aXi+bYiN (19)进而求得A,B,R的估计拟合值：A=a-2 B=b-2 R=12a2+b2-4c C=ZiN编写matlab函数function R,A,B,C=fitted_model(x,y,z,N)计算得到各次测量的古塔各层中心坐标如表1：观测年份塔层中心坐标半径1986年第1层(566.6650,522.7090，1.7874)5.4276第2层 (566.7225,522.6713, 7.3202)05.2263第3层 (566.7786,522.6346, 1

12、2.7552)5.0288第4层(566.8226,522.6056, 17.0783)4.8724第5层(566.8698,522.5746, 21.7205)4.7046第6层(566.9185,522.5443, 26.2351)4.5416第7层(566.9521,522.5271, 29.8369)4.2731第8层(566.9846,522.5102, 33.3509)4.0118第9层(567.0176,522.4933, 36.8549)3.7511第10层(567.0476,522.4789, 40.1721)3.5046第11层(567.1013,522.4385, 44.

13、4409)3.2772第12层(567.1552,522.3981, 48.7119)3.0496第13层(567.2025,522.3841, 52.8400)2.8176塔顶(567.2434,522.2393, 55.1232)0.01071996年第1层(566.6652,522.7087, 1.7830)5.4276第2层(566.7234,522.6703, 7.3146)5.2263第3层(566.7803,522.6330, 12.7508)5.0288第4层(566.8248,522.6034, 17.0751)4.8724第5层(566.8726,522.5717, 21.

14、7160)4.7046第6层(566.9219,522.5409, 26.2295)4.5416第7层(566.9559,522.5232, 29.8322)4.2730第8层(566.9888,522.5060, 33.3454)4.0118第9层(567.0223,522.4883, 36.8482)3.7511第10层(567.0527,522.4737, 40.1676)3.5045第11层(567.1070,522.4327, 44.4354)3.2772第12层(567.1612,522.3919, 48.7074)3.0497第13层(567.2086,522.3778, 52.

15、8360)2.8176塔顶(567.2494,522.2286, 55.1198)0.01302009年第1层(566.7447,522.7004, 1.7645)5.4269第2层(566.7790,522.6720, 7.3090)5.2162第3层(566.8118,522.6455, 12.7323)5.0104第4层(566.8388,522.6232, 17.0697)4.8460第5层(566.8671,522.6003, 21.7094)4.6704第6层(566.9561,522.5512, 26.2110)4.2944第7层(566.9892,522.5289, 29.82

16、46)4.1544第8层(567.0417,522.4967, 33.3399)3.9316第9层(567.0940,522.4640, 36.8438)3.7100第10层(567.1484,522.4100, 40.1611)3.5274第11层(567.1908,522.3695, 44.4326)3.2917第12层(567.2329,522.3298, 48.6998)3.0568第13层(567.2814,522.2844, 52.8184)2.8309塔顶(567.3165,522.2518, 55.1067)0.01102011年第1层(566.7448,522.7003, 1

17、.7632)5.4270第2层(566.7792,522.6718, 7.2905)5.2163第3层(566.8122,522.6458, 12.7269)5.0104第4层(566.8392,522.6228, 17.0520)4.8460第5层(566.8677,522.5997, 21.7039)4.6704第6层(566.9567,522.5504, 26.2045)0.0133第7层(566.9900,522.5281, 29.8170)4.1544第8层(567.0425,522.4959, 33.3366)3.9316第9层(567.0950,522.4630, 36.8222

18、)3.7100第10层(567.1495,522.4089, 40.1441)3.5274第11层(567.1919,522.3684, 44.4249)3.2917第12层(567.2342,522.3285, 48.6839)3.0568第13层(567.2827,522.2829, 52.8131)2.8309塔顶(567.3375,522.2135, 55.087)0.0133表1 古塔各层中心坐标模型可视化结果如图3：各年份拟合曲线俯瞰图，由图可以看出，拟合曲线完美的拟合了图1中的观测数据，结果准确。图3 古塔拟合曲线俯瞰图5.2问题2 模型建立与求解（1） 古塔倾斜分析思路：线性回

19、归是一种有效的机器学习算法，根据问题1中求得的古塔各年中心的坐标，通过线性回归的方法对中心坐标进行拟合，求得三维空间的一条直线，根据垂直投影的方法求得直线的斜率，从而求出古塔的倾斜度。（2）数学模型：对问题1中求得的数据进行从新整理，例如1986年的数据如下表：Z(即R)x0x1(即A)x2（即B）1.78741566.6650522.70907.32021566.7225522.671312.75521566.7786522.634617.07831566.8226522.605621.72051566.8698522.574626.23511566.9185522.544329.83691

20、566.9521522.527133.35091566.9846522.510236.85491567.0176522.493340.17211567.0476522.4789假设目标拟合直线方程为：hx=0+1x1+2x2=Tx 其中x=1x1x2， =012 定义损耗函数：J(）= 12mi=1mhxi-zi2则满足min（J)的即为所求的值。（3）模型求解：运用梯度下降法求的值。梯度下降法：Repeat0：0- 1m i=1m(h(xi)-zi)x0(i)1: 1- 1m i=1m(h(xi)-zi)x1(i)2: 2- 1m i=1m(h(xi)-zi) x2(i)重复迭代，直到找到使

21、J(）最小的。使用matlab仿真（程序见附表2）求得的值如下表：观测年份值1986年30.-7.1996年30.-7.2009年30.-8.2011年30.-8.进而求得拟合直线方程如下表：观测年份所求得的空间拟合直线1986年hx=30.+9.40016x1-7.x21996年hx=30.+9.415x1-7.x22009年hx=30.01871+8.x1-8.26129x22011年hx=30.00798+8.x1-8.x2（4）结果分析：O(x2,y2)(x1,y1)(x0,y0)xyz取所求直线和xoy面的交点x0,y0，直线上任意一点x1,y1，以及点x1,y1在xoy面上的投影点

22、x2,y2。则可求得古塔倾斜角=arctan(y1-y2)/(x2-x0)2+(y2-y0)2)5.2.1古塔弯曲、扭曲度建模求解思路通过将问题1中求得的各层中心坐标投影到xoz和yoz平面内，然后利用曲线拟合的方法，对投影点进行曲线拟合，最后利用曲率的相关数学方法对拟合曲线进行曲率计算。进而求得古塔的弯曲、扭曲程度。（1）数学模型对问题1中求得的数据进行从新整理，例如1986年的数据如下表：Zx0x1(即A)1.78741566.66507.32021566.722512.75521566.778617.07831566.822621.72051566.869826.23511566.918

23、529.83691566.952133.35091566.984636.85491567.017640.17211567.0476以投影到xoy坐标平面为例，假设目标拟合曲线方程为：h(x)=0+1x+2x2+3x3定义损耗函数：J(）= 12mi=1mhxi-zi2则满足min（J)的即为所求的值。（2）模型求解：运用梯度下降法求的值:梯度下降法：Repeat jj-jJ()重复迭代，直到找到使J(）最小的。利用曲率知识定量分析古塔弯曲、扭曲程度：曲率知识介绍如图1所示，设曲线上一段弧长为，在点作切线，在点作切线 ，记切线转过的角度。对于同样的弧长，切线转过的角度。显然，比弯曲程度大。因此弧

24、长相等时，转角越大，曲率越大。图1另一方面，如图2所示，若转角相等时，弧长越短，弯曲程度越大。图2如上图所示，如果曲线的长度为，两点切线正向所夹的角为，将转角与弧长的比值定义为曲率。下面给出曲率的计算公式设曲线的方程为，且具有二阶导数，则曲线的曲率为：5.3问题3 形变趋势分析由问题2的建模分析、理论验证以及matlab求解和问题3的仿真作图可得塔的上半部分呈下降趋势，而塔的底部呈上升趋势，由此可知塔身呈一定程度的倾斜。通过塔中心坐标的研究，发现塔身呈一定程度的弯曲和扭曲，需要引起有关部门的重视，加强古塔文物的保护和监测，为古塔提供有力的保障。八、模型改进及评价8.1 模型的改进8.2 模型的

25、评价优点：模型方法比较具有通用合理性，结果准确，计算精确，科学合理缺点（结合模型假设）九、模型推广 参考文献：1黄强,古塔变形监测的探讨，上海，测绘与空间地理信息，2013.62 刘珂，周富强，张广军，半径约束最小乘圆拟合方法及误差分析，光电子激光 ，2006.53 梁海奎，古塔变形测量方法探讨，城市勘测，2011.64何凭宗，圆形立柱倾斜度测定的解析方法，矿山测量，2001.12附件附件1：function R,A,B,C=fitted_model(x,y,z,N)%fitted_model.m 最小二乘圆拟合% x 测量点x坐标值% y 测量点y坐标值% N 测量数据集中数据点个数% R

26、拟合得到的圆半径% A 拟合得到的圆心横坐标% B 拟合得到的圆心纵坐标x1 = 0;x2 = 0;x3 = 0;y1 = 0;y2 = 0;y3 = 0;x1y1 = 0;x1y2 = 0;x2y1 = 0;for i = 1 : Nx1 = 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 = x2

27、y1 + x(i)*x(i)*y(i);endC = 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 坐标C = mean(z); % z 坐标R = sqrt(a * a + b * b - 4 * c)/2;附件2（线性回归matlab代码）：function liner_regressioncle