古塔的变形
摘要
本文对古塔的各层中心坐标，倾斜、弯曲、变形问题，采用数据插值、最小二乘法、回归分析法及Excel软件与Matlab软件，进行求解，具体步骤如下：
1．第13层中第5个观测点空白数据处理
对1986年和1996年，第13层中第5个观测点进行了估计，采用数据插值方法中的3次样条插值，算出1986年(xi,yi,zi)数据分别为566.308，519.7624， 52.7686及1996年(xi,yi,zi)数据568.0575，519.7562，52.7657 。
2．模型建立及求解各层中心坐标
将空间中每层观测的点(xi,yi,zi)分别投影到xoy,yoz,xoz各平面，通过最小二乘法确定出以投影的点为中心的中心直线，以三条直线为母线做平行于三角坐标的平面，求解三个平面的交点即为中点(x0,y0,z0)的坐标。 利用Excel及Matlab软件求解出各层中心坐标。
3．古塔倾斜分析
将古塔各层相邻的中心坐标连接起来，构成空间向量，通过求解向量的夹角
?、?、?来确定古塔的倾斜情况。计算与仿真结果表明，古塔中心各层间都有
倾斜，但在8-9层，12-13层时倾斜更为严重。
4．古塔扭曲，变形分析
将中心坐标投影到xoy面，计算相邻层上投影点(xi,yi)构成的直线的斜率
k?tan?，得到夹角?。该角度说明了中心发生倾斜时的方向、角度大小描述了
扭曲程度大小。结果表明古塔朝着xoy平面的第四象限偏移大。
5．古塔变形预测
采用多元回归分析法，根据最小二乘法原理，拟合出古塔最终变形的趋势。从图中可以看出古塔是各层不断的发生变化。并且在8-11层变形严重。
关键词：中心直线 插值 曲线拟合 回归分析
1
一、问题重述
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。
某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。
请你们根据附件1提供的4次观测数据，讨论以下问题：
1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。
2. 分析该塔倾斜、弯曲、扭曲等变形情况。 3. 分析该塔的变形趋势。
二、模型假设
1．假设对模型计算结果均保留四位有效数字，不影响其变形趋势； 2. 假设任意两层中心的连线为直线时，考虑其倾斜程度； 3. 假设任意两层中心的连线为曲线时，考虑其弯曲、扭曲程度。 4.假设古塔是规则图形；
5.假设古塔刚建好时地面平坦没有倾斜和扭曲。
三、符号说明
1．(xi,yi,zi)，i?1,2,,n表示空间中点的坐标；
2.a,b,c,d,e,f为投影平面直线方程的参数；
3.M表示空间投影到平面上的各点距离之和的最小值； 4.(x0,y0,z0)三个平面的交点；
5.?、?、?（0????，0????，0????）表示向量的方向角； 6.?两点构成的直线的倾斜角； 7.b0,b1,b2回归中未知参数。
2
四、模型建立及求解
1. 中心位置模型建立
以古塔其中一层为研究对象，建立中心位置的数学模型。
设(xi,yi,zi)，i?1,2,,8是古塔第i?1,2,,n层的观测的数据，为空间中
的点。首先，将(xi,yi,zi)，i?1,2,,n投影到xoy平面上，得到观测数据记为
(xi,yi)，i?1,2,,n。
在平面直角坐标系中，建立以这些点构成的中心直线方程，记为
y?ax?b，要求使各点(xi,yi)，i?1,2,n,n到直线的距离之和最小，即用最小二
乘法。将问题转化为求解M=?[yi?(axi?b)]2的最小值。由二元函数极值原理，令
i?1n?M=?2?[yi?(axi?b)]xi?0 ?ai?1n?M=?2?[yi?(axi?b)]?0 ?bi?1即可求得a、b，确定出直线方程
y?ax?b （1）
i?1,2,同理，将(xi,yi,zi)，i?1,2,,n投影到yoz平面上，得到观测数据记为(yi,zi)，
,n。采用二元函数极值原理，可确定出直线方程
z?cy?d （2）
将(xi,yi,zi)，i?1,2,i?1,2,,n投影到xoz平面上，得到观测数据记为(zi,xi)，
,n。可确定出直线方程
x?ez?f （3）
确定出的这三条直线（1）、（2）、（3）在空间分别表示三个母线依次平行于三个坐标轴的平面。求解出三个平面的交点，就是我们要找的每层的中心，记作
(x0,y0,z0)。
2.模型求解
（1）题目中给出的1986年和1996年观测数据中，第13层中第5个观测点没有数据。
利用数据插值的方法，将这个数据估计出来。在此，利用Matlab中3次样条插
3
值，给出缺失数据。调用函数spline(X,Y,X1)，数据计算如下所示：（计算的缺失数据加黑显示）
表一 1986年第13层观测值 表二 1996年第13层观测值
y/m z/m x/m y/m z/m 13层 x/m 13层 1 566.308 525.092 52.866 1 566.3142 525.0857 52.864 2 564.716 523.616 52.878 2 564.7222 523.6098 52.877 3 564.418 521.521 52.897 3 564.4872 521.5148 52.889 4 565.91 519.893 52.88 4 565.9161 519.8868 52.878 5 566.308 519.7624 52.7687 5 568.0575 519.7562 52.7657 6 569.701 521.05 52.703 6 569.7072 521.0438 52.696 7 569.897 523.188 52.794 7 569.9032 523.1817 52.789 8 568.582 524.822 52.822 8 568.588 524.8157 52.817 （2）将每一层都采取上述的建模方法，计算每层的中心位置，其中n?8。计算时采用Excel、Matlab软件对数据进行求解，求解过程见附件一。解得每年各层中心位置坐标如下表所示：
表三 每年各层中心位置坐标
层层X Y Z 数 数 1 566.6645 522.7131 1.7794 1 2 566.7196 522.6549 7.3046 2 3 566.7729 522.6038 12.7386 3 4 566.8161 522.6068 17.0532 4 5 566.8621 522.5907 21.7445 5 6 566.9084 522.5435 26.2351 6 7 566.9408 522.5066 29.8322 7 8 566.9468 522.4194 33.3356 8 9 567.0218 522.5105 36.8667 9 10 567.0569 522.5449 40.194 10 11 567.1045 522.4644 44.4341 11 12 567.1518 522.4153 48.7354 12 13 567.2011 522.738 52.83 13 14 567.2471 522.2437 55.12335 14 2009年观测数据 中心点坐标 层层X Y Z 数 数 1986年观测数据 中心点坐标 1996年观测数据 中心点坐标 X 566.665 566.7205 566.7751 566.8183 566.8641 566.9118 566.9506 566.9884 567.0265 567.0169 567.1102 567.1578 567.2012 567.2544 Y 522.7097 522.6598 522.5937 522.6247 522.5865 522.5699 522.4961 522.4591 522.4542 522.5577 522.3371 522.363 522.3684 522.2367 2011年观测数据 中心点坐标 Z 1.7773 7.3235 12.7278 17.056 21.7367 26.2518 29.8261 33.3284 36.8404 40.1872 44.392 48.6833 52.8266 55.1198 x y z 4
1 2 3 4 5 6 7 8 9 10 11 12 13 14 566.7267 522.6787 566.764 522.6517 566.8002 522.6326 566.8293 522.6228 566.8603 522.4884 566.9469 522.5532 566.9793 522.5376 567.0305 522.4588 567.0815 522.4465 567.1368 522.3662 567.1798 522.3755 567.2224 522.3348 567.2712 522.249 567.336 522.2146 1.7437 7.2881 12.7147 17.0442 21.6853 26.2178 29.8348 33.3533 36.865 40.1744 44.4228 48.7135 52.8053 55.089 1 2 3 4 5 6 7 8 9 10 11 12 13 14 566.72695 522.7234 566.7642 522.6899 566.80044 522.6553 566.82971 522.5793 566.86096 522.4982 566.94784 522.5478 566.98 567.0313 567.0823 522.5181 522.4531 522.4471 1.7792 7.2951 12.7128 17.0387 21.6821 26.1841 29.8301 33.3284 36.802 40.1375 44.4024 48.7088 52.8037 55.087 567.13806 522.3649 567.18096 522.3752 567.22383 522.3329 567.27251 522.2473 567.3374 522.2132 通过数据，可以看出古塔各层中心坐标都有不同幅度的变动，为了更加清楚的描述数据间的关系，通过Matlab软件仿真，得到下图：
图1 1986年中心坐标仿真图 图2 1996年中心坐标仿真图
图3 2009年中心坐标仿真图 图4 2011年中心坐标仿真图
5