高教社杯全国大学生数学建模竞赛C题论文穆志勇饶兵.docx
《高教社杯全国大学生数学建模竞赛C题论文穆志勇饶兵.docx》由会员分享，可在线阅读，更多相关《高教社杯全国大学生数学建模竞赛C题论文穆志勇饶兵.docx（27页珍藏版）》请在冰豆网上搜索。
高教社杯全国大学生数学建模竞赛C题论文穆志勇饶兵
2012高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。
如有违反竞赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）：
C
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：
临沧师范高等专科学校
参赛队员（打印并签名）：
1.饶兵
2.田兴平
3.熊生
指导教师或指导教师组负责人（打印并签名）：
穆志勇
日期：
2013年9月15日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2012高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
评
阅
人
评
分
备
注
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
2013高教社杯全国大学生数学建模竞赛C题论文
摘要：
文物是人类的重要历史文化遗产，它不仅反映过去，还可以借鉴未来，文物保护是一项重要的工作，不仅要有高度的文物保护意识，还应该有强有力的措施和先进的科学技术手段。
本题给出了某古塔在1986年7月、1996年8月、2009年3月和2011年3月进行的4次观测数据，数据包含了古塔每一层以及塔尖上的采样点的三维坐标值（x,y,z）。
要求确定古塔各层的中心坐标，分析古塔的倾斜和变形情况，最后再分析古塔的变形趋势，以便于加强防护措施，保护古塔。
本文根据题目附件给出的三维坐标数据，做出了各次测量数据采样点的三维散点图，画出各条相连边的连线，求出该空间多边形的中心点坐标值。
算出了各楼层的中心坐标值之后，可以画出4次测量的各楼层中心连线图，再连接上塔的顶点，就刻画出古塔的外形。
根据该连线图，我们可以分析古塔的倾斜、扭曲、弯曲等变形情况。
再将4次测量画出的中心连线图作比较，就可以知道古塔的变形趋势。
关键词：
文物保护、古塔、垂径定理、中心、线性回归
一、问题重述
古塔的变形
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。
为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。
某古塔已有上千年历史，是我国重点保护文物。
管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。
请你们根据附件1提供的4次观测数据，讨论以下问题：
1.给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。
2.分析该塔倾斜、弯曲、扭曲等变形情况。
3.分析该塔的变形趋势。
二、问题分析
根据数据的散点图初步判断古塔是至少具有八个侧面的多边形棱锥体（或者棱台体），也有可能是圆锥体（或者圆台体）。
每层楼具有八个采样点，将八个点投影在水平面（x,y平面）上，采样点的在水平面的坐标值为（x,y），将相连的两个采样点用直线连接，连线就是外接圆的玄，做出八条边（亦及八条外接圆的玄）的中垂线，中垂线的交点（或者多个交点的平均值）就是圆心的水平坐标（xo,yo）,而（外接圆的圆心）是该楼层的中心在水平面上的投影，也就说楼层中心的水平坐标值也是（xo,yo）。
再根据该楼层八个采样点的垂直坐标（z轴）的平均值，来近似的估算中心点的垂直坐标zo的值。
最后就得出了该楼层的中心三维坐标值
（xo,yo,zo）。
三、模型的假设
3.1、假设古塔每一层的水平横截都是中心对称图形。
（保证中心点坐标的平均值居中）
3.2、古塔每层的重心就是它的几何中心。
四、定义与符号说明
直线AB的斜率
直线AB的中垂线的斜率
AB与BC的中垂线的交点
xy平面内直线和y轴的夹角
θ空间直线和xy水平面的夹角
五、模型建立与求解
5.1画散点图，对塔的形体做初步判断
在第u次测量中，根据数据附件中给出的每次测量的数据，画出了各次采样点的散点图，并将古塔同一楼层的点用线连起来，例如2011年的测量数据采样点连线图如下：
根据图像初步判断塔体是至少具有八个侧面的多边形棱锥体（或者棱台体），也有可能是圆锥体（或者圆台体）。
并且假设其每个水平截面都是中心对称图形。
5.2求出古塔每层楼的中心坐标值并画出连线图
在第u次测量过程中，第j楼的八个观测点的坐标值分别为：
（xj,yj,zj）,楼层序数j=1,2,...8。
提取八个采样点在水平面（x,y）上的投影坐标值：
（xj,yj）,j=1,2,...8，连接相连的两个点，画得一个八边形，如图所示：
做出第一条变AB的中垂线，中垂线与AB交于PAB点，坐标值为（x12,y12）,根据中点坐标公式可得：
，
。
直线AB的斜率为：
，那么直线AB的中垂线的斜率为：
根据直线的点斜式方程，已知直线上的一个点，又已知直线的斜率，则AB的中垂线的方程可以确定。
方程为：
，带入
的值得：
。
同理可算得BC的中点PCB的坐标值和BC边中垂线的斜率。
PCB（x23,y23）中，
。
也可以确定了BC的中垂线的方程
。
AB和BC直线的中垂线交于
点，已知两直线的方程，联立两条中垂线的方程，可以解得它们的交点坐标值：
。
以此类推，可以算得所有边（八条边）的中垂线的方程，然后出它们两两相交的交点坐标。
如果八边形的外接圆存在，那么所有中垂线的交点都应该重合于圆心，但是由于测量误差和形变，这些中垂线的交点并不重合，因此，取所有交点坐标的平均值，作为外接圆的圆心值，也就是八边形的中心坐标值。
实现的方法是用组合，先算出第1条中垂线（即AB边的中垂线）与其它7条边的中垂线的交点，再算第2条中垂线（即BC边的中垂线）与第3条到第8条中垂线的交点，依次类推，算出所有中垂线之间的交点坐标。
如果是八条边则产生
个交点，如果是L条边则产生
个交点。
总的交点的坐标值为：
其中i,j为边的序号，并且要求
另外
是等价的，只能出现一次。
对于中心点的垂直坐标，用
其中L是边的条数，也是顶点的个数。
根据上述方法，算得每层楼的中心点坐标
，然后用直线连接，如图所示：
5.3数据和分析
从图中可以发现，第1、2次（1986年和1996年）测量的数据非常接近，图像几乎重合；第3、4次（2009年和2011年）测量的数据非常接近，图像几乎重合。
而在第2次（1996年）测量值和第3次（2009年）测量值之间，存在着较大差异。
第1次1986年测量各楼层中心坐标值列表
楼层
坐标轴
1
2
3
4
5
6
7
8
9
10
11
12
13
楼顶
x
566.68
566.73
566.77
566.81
566.85
566.89
566.93
566.96
567
567.04
567.11
567.17
567.21
567.25
y
522.7
522.66
522.63
522.59
522.56
522.53
522.5
522.47
522.44
522.42
522.39
522.36
522.32
522.24
z
1.7874
7.3202
12.755
17.078
21.721
26.235
29.837
33.351
36.855
40.172
44.441
48.712
52.834
55.123
第2次1996年测量各楼层中心坐标值列表
楼层
坐标轴
1
2
3
4
5
6
7
8
9
10
11
12
13
楼顶
x
566.68
566.73
566.77
566.81
566.85
566.9
566.93
566.97
567.01
567.04
567.11
567.18
567.22
567.25
y
522.7
522.66
522.62
522.59
522.56
522.52
522.5
522.47
522.44
522.42
522.39
522.35
522.31
522.24
z
1.783
7.3146
12.751
17.075
21.716
26.23
29.832
33.345
36.848
40.168
44.435
48.707
52.83
55.12
第3次2009年测量各楼层中心坐标值列表
楼层
坐标轴
1
2
3
4
5
6
7
8
9
10
11
12
13
楼顶
x
566.76
566.81
566.84
566.88
566.91
567.02
567.06
567.12
567.19
567.24
567.28
567.32
567.37
567.34
y
522.69
522.65
522.61
522.58
522.54
522.52
522.51
522.5
522.49
522.43
522.38
522.33
522.28
522.21
z
1.7645
7.309
12.732
17.07
21.709
26.211
29.825
33.34
36.844
40.161
44.433
48.7
52.818
55.091
第4次2011年测量各楼层中心坐标值列表
楼层
坐标轴
1
2
3
4
5
6
7
8
9
10
11
12
13
楼顶
x
566.77
566.81
566.85
566.88
566.91
567.02
567.06
567.12
567.19
567.24
567.28
567.32
567.37
567.34
y
522.69
522.65
522.61
522.58
522.54
522.52
522.51
522.5
522.49
522.43
522.38
522.33
522.28
522.21
z
1.7632
7.2905
12.727
17.052
21.704
26.205
29.817
33.337
36.822
40.144
44.425
48.684
52.813
55.087
5.3.1静态分析
由于第1、2次（1986年和1996年）测量的数据非常接近，所以选取第2次（1996年）测量的数据为例进行分析。
根据1996年数据画得古塔各楼层中心连线图，如图红色所示，为了便于观察，相对前一个图而言，将图像旋转了大约180度。
在xy水平面（z=0）、xz侧面（y=552.8时的平面）和yz侧面（x=566.5时的平面）的投影都接近直线，因此可以用线性拟合的方法，将数据拟合为空间直线。
因此对连线做了空间直线拟合，拟合后的图像如图所示。
另外因为楼顶坐标（上图中的五角星）不是某一楼层的中心，因此我们不放在拟合直线的参考范围之内。
算出了该拟合直线的空间直线方程为x=Az+C1,y=Bz+C2,据此我们可以算出直线在xy水平面内的投影方程为
，根据此方程可以算得古塔中心连线的倾斜方向与x轴和y轴的夹角。
与y轴的夹角
，
=-54.439°；那么与x轴的夹角等于
144.439°。
为了算得空间直线与xy水平面的夹角，我们在空间直线上任取两点（x1,y1,z1）,（x2,y2,z2）,课算得空间直线和它在xy平面内的夹角，可证明这个夹角就是空间直线和xy平面的夹角。
设夹角为：
θ，则
，就得
，任取两点带入，可算得
89.276°，也就说该直线偏离了垂直（z轴）方向0.72353°。
第3、4次（2009年和2011年）测量的数据非常接近，因此取第3次（2009年）测量的数据为例进行分析。
根据数据画得如下的图像：
从图中可以发现，各个楼层的中心连线图是扭曲的，在xy水平面（z=0）、xz侧面（y=552.8时的平面）和yz侧面（x=566.5时的平面）的投影都是弯曲的，说明各楼层的中心分布不是线性的，不能简单的用直线来代替。
但我们还是先用线性拟合直线的算法来算一个平均值，后面再具体分析扭曲的情况。
用前面的算法，可以算得xy水平投影与y轴夹角为
-60.120°，与x轴夹角为：
150.12°根据拟合直线的方程，可算得拟合直线和水平面的夹角为
89.180°，也就说该直线偏离了垂直方向（z轴方向）0.81989°。
根据中心连线及其投影图，我们将它分为三段来进行计算：
1-5楼、5-9楼、9-13楼。
分别用直线代替了这三条线段，如图所示：
根据拟合直线的方程算得：
夹角名称
楼层段
水平投影和y轴的夹角（度）
水平投影和x轴的夹角（度）
所在直线和xy水平面的夹角（度）
偏离垂直线的角度（与z轴的夹角）
1-5层
-44.0669
134.0669
89.3986
0.6013
5-9层
-79.6511
169.6511
88.9719
1.0280
9-13层
-40.2861
130.2861
89.0282
0.9717
根据图表数据分析，中间段5-9层的角度与上下两层的角度相差较大，偏离垂直线的角度也是最大，倾斜严重，最有可能在此段发生坍塌。
5.3.2动态分析。
为了比较1996年和2009年的数据，从总体拟合直线上来比较两组数据，列表如下：
数据名
测量序数
水平投影和y轴的夹角（度）
水平投影和x轴的夹角（度）
所在直线和xy水平面的夹角（度）
偏离垂直线的角度（与z轴的夹角）
2009年（第3次）
-60.120
150.12
89.180
0.81989
1996年（第2次）
-54.439
144.43
89.276
0.72353
差值
-5.681
-5.690
-0.096
0.09636
可以发现总体上它们的倾斜的角度的差值并不明显。
将1996年和2009年的数据拟合图的起点重合后，两次拟合连线的对图比如下：
将2009年（第3次）分为三段拟合直线的夹角，分别和1996年（第2次）的拟合直线的夹角值进行比较。
古塔2009年测量所得角度和1996年测量所得角度分层比较表：
角度
量分段数
水平投影和y轴的夹角之差
水平投影和x轴的夹角之差
所在直线和xy水平面的夹角之差
偏离垂直线（z轴）的角度差
1-5层
10.372
-10.363
0.1226
-0.1222
5-9层
-25.212
25.221
-0.3041
0.3044
9-13层
14.153
-14.143
-0.2478
0.2481
由表格可以看出，5-9层在1996年至于2009年期间，发生了较大的形变，角度变化最大，致使整个塔身发生了扭曲。
六、对模型的评价
6.1模型优点
简化了塔的模型，确定了每层塔的中心（重心），以连线代替塔体。
可以充分展示塔身的倾斜、变形、弯曲和扭曲情况。
6.2模型缺点
实际的塔体情况可能很复杂，其切面不一定是中心对称图形。
只能通过平均值来尽量接近中心值。
6.3改进
如果能增加测量点的数据个数是最好的，在限制现有数据个数的情况下，可以进一步用数据筛选方法，淘汰偏差较大的数据。
七、参考文献
[1]李柏年.MATLAB数据分析方法[M].北京:
机械工业出版社,2012:
12-19.
[2]王能超.数值分析简明教程[M].北京：
高等教育出版社，2003：
36-37.
[3]于晓秋，李欣，张宏礼.MATLAB数值实验[M].沈阳：
辽宁科学技术出版社，2004：
34-40.
八、附录
附录一，主程序（保存为文件C2013.m）：
gcsj1=xlsread（'D:
\gcsj.xls','sheet1','C4:
E111'）;%从电子表格中读取对应测量序数的数据
gcsj2=xlsread（'D:
\gcsj.xls','sheet1','I4:
K111'）;
gcsj3=xlsread（'D:
\gcsj.xls','sheet1','O4:
Q111'）;
gcsj4=xlsread（'D:
\gcsj.xls','sheet1','U4:
W111'）;
symsYoXoZoyoxoXYyo2xo2
%画散点连线图
foru=1:
4
figure;
fork=0:
12
x=gcsj（（k*8+1）:
（8*k+8）,1）;%读取第k+1楼层的坐标数据值
y=gcsj（（k*8+1）:
（8*k+8）,2）;
z=gcsj（（k*8+1）:
（8*k+8）,3）;
x=x（~isnan（x））;%去除空白数据NaN
y=y（~isnan（y））;
z=z（~isnan（z））;
r1=size（x）;
r=r1
（1）;
plot3（x,y,z,'r-*'）;%画采样坐标点连线图
gridon
holdon
plot3（[x
（1）,x（r）],[y
（1）,y（r）],[z
（1）,z（r）],'r-*'）;%画第1点和最后1个的连线图
xlabel（'x轴'）;ylabel（'y轴'）;zlabel（'z轴'）;
end
end
%计算古塔每个楼层的中心坐标值
figure;
foru=1:
4
switchu%读取第u次测量数据
case1
gcsj=gcsj1;
case2
gcsj=gcsj2;
case3
gcsj=gcsj3;
case4
gcsj=gcsj4;
end
Yo=0;
Xo=0;
fork=0:
12
xq=gcsj（（k*8+1）:
（8*k+8）,1）;%读取第k+1楼层坐标数据
yq=gcsj（（k*8+1）:
（8*k+8）,2）;
zq=gcsj（（k*8+1）:
（8*k+8）,3）;
Yo=0;
Xo=0;
ll=0;
x=xq（~isnan（xq））;%去除空白数据NaN
y=yq（~isnan（yq））;
z=zq（~isnan（zq））;
r1=size（x）;
r=r1
（1）;
forj=1:
r-2
forn=j+1:
r-1
K1=（x（j）-x（j+1））/（y（j+1）-y（j））;%算中垂线斜率
K2=（x（j+2）-x（j+1））/（y（j+1）-y（j+2））;
equation1=sym（'Y-（y（j）+y（j+1））/2-K1*（X-（x（j+1）+x（j））/2）=0'）;%解中垂线交点方程
equation2=sym（'Y-（y（j+2）+y（j+1））/2-K2*（X-（x（j+1）+x（j+2））/2）=0'）;
[xo,yo]=solve（equation1,equation2,'X','Y'）;
yo2=eval（yo）;xo2=eval（xo）;
Yo=Yo+yo2;%累加中心点坐标值
Xo=Xo+xo2;
ll=ll+1;%累加中心点个数
end
K1=（x（r-1）-x（r））/（y（r）-y（r-1））;
K2=（x
（1）-x（r））/（y（r）-y
（1））;
equation1=sym（'Y-（y（r-1）+y（r））/2-K1*（X-（x（r）+x（r-1））/2）=0'）;
equation2=sym（'Y-（y
（1）+y（r））/2-K2*（X-（x
（1）+x（r））/2）=0'）;
[xo,yo]=solve（equation1,equation2,'X','Y'）;
yo2=eval（yo）;xo2=eval（xo）;
Yo=Yo+yo2;
Xo=Xo+xo2;
ll=ll+1;%累加中心点个数
end
Yok（u,k+1）=Yo/ll;%计算中心点坐标平均值
Xok（u,k+1）=Xo/ll;
Zok（u,k+1）=sum（z）/r;
end
s=size（gcsj）;%开始计算顶点坐标
xd=gcsj（105:
s
（1）,1）;
yd=gcsj（105:
s
（1）,2）;
zd=gcsj（105:
s
（1）,3）;
xd2=xd（~isnan（xd））;%去除空白数据NaN
yd2=yd（~isnan（yd））;
zd2=zd（~isnan（zd））;
p1=size（xd2）;
p=p1
（1）;
xt（u）=sum（xd2）/p;%算得顶点的坐标值（xt,yt,zt）
yt（u）=sum（yd2）/p;
zt（u）=sum（zd2）/p;
huatu（u,Xok,Yok,Zok,xt,yt,zt）;%画各楼层中心连线图
end
%1996年数据分析
figure;
title（'1996年测量数据楼层中心连线分析图'）;
plot3（Xok（2,:
）,Yok（2,:
）,Zok（2,:
）,'r-o'）%画中心坐标连线图
holdon;
gridon;
xlabel（'x轴'）;ylabel（'y轴'）;zlabel（'z轴'）;
plot3（[Xok（2,13）,xt
（2）],[Yok（2,13）,yt
（2）],[Zok（2,13）,zt
（2）],'r-pentagram'）;%最后和顶点连线
plot（Xok（2,:
）,Yok（2,:
）,'k-o'）%中心坐标连线在xy平面上的投影图
plot（[Xok（2,13）,xt
（2）],[Yok（2,13）,yt
（2）],'k-pentagram'）;%和顶点的连线在xy平面上的投影
plot3（566.5*ones（1,13）,Yok（2,:
）,Zok（2,:
）,'g-o'）%画连线在zy平面上的投影
plot3（[566.5,566.5],[Yok（2,13）,yt
（2）],[Zok（2,13）,zt
（2）],'g-pentagram'）;%最后和顶点连线
plot3（Xok（2,:
）,522.8*ones（1,13）,Zok（2
升级会员 