古塔变形分析.doc

古塔变形的分析摘要：文物部门需要适时的对古塔进行观测，了解各种变形量，以制订必要的保护措施。附件1给出了对古塔的四次测量数据，是每层八个测量点的三个坐标值，但1986年和1996年对塔的第13层的测量点5的数据值是缺失的，为了对塔的变形进行准确的分析研究，我们首先应用MATLAB7.0工具箱的数据拟合功能，取1986年和1996年每层观测点5的测量数据和坐标值，用这些已知数据值对第13层的测量点5的和坐标值进行二次曲线拟合，假设塔的每层高度是相同的，用第13层的z的平均值来预测点5的z坐标，则1986年和1996年第13层点5的坐标分别是x=567.992，y=519.727，z=52.83429和x=567.992，y=519.737，z=52.83对所给测量点坐标应用MATLAB7.0做3维线图，直观看出测量点在每层为均匀分布的8个点，各层为正八边形，且随着楼层的增高每层边长逐渐减小，可以推测该古塔为上小下大的八角塔。进一步通过对测量得到的各层的相邻两点间的坐标值应用空间两点间的距离公式：计算边长，每层边长近似相等，可以进一步看出该古塔是一个八角形塔。对第一问：因为该塔是一个正八边形塔，对于中心坐标我们认为是每层的中心，就是该层的坐标中心位置，我们取各层8个测量点坐标的平均值作为该层的平均值，并计算出了该塔每层的中心坐标。对第二问：以第一层中心坐标为起始点向上一层的中心坐标做向量计算各向量与X 、Y 、Z 轴的夹角来判断古塔的倾斜、弯曲、扭曲，根据4次各层中心向量与X、Y 、Z 轴的夹角的变化我们得出1986和1996年塔的倾斜、弯曲、扭曲等变形情况不明显，2009和2011年情况基本相同，但是在第5层到第7层发生扭曲，第8到11层发生弯曲、扭曲，第12到13层发生倾斜，最后倾向于X轴。对于第三问：由第二问的分析情况可以得出该塔发生扭曲、弯曲、倾斜等情况，所以该塔会随着时间的推移第12到13层会发生倾倒。关键词：曲线拟合；空间两点间距离公式；空间物体几何中心；MATLAB7.0；空间向量方向角一、问题的重述 某古塔已有上千年的历史，是我国的重点保护文物。由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响、古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月、2011年3月对该塔进行了4此观测。文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。 希望根据附件一提供的4次观测数据讨论以下问题：1.给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。2. 分析该塔倾斜、弯曲、扭曲等变形情况。3. 分析该塔的变形趋势。二、模型准备（1）测量缺失值的预测由于1986年和1996年第13层测量点5的测量坐标值缺失，为了方便研究首先对缺失值进行预测。把测量点5作为观测点，取第1层到第12层的观测点5的横坐标x的测量值为观测数据，应用MATLAB7.0工具箱的曲线预测功能，对1986年第13层测量点5缺失值的横坐标进行预测，输出下结果：Linear model Poly2: f(x) = p1*x2 + p2*x + p3Coefficients (with 95% confidence bounds): p1 = -0.005984 (-0.006677, -0.00529) p2 = 0.08906 (0.07979, 0.09832) p3 = 567.8 (567.8, 567.9)Goodness of fit: SSE: 0.001129 R-square: 0.9832 Adjusted R-square: 0.9794 RMSE: 0.0112即预测方程：f(x) = -0.005984*x2 + 0.08906*x + 567.8从而可预测出第13层测量点5的横坐标x=567.992对纵坐标y的预测，输出结果：Linear model Poly2: f(x) = p1*x2 + p2*x + p3Coefficients (with 95% confidence bounds): p1 = 0.005767 (0.002656, 0.008877) p2 = 0.1125 (0.07092, 0.154) p3 = 517.3 (517.2, 517.4)Goodness of fit: SSE: 0.02271 R-square: 0.9955 Adjusted R-square: 0.9945 RMSE: 0.05023可得纵坐标y的预测方程：f(x) = 0.005767*x2 + 0.1125*x + 517.3，从而可预测出第13层测量点5的纵坐标y=519.727图1： 1986年第13层测量点5的横坐标x的拟合曲线 图2:1986年第13层测量点5的纵坐标y的拟合曲线对缺失值竖坐标z的预测，由于z表示层高，则用该层其它7个测量点高度的平均值来确定：z=52.83429即1986年第13层观测点5的坐标为：x=567.992，y=519.727，z=52.83429同样的方法对1996年第13层测量坐标的缺失值进行预测，得到预测曲线方程如下： f(x) = -0.005983*x2 +0.08906*x + 567.8 - x的预测曲线 f(x) = 0.005767*x2 +0.1125*x + 517.3 y的预测曲线则1996年第13层观测点5的坐标：x=567.992，y=519.737，z=52.83图3： 1996年第13层测量点5的横坐标x的拟合曲线图4： 1996年第13层测量点5的纵坐标y的拟合曲线（2）古塔形状的推测为了确定古塔的形状，由附件1提供的4次测量点的坐标数据，用MATLAB7.0做出古塔的3维线图如下：图5图6:由测量点坐标的3维线图，可以认为古塔每层的边长为正八边形，古塔形如上小下大的八角楼。用空间两点的距离公式：计算每层相邻测量点的距离如表1，可以看出测量点的选取是均匀分布的，8个测量点构成了一个正八边形。表1 1986年各楼层测量点相邻距离（单位：m)测量点第1层第2层第3层第4层第5层第6层第7层第8层第9层第10层第11层第12层第13层第14层1点4.20 4.03 3.87 3.75 3.61 3.47 3.27 3.08 2.89 2.70 2.52 2.34 2.17 0.03 2点4.15 3.99 3.82 3.70 3.56 3.43 3.24 3.05 2.87 2.70 2.50 2.30 2.11 0.03 3点4.11 3.96 3.81 3.70 3.57 3.45 3.25 3.06 2.87 2.68 2.51 2.33 2.17 0.01 4点4.22 4.05 3.88 3.75 3.61 3.47 3.26 3.06 2.86 2.66 2.47 2.28 1.98 0.01 5点4.12 3.97 3.83 3.72 3.59 3.48 3.27 3.06 2.86 2.67 2.53 2.39 1.99 6点4.15 4.00 3.85 3.73 3.60 3.47 3.26 3.05 2.84 2.64 2.48 2.31 2.15 7点4.19 4.02 3.85 3.72 3.59 3.45 3.26 3.08 2.89 2.71 2.50 2.30 2.10 8点4.10 3.98 3.87 3.78 3.68 3.59 3.36 3.13 2.91 2.70 2.56 2.42 2.29 平均边长4.15 4.00 3.85 3.73 3.60 3.48 3.27 3.07 2.87 2.68 2.51 2.34 2.12 0.02 由表1看出，随着楼层的增加，边长逐渐减小，可以认为古塔为上小下大的八角形塔。三、模型的基本假设 1.假设不考虑测绘工作人员产生的数据误差； 2.附件1给出每层八个点的坐标，假设该古塔为八边形塔，每次测量位置在八个角点上面； 3.由于古塔有上千年的历史，假设排除地球引力的因素； 4.假设在几千年前塔在完工之后没有任何倾斜、弯曲、扭曲；5.假设本题目所给出的数据都是真实可靠的，不存在虚假数据。四、符号说明：第层中心到第层中心的向量；、：表示向量与三个坐标轴的夹角、：表示向量的三个坐标；、：表示中心坐标；五、模型的建立对于问题1：古塔每层中心坐标的确定用公式： ， ， 对于问题2：分析该塔倾斜、弯曲、扭曲等变形情况。要研究该塔的倾斜、弯曲、扭曲等变形情况，为了简化问题我们把塔的每层的中心坐标作为研究对象，用表示第层的中心坐标（）定义=为第层中心到第层中心的向量，我们可以用它与三个坐标轴的夹角、来表示它的方向， 、与三个坐标轴夹角的的方向余弦： 根据、的变化情况来研究塔的倾斜、弯曲、扭曲等。对于问题3：根据问题二得出的各层中心向量与坐标轴夹角的变化分析该塔的变形趋势。对各层中心坐标数据值进行曲线拟合，由拟合曲线给出该塔的变化趋势。 六、模型的求解问题一的求解：利用 ， ， ， excel的函数命令，计算出四次测量各楼层的中心坐标如下表： 表2 1986年古塔各层中心坐标 单位：m楼层坐标横坐标x纵坐标y竖坐标z1566.6648522.71051.7873752566.7196522.66847.320253566.7735522.627312.755254566.8161522.594417.078255566.8621522.559121.72056566.9084522.524426.235137566.9468522.508129.836888566.9843522.492433.350889567.0218522.476436.8548810567.0569522.462440.1721311567.1045522.42344.4408812567.1518522.383648.7118813567.1984522.364952.83429塔顶567.2473522.243855.12325表3 1996年古塔各层中心坐标 单位： m楼层坐标横坐标x纵坐标y竖坐标z1566.665522.71021.7832566.7205522.66747.3146253566.7751522.625612.750754566.8183522.592217.075135566.8649522.556321.7166566.9118522.52126.22957566.9506522.504229.832258566.9884522.488133.345389567.0265522.471436.8482510567.062522.457240.1676311567.1102522.417344.4353812567.1578522.377548.7073813567.2038522.359452.83塔顶567.2591522.236755.11975表4 2009古塔各层中心坐标 单位：m楼层坐标横坐标x纵坐标y竖坐标z1566.7268522.70151.76452566.764522.66937.3093566.8001522.638412.732254566.8293522.613217.069755566.8604522.586621.709386566.9471522.534226.2117566.9792522.512329.824638567.0305522.479733.339889567.0816522.446636.8437510567.137522.393740.1611311567.1799522.354744.4326312567.2225522.31649.1618913567.2712522.271553.07089塔顶567.336522.214855.091表5 2011年古塔各层中心坐标 单位：m楼层坐标横坐标x纵坐标y竖坐标z1566.727522.70141.7582566.7642522.6697.29053566.8004522.638712.726884566.8297522.612717.0525566.861522.58621.703886566.9478522.533526.20457566.98522.511529.8178567.0313522.478833.336639567.0825522.445736.8222510567.1381522.392640.1441311567.181522.353544.4248812567.2238522.314748.6838813567.2725522.270152.81313塔尖567.3375522.213555.087问题二求解： 用a,b,c表示中心的三个坐标值，dx、dy、dz分别表示向量与x轴、y轴、z轴的夹角（弧度），tx、ty、tz分别表示向量与x轴、y轴、z轴方向夹角（度），计算用MATLAB命令实现：a= ;b= ;c= ;dx=acos(a/(sqrt(a2+b2+c2)dy=acos(b/(sqrt(a2+b2+c2)dz=acos(c/(sqrt(a2+b2+c2)tx=dx*180/pity=dy*180/pitz=dz*180/pi表格6 1986年1996年2001年2011年各层方向角 时间方向第1层第2层第3层第4层第5层第6层第7层第8层第9层第10层第11层第12层1986年x-轴tx89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.3 y-轴ty90.4 90.4 90.4 90.4 90.4 90.3 90.3 90.3 90.2 90.5 90.5 90.3 z-轴tz0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 0.7 1996年x-轴tx89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 89.4 y-轴ty90.4 90.4 90.4 90.4 90.4 90.3 90.3 90.3 90.2 90.5 90.5 90.3 z-轴tz0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.8 0.7 2009年x-轴tx89.6 89.6 89.6 89.6 88.9 89.5 89.2 89.2 89.1 89.4 89.5 90.3 y-轴ty90.3 90.3 90.3 90.3 90.7 90.3 90.5 90.5 90.9 90.5 90.5 89.3 z-轴tz0.5 0.5 0.5 0.5 1.3 0.6 1.0 1.0 1.3 0.8 0.7 0.8 2011年x-轴tx89.6 89.6 89.6 89.6 88.9 89.5 89.2 89.2 89.0 89.4 89.4 89.3 y-轴ty90.3 90.3 90.3 90.3 90.7 90.3 90.5 90.5 90.9 90.5 90.5 90.6 z-轴tz0.5 0.5 0.5 0.5 1.3 0.6 1.0 1.0 1.3 0.8 0.8 0.9 对四次古塔每层间的中心向量与X、Y、Z轴的夹角做折线图得：图7说明：1986年与1996年古塔中心向量与x轴夹角不是很明显，而2009年和2011年5-6层中心点向量值与x轴夹角明显增大，7-11层中心点向量与x轴也有明显变化，第12-13中心点向量与x轴夹角更加明显。图8 说明：在四次观测中，在第一到第六层，古塔中心与y轴的夹角大致相同，则古塔在第1到第6层在y轴方向没有扭曲变化，而从第6层以后与y轴的夹角有了明显的波动，到第12层2009年的变动比较突出。图9说明：从图中可知1986年和1996年在6到11层有扭曲情况，12到13层发生倾斜。而2009年和2011年从第5层开始变化比较大，扭曲程度也很明显。图10 从这幅图可以看出古塔1986年和1996年中心位置基本一致，2009年和2011年的中心位置基本一致，然而1996年和2009年古塔的中心位置差别有点大，说明倾斜度大了，这点可以推导出1996年到2009年之间也许发生过地震等一些自然灾害，导致了古塔的倾斜增大。问题三求解： 根据问题一和二的分析说明该古塔的第12到第13层随着时间的增长，倾斜度会越来越大，最后很有可能会导致倾倒。七、模型的评价和改进方向模型的优点：1：利用MATLAB7.0工具箱的曲线预测功能，对1986年第13层测量点5缺失值的横坐标进行预测。2：利用Excel的函数功能计算,方便快捷。3：充分利用了MATLAB7.0 的绘图与计算功能。 模型的不足：1：在考虑问题方面不是很全面。2：为了方便问题的解决没有考虑塔尖。参考文献 1姜启源,谢金星,叶俊，数学模型（第三版），北京：高等教育出版社，20032于洪彦.Excel统计分析与决策M.北京:高等教育出版社.2006年4月；3胡守信,李柏年.基于MATLAB的数学实验M.北京:科学出版社.2004年6月。附录：1.图5：1986年古塔测量坐标的3维线图绘制的MATLAB程序：x=565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5;y=528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356;z=1.792 1.818 1.783 1.769 1.772 1.77 1.794 1.801;x1=565.48 562.238 561.663 564.001 567.995 571.165 571.801 569.414;y1=527.764 525.364 521.42 518.226 517.563 519.961 523.908 527.141;z1=7.326 7.351 7.314 7.301 7.306 7.304 7.324 7.336;x2=565.506 562.415 561.931 564.216 568.048 571.076 571.666 569.33;y2=527.52 525.188 521.394 518.343 517.716 520.063 523.864 526.93;z2=12.761 12.786 12.749 12.736 12.741 12.74 12.758 12.771;x3=565.526 562.555 562.144 564.387 568.091 571.005 571.558 569.263;y3=527.327 525.047 521.373 518.435 517.838 520.144 523.829 526.762;z3=17.084 17.109 17.072 17.059 17.064 17.063 17.081 17.094;x4=565.548 562.706 562.373 564.571 568.136 570.929 571.443 569.191;y4=527.119 524.896 521.351 518.534 517.969 520.232 523.791 526.581;z4=21.726 21.751 21.714 21.701 21.705 21.708 21.723 21.736;x5=565.57 562.854 562.6 564.752 568.18 570.857 571.333 569.121;y5=526.915 524.748 521.329 518.632 518.095 520.315 523.755 526.406;z5=26.267 26.309 26.308 26.264 26.189 26.136 26.164 26.244;x6=565.671 563.132 562.883 564.949 568.172 570.679 571.094 568.994;y6=526.652 524.585 521.356 518.846 518.346 520.441 523.672 526.167;z6=29.869 29.911 29.91 29.866 29.791 29.737 29.765 29.846;x7=565.77 563.403 563.158 565.141 568.164 570.506 570.862 568.87;y7=526.397 524.427 521.382 519.055 518.59 520.564 523.591 525.933;z7=33.383 33.425 33.424 33.38 33.305 33.251 33.279 33.36;x8=565.868 563.674 563.433 565.333 568.156 570.333 570.63 568.747;y8=526.141 524.268 521.408 519.263 518.834 520.686 523.51 525.701;z8=36.887 36.929 36.928 36.884 36.809 36.755 36.783 36.864;x9=565.961 563.927 563.693 565.516 568.148 570.171 570.408 568.631;y9=525.9 524.12 521.433 519.462 519.068 520.801 523.433 525.482;z9=40.201 40.214 40.244 40.223 40.171 40.038 40.129 40.157;x10=566.078 564.193 563.958 565.649 568.094 570.013 570.236 568.615;y10=525.628 523.95 521.463 519.607 519.242 520.885 523.35 525.259;z10=44.472 44.485 44.505 44.486 44.442 44.309 44.4 44.428;x11=566.195 564.459 564.224 565.782 568.039 569.854 570.063 568.598;y11=525.355 523.78 521.492 519.753 519.415 520.969 523.268 525.037;z11=48.743 48.756 48.776 48.757 48.713 48.58 48.671 48.699;x12=566.308 564.716 564.481 565.91 567.8055 569.701 569.897 568.582;y12=525.092 523.616 521.521 519.893 520.4715 521.05 523.188 524.822;z12=52.866 52.878 52.897 52.88 52.834 52.703 52.794 52.822;plot3(x,y,z,r)hold onplot3(x1,y2,z3,r)hold onplot3(x2,y2,z2,r)hold onplot3(x3,y3,z3,r)hold onplot3(x4,y4,z4,r)hold onplot3(x5,y5,z5,r)hold onplot3(x6,y6,z6,r)hold onplot3(x7,y7,z7,r)hold onplot3(x8,y8,z8,r)hold onplot3(x9,y9,z9,r)hold onplot3(x10,y10,z10,r)hold onplot3(x11,y11,z11,r)hold onplot3(x12,y12,z12,r)grid on2. 楼层边1986年1996年2009年2011年11-24.204.204.124.122-34.154.154.154.153-44.114.114.194.194-54.224.224.164.165-64.124.124.194.196-74.154.154.164.167-84.194.194.124.128-14.104.104.154.1521-24.034.033.963.962-33.993.993.973.973-43.963.964.024.024-54.054.053.993.995-63.973.974.044.046-74.004.003.993.997-84.024.023.953.958-13.983.984.014.0131-23.873.873.813.812-33.823.833.803.803-43.813.813.863.864-53.883.883.823.825-63.833.833.903.906-73.853.853.833.837-83.853.853.793.798-13.873.873.873.8741-23.753.753.693.692-33.703.703.663.663-43.703.703.733.734-53.753.753.693.695-63.723.723.793.796-73.733.733.703.707-83.723.723.663.668-13.783.783.763.7651-23.613.613.563.562-33.563.563.523.523-43.573.573.593.594-53.613.613.543.555-63.593.593.673.676-73.603.603.563.567-83.593.593.523.528-13.683.683.643.6461-23.473.473.243.252-33.433.433.243.243-43.453.453.303.304-53.473.473.263.265-63.483.483.363.366-73.473.473.303.307-83.453.453.253.258-13.593.593.343.3471-23.273.273.133.132-33.243.243.143.143-43.253.253.203.204-53.263.263.163.165-63.273.273.243.246-73.263.263.203.207-83.263.263.153.158-13.363.363.233.2381-23.083.082.952.952-33.053.052.982.983-43.063.063.033.034-53.063.062.992.995-63.063.063.053.056-73.053.053.053.057-83.083.082.982.988-13.133.133.053.0591-22.892.882.772.772-32.872.872.822.823-42.872.872.862.864-52.862.862.832.835-62.862.862.862.876-72.842.842.892.897-82.892.892.822.828-12.912.912.874.42101-22.702.702.642.642-32.702.702.672.673-42.682.682.742.744-52.662.662.692.695-62.672.672.712.716-72.642.642.772.777-82.712.712.682.688-12.702.702.712.71111-22.522.522.482.482-32.502.502.482.483-42.512.512.582.584-52.472.472.502.505-62.532.532.502.506-72.482.482.612.617-82.502.502.512.518-12.562.562.512.51121-22.342.342.322.322-32.302.302.292.293-42.332.342.412.414-52.282.282.322.325-62.392.392.302.306-72.312.312.452.457-82.302.302.342.348-12.422.422.302.30131-22.172.172.172.172-32.112.112.102.103-42.172.172.262.264-52.