古塔的变形
古塔的变形
摘要：
在本文中，我们将对管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测的数据进行分析处理，来获得塔的变形情况，并进行拟合。
在这里，我先对各个变形进行描述：
倾斜：建筑中心线或其墙、柱等，在不同高度的点对其相应底部点的偏移现象。在本文中定义α，来表述塔的倾斜度。通过对中心线的直线拟合来确定倾斜度。
弯曲：按照材料力学里的定义，当杆件受到与杆轴线垂直的外力或在轴线平面内的力偶作用时，杆的轴线由原来的直线变成弯曲，这种变形叫弯曲变形。当发生弯曲变形之后，轴线将变成一条曲线，而这条曲线，我们称之为挠曲线。杆件对其原来位置转过的角度θ，来表征弯曲程度，我们可以对挠曲线进行求导，
dw dw
=K 。但在本文中采可求得各个点的导数，来进行对θ的求解。即tan θ=
dx dx
用的是将各个层的平面拟合出来，那么转过的角度θ就可以通过平面的法向量来求解。
扭曲：建筑产生的非竖向变形。由于扭曲为非竖向的变形，讨论古塔扭曲时只需考虑水平方向的坐标变化，即x ，y 坐标的水平旋转，因此我们用古塔水平旋转角度的扭曲度γ来描述。
这里需要强调的是，每个变形都需要一个相对量，在倾斜度上有竖直线；在弯曲时，也可以选择竖直线来进行参照；而对于扭曲程度，则没法选择一个原有的参照量，所以在本文当中，对于扭曲程度，我选择1986年的数据来做为参照量来进行分析。
最后通过多项式拟合得到的结果有
α=0.56865-0.031985t +0.004411t 2-0.000103t 3θ=1.577459+0.017723t -0.0019574t 2+0.000055t 3
，，
γ=-0.444206t +0.098426t 2-0.00247t 3，
这里的时间是经过变换的，从1986年开始t =0，然后以年份为单位。
问题陈述：
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。
某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后
于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。
请你们根据附件1提供的4次观测数据，讨论以下问题：
1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。
2. 分析该塔倾斜、弯曲、扭曲等变形情况。 3. 分析该塔的变形趋势。
对于问题1表述：
将八边形分成6个三角形，并分别将各个三角形的形心位置球出来，然后再根据比重，求出八边形的中心位置，即形心位置。
并且在数据当中有些数据缺少了，我将通过插值的方式来将这些数据补上。 对于问题2表述：
根据摘要里的α、θ、γ来分析该塔的倾斜、弯曲、扭曲等变形情况。 对于问题3表述：
θ、γ的数据情形来分析该塔的变形趋势。 结合Matlab 所画出的图形，和α、最后，我想通过建模对α、θ、γ给出确定的表达式，并进行误差分析。
模型准备：
画出1986年塔的大概样子：
据此图片可以分析，并得出几个假设。
模型假设：
1、假设每层的截面的形心位置就为重心位置，即此塔质量均匀分布。 2、假设该塔刚建好时，轴心线是竖直的。
3、假设在1986年时，该塔的扭曲情形较小，可以忽略不计。 4、假设该塔的截面原本就是个正八边形。
填补数据：
根据测来的数据，我发现在1986年和1996年的数据当中的第13层的第五个点缺少了。所以我将用插值的方法来进行填补。（对于x 和y 、z 则用）拉格朗日插值函数来预测得到）
拉格朗日插值函数：
L n (x ) =∑y k l k (x )
k =0
n
其中：
(x -x 0) ∙... ∙(x -x k -1)(x -x k +1) ∙... ∙(x -x n )
l k (x ) =
(x k -x 0) ∙... ∙(x k -x k -1)(x k -x k +1) ∙... ∙(x k -x n )
带入数据求得：
由于本来在此层的数据并不多，所以我决定用上这一层的所有数据，我们设
x 和y 、z 都是测量点[1,2,3,4,5,6,7,8]的函数。
所以得到1986年的缺少的观测数据：x =568.3565;y=519.6883;z=52.7791 1996年的缺少的观测数据：x =568.0712;y= 519.6821;z=52.7785
确定各层之间的中心：
首先，在原数据当中，八个点肯定不可能都在同一个平面上，所以我这里先进行平面拟合，并将八个点投影到此平面上：
设平面方程为A x +B y +C +z
D 0=，则各个点到此平面的距离
为
8
d i =
⎧∂T ⎪∂A ⎪
⎪∂T ⎪∂B ⎨
⎪∂T ⎪∂C ⎪∂T ⎪
⎩∂D
=0=0
，所以令T
D
=∑d i 2。
i =1
根据最小二乘法得
=0=0
通过次式子来得到A , B , C , D 。接下来我将八个点投影到此平面上，设投影点为(x i , y i , z i ) ，而原来的点为(x j , y j , z j ) 。则得到以下方程：
x i -x j y i -y j
=
A x i -x j A ; =; Ax +By +Cz +D =0。联立这三个方程可以解出B z i -z j C
(x i , y i , z i )
接下来我将通过新获得的点通过以下方法来计算出中心： 设中心坐标为(x , y , z ) 。
设从红色起的第一个点为点1，顺时针依次为点2，点3，点4，点5，点6，点7，点8。
(x 1, y 1, z 1);(x 2, y 2, z 2);(x 3, y 3, z 3);(x 4, y 4, z 4) (x 5, y 5, z 5);(x 6, y 6, z 6);(x 7, y 7, z 7);(x 8, y 8, z 8)
将每个三角形的形心(x , , y , , z , ) ，面积S 求出来，从右到左依次为三角形1，三角形2，三角形3，三角形4，三角形5，三角形6。
x 1+x 2+x 3y 1+y 2+y 3z 1+z 2+z 3⎧, , ,
(x , y , z ) =(, , ) ⎪111
333
⎪
⎪(x , , y , , z , ) =(x 1+x 4+x 3, y 1+y 4+y 3, z 1+z 4+z 3) ⎪222333⎪
x +x +x y +y +y 5z 1+z 4+z 5, , , ⎪(x 3, y 3, z 3) =(145, 14, )
⎪333
⎨
x +x +x y +y +y z +z +z 6⎪(x , , y , , z , ) =(156, 15, 156) 444
⎪333⎪x +x +x y +y +y 7z 1+z 6+z 7, , , ⎪(x 5, y 5, z 5) =(167, 16, )
333⎪
⎪, , , x 1+x 7+x 8y 1+y 7+y 8z 1+z 7+z 8(x , y , z ) =(, , ) ⎪666
333⎩
利用向量叉乘来计算面积S ，这里只给出三角形1的面积计算公式。
，其中
S 1=
a =x 2-x 1; b =y 2-y 1; c =z 2-z 1d =x 3-x 3; e =y 3-y 1; f =z 3-z 1
其他的与三角形1的算法一样得到S 2; S 3; S 4; S 5; S 6。 所以根据重心的计算公式可以得到重心坐标(x , y , z ) 。
, , , , ,
x 1, S 1+x 2S 2+x 3S 3+x 4S 4+x 5S 5+x 6S 6
x =
S 1+S 2+S 3+S 4+S 5+S 6, , , , ,
y 1, S 1+y 2S 2+y 3S 3+y 4S 4+y 5S 5+y 6S 6
y =
S 1+S 2+S 3+S 4+S 5+S 6, , , , , , z 1S 1+z 2S 2+z 3S 3+z 4S 4+z 5S 5+z 6S 6
z =
S 1+S 2+S 3+S 4+S 5+S 6
将数据代入得到以下表格：（对于1986、1996年的塔顶的中心点取其各项的
倾斜：
在求得中心坐标后，我将对其进行直线拟合：
这里是进行空间三维的直线拟合，所以我将其转化为平面二维的直线拟合。 在xOy 平面上：
设方程为y =h 2x +h 1，则通过最小二乘法得到以下方程组：
⎡n ⎤⎡n ⎤n , x y ∑∑i i ⎥⎢⎥⎡h ⎤⎢n n n n
i =1i =112⎢⎥⎢⎥=⎢⎥，令q 1=∑x i ; q 2=∑x i ; q 3=∑x i y i ; q 4=∑y i ；n n n ⎢⎢⎥i =1i =1i =1i =12⎥⎣h 2⎦x , x x y ⎢∑i ∑i ⎥⎢∑i i ⎥
i =1⎣i =1⎦⎣i =1⎦
q 2q 4-q 1q 3nq 3-q 1q 3
。 ; h =2
nq 2-q 12nq 2-q 12
得到：h 1=
在yOz 平面上：
,
1
,
q 2, q 4, -q 1, q 3, , nq 3, -q 1, q 3,
设方程为z =h +h y ，则同理可得：h =。 ; h 2=
nq 2, -q 1,2nq 2, -q 1,2
,
1
, 2
, 1
,,
在xOz 平面上：设方程为z =h 1,, +h 2x ，则同理可得：
, ,
1
, 3
,
, ,
, ,
q 2, q , 4-, q , q 1, 3, , nq -q 3q ,
h =; h =2,,2,,2
nq 2,, -q 1nq 2,, -q 1
这些方程在三维坐标中都是垂直于相对应的坐标平面的平面，则两两平面相交的直线有：
h 1,, h 1, x +,, y +,
h 1y z -h 1, y -h 1z -h 1,, h 2h 2z
。 1. x +==;2. x ==;3. ==, ,, , ,, , ,,
h 2h 2h 2h 2h 2h 2h 2h 2h 2h 2通过各个点与各个直线之间的距离，来选择哪条直线。 通过计算得到以下结果：
α1=0.56865; α2=0.58673; α3=0.91080; α4=0.91340
弯曲：
一、在上述拟合平面的时候，可以获得每层平面的法向量n =(N 1, N 2, N 3)
→
→
→
弯曲转过的角度θ，可以通过n =(N 1, N 2, N 3) 与n 0=(0,0,1)来计算，其实这两个向量之间的夹角就是弯曲转过的角度θ。所以就有：
→
→
cos θ=
n ∙n 0
→
n n 0
→
鉴于上述方法并不能求得塔顶的弯曲程度，所以这里又采用了三次样条插值函数来求得塔顶的弯曲程度，并且以塔顶的弯曲程度为整体的弯曲程度。原因是在塔顶上的弯曲程度最大。下面先介绍下三次样条插值函数：
这里采用的边界条件为自然边界条件。设各个点的二次导数为
M =[M 0, M 1... M n -1, M n ], ，因为是采用自然边界条件，所以M 0=0; M n =0。
有(x 0, y 0)...(x n , y n ) ，根据三次样条插值可以得到下列方程组:
⎧h i =x i -x i -1⎪
⎧2M 1+λM 2=g 1⎪μi =h i ⎪μM +2M +λM =g ⎪h i +h i +1⎪21232，其中(i =1, 2,..., n -1) 。⎨⎨
⎪λi =1-μi ⎪.......
⎪⎪y -y y -y i -16⎩μn -1M n -2+2M n -1=g n -1
(i +1i -i ) ⎪g i =
h +h h h i i +1i +1i ⎩
(x 13-x 12) 2(y 13-y 12) h 13(M 12-M 13)
所以获得的端点的弯曲为S (x 13) =M 13，++
3h 13h 136
,
13
再通过导数，获得此点的方向向量。
在本文中，这里是三维的三次样条插值函数，所以我先将其投影到xOy , yOz 的平面上做二维的三次样条插值函数，再通过过这两个平面上的方向向量，来获得最终要获得的三维上的方向向量。
带入数据得到：
由于时间关系，和知识所限，所以这里并没有给出结果。我以第一种方法的结果作为结果。
扭曲：
由于要知道变形程度的大小，则必须知道原本的样子，所以我在这里将以1986年的数据为参照量来进行分析。在这里我们就先把每层各个相对应的点，进行最小二乘法拟合成一条直线，再通过与1986年的进行比较，求出两条直线之间的夹角，则我定义这两条直线之间的夹角就为扭曲度γ。在讨论倾斜时就已经讨论过空间三维直线的拟合，所以在这里不多赘述。所以每年的拟合出来的直线有8条，所以我就将这八条算出来的扭曲度进行平均，来作为最终的结果。
假设1986年的每条拟合直线的方向向量为n i =(l i , m i , n i ) ，而其他年份的为
→
→
N i =(a i , b i , c i ) 。所以cos γi =
n i ∙N i n i N i
→
→
。所以=
∑γ
i =1
8
i
8
带入数据得到1986年的八条直线的方向向量分别为：
n 1=(0.0001,0.2226,-1.6473); n 2=(0.0010,0.8786,-6.1547) n 3=(0.0010,0.4503,-5.6720); n 4=(0.0010,-0.1136, -9.8889) n 5=(0.0010,2.2330,5.5791);n 6=(0.0010,0.8563,7.5760)n 7=(0.0010,0.3137,6.7957);n 8=(0.0001,-0.0777,1.4506) 而根据其他测出的数据得到几个扭曲度：
1=2.9308; γ2=11.8015; γ3=11.8217
拟合：
进行多项式拟合得到：α=
0.56865-0.031985t +0.004411t 2-0.000103t 3 进行多项式拟合得到：θ=
1.577459+0.017723t -0.0019574t 2+0.000055t 3 进行多项式拟合得到：γ=-0.444206t +0.098426t 2-0.00247t 3
参考文献: 1、《数值计算方法》 朱建新 李有法编 高等教育出版社 2、《详解Matlab 在科学计算中的应用》 陈泽 占海明编 电子工业出版社 3、《空间直线拟合研究》 王伟峰 温耐著 许昌学院报第29卷第五期 2010年9月