2013数学建模——古塔的变形

承诺书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

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

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

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

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C
我们的参赛报名号为（如果赛区设置报名号的话）：5339
所属学校（请填写完整的全名）：
参赛队员(打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人(打印并签名)：
（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。

以上内容请仔细核对，提交后将不再允许做任何修改。

如填写错误，论文可能被取消评奖资格。

）
日期： 2013 年 09 月 16 日赛区评阅编号（由赛区组委会评阅前进行编号）：
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：
古塔的变形数学模型
摘要:
本文是研究关于古塔变形类型以及变形分析的模型，用Matlab画出古塔的三维结构可以看出它是近似于正八边形的形状。

因此，问题一我们用每层各个测量点坐标的平均值作为塔每层的中心坐标，再用中心坐标的三个坐标值分别对时间t做回归来得到确定古塔各层中心位置的通用方法。

对于问题二，我们分别研究该塔倾斜、弯曲、扭曲等变形情况，通过建立数学模型来确定变形的程度。

首先，用各层中心点坐标的z坐标值与其相应点的x,y坐标值做多元线性回归。

然后得到的回归方程所表示的回归平面与z轴正方向的夹角就可以表示古塔的倾斜程度大小。

最后根据各层中心的分布和变化趋势方向，确定古塔的倾斜方向。

用古塔各层中心点进行平面拟合，从效果上观察，较为精确地反映了实例中的问题，由此也说明了我们所建模型的合理性。

古塔的倾斜变形必然会导致在同一层中，测点存在高程的绝对差h，如果古塔只存在倾斜变形的话，每层的h值会相等；如果古塔存在倾斜变形的同时也存在弯曲变形的话，则每层的h值会发生改变。

所以相邻两层的高程绝对差的变化量，表示古塔每层弯曲程度大小。

再根据每层出现高程绝对差h的两个测量点的连线，确定每层弯曲方向。

古塔的扭曲变形，首先每层选取两对相同的对测量点，并做连线。

然后通过每层对测量点的连线，分别与第一层相同对测量点的连线所成的角度的平均值来衡量古塔的扭曲情况。

对于该塔的变形趋势的研究，将倾斜指标、弯曲指标、扭曲指标对时间的回归。

再用得到的回归方程预测未来几年的数据，结合用excel画出的图来预测古塔在未来时间里的变形趋势。

关键字：线性回归变化趋势拟合预测
一、问题重述
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古
塔会产生各种变形，诸如倾斜、弯曲、扭曲等。

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

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

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

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

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

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

二、问题分析
（一）、对问题一的分析
问题一中确定古塔各层中心位置的通用方法。

因为古塔各层为近似正八边形，根据正八边形图形特征，可以用每次测量时，古塔各层测量点坐标的平均值作为各层中心点坐标。

然后将各层中心点坐标对时间回归，可得到各层中心点坐标对时间的回归方程。

根据方程就可以确任意时间各层中心点坐标。

（二）、对问题二的分析
问题二要求我们确定塔的倾斜、弯曲、扭曲等变形情况。

为了简化模型，我们分别对古塔的倾斜、弯曲、扭曲进行讨论。

对于倾斜，首先根据不同年份，各层中心点坐标的z坐标值与其相应点的x,y坐标值做多元线性回归。

然后得到的回归方程所表示的回归平面与z轴正方向的夹角就可以表示古塔的倾斜程度大小。

最后根据各层中心的分布和变化趋势方向，确定古塔的倾斜方向。

对于古塔的弯曲，首先求出每层高程绝对差t i h,，然后相邻两层的高程绝对差的变化量，表示古塔每层弯曲程度大小。

再根据每层出现高程绝对差t i h,的两个测量点的连
线，确定每层弯曲方向。

对于古塔的扭曲变形，首先每层选取两对相同的对测量点，并做连线。

然后通过每层对测量点的连线，分别与第一层相同对测量点的连线所成的角度的平均值来衡量古塔的扭曲情况。

（三）、对问题三的分析
问题三要求我们分析该塔的变形趋势，这个问题属于预测的数学问题。

对于这个问题我们一般用回归的方法来求解，得出倾斜指标、弯曲指标、扭曲指标对时间的回归方程，并作出各自的图像，观察趋势。

三、模型假设
1.假设古塔只存在倾斜，弯曲，扭曲的三种变形情况；
2.假设在1986年到2011年没有对古塔进行人为的保护，如加固或修补；
3.忽略1986年与1996年观测的第13层第5个测量点所少数据；
4.假设古塔的变形是连续的；
四、符号说明
五、模型的建立与求解
5.1问题一的求解：
为观察同层各观测点的大概位置，做出1986年古塔内同层观测点连线的俯视图进行分析，做出下图：
图1-1
图1-1是通过1986年每层各测量点的坐标点连起来的（用CAD 制）图，每层所测的点相交构成一个多边形，得到每层的近似平面图，可以近似地把每层当作正八边形。

根据正八边形图形特征，古塔各层测量点坐标（j x ，j y ，j z ）的平均值作为各层中心点坐标。

即：
8
1
8
j
j i x
X ==
∑， 8
1
8j
j i y
Y ==
∑ ， 8
1
8
j
j i z
Z ==
∑
算出的各层中心坐标如下：
表1-1.所测年数各层中心坐标表
1986年
1996年
楼层i X Y
Z
楼层i X
Y Z 1 566.6648 522.7105 1.7874 1 566.665
522.7102
1.7102 2 566.7196 52
2.6684 7.3203 2 566.7205 522.6674 7.3146 3 566.7735 522.6273 12.7553 3 566.7751 522.6256 12.7508 4 566.8161 522.5944 17.0783 4 566.8183 522.5922 17.0751 5 566.8621 522.5591 21.7205 5 566.8649 522.5563 21.716 6 566.9084 522.5244 26.2351 6 566.9118
522.521
26.2295 7
566.9468 522.5081 29.8369
7
566.9506 522.5042
29.8323
用上表得到的数据，把每层中心点的X，Y，Z坐标分别对时间t（设1986年为第一年，即t=1）做回归，得到下表的一系列回归方程，用以下的方程就能算出古塔任意一年任意一层的中心坐标,即为确定古塔各层中心位置的通用方法。

5.2问题二的求解：
问题二要求我们确定塔的倾斜、弯曲、扭曲等变形情况。

为了简化模型，我们分别对古塔的倾斜、弯曲、扭曲进行讨论。

5.2.1倾斜变形求解：
1) 倾斜度大小的求解 首先用第t 年中,各层中心点坐标的t i Z ,对t i X ,，t i Y ,做多元线性回归。

t i Z ,=t i t i Y a X a a ,2,10++ （1）
该回归方程在空间直角坐标系中是一个平面，表示各层中心近似所处的平面。

运用Excel 软件根据（1）式求解各年回归方程。

统计各年回归方程系数可得下表：
表2-1 回归方程系数
运用Matlab 做出2009年回归方程对应的回归平面（过程见附录1）。

图2-1
由图2-1可以直观地看出各层中心点贴近回归平面，证明上面所建立模型的准确性。

上述方法所得的回归平面与z 轴正方向的夹角可以表示塔的倾斜角。

原理解释如下：
图2-2
如图，空间直角坐标系（由CAD 制作）中有下列关系： 其中，平面ABC 是同一年塔内各层中心点的回归平面。

AB 垂直于OD ，AB 垂直于OC ，即AB 垂直于CD ，即∠ODC 为平面OAB 与平面ABC 夹角。

所以平面ABC 与z 轴夹角α为： ODC ∠-=︒90α
即角α为回归平面与z 轴正方向的夹角可以表示塔的倾斜角。

令t i Z ,=0，即方程为t i t i Y a X a a ,2,10++=0为AB 所在直线方程，所以OD 为点O 到AB 的距离，根据点到直线的距离公式可得：
2
2
210
22
21
2100
0a a a a
a a a a OD +=
+⨯+⨯+=
令x 和y 等于0，可以得出 0a OC = 根据正切性质得：
22
21)90tan(a a OD
OC +==-︒
α 根据反三角函数，可知：
)arctan()90(2
221a a +=-︒α (2)
根据（2）式可以得算出古塔每年的倾斜角α，列表如下：
表2-2 塔的倾斜角度，单位：（°）
夹角α的值可以表示古塔的倾斜程度大小。

2) 倾斜方向的求解
根据古塔各层中心点在水平面xoy 中的投影的分布和变化趋势，来确定古塔的倾斜方向。

下面以2009年数据为例。

用Matlab 作2009年各测点与中心点的平面图（过程见附录2）：
图2-3
年份 1986 1996 2009 2011 塔的倾斜角
α
0.4739
0.4782
0.6153
0.6208
1
2
3
4
5
6
7
8
由图2-3可以看出，各层中心点都大致分布在第2，6个测量点的对角线上。

再根据中心点投影位置随楼层的增加而自测量点2向测量点6移动。

可以知道古塔的倾斜方向大致是沿测量点2向测量点6方向倾斜。

5.2.2弯曲变形求解
1) 弯曲程度大小的求解
对于古塔的弯曲情况，我们通过每层平面倾斜的变化程度初步分析，然后再结合整栋古塔，得出古塔的大概外形，从局部到整体分析古塔的变形。

首先，计算第t 年，i 层测量点高程的绝对差t i h ,，为：
min max ,j
j t i z z h -= 它能直观地反映在各层内最大倾斜程度，但绝对差不能全面的表现出弯曲的情况。

而第t 年i 层到i+1层高程的绝对差的变化量t i i u ,1~+∆，可以反映相邻两层的弯曲的大小程度。

公式如下：
t i t i t i i h h u ,,1,1~-=∆++
运用excel 计算与排列每年各层的测量点绝对差，得下表：
通过上表中的绝对差计算第t 年i 层到i+1层高程的绝对差的变化量t i i u ,1~+∆。

得出1986年，1996年，2009年，2011年各个t i i u ,1~+∆值，得出下表：
从所得各年t i i u ,1~+∆值的统计可得出结论：在各年中都是5~6层间的t i i u ,1~+∆值为最大，从1986年到2011年，t i i u ,1~+∆值呈某部分增大，而且总体t i i u ,1~+∆的总值都在变大。

对统计数据进行分析：在5~6层间的倾斜程度变化很大，表现为5~6层间产生大的弯曲，在1986年到2011年塔的相邻层之间的倾斜变化程度越来越大。

表现为随着时间的推移塔的弯曲程度越来越严重。

t i i u ,1~+∆的值主要在第五与第六层之间有较大值，其他值都较小。

所以古塔主要弯曲的地方是在第5层与第6层。

所以用第5层到第6层间的高程绝对差的增量t u ,6~5∆的值来表示古塔的弯曲大小程度。

由上表可以知在每年中，每层的高程的绝对差都是较小的数，即表示该塔同层内高程起伏较小，基本处于同一水平面；在第5层至第6层高程的绝对差异常增大，且使在第6层以后每层的绝对差都稳定在较大数值；在2011年中，在相邻两层的高程的绝对差的数值变化较大。

图2-4
图2-4是用Matlab作出古塔的三维结构图（过程在附件2），可以看出古塔中间的地方比较弯曲，由此也说明了我们所建模型的合理性。

2)弯曲方向的求解
列出在第5、6层各点高程从高到低的测量点序号：
表2-5 第5、6层各年各测量点的排序
1986年1996年2009年2011年
第5层
2 2 2 2
8 8 8 1
1 1 1 7
7 7 7 5
3 3 6 6
6 5 3 4
5 6 5 3
4 4 4 8 1986年1996年2009年2011年
第6层2 2 3 3
3 3 2 2 1 1 1 1
4 4 4 4 8 8 8 8
5 5 5 5 7 7 7 7
从上表可得出结论：每年的第5层楼到第6层楼，测量点6相对其他点降低的程度最大，即可知在第5层到第6层楼中古塔向测量点6方向弯曲。

5.2.3扭曲变形求解
对于古塔的扭曲变形，首先每层选取两对相同的对测量点，并做连线。

然后通过每层对测量点的连线，分别与第一层相同对测量点的连线所成的角度的平均值来衡量古塔的扭曲情况。

由于古塔首层是与稳定的地基连接，所以首层的扭曲变形几乎可以忽略不计，所以取每层相同对角线在平面xoy 投影的夹角作为扭曲变形的量度：夹角越大，则扭曲变形越严重。

具体原理如下：
如图2-5所示，以第1与第13层为例，首先选取两对测量点，分别为2-6、4-8。

取第一层的2-6对角线为L1，第13层的2-6对角线为L2，L1与L2的夹角为62,113--λ。

取第一层的4-8对角线为L3，第13层的4-8对角线为L4，L3与L4的夹角为84,113--λ。

用62,113--λ与84,113--λ的平均值113-λ作为古塔第13层相对于第一层扭曲变形的量度。

图2-5
同理，其他层2-6测量点对角线与第一层2-6测量点对角线L1的夹角为62,1--i λ，其他层
4-8测量点对角线与第一层4-8测量点对角线L3的夹角为6 6 6 6
62,113--λ
84,113--λ
1 2
3
4 5
6
7 8
L3
L4
84,1--i λ。

它们的平均值1-i λ反映其他层相对第一层的扭曲度。

即： 28
4,162,11-----+=i i i λλλ
计算结果如下表
表2-7 每层相对于第一层扭曲的度，单位：（°）
得出结论：第2层到第10层随着层数的增高1-i λ的数值越来越大，第11层到第13层随着层数的增高1-i λ的数值越来越小。

表现为在第2层到第10层的各层塔向逆时针扭曲，第11层到第13层的各层塔向顺时针扭曲。

5.3问题三的求解：
5.3.1塔倾斜变形的趋势：
要确定塔的变形趋势，对于倾斜变形，我们利用问题二中所求出的塔的倾斜角α来对时间做回归，得到4478.0,006598.021==p p ，可以知道回归方程为： t 006598.04478.0+=α
利用上述公式可以得到2011年之后每年的倾斜角，列表如下：
图3-1 塔的倾斜角与时间的变化曲线
由上图所得出来的拟合曲线图我们可以得出结论：随着时间的增长，塔倾斜所成的倾斜角度值会越来越大。

预测在2025年塔的倾斜角度值达到0.7º，2039年塔的倾斜角度值达到0.8°。

5.3.2塔弯曲变形的趋势：
讨论该塔弯曲程度的变化趋势，即可变为讨论不同时间第5层到第6层间的高程绝对差的增量t u ,6~5∆随时间变化的变化。

因为每年都是在第5层到第6层为最大弯曲，也为整个塔的主要弯曲的表现，所以可用t u ,6~5∆对t 进行回归来说明塔随时间的变化导致其弯曲程度的变化，对回归方程求解；
t u ,6~5∆=t q q 12+
运用Excel 软件的LINEST 函数求解该回归方程系数为：
2q =0.11461 1q =0.001186
可以得出：
t u t 001186.011461.0,6~5+=∆ （4）
用（4）式在excel 软件中画出弯曲程度的拟合曲线，如下：
图3-2 塔的弯曲程度与时间的变化曲线
预测在2011年以后，每年塔的弯曲程度都会增大。

预测在2020年塔弯曲程度为0.15612，即为塔在第5层到第6层的测量点最大高程差的变化增加了0.15612米，在2050年他弯曲程度为0.17391，即为塔在第5层到第6层的测量点最大高程的变化增加了0.17391米。

5.3.3塔扭曲的变形趋势：
由于每层相对于第一层都存在一定的扭曲，因此取各层对第一层的扭曲度数的平均值作为古塔的扭曲程度λ。

用λ对时间t 的回归求回归方程，这样就得到古塔扭曲变形的趋势。

回归方程为：
t p p 12+=λ
运算得到1p 与2p 的值(过程见附录.3)，则回归方程解析式为;
t ⨯=0.02031-2.862887λ
图3-3
由图3-3可以知道用excel 作出来的曲线的斜率为负值，因此古塔的扭曲程度随时
间的增加越来越小。

六、模型的评价与推广
6.1 模型的评价
本模型从实际出发，分析了多种情况，将古塔的变形分成倾斜、弯曲、扭曲三种类型来分别讨论。

先定性地分析古塔的变形，再定量解出古塔变形了多少。

在通过结合回归思想，将变形的数据转化为回归方程，建立了古塔变形趋势的模型，使得对古塔的变形了解得更全面，能为古塔有效地制出保护措施提供依据。

同时模型也存在着不足的地方，由于题目给出的数据不是很全面，所以很难算出古塔准确的变形，只能大概地算出一些近似值。

但这模型也有很大的意义，不但对该古塔适用，也能把它用于其他古塔的变形分析，从而保证了对古塔进行维护的快速性和高效性。

6.2 模型的推广
由图1-1中同层视为平面时可看出古塔每层的相邻连线图像相似于是正八边形，但相对于正八边形确有一定差距，但从古塔的设计理念应该使古塔完全成为正八边形从而增强整个塔的稳固性，但是实际塔的形状确实产生了误差从而使整个塔的形状产生的变化。

因为原本的形状为8条等边的正八边形，而古塔各层形状的变化必引起8条边的长度变化，使8条边的长度变为不等，即可用八边形边长的长度的标准差来判断该层形状的变化程度，用各边长长度大小的比较来判断该层变化大概形状。

七、参考文献
【1】数学建模论文格式要求 360问答 /q/65170
【2】李大潜主编中国大学生数学建模竞赛高等教育出版社 (1998)
附件
附录一：
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-')
plot3(mean(x1(i,:)),mean(y1(i,:)),mean(z1(i,:)),'k.')
cen(i,:)=[mean(x1(i,:)),mean(y1(i,:)),mean(z1(i,:))]
end
plot3(d(1,1),d(1,2),d(1,3),'k.')
[X,Y]=meshgrid(cen(:,1),cen(:,2))
Z=-58296.37472+92.78414147.*X+10.94248841.*Y
mesh(X,Y,Z)
axis([560 575 515 530 0 60])
xlabel('x', 'FontWeight', 'bold');
ylabel('y', 'FontWeight', 'bold');
zlabel('z', 'FontWeight', 'bold');
zoom on
axis square
view(45,45)
附录二：
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
2013数学建模——古塔的变形
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]
x=[x1,x1(:,1)]
y=[y1,y1(:,1)]
hold on
for i=1:1:13 plot(x(i,:),y(i,:),'k.')
plot(x(i,:),y(i,:),'k-')
plot(mean(x1(i,:)),mean(y1(i,:)),'k.')
end
axis([560 575 515 530])
xlabel('x', 'FontWeight', 'bold');
ylabel('y', 'FontWeight', 'bold');
zoom on。