2013数学建模C题古塔问题探究.doc

古塔的变形监测与预测摘要：本文针对古塔长时间承受自重和外力作用引起的组合变形问题，采用数据处理、几何分析及曲线拟合的方法，利用matlab、lingo等数学软件编程、计算，给出了确定古塔各层中心位置的通用方法，根据所给的各次测量数据计算出古塔各层的中心坐标（见表9），分析了该塔的倾斜、弯曲、扭曲等变形情况并对变形趋势做了预测。得出如下结论：1、随着时间的增长，古塔的倾斜弧度与时间呈二次抛物线分布，总的变形是先减小后增大趋势。2、随着时间的增长，古塔的弯曲曲率与时间呈三次函数分布，曲率随时间先增大，后减小，再增大的趋势。3、随着时间的增长，古塔的扭曲角弧度与时间呈二次抛物线分布，总的扭曲变形是先减小后增大趋势。关键词：中心坐标 投影 组合变形 倾斜 弯曲 扭曲 1、 问题重述由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，防止其变形而导致破坏，文物部门需适时对古塔进行观测，通过测量所得数据可以了解各种变形量，以制定必要的保护措施。某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。请根据提供的4次观测数据，讨论以下问题：1. 列表给出各次测量的古塔各层中心坐标，给出确定古塔各层中心位置的通用方法。2. 分析该塔倾斜、弯曲、扭曲等变形情况。3. 分析该塔的变形趋势。 二、模型假设1、 假设塔的每层是一个平面，不考虑楼层厚度；2、 假设古塔在2011年3月前没有进行过大的维修工程，如地基处理；古塔加固；3、假设观测过程中测的数据精确度高，误差很小，可以忽略不计；4、假设把各种外界因素对塔的影响看做一个整体来研究；5、假设古塔在受到外界因素干扰时，未受到严重破坏；三、符号说明：第i层观测点的横坐标；：第i层观测点的纵坐标；：第i层观测点的竖坐标；G：表示的各层中心坐标；：第1层塔中心塔尖中心连线的倾斜角；：塔尖到第1层竖坐标差；：曲线弯曲曲率；：扭转角；四、模型准备 观察所给的数据，发现 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层的结构相似，而塔尖和塔层观测点数目不同，而且结构各异，由此我们分别讨论塔层和塔尖，找到它们的中心。图11 第1层到13层1.1 问题分析每层塔的8个观测点的z不全相同，它们构成一个不规则的曲面。首先，将8个观测点投影到xoy面，投影点分别为。利用matlab拟合一个斜平面，满足：8个观测点与斜平面的交点与对应观测点的距离之和最小。依次连接这8个交点，所得图形的中心就是要找的塔层中心。再将xoy面的投影点依次连接，利用lingo软件拟合出一个到8个顶点的距离之和最小的圆，圆心G可以近似看做该八边形的中心。最后，过圆心作平行于z轴的直线，它与斜平面的交点就是我们所要找的塔层的中心点。1.2模型建立设拟合的斜平面方程为： （1）建立模型： min （）设xoy面内拟合的圆的方程为： （2） 建立模型：min 1.3模型求解根据已知数据，利用matlab拟合可以算出A，B，C，得到斜平面的方程，做出图像（如图2）；利用lingo软件拟合，可以得到圆的方程，如图3。将圆心坐标，代入斜平面的方程可解出，即得到中心坐标的通式 图2图3将四年测量的数据代入之后得到1到13层的中心坐标（如表1）表1年份19861996坐标层数XYZXYZ1566.665522.7091.7874566.6652522.70871.7832566.7225522.67147.3203566.7234522.67047.31463566.7787522.634812.7553566.7803522.633112.75084566.8227522.605817.7553566.8248522.603617.07525566.87522.574621.7206566.8727522.571821.71616566.9188522.544326.2351566.9222522.540826.22947566.9522522.527129.8369566.9561522.523229.83238566.9847522.510333.351566.9888522.50633.34559567.0176522.493436.8551567.0223522.488436.848410567.0476522.479240.1723567.0527522.473940.167811567.1014522.438844.441567.1071522.43344.435512567.1554522.398448.7118567.1614522.392248.707313567.202522.3852.8288567.2082522.373752.8243塔尖567.3785522.113455.1818567.3375522.266355.244年份20092011 坐标层数XYZXYZ1566.7447522.70031.7645566.7448522.70021.76322566.7790522.67197.3089566.7792522.67177.29053566.8117522.645412.7322566.8121522.645712.72684566.8386522.623217.0697566.839522.622717.05195566.8666522.600221.7093566.8672522.599621.70386567.0103522.639626.2085566.8967522.7426.20037567.2446 522.791929.8176566.9164522.746629.81328566.9645522.712833.3353566.9653522.71233.3329567.0128522.677236.8391567.0138522.676336.816710567.1112522.521340.1583567.0686522.622540.138111567.1909522.3743.4323567.192522.368944.424512567.2331522.330248.6994567.2344522.328948.683613567.2817522.284652.818567.283522.283252.8127塔尖567.336522.214855.091567.3375522.213555.0871.4模型分析与评价该模型的建立，采用投影、拟合的方法将不规则转化为规则，将空间转化到平面，近似找到塔层的中心坐标，由于第1到13层的结构相似，因此方法具有通用性。即使1986年和1996年的测量数据中缺失了第13层观测点5的数据，但是我们根据其他各层的数据猜测得到，不影响方法的使用。但是也因为猜测和两次拟合，塔层中心坐标的计算产生误差。2 塔尖2.1 问题分析考虑到塔尖的横截面面积相对塔层横截面较小，且已知的测点数据表明，塔尖各测点变形相对较小，因此我们采用求平均值的方法来计算塔尖的中心坐标。2.2 模型建立将塔尖的4个观测点坐标设为，中心坐标的计算式为： 2.3 模型求解将四次观测数据分别代入平均值计算的中心坐标公式： 算出塔尖的中心坐标（如表2）表2年份19861996坐标层数XYZXYZ塔尖567.3785522.113455.1818567.3375522.266355.244年份20092011 坐标层数XYZXYZ塔尖567.336522.214855.091567.3375522.213555.0872.4 模型分析与评价该模型可以简单、方便地计算出塔尖的中心坐标，但是2009年和2011年观测点只有1个，我们只能将此观测数据作为塔尖的中心坐标。而且平均值计算的方法有它的缺点，观测点的数目越少，结果的误差就会越大。这是该模型有待改进的地方。（二）问题二1 问题分析在古塔刚修建成功时，不考虑修建过程中的误差，各层中心的连线（中性轴）应该是垂直于底平面的一条垂线。但是由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生组合变形，包括倾斜、弯曲、扭曲等，此时其中性轴也会发生相应的变形（见图4）。因此，我们可以通过对古塔中性轴的变形来讨论古塔的变形，分别单独讨论倾斜、弯曲和扭曲。图4（中心点连线的图）1.1 倾斜变形1.1.1问题分析塔的倾斜变形是一种塑性变形，是塔基固定，塔身和塔尖产生倾斜的一种变形。这种变形反映了整座塔的变形，因此我们考虑塔尖中心相对于第1层中心的倾斜。连接第1层的中心和塔尖的中心，将其投影到过且平行于z轴的直线上，的投影点到的距离记作，到的距离记作，倾斜角。（如图5） 图51.1.2模型建立1.1.3 模型求解将第一题中计算出来的数据带到以上公式，分别计算出四年古塔的倾斜角度（表3）表3年份1986199620092011倾斜位移值0.92940.80480.76510.7669倾斜度0.01680.01460.01390.0139倾斜角5949.995144.874919.174926.281.1.4模型分析与评价根据上述表格所得的数据可知，古塔的倾斜角度是很微小的。在古塔开始发生倾斜时，由于地基基础下的地质分布不同及其他外界不利因素，导致古塔刚开始是向一边倾斜。但是随着时间的变化，以及各种作用下，古塔开始向其另一边发生倾斜，从而抑制塔向一边倾斜的趋势。由二力平衡可知，两边对塔产生力的作用，这样使塔处于相对平衡的状态。由此我们可以从上述计算出的结果看出，在没有大的事故作用下，塔在一定的范围时间内不会被外界作用所毁坏。2.1 弯曲变形2.1.1问题分析曲线的弯曲程度可以用曲率来度量，依次连接塔层中心点，得到的是空间折线段，将它拟合成一条曲线，参数方程为：。代入曲率公式可以计算出各点的曲率，由于曲线各点的弯曲程度不同，我们用曲线的最大曲率来度量古塔的弯曲变形程度。2.1.2模型建立曲线方程：max 2.1.3模型求解 利用第一题计算的中心坐标，运用matlab软件拟合出曲线（见图6）的方程，计算出参数方程的一阶导数、二阶导数代入曲率公式，得到再利用matlab对k求导数，计算出当时，k取得最大值。图6将第一题求得的中心坐标带入，按照上述方法计算，可以分别计算不同年古塔的最大曲率（见表4）。表4年份1986199620092011曲率1.030310-97.873610-103.616010-78.00310-73.1扭曲变形3.1.1问题分析 扭曲变形是指各塔层横截面绕轴线的相对转动，转动后与其原来的位置形成一定的角度，称这个角为相对扭转角。考虑古塔的扭曲是绕着中性轴的转动，而每层扭转的程度不一样，因此分层讨论。假设第1层为基准面，不发生扭曲，考虑其他各层相对于它绕着中性轴的扭曲程度。而塔尖的数据不全，且近似看作在中性轴上，因此只考虑第2层到第13层的扭曲。每一个塔层有8个观测点，任选一个作为扭曲变形的研究对象，不妨假设观测点1。将第1层的观测点1投影到第i层中心点所在的水平面，投影点记作，分别连接与以及第i层观测点1（）与，它们之间的夹角就是要找的相对扭转角（如图7）。图73.1.2模型建立已知坐标，, ，利用余弦定理列出计算式(i =2、3.13) 。 ， 3.1.3模型求解将题中所给的数据以及第一题计算的中心坐标代入上式，利用matlab可以计算出每一层相对于第一层的扭转角（见表5）。表5 1986年的值1996年的值2009年的值2011年的值 0 0.0469 0.0945 0.1334 0.1766 0.2205 0.2810 0.3444 0.4131 0.4832 0.5666 0.6581 0.7597 0 0.0470 0.0962 0.1334 0.1768 0.2223 0.2811 0.3451 0.4133 0.4852 0.5673 0.6616 0.7613 0 0.0529 0.1054 0.1498 0.2012 0.2950 0.3227 0.4093 0.4807 0.5412 0.6106 0.6957 0.7886 0 0.0493 0.1053 0.1497 0.2011 0.3021 0.3460 0.4134 0.4806 0.5461 0.6103 0.6954 0.7883 （三）问题三3.1问题分析 变形趋势就是一种惯性，没有外力不会轻易改变状态。由于古塔的变形是多种影响的组合变形，如倾斜、弯曲、扭曲等等。这些影响来自各个方面，不便于综合确定各种影响的程度，这样就造成了无法整体去研究它的变形趋势。因此我们还是分开考虑，分别考虑倾斜、弯曲、扭曲对塔的变形趋势。这些变形趋势都是建立在问题二的基础之上，模型相同，因此在以下变形趋势中同样适用。3.1.1倾斜变形趋势分析在问题二倾斜变形模型求解中，我们计算出了塔的整体倾斜角（见表3），我们再用时间（单位是年）与倾斜角（单位是弧度）建立平面直角坐标系，以O为坐标原点，时间（单位是年）为X轴,倾斜角（单位是弧度）为Y轴，坐标系表示倾斜角与时间变化的关系。所以我们可以得出四个坐标点，再通过Excel处理拟合曲线作图（见表6、图8）。表6年份1986199620092011倾斜角(弧度)0.99720.86250.8220.824图8根据上述图表预测，随着时间的增长，古塔的倾斜弧度与时间呈二次抛物线分布，当（年）时塔的倾斜角最小，即倾斜程度最小，塔的变形最小。总的变形是先减小后增大趋势。3.1.2 弯曲变形趋势分析在问题二弯曲变形模型中，我们求解出了曲线的最大曲率，曲线的最大曲率能够帮助我们分析曲线的弯曲程度。基于此情况，我们在问题二的弯曲变形模型基础上去分析弯曲变形趋势。首先，建立平面直角坐标系，O为坐标原点，OX轴代表时间（年），OY轴代表曲率，此坐标系是表示曲率与时间的变化关系。 所以我们可以得出四个坐标点，再通过Excel处理拟合曲线作图（见表7、图9）。 表7 时间1986199620092011曲率（10-10）1.03037.873636168003图9根据上述图表预测，随着时间的增长，古塔的弯曲曲率与时间呈三次函数分布，曲率随时间先增大，后减小，再增大的趋势。3.1.3扭曲变形趋势分析在问题二扭曲变形模型中，已经求解出了各层相对于第1层的扭曲角（见表5），观察表的数据，可以找出一规律：随着塔层的增高，扭曲角增大，即扭曲程度增大。最大扭曲角来最能够代表整座塔的变形，因此我们还是选取它来研究扭曲变形趋势。同样建立平面直角坐标系，O为坐标原点，OX轴代表时间（年），OY轴代表最大扭曲角（弧度）。此坐标系是表示扭曲角随时间的变化关系。所以我们可以得出四个坐标点，再通过Excel处理拟合曲线作图（见表8、）。 表8年份1986199620092011扭曲角（弧度）0.75970.76130.78860.7883图10根据上述图表预测，随着时间的增长，古塔的扭曲角弧度与时间呈二次抛物线分布，当（年）时塔的扭曲角最小，即是扭曲程度最小。总的扭曲变形是先减小后增大趋势。六、模型推广与评价一、优点：本文通过对数据的处理、分析，利用MATLAB和LINGO软件对古塔每年每层的测量点数据进行拟合，有效地将不规则立体转化为规则的平面图形，利用规则平面图形的中心来找到塔层的中心。这样的处理很方便、明了，便于实现。此模型可用于大多数文物建筑，通用性较强。二、缺点：在数据的处理中，由于对数据采用了“四舍五入”方法，并且预测补充缺失数据，多次拟合，从而使数据误差增大，降低了其精确性。 7、 参考文献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)*y2、运用lingo软件对数据处理的程序：model:sets:dianji/1.8/:x,y,d;endsetsdata: 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;enddatamin=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、 Variable Value Reduced Cost X0 567.0476 0.000000 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）：年份19861996坐标层数XYZXYZ1566.665522.7091.7874566.6652522.70871.7832566.7225522.67147.3203566.7234522.67047.31463566.7787522.634812.7553566.7803522.633112.75084566.8227522.605817.7553566.8248522.603617.07525566.87522.574621.7206566.8727522.571821.71616566.9188522.544326.2351566.9222522.540826.22947566.9522522.527129.8369566.9561522.523229.83238566.9847522.510333.351566.9888522.50633.34559567.0176522.493436.8551567.0223522.488436.848410567.0476522.479240.1723567.0527522.473940.167811567.1014522.438844.441567.1071522.43344.435512567.1554522.398448.7118567.1614522.392248.707313567.202522.3852.8288567.2082522.373752.8243塔尖567.3785522.113455.1818567.3375522.266355.244表9年份20092011 坐标层数XYZXYZ1566.7447522.70031.7645566.7448522.70021.76322566.7790522.67197.3089566.7792522.67177.29053566.8117522.645412.7322566.8121522.645712.72684566.8386522.623217.0697566.839522.622717.05195566.8666522.600221.7093566.8672522.599621.70386567.0103522.639626.2085566.8967522.7426.20037567.2446 522.791929.8176566.9164522.746629.81328566.9645522.712833.3353566.9653522.71233.3329567.0128522.677236.8391567.0138522.676336.816710567.1112522.521340.1583567.0686522.622540.138111567.1909522.3743.4323567.192522.368944.424512567.2331522.330248.6994567.2344522.328948.683613567.2817522.284652.818567.283522.283252.8127塔尖567.336522.214855.091567.3375522.213555.0873、根据空间参数方程： 通过matlab软件对上述数据进行拟合： x=data(:,1); z=data(:,3); polyfit(z,x,2); y=data(:,2); z=data(:,3); polyfit(z,y,2);得到具体的空间曲线方程（以1986年的为例）：通过下面这个程序，画出1986年的空间曲线： t=0:100; x=1.556978433478483e-05*t.2+0.0095*t+566.6518; y=-3.775091994799030e-06*t.2 -0.0061*t+ 522.7157; z=t; plot3(x,y,z,r-);然后用此方法分别对每年的数据进行拟合，这样我们就可以分别得到每年空间曲线方程和空间曲线。用diff命令对空间曲线方程求导：syms tdiff(1.556978433478483e-05*t2+0.0095*t+566.6518)diff(1.556978433478483e-05*t2+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