2013高社杯全国大学生数学建模竞赛C题

2013高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

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

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

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

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：
参赛队员 (打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人 (打印并签名)：
日期：年月日赛区评阅编号（由赛区组委会评阅前进行编号）：
2013高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：
全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：
古塔的变形
摘要：
本文研究的古塔的变形问题，通过对问题背景及附件资料进行深入地分析，采用数据拟合、求平均值等方法整理出具有科学性的分析数据。

通过对建筑物位移监测数据处理方法的研究, 采用自回归模型对位移监测数据进行处理, 根据建立的模型对具体建筑物的监测点
的位移变化量进行预报。

经过计算分析, 根据位移量之间变化的关系而建立的自回归预测模型具备较高的拟合及预测精度,运用三维坐标系和数学软件将古塔的模型以空间模型的形式表现出来，直观且科学，对于研究古塔的变形具有较高的科学性和说服性。

再通过三维坐标之间的回归和三维坐标与时间的回归而分析出古塔的倾斜，弯曲，扭曲等变形状况，通过数学软件的计算及列表列图的方法将结果直观体现，通过大量的计算与分析，运用几何和代数方法将古塔的变形量以数学的方式说明。

对于分析古塔变形趋势中，运用了位移差和位移残差平方公式等量及与时间的关系来说明其变形趋势。

对于问题一，通过对监测数据的分析，得出此塔为八边形的塔，并通过平均值法求出古塔各层的中心坐标，具体见表（一）。

对于问题二，通过问题一对变形监测数据的研究和处理，我们组运用了自回归模型的方法，利用Z和X,Y之间的回归关系，我们运用数学软件计算求出a1和a2，并通过代数及其几何关系，求出每年
监测出的古塔倾斜角度。

具体见图（2）和表（2）。

对于问题中的弯曲问题，我们用古塔的中点高度发生变化的多少来表示弯曲程度，由于弯曲程度主要是随着时间的变化而变化，所以我们用时间和变形的监测数据进行回归拟合得出△Z（表示古塔的弯曲程度），经过数学软件的大量计算，用列表的方式将每次监测所得每层的弯曲程度表示出来，具体见表（3）。

对于问题中的扭曲问题，我们用第一层作为基层，即不扭动层，其他层相对于一层扭动了多少度来说明古塔的扭曲程度。

为了这一说明，我们取每一层第一点和第五点作直线。

并将每一层的直线与第一层的直线做对比，求出两直线之间的夹角，并用此夹角来说明除第一层以外的每一层相对于第一层扭动了多少度。

以此角度大小来表现古塔扭曲程度。

具体见表（4）
对于问题三，是分析古塔的变形趋势，根据监测的变形数据和位移与时间及波动稳定性的关系，列出了一个时期对上一时期的位移差，看出古塔变形的趋势是向哪个方向，具体见图（5）。

再通过位移平方差公式，对古塔整体变形的趋势波动进行分析，并结合时间等因素推算出古塔变形趋势在增加。

具体见表（5）。

关键词：变形数据拟合平均值自回归模型位移变化
三维坐标系数学软件几何代数方法位移差位移
平方差
一、问题重述
1.1基本情况：
某古塔在我国已有上千年的历史，是我国的重点保护文物。

但由于长时间受到自重、气温、风力、地震、飓风的影响，古塔出现了各种变形，如倾斜、弯曲、扭曲等。

为保护古塔，文物部门需要适时对古塔进行观测，了解各种变形量，以制定必要的保护措施，因此管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该古塔进行了4次观测。

1.2需解决的问题:
问题一：根据附件1提供的4次观测数据，给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。

问题二：利用数学模型和所得数据分析该塔倾斜、弯曲、扭曲等变形情况。

问题三：分析该塔的变形趋势。

二、问题分析
2.1问题一的分析：
对于问题一，通过观察所给的数据，每一层的8个观测点都位
于古塔的每层的八个角落。

据此，将古塔的每一层类似看做一个正八边形，可用求正八边形的中点的方法来确定古塔每一层的中点坐标。

通过所给每层的观测数据，分别求X，Y，Z坐标的平均值，近似作为古塔各层中心坐标。

2.2问题二的分析：
对于问题二，主要是研究古塔的倾斜，弯曲，扭曲等变形量，据此分三点来分析问题二。

对于倾斜，主要是三维坐标中X，Y轴对于Z轴的倾斜角度，即是对中心轴的倾斜角度，而第一题已经对每年的监测数据进行了总结，得出每次监测各层的中心坐标，故运用X，Y与Z的自回归模型，再通过代数和几何关系，求出每年古塔X，Y轴的倾斜角度，并加以说明，具体结果见表（2）。

对于弯曲，主要是要找出能表示弯曲程度的量，经过查找变形量具体分析的资料，知道弯曲主要是中点的高度（Z）发生了变化，列出时间和每层的中心坐标的回归拟合，并且通过数学软件计算，得出弯曲程度△Z，并列出表直观说明，易于比较。

具体结果见表（3）。

对于扭曲变量，通过matlab软件做出的古塔空间图形，及查找了扭曲问题的分析的资料，我们最终采用以第一层为基层，其他层相对于古塔扭曲的度数，主要是通过每层的两个点算出每层的斜率，再通过斜率公式算出角度，以此说明古塔每次扭曲程度的变化及大小。

具体
见表（4）
2.3问题三的分析：
对于该塔的变形趋势，经过第一题和第二题，已将古塔的变形位置的中心和变形量进行了具体分析。

本题要分析变形趋势，具体是从古塔各中心点在各个时期的位移变化量，通过列位移残差公式，得出该塔的变形趋势。

见图（5）。

再通过整体的位移残差平方公式，得到古塔每层在某个时期的变形趋势波动，具体见表（5）。

三、问题假设
１．假设该塔为正八边形的塔
２．假设该塔的底层不扭曲变形
３．假设监测该塔时后两次改变了监测位置
４．假设地质的变动对古塔无影响
四、建立模型及求解
（一）符号说明：
为古塔与X轴构成的角度
1∠A
1
2∠A
为古塔与Y轴构成的角度
2
3 Z为中点的高度
4 △Z 来表示古塔的弯曲程度
5 Bi 为表示两直线的夹角
6△X 为X 轴的位置偏移量
7△Y 为Y 轴的位置偏移量
8△Z1为Z 轴的位置偏移量
9△S 2为古塔整体的位置偏移平方差
10 i 为古塔层数
11 t 为从1886年开始所经过的时间
（二）模型建立及求解
3.1 对问题1的求解
根据数据中给出的各层各个点坐标的值，可把古塔的每一层类似的看成一个正八边形。

根据数学逻辑和模型，正八边形的中心坐标可用平均值法求出，故塔的中心坐标为（-
--Z Y X ，，）。

附件一是每次测量出的古塔的各层的各个点的坐标，通过平均值法整理数据，得出以下结果：
表（一）
其中的X=(x1+x2+……+x8)/8得来，Y，Z值同理。

通过matlab做出古塔的模型图片：
下图（1）
图（1）
3.2对问题二的求解：
对于古塔的倾斜，弯曲、扭曲等变形情况的分析，第一题通过科学分析和数学逻辑思维，用平均值法对建筑物位置的监测数据进行总结，得出各个时期各个层的中心坐标值，通过这些数据，运用过对建筑物位移监测数据处理方法的研究, 采用自回归模型对位移监测数据进行处理, 根据建立的模型对具体建筑物的监测点的位移变化量进行预报。

经过计算分析, 对古塔的变形情况进行说明。

3.2.1关于古塔的倾斜问题，我们建立了空间直角坐标（如上图1）对此进行了讨论，发现古塔与X轴构成的角度为∠A
1
,与Y轴构成的
角为∠A
2。

因此建立了-Z i,t与-X i,t、-Y i,t的自回归模型如下式○1,
且∠A
1=arctan a1 ，∠A
2
=arctan a2 。

-
Z i,t=a0+a1*-X i,t+a2*-Y i,t ○1从问题一中所求得的平均值可得出下图2：
图2
因此，在所观测的4年里古塔对于X轴和Y轴的倾斜角度如下表（2）
表（2）
由分析和计算得出每次监测时对中心轴的倾斜角度都有变大的趋势，
古塔的变形不容小觑。

3.2.2 弯曲
关于古塔的弯曲，我们考虑到古塔弯曲后，中点的高度会发生改变，即Z发生了变化，故我们用△Z来表示古塔的弯曲程度，即
△z-= z'-Z
z'=a t i,+b t i,*t
具体计算值如下表（3）
表（3）
表中表示的是古塔每次每层的中心点的变化及弯曲程度，对于1986和1996及2009和2011的比较，前者（即1986和1996）弯曲程度明显加大，这不仅由于时间长达十年的关系，还有各个方面的影响，比如倾斜，风化，地壳运动等对变形位移数据的影响，肯定的是古塔的弯曲程度在很小很小范围内变化，后者（2009与2011）虽没有前者那么明显，但是考虑到时间因素的话，平均每年的弯曲程度的变化还是高于前者的，所以古塔的安全维护必须尽快落实做好。

3.2.3 扭曲
对于古塔的扭曲,我们取每一层第一点和第五点作直线。

并将每一层的直线与第一层的直线做对比，求出两直线之间的夹角Bi(下图为第一层和第十三层的事例)，并用此夹角来说明除第一层以外的每一层相对于第一层扭动了多少度，计算公式如下：
B i=arctan(Y2-Y1)/(X2-X1)
L2
L1
具体数值如下所示：
1986(1)
k Bi 565.454528.012567.941517.407-4.26417
565.48527.764567.995517.563-4.05606-0.732
表（4）
此问题通过最简单的斜率计算得出古塔每次每层相对于基层的扭曲
度数，用数字说明了扭曲程度的变化。

3.3对问题三的求解：
问题三是分析古塔的变形趋势，我们通过每个时期与上一时期古塔每层的中心位置的偏移，即偏移量△X，△Y，△Z1；又△X=X2-X1; △Y，△Z1亦同理求出其值，具体见图（5）。

再通过算位移残差平方公式△S2=△X2+△Y2+△Z12,可以求出古塔每层每两个时期的变形波动，对于变形趋势的稳定性和变化有重要作用。

具体见表（5）。

2,
1986与1996的比较
2009和2011的比较
图（5）
通过△X，△Y，△Z1的值大小，能看出古塔变形的趋势向哪个轴的方向偏移多一点，由两次比较得出古塔的变形趋势向Z轴，即向中心线的斜方向变形。

下面是位移平方差的结果：
1986与1996年的位移差平方
2009和2011的位移平方差
表（5）
通过数学软件对四次测量结果的位移平方差的对比，其变形波动率虽然小，但是其波动已经能通过数字表现出来了，说明古塔的变形随时间变化会越来越明显，而通过观察计算出来的结果，虽然1996年和1986年的位移平方差比2011和2009的大，可是从年平均变形稳定的角度，2011和2009年的还是比之前明显的，由此可见古塔变形趋势的波动在增大，望有关部门能重视。

五模型检验
1、问题一我们采用求平均值的方法，求古塔每层中心点坐标，由数据中给出的8个基点的位置坐标去求，又假设其为正八边形的塔，所以用代数方法直接求出，给出其通用算法并将结果以列表形式表达出来，科学合理。

2、对于问题二，通过查找变形资料的数据处理，还有运用自回归方程及各个潜在变量的回归关系，再通过代数和几何关系将其各个变形量用其他关系量形象的表达出来，具有很强的数学逻辑思维和建模思维。

3对于问题三，通过时间及位移差经典公式，将已经处理好的具有科学性分析的数据代入，通过计算位移差和整体位移平方量，得出其古塔的变形趋势的方向和程度，其应用价值在防止古塔倒塌方面及安全保护预测具有实质性的作用。

六、模型评价
（一）模型的优点：
1.建立的模型能与实际紧密相联系，综合实际情况对问
题进行求解，使模型具有很好的通用性和推广性；
2.模型的计算采用专业的数学软件，可信度高；
3.对附件中的表格进行处理，找出了许多变量之间的潜在关系；
4.对模型中涉及的众多影响因素进行量化分析，使得论文有说服力。

（二）模型的缺点
1.在编程中没有加入题目中的约束条件，导致最终的运算结果出
现小数，最后我们采用人工方法进行较好的弥补。

2、在计算倾斜角度及等变形量时肯定会有出入，但在模型计算
中，我们去的预计值做为近似值为计算，这与实际值有出入
七、意见
为确保古塔安全，应紧密结合经常性的目测检查与专门的全面检查，以便发现问题，及时处理。

目前塔体的变形，有些是可逆的，如由温、湿度造成的物理变化；有些是不可逆的，如转表层的风化、墙面裂缝、砖间黄泥层的压缩等。

其中由于自然环境诸如温度、湿度、风雨、虫鸟草树等多种因素的影响所造成的塔体表面风化日趋严重，这些不可逆变化累积到一定数量后就会产生质变，导致塔体整体变形。

例如该塔应定期进行检测，防范雨水冲刷、虫鸟草树侵扰，封堵土蜂洞穴，堵塞墙面裂缝，进行经常性的维护。

八、参考文献
1] ?? 黄声享, 尹晖, 蒋征. 变形监测数据处理[M ] . 武汉: 武汉大学出版社, 2007
[ 2] ?? 王佳璆. 时空序列数据分析建模[ D] . 广州: 中山大学硕士学位论文, 2008
[ 4] ?? 金彪, 吴北平, 李艳芳. 曲线拟合与自回归模型在地铁变
形监测中的运用[ J]. 地矿测绘, 2009, 25( 1): 35~ 37
[ 5] 自回归模型在沉降预报中的应用
附录:
close all
clc
clear
x1=[561.4782 563.9092 567.9981 571.3812
572.0142 569.5162 565.4062 562.1112; 561.7463 564.1008 568.0171 571.243 571.8122 569.4066 565.4625 562.3235; 562.0086 564.289 568.0358 571.107 571.6131 569.3 565.5175 562.5296;
562.2183 564.4385 568.0507 570.9988 571.4555 569.2151 565.562 562.6956; 562.4427 564.5997 568.0667 570.8829 571.2857 569.1237 565.6087 562.8727; 562.8807 564.9208 568.1219 570.6945 571.0345 568.984 565.7192 563.2215; 563.0424 565.0389 568.1419 570.6257 570.941 568.93
565.7612 563.3528; 563.3012 565.2289 568.1741 570.5146 570.7907 568.8459 565.8278 563.5608; 563.5593 565.4173 568.2063 570.4043 570.6414 568.7614 565.8943 563.7683; 563.7822 565.5507 568.1927 570.2892 570.5166 568.7563 566.0125 563.9956; 564.0405 565.6943 568.1455 570.1131 570.3292 568.7197 566.135 564.2616; 564.2984 565.8377 568.0986 569.9368 570.1419 568.6831
566.2569 564.5268; 564.5544 565.9834 568.0609 569.7744 569.9704 568.6554 566.3814 564.7894]
y1=[521.4177 518.0897 517.4028 519.8707
523.9807 527.3387 527.9887 525.5227; 521.4074 518.2185 517.5595 519.9587 523.9056 527.1539 527.7745 525.3761; 521.3974 518.345 517.7137 520.0452 523.8316 526.975 527.5655 525.2338; 521.3894 518.4456 517.8365 520.114 523.773 526.831 527.397 525.119; 521.3808 518.5539 517.9679 520.1879 523.7099 526.6769 527.2189 524.9969; 521.363 518.8401 518.318 520.3849 523.6306 526.288 526.7638 524.6852; 521.3565 518.9453 518.4445 520.456 523.6 526.1379 526.5906 524.5678; 521.346 519.1146 518.6487 520.573 523.5537 525.9031 526.3165 524.382; 521.3356 519.2826 518.8526 520.6886 523.5066 525.6676 526.0426 524.1966; 521.3422 519.3767 518.9714 520.7323 523.4077 525.462 525.8122 524.045; 521.3802 519.5287 519.1543 520.8171 523.3109 525.2229 525.5456 523.8776; 521.4182 519.6805 519.3369 520.9021 523.2143 524.9848 525.2806 523.7107; 521.4476 519.8195 519.5048 520.9766 523.1144 524.7485
525.0184 523.5424]
z1=[1.777 1.802 1.773 1.757 1.765 1.723 1.76 1.759; 7.322 7.335 7.298 7.287 7.292 7.302
7.31 7.326; 12.745 12.769 12.734 12.725 12.739 12.693 12.724 12.729; 17.081 17.088 17.065 17.052 17.049 17.048 17.092 17.083; 21.721 21.741 21.697 21.691 21.695 21.699 21.706 21.725; 26.251 26.285 26.293 26.248 26.182 26.095 26.132 26.202; 29.864 29.894 29.893 29.851 29.776 29.734 29.75 29.835; 33.377 33.421 33.406 33.374 33.299 33.235 33.259 33.348; 36.882 36.919 36.911 36.874 36.799 36.746 36.766 36.853; 40.195 40.21 40.226 40.217 40.165 40.022 40.109 40.145; 44.467 44.477 44.492 44.477 44.438 44.3 44.389 44.421; 48.737 48.745 48.758 48.746 48.702 48.57 48.653 48.687; 52.859 52.873 52.878 52.873 52.796 52.686 52.773 52.809]
x=[x1,x1(:,1)]
y=[y1,y1(:,1)]
z=[z1,z1(:,1)]
d=[567.336 522.2148 55.091]
hold on
for i=1:13
plot3(x(i,:),y(i,:),z(i,:),'k.') plot3(x(i,:),y(i,:),z(i,:),'k-') end
plot3(d(1,1),d(1,2),d(1,3),'k.') xlabel('x', 'FontWeight', 'bold'); ylabel('y', 'FontWeight', 'bold'); zlabel('z', 'FontWeight', 'bold'); zoom on
axis square
view(45,45)
2020-2-8。