1、古塔变形1古塔变形一、问题分析：通过对数据的初步观察，知道个问题其实就是一个通过古塔的观测数据找到古塔变形的规律进而推测其未来发展的问题。我认为题目的关键是对这些数据进行拟合，然后根据函数的变化规律来进行预测。先对古塔大致形态进行分析，我先用 Matlab 将原始数据用三维图的形式描绘出来，这样有助于我对问题进一步分析与求解：由于前两次数据有缺失，故选择 2009 年的数据：古塔侧视图：古塔变形2古塔俯视图：通过对这个大致图形的观察，我掌握了古塔的一些基本信息，这是一个十三层的八角古塔，每一层都是由正八边形（由于测量误差或变形可能不那么正，但古人肯定是这么设计的） ，于是可以大胆假设古塔的各层
2、都是正八边形。二、模型假设：1）假设古塔的各层质量分布均匀，即重心即为古塔几何中心2）假设古塔的各层都是正八边形，即便有些数据误差，不过可以忽略二、符号申明：古塔变形3三、模型求解：由于 1986 年和 1996 年的数据有缺失（如下图 13 层第 5 点）：这对于利用几何法求中心是有阻碍的，于是先将它们补齐。根据假设 2：古塔各层都是正八边形，故这里用几何法进行推算:由于是正八边形，如图，E 为所求的点：由几何关系应该有： ahefabedkk用坐标表示出来：即为：古塔变形4efabefahedeyyxxyyxx纵坐标可直接利用其余各点平均值求得。利用 Matlab 求解这个方程，分别将 1
3、986 和 1996 的数据带入方程，解得缺失的中心点的坐标分别为：年份 X Y Z1986 568.1644 519.6253 52.834291996 568.1704 519.6191 52.83于是现在每一次数据都是齐全的了，下面来求各层的中心坐标点：首先给出各层中心点求法的通用方法：首先，先将各层都拟合到一个平面，由于竖直面上各层差异不大，于是采用对各层各点 z 坐标取均值的方法将古塔的每一层拟合到同一个平面上,然后再求中心坐标，根据假设 1：塔的各层分布均匀，故几何中心即为重心，在网上找到一篇论文就是研究关于求多边形几何中心的方法，这里不妨直接拿来用一下：古塔变形5在这里，只需将
4、n 设置为 8，利用 Matlab 求解公式即可得到各层的中心坐标。得到的结果，如下：1986 年：566.6648 522.7105 1.787375566.7196 522.6684 7.32025566.7735 522.6273 12.75525566.8161 522.5944 17.07825566.8621 522.5591 21.7205566.9084 522.5244 26.23513566.9468 522.5081 29.83688566.9843 522.4924 33.35088567.0218 522.4764 36.85488567.0569 522.4624
5、40.172131996 年：566.665 522.7102 1.783566.7205 522.6674 7.314625566.7751 522.6256 12.75075566.8183 522.5922 17.07513566.8649 522.5563 21.716566.9118 522.521 26.2295566.9506 522.5042 29.83225566.9884 522.4881 33.34538567.0265 522.4714 36.84825567.062 522.4572 40.167632009 年：古塔变形6566.7268 522.7015 1.76
6、45566.764 522.6693 7.309566.8001 522.6384 12.73225566.8293 522.6132 17.06975566.8604 522.5866 21.70938566.9471 522.5342 26.211566.9792 522.5123 29.82463567.0305 522.4797 33.33988567.0816 522.4466 36.84375567.137 522.3937 40.161132011 年：566.727 522.7014 1.76325566.7642 522.669 7.2905566.8004 522.6387
7、 12.72688566.8297 522.6127 17.052566.861 522.586 21.70388566.9478 522.5335 26.2045566.98 522.5115 29.817567.0313 522.4788 33.33663567.0825 522.4457 36.82225567.1381 522.3926 40.14413散点连线图：古塔变形7将四年的放在一起看：古塔变形8现在大致可以知道 1996 到 2009 年这个时间段古塔变形较大，下面进行进一步量化分析：（一） 、倾斜程度：对中心点进行线性拟合，求出拟合直线与竖坐标轴的夹角，夹角的余弦值可作为倾
8、斜程度的标识：古塔变形9即： 222cos| llllxnlayzA对各个年份的中心点进行拟合，采用线性插值法得：古塔变形10放大后：解得的数据为：年份 余弦值 角度值 正切值 余弦值1986 0.010654 1.560142 93.85223 0.0106541996 0.010775 1.560021 92.80484 0.0107752009 0.011579 1.559217 86.35388 0.0115792011 0.011604 1.559192 86.17069 0.011604有表格可以看出，塔身与地面的夹角从 1.560142 变化到1.559192，逐渐靠近地面，下面
9、利用 Matlab 对其拟合分析：古塔变形11原始散点图：年份-角度用三次样条插值后：古塔变形12由模拟的函数图可知，古塔的倾斜度每年都有所增加，但是在 2000年到 2010 年这个时间段内变化特别明显，这说明在这个时间段内自然环境可能有比较剧烈的变化，联想到 2008 年我国的汶川地震，可以知道那段时间我国的地壳运动是比较活跃的，由此猜测也许与地壳活动剧烈有关，到了 2010 年以后，变化速度急剧变慢，回到了90 年代的速度，由此我大胆的做出一个与本题无关的推测，趋势变缓意味着地壳运动变缓，也就是说在未来的一段日子里我国应该不会发生强烈的地震灾害，所以大家可以安稳的生活。同时，我对古塔未来
10、一段时间的倾斜度变化做出一个预测：变化存在，但很缓慢，与地面的夹角逐渐在变小，由图形走向来看，可能趋向停止（如果没有不可抗力因素，如地震狂风的影响） 。（二） 、弯曲程度：由于弯曲程度是指中心轴线的弯曲度，将各个中心点连接起来，依次求出相邻点的之间的夹角，然后取平均值，作为弯曲度的度量：用 MatLab 解得为：年份 余弦 角度（弧度）1986 -9.2E-05 1.5708891996 -0.00019 1.5709832009 -0.00029 1.5710872011 -0.0004 1.571193古塔变形13用三次样条拟合后得到的图像：弯曲程度在 1986 到 2007 年之前都基本
11、保持线性，但是 2008 到2010 突然变为指数级增长，估计为自然灾害影响。所以可以预测，2011 年以后，若无自然灾害或其他不可抗力的因素的影响，弯曲程度将与 1986 年左右持平，每年大约弯曲 ,其实这个只是很76.10rad小的，至少要等上千年才能看到较明显的变化。分析这个问题我还采用了另外一种做法，就是将各层坐标投影到xoz 平面上，然后通过观测曲线曲率的变化来分析弯曲程度的变化：具体做法是，分别对每年都做一次三次样条插值，然后在曲线上抽样一些点，取这些点的曲率的平均值作为该年度塔的弯曲度的特征值，然后通过对改特征值的分析来找到弯曲度变化的规律。古塔变形14对曲率拟合后的到的结果得到
12、的结果与前面基本相同。（三） 、扭曲程度：平面间的旋转角度可作为扭曲的度量标识，先将各层中心与塔角的向量与底部的同位置向量求余弦，然后将所有的点都求一次，取平均值作为该年度塔的弯曲程度。古塔变形15原理示意图:余弦公式： 12121122cosxyx将这些数据利用 MatLab 求解得到四年的数据如下：年份 余弦值1986 0.9985331996 0.9985332009 0.9989162011 0.998916古塔变形16现在对余弦值进行分析,进行三次样插值后得到的结果：由图中可以看出，其实在 1996 年以前几乎是没有多少扭曲的，但是在 1996 到 2009 年这段时间里变化十分巨大
13、，由此估计在这段时间发生过一些比较特别的事，比如地震，这和之前的论断不谋而合，但到了 2011 年以后，角度变化又趋于平缓，可以预测，在未来的几年内如果没有不可抗力因素的影响，古塔的是不会发生太大的扭曲的。现在来观察一下古塔扭曲的方向：我的做法是先将各个平面都投影到地面，然后平移他们，使他们几何中心重合到一点，然后将他们边角上的 按顺序相连和拟合得到古塔的扭曲形状。坐标变换公式为：古塔变形170(,)(,)(,)xyxyxy为原坐标,为中心点坐标0()xy为平移后的坐标,对边沿点的拟合曲线为：理论上如果古塔没有扭曲，那么边沿拟合线应该是直线，但现在边沿线呈现出了螺旋的形状，并且还呈现出了逆时针
14、扭曲的趋势，至于为什么会出现逆时针扭曲趋势，其实可以用地转偏向力来解释，地转偏向力理论根据大地的磁场提出，北半球的物体在运动时会产生地转偏向力，这个里的方向可以有左手理论提出：古塔变形18地转偏向力示意：即古塔（肯定在北半球）在偏移过程收到地转偏向力的作用，他的偏转方向会偏向运动方向的左边，在整体上表现就是扭曲方向（俯视）呈现逆时针。至此，三种情况个数据都已分析完全，对古塔 1986 到 2011 这段时间内的变化分析可以得到如下结论，1986 到 1996 年间无论是扭曲，倾斜还是弯曲表现得都很微弱，但这以后，三者的变动极大，估计古塔变形19这段时间自然环境发生过一些大的变化或者是遭遇过自然
15、的灾害（初步估计为地震） 。到了 2009 年以后，变化速度急剧减小，特别是扭曲速度，基本减到了 0，倾斜速率也减小到即为微弱，但是古塔的趋势是与地面夹角在逐渐减小。五、附件：源代码：限于篇幅，值附上部分代码1.Tempshape.m %对古塔形状进行大致描绘point=xlsread(C:Usersdsd-dellDesktopC2009.xlsx);grid on;hold on;for i=0 : 1 : 13if i 13x=point(i*8+1 : i*8+8 , 1);y=point(i*8+1 : i*8+8 , 2);z=point(i*8+1 : i*8+8 , 3);pl
16、ot3(x,y,z,red-*);%meshz(x,y,z);plot3(x(1);x(8),y(1);y(8),z(1);z(8),red-*);else if i = 13x=point(i*8+1 : i*8+4 , 1);y=point(i*8+1 : i*8+4 , 2);z=point(i*8+1 : i*8+4 , 3);plot3(x,y,z,red-*);plot3(x(1);x(4),y(1);y(4),z(1);z(4),red-*);endendend%cen=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);plot3(c
17、en(1:13,1),cen(1:13,2),cen(1:13,3),red-*);2.Center.m %几何法补充缺失坐标古塔变形20point=xlsread(C:Usersdsd-dellDesktopC2011nihe.xlsx);syms cen xg ygcen=zeros(13,3);for i=0:1:12x=point(i*8+1 : i*8+8 , 1);y=point(i*8+1 : i*8+8 , 2);z=point(i*8+1 : i*8+8 , 3);xg=0;yg=0;for j=1:1:8xg=xg+x(j);yg=yg+y(j);endcen(i+1,1)
18、=xg/8;cen(i+1,2)=yg/8;cen(i+1,3)=z(1);endxlswrite(C:Usersdsd-dellDesktopC2011center.xlsx,cen);3.QinxieNihe.m %对倾斜程度进行分析%point=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);p1=xlsread(C:Usersdsd-dellDesktopC1986center.xlsx);p2=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);p3=xlsread(C:Usersdsd-de
19、llDesktopC2009center.xlsx);p4=xlsread(C:Usersdsd-dellDesktopC2011center.xlsx);%xi=522.4:0.001:522.6;%yi=interp1(C,xi,linear);%mesh(xi,yi,zi);%plot3(point(1:13,1),point(1:13,2),point(1:13,3),r-*);grid on;hold on;%plot(xi,yi,black-*);%xi=522:0.01:523;%yi=interp1(point(1:13,2),point(1:13,3),xi,linear);
20、%plot(xi,yi,red-*);x=zeros(4,4);p=polyfit(p1(1:13,1),p1(1:13,3),13);f=polyval(p,p1(1:13,1);%plot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(1,1)=1986;x(1,2)=(f(13)-f(1)/(p1(13,1)-p1(1,1);x(1,3)=atan(x(1,2);古塔变形21x(1,4)=cos(x(1,3);p=polyfit(p2(1:13,1),p2(1:13,3),13);f=polyval(p,p2(1:13,1);%pl
21、ot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(2,1)=1996;x(2,2)=(f(13)-f(1)/(p2(13,1)-p2(1,1);x(2,3)=atan(x(2,2);x(2,4)=cos(x(2,3);p=polyfit(p3(1:13,1),p3(1:13,3),13);f=polyval(p,p3(1:13,1);%plot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(3,1)=2009;x(3,2)=(f(13)-f(1)/(p3(13,1)-p3(1,1)
22、;x(3,3)=atan(x(3,2);x(3,4)=cos(x(3,3);p=polyfit(p4(1:13,1),p4(1:13,3),13);f=polyval(p,p4(1:13,1);%plot(point(1:13,1),point(1:13,3),*,point(1:13,1),f,-);x(4,1)=2011;x(4,2)=(f(13)-f(1)/(p4(13,1)-p4(1,1);x(4,3)=atan(x(4,2);x(4,4)=cos(x(4,3);xlswrite(C:Usersdsd-dellDesktopCjiaoduqinxie.xlsx,x);4.NiuquNi
23、he.m %对扭曲程度进行分析cen1=xlsread(C:Usersdsd-dellDesktopC1986center.xlsx);cen2=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);cen3=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);cen4=xlsread(C:Usersdsd-dellDesktopC2011center.xlsx);p1=xlsread(C:Usersdsd-dellDesktopC1986nihe.xlsx);p2=xlsread(C:Usersdsd-de
24、llDesktopC1996nihe.xlsx);p3=xlsread(C:Usersdsd-dellDesktopC2009nihe.xlsx);p4=xlsread(C:Usersdsd-dellDesktopC2011nihe.xlsx);syms c t%古塔变形22for i=1:1:13a=cen1(i,1),cen1(i,2);for j=1:1:8p1(j+(i-1)*8,1)=p1(j+(i-1)*8,1)-a(1);p1(j+(i-1)*8,2)=p1(j+(i-1)*8,2)-a(2);endendfor i=1:1:13a=cen2(i,1),cen2(i,2);for
25、 j=1:1:8p2(j+(i-1)*8,1)=p2(j+(i-1)*8,1)-a(1);p2(j+(i-1)*8,2)=p2(j+(i-1)*8,2)-a(2);endendfor i=1:1:13a=cen3(i,1),cen3(i,2);for j=1:1:8p3(j+(i-1)*8,1)=p3(j+(i-1)*8,1)-a(1);p3(j+(i-1)*8,2)=p3(j+(i-1)*8,2)-a(2);endendfor i=1:1:13a=cen4(i,1),cen4(i,2);for j=1:1:8p4(j+(i-1)*8,1)=p4(j+(i-1)*8,1)-a(1);p4(j+
26、(i-1)*8,2)=p4(j+(i-1)*8,2)-a(2);endend%theta=zeros(4,2);theta(1:1:4,1)=1986,1996,2009,2011;c=0;for i=1:1:12t=0;for j=1:1:8t=t+(p1(j,1)*p1(i*8+j,1)+(p1(j,2)*p1(i*8+j,2)/(p1(j,1)2+p1(j,2)2)0.5)*(p1(i*8+j,1)2+p1(i*8+j,2)2)0.5)/8;end古塔变形23c=c+t/12;endtheta(1,2)=c;c=0;for i=1:1:12t=0;for j=1:1:8t=t+(p2(j
27、,1)*p2(i*8+j,1)+(p2(j,2)*p2(i*8+j,2)/(p2(j,1)2+p2(j,2)2)0.5)*(p2(i*8+j,1)2+p2(i*8+j,2)2)0.5)/8;endc=c+t/12;endtheta(2,2)=c;c=0;for i=1:1:12t=0;for j=1:1:8t=t+(p3(j,1)*p3(i*8+j,1)+(p3(j,2)*p3(i*8+j,2)/(p3(j,1)2+p3(j,2)2)0.5)*(p3(i*8+j,1)2+p3(i*8+j,2)2)0.5)/8;endc=c+t/12;endtheta(3,2)=c;c=0;for i=1:1:
28、12t=0;for j=1:1:8t=t+(p4(j,1)*p4(i*8+j,1)+(p4(j,2)*p4(i*8+j,2)/(p4(j,1)2+p4(j,2)2)0.5)*(p4(i*8+j,1)2+p4(i*8+j,2)2)0.5)/8;endc=c+t/12;endtheta(4,2)=c;%xlswrite(C:Usersdsd-dellDesktopCniuqu1.xlsx,theta);hold on;古塔变形24grid on;p=zeros(13,2);for i=1:1:8for j=0:1:12p(j+1,1:2)=p1(8*j+i,1),p1(8*j+i,2);endpl
29、ot(p(1:13,1),p(1:13,2),black-*);end plot(0,0,black-*);5.WanquNihe.m %对弯曲程度进行拟合p1=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);p2=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);p3=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);p4=xlsread(C:Usersdsd-dellDesktopC2011center.xlsx);p=xlsread(C:Usersd
30、sd-dellDesktopCniuqu1.xlsx);xi=522:0.02:523;yi=interp1(p1(1:13,2),p1(1:13,3),xi,cubic);plot(xi,yi,black-*);%plot(point(1:4,1),point(1:4,3),-*);xi=1986:1:2011;yi=interp1(p(1:4,1),p(1:4,2),xi,cubic);plot(xi,yi,black-*);grid on;hold on;plot(p(1:4,1),p(1:4,2),r-*);6.QingxieDu.m 对倾斜度的计算syms c p1jiao=zero
31、s(4,3);jiao(1:4,1)=1986,1996,2009,2011;c=0;for j=1:1:4if j=1p1=xlsread(C:Usersdsd-dellDesktopC1986center.xlsx);elseif j=2p1=xlsread(C:Usersdsd-dellDesktopC1996center.xlsx);elseif j=3p1=xlsread(C:Usersdsd-dellDesktopC2009center.xlsx);elsep1=xlsread(C:Usersdsd-dellDesktopC2011center.xlsx);endfor i=2:1
32、:12古塔变形25x0=p1(i-1,1);x1=p1(i,1);x2=p1(i+1,1);y0=p1(i-1,2);y1=p1(i,2);y2=p1(i+1,2);z0=p1(i-1,3);z1=p1(i,3);z2=p1(i+1,3);t1=x0-x1,y0-y1,z0-z1;t2=x2-x1,y2-y1,z2-z1;co=(t1(1)*t2(1)+t1(2)*t2(2)+t1(1)*t2(2)/(t1(1)2+t1(2)2+t1(3)2)0.5*(t2(1)2+t2(2)2+t2(3)2)0.5);co=co/11;c=c+co;endjiao(j,2)=c;endfor i=1:1:4jiao(i,3)=acos(jiao(i,2);endxlswrite(C:Usersdsd-dellDesktopCqingxiedu.xlsx,jiao);grid on;hold on;plot(jiao(1:4,1),jiao(1:4,3),r-*);xi=1986:1:2011;yi=interp1(jiao(1:4,1),jiao(1:4,3),xi,cubic);plot(xi,yi,black-*);
