古塔的变形模型

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

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

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

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

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C
我们的参赛报名号为（如果赛区设置报名号的话）：sm063
所属学校（请填写完整的全名）：山东现代职业学院
参赛队员(打印并签名) ：1. 马昱轩
2. 李倩
3. 梁帅健
指导教师或指导教师组负责人(打印并签名)：宋祖芳
日期：2013年 9月 16 日赛区评阅编号（由赛区组委会评阅前进行编号）：sm06303
2013高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
评
阅
人
评
分
备
注
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
古塔的变形模型
摘要
某古塔是我国重点保护文物,已有上千年历史。

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

为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。

管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。

基于附件1提供的4次观测数据：对于问题1，1986年和1996年的观测数据中都缺少13层一个点的数据（因MATLA B程序中用的是循环语句，所以计算时赋予0值），其它各层均给出8个点的观测值，为使所得数据更具真实性，确定古塔各层中心位置的方法更具适用性，本文将每层所给8点构成的图形看做不规则八边形，用中垂线求交点法求得古塔各层中心坐标。

对于问题2，结合问题1的分析，采用垂直投影法[1]求古塔的倾斜度，根据所得数据，分析古塔的倾斜程度（因1986年和1996年13层赋予值后所得数据偏差较大，为使所得数据更具真实性，所以本问题起1986年和1996年13层数据予以舍弃）；弯曲是建立在二维平面上的一条曲线，通过截取古塔过x轴、z轴的界面，求出每层古塔的倾斜度，从而分析得出古塔塔身在x轴、z轴的界面的弯曲程度。

同理，也可分析y轴、z轴的界面的弯曲程度；扭曲同样是采用垂直投影法[1]求古塔每层的倾斜度，建立三维立体空间，根据所得数据，分析古塔塔身的扭曲程度。

对于问题3，利用问题1、2所得数据，进行合理的分析与猜想，进而分析出古塔塔身的倾斜、弯曲、扭曲等变化趋势。

本文的模型解决了题目给出的问题，计算过程中充分尊重观测数据，给出更符合实际的结果。

本文所得结果大部分由图表给出，结合图像，较为直观地表现出古塔变形情况。

结果表明，采用MATLAB数学软件可得出可靠结论。

关键词：古塔变形监测MATLAB垂直投影法倾斜度
一、问题的重述
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。

为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。

某古塔已有上千年历史，是我国重点保护文物。

管理部门委托测绘公司先后于1 986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。

请你们根据附件1提供的4次观测数据，讨论以下问题：
1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。

2. 分析该塔倾斜、弯曲、扭曲等变形情况。

3. 分析该塔的变形趋势。

二、问题的分析
2.1问题的背景
某古塔是我国重点保护文物。

已有上千年历史，由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生诸如倾斜、弯曲、扭曲等各种变形。

对现存历史文物古塔的保护，掌握古塔的形变，显得尤为重要，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。

管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。

2.2问题的分析
基于附件1提供的4次观测数据：
2.2.1问题一的分析
本题建立模型是利用中垂线（即中垂线相交于一点为中点)法，求解古塔各层中心的坐标，以确保古塔中心坐标的准确性。

2.2.2问题二的分析
根据层与层之间的中心点，塔体倾斜变化采用垂直投影法[1](即垂线法，利用相对位移量进行监测)；弯曲是建立在二维平面上的一条曲线，通过截取古塔过x轴、
z 轴的界面，求出每层古塔的倾斜度，从而分析得出古塔塔身在x 轴、z 轴的界面的
弯曲程度。

同理，也可分析y 轴、z 轴的界面的弯曲程度；扭曲同样是采用垂直投影法求古塔每层的倾斜度，建立三维立体空间，根据所得数据，分析古塔塔身的扭曲程度。

2.2.3问题三的分析
我们在对古塔进行安全性监测方面要满足两个方面：①塔体的倾斜变化；②塔身各段的变化[1]。

针对上述几个方面就采取对应的监测方法。

三、问题的假设
1.假设求解古塔各层中心坐标时纵坐标先不予以考虑；
2.假设古塔的地基土抗压能力强、韧性大、不易收缩，对古塔的变形不产生任何影响；
3.假设对古塔观测数据的缺失可忽略不计；
四、符号及文字说明
x 表示建立古塔坐标轴的横坐标
i x 表示古塔各层的第i (i =1,2,...,8）个点的横坐标
y 表示建立古塔坐标轴的纵坐标
i y 表示古塔各层的第i （i =1,2，...,8)个点的纵坐标
z 表示建立古塔坐标轴竖坐标
i z 表示古塔各层的第i (i =1,2,...,8)个点的竖坐标 i k 表示古塔的第i l 条直线的斜率
'x 表示两个x 点之间的中点
'y 表示两个y 点之间的中点 H 表示古塔中心点到水平面的距离 D 表示古塔中心点的位移量
Q表示古塔的倾斜度
j表示古塔的第j(j=1,2,...,13)层数（MATLAB程序中）
i表示古塔各层的第i(i=1,2,...,8）个点（MATLAB程序中）
X表示关于横坐标的8
13⨯矩阵（MATLAB程序中）
Y表示关于纵坐标的8
13⨯矩阵（MATLAB程序中）
K表示各层每两点之间中垂线的斜率（MATLAB程序中）
xx表示每层8条中垂线两两相交所得交点的关于横坐标的8
13⨯矩阵（MATLAB程序中）
yy表示每层8条中垂线两两相交所得交点的关于纵坐标的8
13⨯矩阵（MATLAB程序中）
五、模型的建立与求解
（一）问题一的模型建立与求解
1.通过中垂线的方法对古塔各层中心坐标进行求解，根据附件1表格数据，1986年和1996年的观测数据中都缺少13层一个点的数据计算时赋予0值。

然后观察数据，古塔的每i层八个点的高度变化不大，假设把古塔的每一层的点看作是在一个平面上，求解古塔各层中心坐标时纵坐标先不予以考虑，之后用平均值的方法计算出。

分析过程如下图1-1：（用MATLAB软件作图（附件1，程序1））
1-1 根据上图，做出任意两边的中垂线，交与一点O，如下图所示1-2（在上图上用画图工具作图）：
1-2
由图可知，'A 点坐标为()','y x ，'B 点坐标)","(y x 。

根据公式2'21x x x +=，得:
2'21x x x +=
，2'2
1y y y +=；2
"32x x x +=，2"32y y y += ① 根据直线的斜率公式，求出直线AB l 的斜率A k 。

同理，求出直线BC l 的斜率B k 。

121
2x x y y k A --=
2
323x x y y k B --= 再根据中垂线定理1'-=kk ，求出中垂线1l ，2l 的斜率1k ，2k ； 122111y y x x k k A --=-=
2
33
221y y x x k k B --=
-= ② 设中垂线的直线方程1l 为：
)'('1x x k y y -=- ③ 设中垂线的直线方程2l 为：
)"("2x x k y y -=- ④ 通过方程③、④建立方程组:
⎩
⎨⎧-=--=-)"(")
'('21x x k y y x x k y y
解得： 2
121"
''"k k x k x k y y x -+--=
； ''21y x k x k y +-= ⑤
再把①、②式中所求的数据代入⑤中，即可解出两中垂线交点O 坐标。

同理，通过数学软件MATLAB 编程求解（附录1，程序2）可得出每层八条中垂线两两相交的八个交点坐标，得出数据并对这八个中垂线的交点求平均值，因1986年和1996年的观测数据中都缺少13层一个点的数据计算时赋予0值得出的数据偏差太大，此时予以舍弃。

塔顶数据是塔顶的四个点求平均值而得。

各层中心坐标如下表：
表一： 古塔层次
各层中心坐标),(i i i y x O
i x
i y 1 566.664514 522.7091368 2 566.7215137 522.6721729 3 566.7773312 522.6360299 4 566.821135 522.607225 5 566.8683524 522.5760928 6 566.9173327 522.5455899 7 566.9505701 522.5283882 8 566.9835974 522.5115552 9 567.0181737 522.4949434 10 567.0509921 522.4813031 11 567.1072922 522.4414325 12 567.1630598
522.4013267
13 --- --- 塔顶 567.24725
522.24375
而每层古塔都有一定的高度差，因此，假设每层古塔的中心高度为每层八点的平均高度，即：),,(i i i i z y x O ，因此，古塔每层的中心坐标是： 表二：1986年 古塔层次
古塔各层中心坐标()i i i i z y x O ,,
i x
i y i z 1 566.664514 522.7091368 1.787375 2 566.7215137 522.6721729 7.32025 3 566.7773312 522.6360299 12.75525 4 566.821135 522.607225 17.07825 5 566.8683524 522.5760928 21.7205 6 566.9173327 522.5455899 26.235125 7 566.9505701 522.5283882 29.836875 8 566.9835974 522.5115552 33.350875 9 567.0181737 522.4949434 36.854875 10 567.0509921 522.4813031 40.172125 11 567.1072922 522.4414325 44.440875 12
567.1630598
522.4013267
48.711857
13 --- --- --- 塔顶 567.24725 522.24375 55.123525
表三：1996年 古塔层次
古塔各层中心坐标()i i i i z y x O ,,
i x
i y i z 1 566.664714 522.7088368 1.783 2 566.7224008 522.6712321 7.314625 3 566.7789381 522.634408 12.75075 4 566.8233162 522.6050197 17.075125 5 566.8710945 522.5732208 21.76 6 566.9207197 522.5421644 26.2295 7 566.9543744 522.5245028 29.836875 8 566.9877018 522.5073211 33.350875 9 567.0229012 522.4899727 36.854875 10 567.0561171 522.47609 40.167625 11 567.1130185 522.4356821 44.440875 12 567.1690607
522.395182
48.707375 13 --- --- --- 塔顶 567.25435
522.23665
55.11975
表四：2009年 古塔层次
古塔各层中心坐标()i i i i z y x O ,,
i x
i y i z 1 566.7448287 522.7001717 1.7645 2 566.7794913 522.6717696 7.309 3 566.8123971 522.6450678 12.73225 4 566.839279 522.6225762 17.06975 5 566.8672333 522.5991397 21.709375 6 566.9573958 522.5522863 26.211 7 566.991 522.5309222 29.824625 8 567.0440874 522.5004085 33.339875 9 567.0968128 522.469119 36.84375 10 567.1524056 522.4153501 40.161125 11 567.1963122 522.3745082 44.432625 12
567.2399861
522.3339371
48.69975
13 567.2901509 522.2874247 52.818375 塔顶 567.336 522.2148 55.091
表五：2011年 古塔层次
古塔各层中心坐标()i i i i z y x O ,,
i x
i y i z 1 566.7449676 522.7000432 1.76325 2 566.7796914 522.6715196 7.2905 3 566.81276 522.6453927 12.726975 4 566.8396637 522.6220928 17.025 5 566.8678429 522.5985264 21.703875 6 566.9580489 522.5514853 26.2045 7 566.9918 522.5301222 29.817 8 567.044881 522.4995201 33.336625 9 567.0977815 522.4681716 36.82225 10 567.1534953 522.4142346 40.144125 11 567.1974118 522.3733937 44.424875 12 567.2412861 522.3326371 48.683875 13 567.2914636 522.2859839 52.813125 塔顶 567.3375
522.2135
55.087
（二）问题二的模型建立与求解
1，根据问题一求出的每层的中心坐标数据，通过分析古塔的倾斜、弯曲、扭曲及层与层之间的中心点，确定了塔体倾斜变化采用垂直投影法[1](即垂线法，利用相对位移量进行监测);弯曲、扭曲即为塔身各段的变化采用求倾斜度的方法。

用MATLAB 数学软件编程（附件1，程序3）和画图工具画图，如下图2-1：
2-1 通过对上图的观察与分析，每层的中心大致在一条倾斜的直线上，采用垂直投影的方法，对古塔的倾斜角进行求解，如下图通过MATLAB数学软件编程和画图工具
并用作图（附件1，程序4），建立的数学模型，如下图2-2；
2-2
通过对图形的观察，()o o o z y x A ,,点为古塔塔顶的中心坐标，()111,,z y x O 点为古塔第一层的中心坐标，()1,,z y x B o o 点是A 点在O 点所在水平面上的投影， 根据倾斜度公式D
H Q =
,[3]
计算出Q 的值： 2
11212
1)()(z z y y x x OB D o o -+-+-==）（
1z z BA H o -== 2
1121211
)()()(z z y y x x z z Q o o o -+-+--=
根据以上公式，得出以下数据，如表六：
倾斜度 1986 1996
2009 2011 Q 71.5191818 70.6072941
69.7171388 69.5506006
比较数据：
2011200919961986Q Q Q Q >>>
分析数据，Q 的值越小，倾斜程度越严重，古塔随着年代的增加而越来越倾斜。

2.对于弯曲的情况，弯曲是建立在二维平面上的一条曲线，在解决问题的过程中，为了对弯曲便以分析，用x 轴、z 轴所在的平面建立二维平面的数学模型，（由于古塔的弯曲方向具有不确定性，此方法具有一定的局限性）通过画图工具画出简易二维图形2-3；
2-3
根据倾斜度公式：D
H Q =
[3]
，求出每层的倾斜度，
每层塔的的高度： 12z z H -= 中心点的位移距离： 12x x D -= 1
21
21x x z z D H Q --=
=
根据上述公式计算出1986、1996、2009、2011年的每层古塔的倾斜度，由于1986、1996年的第十三层丢失数据，因此舍去第十三层的斜率，用MATLAB 数学软件编程（附件1，程序5），计算结果，如下表： 层次 年代 1986 1996 2009 2011 倾斜度
Q 1 97.06849334 95.89065436 159.9562641 159.177567 2 97.37089622 96.15112501 164.8113706 164.3999141 3 98.69006799 97.44389688 161.3539222 159.7559072 4 98.31651044 98.0544515 165.9717825 166.0400224 5 92.17226109 90.0651282 49.92790794 49.89274549 6 108.3643727 107.1878519 107.5349212 107.0335485 7
106.3968293
105.4387681
66.21627731
66.30668224
8 101.341092 99.54715137 66.45516203 65.89020898
9 101.0789679 99.7338624 59.67274539 59.62391723
10 75.82135733 75.0992067 97.28605722 97.47475322
11 76.58536498 76.13013051 97.70423525 97.07277381
12 76.15693988 75.18381555 82.10189216 82.29286035
13 --- --- 49.56749424 49.39298034
根据表格中所得的数据，可以分析出，每层古塔的倾斜度不同，塔身的弯曲程度也不同，倾斜度越小，塔身的弯曲程度越严重，通过这四个年代每层对应的倾斜度的比较，得出古塔的塔身的弯曲程度每年都在加深。

3.扭曲等情况，需要按着上述方法计算出每层古塔的倾斜角Q，通过比较每层的倾斜角Q，可得出古塔的扭曲情况。

下图2-4为古塔每年的中心坐标图（通过MATLAB 数学软件作图（附件1，程序6），旋转选取适当的视角）：
2-4
通过对上图的观察，每层都有一定的扭曲程度，采用第1小题中的方法求出每层古塔的倾斜度Q ，如图（用画图工具作简易图）：
根据倾斜度公式：D
H
Q =
，求出每层的倾斜度， 每层塔的的高度： 12z z H -=
中心点的位移距离： 211212212)()()(z z y y x x D -+-+-= 2112122121
21)
()()(z z y y x x z z D H Q -+-+--=
=
根据上述公式计算出1986、1996、2009、2011年的每层古塔的倾斜度，由于1986、1996年的第十三层丢失数据，因此舍去第十三层的斜率，用MATLAB 数学软件编程（附件1，程序7），计算结果，如下表： 层次 年代
1986
1996 2009 2011 倾斜度Q 1 81.44249931 80.32993855 123.7261777 122.9998537 2 81.73246759 80.56851909 127.9774128 128.9964968 3 82.45893994 81.24437845 123.751475 120.7627543 4 82.08087757 81.62823947 127.1864292 127.3691491 5 78.24063887 76.34695537 44.30313221 44.23875692 6 96.23939862 94.91215913 90.74800462 90.43923785 7 94.79475346 93.71745002 57.40870708 57.44398623 8 91.34568441 89.29121007 57.14949617 56.68476282 9 93.33799073 92.01993259 42.8927531 42.83802791 10 61.87656109 61.23059913 71.2327331 71.37917038 11 62.17642433 61.70395494 71.58332514 71.12121848 12 35.88809682 35.62064489
60.20512701 60.26784157 13 --- ---
26.46080111
26.48113768
观察图表中数据，
（三）问题三的模型建立与求解
古塔存在时间久远，经过长时间大自然的侵袭和破坏，并且还受到人们生产生活的影响。

塔的内部和外部都受到一定的影响。

根据建筑力学得知，当发生地震时，地基结构发生了改变使得古塔发生不均匀沉降，导致古塔发生了变形。

经过四次观测数据得出古塔每次都发生了很微小的变形，如建筑物超过自重，会使得抗外界干扰能力及承载力降低从而造成破坏。

就风压而言，相当于是对古塔产生额外的与重力相同的外力作用，会增加内部古塔柱梁板等构件的内力，可能造成构件内力超过其承载能力，从而造成古塔变形。

就风吸
而言，建筑上一般不考虑其对古塔减轻重力的有利一面，而且需要古塔在风的吸力作用下是否会造成风吸力大于建筑物屋面重力的情况，而造成屋面被掀起。

这样的例子，主要表现在建筑物外围雨棚的设计中，对于一般屋面，其外形一般产生的只是风压。

因此会出现古塔的倾、弯曲、扭曲等各种变形。

根据以上问题得出的数据，分析得出的古塔的倾斜、弯曲、扭曲等变化情况，如果不受外力影响（除重力外），塔身会随着时间的变化会继续倾斜、弯曲、扭曲等变化，如图2-4，根据四次测量的数据及分析，猜想古塔会继续向某些方向倾斜、弯曲、扭曲等变化；如果考虑到外力对塔身的影响，塔身会发生不规则倾斜、弯曲、扭曲等变化。

如果古塔到了一定的年代，塔身的变形程度过大，古塔重心不足以维持塔身的平衡，古塔将会出现倒塌一部分，或者完全倒塌。

六、模型的评价与改进方向
6.1模型的优缺点
（1）本模型综合利用了MATLAB、中垂线法、垂直投影法等数学方法，数学推导严谨，理论性强，所得数据真实可靠。

（2）模型具有一般性，在一定条件下适用于其他古建筑物的变形情况监测。

（3）由于数据较多，计算量过大，对于问题一，算得八条垂直平分线的八个交点后，采用极限思想的计算会使数据更准确一些。

（4）问题二中的弯曲，扭曲在已得的数据上用拟合方法会更优化一些。

6.2模型的改进放向
我们的模型通过一些合理的假设，使得问题的数学描述较为简单直观，但实际上还有很多因素与要讨论的问题密切相关，而且是应该加以考虑的。

对于问题一的求解方法，如果多次计算每次求出的数据，根据极限思想，也就是利用求极限方法求各层中心坐标，最终的结果更无限接近于真实值。

如问题二，在解决弯曲变化时，应多方位考虑弯曲方向等问题，
这些都是在研究最优策略时应该认真考虑的问题，因为时间关系，我们只讨论了比基础、较简单的建模方法，还有待于完善和改进。

七、参考文献
[1]梁海奎.古塔变形测量方法探讨.《城市勘探》2011年第03期
[2]MATLAB软件编程；清华大学出版社
[3]韩煊,李宁.地铁施工引起的建筑物扭曲变形分析.《土木工程学报》2010年1月
八、附录
附件1
程序1：
x=[565.454
562.058
561.39
563.782
567.941
571.255
571.938
569.5
565.454];
y=[528.012
525.544
521.447
518.108
517.407
519.857
523.953
527.356
528.012];
plot(x,y)
程序2：
（1986年）
xa=load('xa.txt');
ya=load('ya.txt');
X=xa';
Y=ya';
K=zeros(13,8);
M=zeros(13,8);
N=zeros(13,8);
x=[X,X(:,1)];
y=[Y,Y(:,1)];
for j=1:13
for i=1:8
M(j,i)=(x(j,i)+x(j,i+1))/2;
N(j,i)=(y(j,i)+y(j,i+1))/2;
K(j,i)=-(x(j,i+1)-x(j,i))./(y(j,i+1)-y(j,i));
end
end
M
N
K
xx=zeros(13,8);
yy=zeros(13,8);
N1=[N,N(:,1)];
M1=[M,M(:,1)];
K1=[K,K(:,1)];
for j=1:13
for i=1:8
xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1))/(K1(j,i)-K1(j,i+1));
yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i);
end
end
xx
yy
xlswrite('data1.xls',xx,1)
xlswrite('data2.xls',yy,1)
（1996年）
xa=load('xb.txt');
ya=load('yb.txt');
X=xa';
Y=ya';
K=zeros(13,8);
M=zeros(13,8);
N=zeros(13,8);
x=[X,X(:,1)];
y=[Y,Y(:,1)];
for j=1:13
for i=1:8
M(j,i)=(x(j,i)+x(j,i+1))/2;
N(j,i)=(y(j,i)+y(j,i+1))/2;
K(j,i)=-(x(j,i+1)-x(j,i))./(y(j,i+1)-y(j,i));
end
end
M
N
K
xx=zeros(13,8);
yy=zeros(13,8);
N1=[N,N(:,1)];
M1=[M,M(:,1)];
K1=[K,K(:,1)];
for j=1:13
for i=1:8
xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1))/(K1(j,i)-K1(j,i+1));
yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i);
end
end
xx
yy
xlswrite('data3.xls',xx,1)
xlswrite('data4.xls',yy,1)
（2009年）
xa=load('xc.txt');
ya=load('yc.txt');
X=xa';
Y=ya';
K=zeros(13,8);
M=zeros(13,8);
N=zeros(13,8);
x=[X,X(:,1)];
y=[Y,Y(:,1)];
for j=1:13
for i=1:8
M(j,i)=(x(j,i)+x(j,i+1))/2;
N(j,i)=(y(j,i)+y(j,i+1))/2;
K(j,i)=-(x(j,i+1)-x(j,i))./(y(j,i+1)-y(j,i));
end
end
M
N
K
xx=zeros(13,8);
yy=zeros(13,8);
N1=[N,N(:,1)];
M1=[M,M(:,1)];
K1=[K,K(:,1)];
for j=1:13
for i=1:8
xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1))/(K1(j,i)-K1(j,i+1));
yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i);
end
end
xx
yy
xlswrite('data5.xls',xx,1)
xlswrite('data6.xls',yy,1)
（2011年）
xa=load('xd.txt');
ya=load('yd.txt');
X=xa';
Y=ya';
K=zeros(13,8);
M=zeros(13,8);
N=zeros(13,8);
x=[X,X(:,1)];
y=[Y,Y(:,1)];
for j=1:13
for i=1:8
M(j,i)=(x(j,i)+x(j,i+1))/2;
N(j,i)=(y(j,i)+y(j,i+1))/2;
K(j,i)=-(x(j,i+1)-x(j,i))./(y(j,i+1)-y(j,i));
end
end
M
N
K
xx=zeros(13,8);
yy=zeros(13,8);
N1=[N,N(:,1)];
M1=[M,M(:,1)];
K1=[K,K(:,1)];
for j=1:13
for i=1:8
xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1))/(K1(j,i)-K1(j,i+1));
yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i);
end
end
xx
yy
xlswrite('data7.xls',xx,1)
xlswrite('data8.xls',yy,1)
程序3：
x=[566.664514
566.7215137
566.7773312
566.821135
566.8683524
566.9173327
566.9505701
566.9835974
567.0181737
567.0509921
567.1072922
567.1630598
567.24725];
y=[522.7091368
522.6721729
522.6360299
522.607225
522.5760928
522.5455899
522.5283882
522.5115552
522.4949434
522.4813031 522.4414325 522.4013267 522.24375]; z=[1.787375 7.32025 12.75525 17.07825 21.7205 26.235125 29.836875 33.350875 36.854875 40.172125 44.440875 48.711857 55.123525]; plot3(x,y,z,'r.')
程序4：（1）
x=[565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5
565.48 562.238 561.663 564.001 567.995 571.165 571.801 569.414 565.506 562.415 561.931 564.216
571.076 571.666 569.33 565.526 562.555 562.144 564.387 568.091 571.005 571.558 569.263 565.548 562.706 562.373 564.571 568.136 570.929 571.443 569.191 565.57 562.854 562.6 564.752 568.18 570.857 571.333 569.121 565.671 563.132 562.883 564.949 568.172 570.679 571.094 568.994 565.77 563.403 563.158 565.141 568.164
570.862 568.87 565.868 563.674 563.433 565.333 568.156 570.333 570.63 568.747 565.961 563.927 563.693 565.516 568.148 570.171 570.408 568.631 566.078 564.193 563.958 565.649 568.094 570.013 570.236 568.615 566.195 564.459 564.224 565.782 568.039 569.854 570.063 568.598 566.308 564.716 564.481 565.91 569.701 569.897
567.255 567.235 567.247 567.252]; y=[528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356 527.764 525.364 521.42 518.226 517.563 519.961 523.908 527.141 527.52 525.188 521.394 518.343 517.716 520.063 523.864 526.93 527.327 525.047 521.373 518.435 517.838 520.144 523.829 526.762 527.119 524.896 521.351 518.534
520.232 523.791 526.581 526.915 524.748 521.329 518.632 518.095 520.315 523.755 526.406 526.652 524.585 521.356 518.846 518.346 520.441 523.672 526.167 526.397 524.427 521.382 519.055 518.59 520.564 523.591 525.933 526.141 524.268 521.408 519.263 518.834 520.686 523.51 525.701 525.9 524.12 521.433 519.462 519.068
523.433 525.482 525.628 523.95 521.463 519.607 519.242 520.885 523.35 525.259 525.355 523.78 521.492 519.753 519.415 520.969 523.268 525.037 525.092 523.616 521.521 519.893 521.05 523.188 524.822 522.238 522.242 522.251 522.244]; z=[1.792 1.818 1.783 1.769 1.772 1.77
1.794 1.801 7.326 7.351 7.314
7.306 7.304 7.324 7.336 12.761 12.786 12.749 12.736 12.741 12.74 12.758 12.771 17.084 17.109 17.072 17.059 17.064 17.063 17.081 17.094 21.726 21.751 21.714 21.701 21.705 21.708 21.723 21.736 26.267 26.309 26.308 26.264 26.189 26.136 26.164 26.244 29.869 29.911 29.91 29.866
29.737 29.765 29.846 33.383 33.425 33.424 33.38 33.305 33.251 33.279 33.36 36.887 36.929 36.928 36.884 36.809 36.755 36.783 36.864 40.201 40.214 40.244 40.223 40.171 40.038 40.129 40.157 44.472 44.485 44.505 44.486 44.442 44.309 44.4 44.428 48.743 48.756 48.776 48.757 48.713
48.671
48.699
52.866
52.878
52.897
52.88
52.703
52.794
52.822
55.128
55.108
55.128
55.129];
plot3(x,y,z,'b.')
（2）
z=-1:1;
x=[565.454 562.058 561.39
563.782 567.941 571.255 571.938 569.5
565.454];
y=[528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356 528.012];
z=[0
0];
plot3(x,y,z,'r')
程序5：
（1986年）
i=1:13;
xa=[566.664514
566.7215137
566.7773312
566.821135
566.8683524
566.9173327
566.9505701
566.9835974
567.0181737
567.0509921
567.1072922
567.1630598
567.24725];
za=[1.787375
7.32025
12.75525
17.07825
21.7205
26.235125
29.836875
33.350875
36.854875
40.172125
44.440875
48.711857
55.123525];
x=xa';
z=za';
Q=zeros(1,13);
for i=1:12
Q(i)=(z(i+1)-z(i))/(x(i+1)-x(i));
Q=Q'
xlswrite('bb.xls',Q,1)
（1996年）
i=1:13;
xa=[566.664714
566.7224008
566.7789381
566.8233162
566.8710945
566.9207197
566.9543744
566.9877018
567.0229012
567.0561171
567.1130185
567.1690607
567.25435];
za=[1.783
7.314625
12.75075
17.075125
21.76
26.2295
29.836875
33.350875
36.854875
40.167625
44.440875
48.707375
55.11975];
x=xa';
z=za';
Q=zeros(1,13);
for i=1:12
Q(i)=(z(i+1)-z(i))/(x(i+1)-x(i)); end
Q=Q'
xlswrite('bb.xls',Q,2)
（2009年）
i=1:13;
xa=[566.7448287
566.7794913
566.8123971
566.839279
566.8672333
566.9573958
566.991
567.0440874
567.0968128
567.1524056
567.1963122
567.2399861
567.2901509
567.336];
za=[1.7645
7.309
12.73225
17.06975
21.709375
26.211
29.824625
33.339875
36.84375
40.161125
44.432625
48.69975
52.818375
55.091];
x=xa';
z=za';
Q=zeros(1,13);
for i=1:13
Q(i)=(z(i+1)-z(i))/(x(i+1)-x(i)); end
Q=Q'
xlswrite('bb.xls',Q,3)
（2011年）
i=1:13;
xa=[566.7449676
566.7796914
566.81276
566.8396637
566.8678429
566.9580489
566.9918
567.044881
567.0977815
567.1534953
567.1974118
567.2412861
567.2914636
567.3375];
za=[1.76325
7.2905
12.726975
17.025
21.703875
26.2045
29.817
33.336625
36.82225
40.144125
44.424875
48.683875
52.813125
55.087];
x=xa';
z=za';
Q=zeros(1,13);
for i=1:13
Q(i)=(z(i+1)-z(i))/(x(i+1)-x(i)); end
Q=Q'
xlswrite('bb.xls',Q,4)
程序6：
（1986年）
x=[566.664514
566.7773312 566.821135 566.8683524 566.9173327 566.9505701 566.9835974 567.0181737 567.0509921 567.1072922 567.1630598 567.24725];
y=[522.7091368 522.6721729 522.6360299 522.607225 522.5760928 522.5455899 522.5283882 522.5115552 522.4949434 522.4813031 522.4414325 522.4013267 522.24375];
z=[1.787375 7.32025
12.75525
17.07825
21.7205
26.235125 29.836875 33.350875 36.854875 40.172125 44.440875 48.711857 55.123525]; plot3(x,y,z)
（1996年）
566.7224008 566.7789381 566.8233162 566.8710945 566.9207197 566.9543744 566.9877018 567.0229012 567.0561171 567.1130185 567.1690607 567.25435];
y=[522.7088368 522.6712321 522.634408 522.6050197 522.5732208 522.5421644 522.5245028 522.5073211 522.4899727 522.47609 522.4356821 522.395182 522.23665];
z=[1.783
7.314625
12.75075
17.075125 21.76
26.2295
29.836875 33.350875 36.854875 40.167625 44.440875 48.707375 55.11975];
plot3(x,y,z)
x=[566.7448287 566.7794913 566.8123971 566.839279 566.8672333 566.9573958 566.991
567.0440874 567.0968128 567.1524056 567.1963122 567.2399861 567.2901509 567.336];
y=[522.7001717 522.6717696 522.6450678 522.6225762 522.5991397 522.5522863 522.5309222 522.5004085 522.469119 522.4153501 522.3745082 522.3339371 522.2874247 522.2148];
z=[1.7645
7.309
12.73225
17.06975
21.709375 26.211
29.824625 33.339875 36.84375
40.161125 44.432625 48.69975
55.091];
plot3(x,y,z)
（2011年）
x=[566.7449676 566.7796914 566.81276 566.8396637 566.8678429 566.9580489 566.9918 567.044881 567.0977815 567.1534953 567.1974118 567.2412861 567.2914636 567.3375];
y=[522.7000432 522.6715196 522.6453927 522.6220928 522.5985264 522.5514853 522.5301222 522.4995201 522.4681716 522.4142346 522.3733937 522.3326371 522.2859839 522.2135];
z=[1.76325
7.2905
12.726975 17.025
21.703875 26.2045
29.817
33.336625
40.144125
44.424875
48.683875
52.813125
55.087];
plot3(x,y,z)
程序7：（1986年）
i=1:12;
xa=[566.664514 566.7215137 566.7773312 566.821135 566.8683524 566.9173327 566.9505701 566.9835974 567.0181737 567.0509921 567.1072922 567.1630598 567.24725];
ya=[522.7091368 522.6721729 522.6360299 522.607225 522.5760928 522.5455899 522.5283882 522.5115552 522.4949434 522.4813031 522.4414325 522.4013267 522.24375];
za=[1.787375 7.32025
12.75525
21.7205
26.235125
29.836875
33.350875
36.854875
40.172125
44.440875
48.711857
55.123525];
x=xa';
y=ya';
z=za';
Q=zeros(1,13);
for i=1:12
Q(i)=(z(i+1)-z(i))/(sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2+(z(i)-z(i))^2)); end
Q=Q'
xlswrite('aa.xls',Q,1)
（1996年）
i=1:12;
xa=[566.664714
566.7224008
566.7789381
566.8233162
566.8710945
566.9207197
566.9543744
566.9877018
567.0229012
567.0561171
567.1130185
567.1690607
567.25435];
ya=[522.7088368
522.6712321
522.634408
522.6050197
522.5732208
522.5421644
522.5073211
522.4899727
522.47609
522.4356821
522.395182
522.23665];
za=[1.783
7.314625
12.75075
17.075125
21.76
26.2295
29.836875
33.350875
36.854875
40.167625
44.440875
48.707375
55.11975];
x=xa';
y=ya';
z=za';
Q=zeros(1,13);
for i=1:12
Q(i)=(z(i+1)-z(i))/(sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2+(z(i)-z(i))^2)); end
Q=Q'
xlswrite('aa.xls',Q,2)
（2009年）
i=1:13;
xa=[566.7448287
566.7794913
566.8123971
566.839279
566.8672333
566.9573958
566.991
567.0440874
567.0968128
567.1963122
567.2399861
567.2901509
567.336];
ya=[522.7001717
522.6717696
522.6450678
522.6225762
522.5991397
522.5522863
522.5309222
522.5004085
522.469119
522.4153501
522.3745082
522.3339371
522.2874247
522.2148];
za=[1.7645
7.309
12.73225
17.06975
21.709375
26.211
29.824625
33.339875
36.84375
40.161125
44.432625
48.69975
52.818375
55.091];
x=xa';
y=ya';
z=za';
Q=zeros(1,13);
for i=1:13
Q(i)=(z(i+1)-z(i))/(sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2+(z(i)-z(i))^2)); end
Q=Q'
xlswrite('aa.xls',Q,3)
（2011年）
i=1:13;
xa=[566.7449676 566.7796914 566.81276
566.8396637 566.8678429 566.9580489 566.9918
567.044881
567.0977815 567.1534953 567.1974118 567.2412861 567.2914636 567.3375];
ya=[522.7000432 522.6715196 522.6453927 522.6220928 522.5985264 522.5514853 522.5301222 522.4995201 522.4681716 522.4142346 522.3733937 522.3326371 522.2859839 522.2135];
za=[1.76325
7.2905
12.726975
17.025
21.703875
26.2045
29.817
33.336625
36.82225
40.144125
44.424875
48.683875
52.813125
55.087];
x=xa';
y=ya';
z=za';
Q=zeros(1,13);
for i=1:13
Q(i)=(z(i+1)-z(i))/(sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2+(z(i)-z(i))^2)); end
Q=Q'
xlswrite('aa.xls',Q,4)。