当前位置:首页 >> 学科竞赛 >>

古塔的变(胡、罗、钟)形








我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参 赛规则》 (以下简称为 “竞赛章程和参赛规则” , 可从全国大学生数学建模竞赛网站下载) 。 我们完全明白,在竞赛开始后参赛队员不能以任何方式(包括电话、电子邮件、网 上咨询等)与队外的任何人(包括指导教师)研究、讨论与赛题有关的

问题。 我们知道,抄袭别人的成果是违反竞赛章程和参赛规则的,如果引用别人的成果或 其他公开的资料(包括网上查到的资料) ,必须按照规定的参考文献的表述方式在正文 引用处和参考文献中明确列出。 我们郑重承诺,严格遵守竞赛章程和参赛规则,以保证竞赛的公正、公平性。如有 违反竞赛章程和参赛规则的行为,我们将受到严肃处理。 我们授权全国大学生数学建模竞赛组委会,可将我们的论文以任何形式进行公开展 示(包括进行网上公示,在书籍、期刊和其他媒体进行正式或非正式发表等) 。 我们参赛选择的题号是(从 A/B/C/D 中选择一项填写) : 我们的参赛报名号为(如果赛区设置报名号的话) : 所属学校(请填写完整的全名) : 参赛队员 (打印并签名) :1. 2. 3. 指导教师或指导教师组负责人 中山职业技术学院 胡学文 罗屯林 钟贵辉 (打印并签名): 杨晶晶 C 5791

(论文纸质版与电子版中的以上信息必须一致,只是电子版中无需签名。以上内容 请仔细核对,提交后将不再允许做任何修改。如填写错误,论文可能被取消评奖资格。) 日期: 2013 年 09 月 16 日

赛区评阅编号(由赛区组委会评阅前进行编号):

2013 高教社杯全国大学生数学建模竞赛

编 号 专 用 页

赛区评阅编号(由赛区组委会评阅前进行编号):

赛区评阅记录(可供赛区评阅时使用): 评 阅 人 评 分 备 注

全国统一编号(由赛区组委会送交全国前编号):

全国评阅编号(由全国组委会评阅前进行编号):

古塔的变形
摘要
本文研究了古塔变形的问题。由于长时间承受自重、气温、风力等各种作用,偶然 还要受地震、飓风的影响,古塔会产生各种变形,诸如倾斜、弯曲、扭曲等。为保护古 塔,文物部门需适时对古塔进行观测,了解各种变形量,以制定必要的保护措施。 为解决问题 1,我们先讨论研究了中古塔各层中心位置的通用方法。为此我们建立 了模型一和优化模型二。通过这四年的数据拟合得出未变形前的古塔的中心轴方程(空 间直线) , 利用中心轴求出层高 H 为 4.254491m, 进而可求出未变形前的古塔各层中心坐 标(见正文表 8) , ,然后再用类似的方法,得用层高可求出测量的这四年各次测量的古 塔各层中心坐标(见正文表 3) 。 为解决问题 2,我们首先分析了该塔的倾斜情况并建立了模型三。用未变形前的古 塔中心轴与这四年各年的中心轴(有变形)之间的夹角来表示各年的塔身倾斜情况。利 用三角函数知识,可求出四年的空间倾斜角度的正切值分别为 0.0119471,0.01211, 0.014978,0.015015。 接着我们分析了该塔的弯曲情况并建立了模型四以及优化模型五。 通过数据拟合以 及曲率公式曲率 k= ,得出 4 个年份的曲率依次为 0.000424,0.00046,0.000143, 0.000145,可以看到古塔弯曲程度先增加,接着有所减少,然后再增加; 对于古塔的扭曲,我们建立了模型六,通过三次多项式拟合以及挠率公式
i? y'' z ''' ? y''' z '' , 得 出 测 量 的 这 四 年 挠 率 依 次 为 0.0000069221 , 2 2 y'' ? z'' +( y' z''- y''z' )
2

0.00000000000002988 ,0.000282427666666667 ,0.000360233333333333 。可见古塔的 扭曲程度先减少,接着再增加; 对于问题 3, 利用塔的倾斜度变化来比较得出塔呈现出递增的倾斜趋势, 倾斜度分 别由 tant1986=0.11947 上升至 tant19966=0.012114 再上升到 tant2009=0.014978 接着 升到 tant2011=0.015015;接着把各个年份的曲率的进行比较、得知塔的弯曲程度先由 K1986=0.000424 上 升 至 K1996=0.00046 下 降 至 K2009=0.000143 再 上 升 为 K2011=0.000145 因此 最后把各年平均挠率进行比较,分析出古塔的扭曲程度先由 i1986=0.0000069221 减 少 到 i1996=0.00000000000002988 再 下 降 为 i2009=0.000282427666666667 再上升为 i2011=0.000360233333333333,因此古塔呈现出 倾斜度,弯曲程度和扭曲程度不规则变化的趋势。我们还画出来图形和表格,来分析这 种不规则变化。

关键词:

数据拟合 斜率

曲率

挠率

MATLAB 软件

1

一、 问题重述
古塔,是中国五千年文明史的载体之一,古塔为祖国城市山林增光添彩,塔被佛 教界人士尊为佛塔。矗立在大江南北的古塔,被誉为中国古代杰出的高层建筑,但由于 长时间受各种因素的影响,古塔会产生各种变形,诸如倾斜、弯曲、扭曲等。为保护古 塔,文物部门需适时对古塔进行观测,了解塔变形的情况,以制定必要的保护措施。现 管理部门委托测绘公司先后于 1986 年 7 月、1996 年 8 月、2009 年 3 月和 2011 年 3 月 对该塔进行了 4 次观测。经过 4 个次的测量,为了解古塔的变形情况,首先需要测量出 古塔各个塔层的中心坐标来判断塔中心轴线(各个塔层中心坐标的连线)是否发生偏移 并且为了方便于以后对古塔的测量,故需要寻求出一个测量的通用方法,以便于以后能 够更好的测出。在通过 4 个年份 4 次分别对古塔每个塔层的多点观测数据,分析出该古 塔在 1986 年到 2011 年间倾斜、弯曲、扭曲的变化情况和各个年之间的趋势分析。

二、 问题假设
1. 古塔各个塔层和塔尖的高度都相同。 2. 古塔每层有 8 个观测点且 8 个点在各层内离层底相对高度相同(注:1986 年和 1996 年 13 层数据中有一组为空白,故把空白数据忽略) 。 3. 每个年份对每一层的观测点都固定 4. 每层观测点在古塔的中心轴周围。

三、 符号说明

( x, y, z ) :表示三维空间中的点的坐标
Mi

( x, y, z ) :表示第 i 层的中心坐标, i =1,2,??,13。

h:表示塔的高度 H:表示古塔各层的层高 N:表示塔的层数 tan t:表示古塔的倾斜程度(斜率) i:表示古塔的扭曲程度(挠率) k:表示古塔的弯曲程度(曲率)

Di : 塔层与塔层间的距离

四、 问题分析
4.1 对于问题一: 先把塔层定义为一层,然后把每年这一层的八个观测进行算术平均处理,接着再 对剩下的塔层进行相同处理。 随后再对每一层 4 个年份进行算术平均, 得出 13 个塔层 4 个年份算术平均的 X Y Z 数据,然后对每个项的 13 组数据进行数据拟合,最终拟合 出一个相似于塔的中心线的方程组, 然后把两次算术平均处理得出的 Z 代入进入方程组 得出 X 拟合,再把 X 拟合代入方程组里得出 Y 拟合如此类推得出 14 组数据,然而每个年份对 每一层的观测点都固定,因此相同位置两两塔层间的距离就等于层的高度 H。但由于两
2

两塔层间的距离不同,因此利用两次算术平均的 Z 把 13 个塔层(塔尖已除外)间的距 离(利用两点间的距离公式)得出 12 组塔层间的距离数据,然后把这 12 组塔层间的距 离数据进行算术平均处理得出平均的层高 H。然而第一层的中心垂直距离为,第二层中
3 1 H N?H ? H 2 。然而各个塔层的中心垂直 心垂直距离为 2 ,所以第 N 层中心垂直距离为

坐标得出来后 ,再次代入拟合后的方程组中 得出各个塔层和塔尖的中心坐标。 4.2 对于问题二: 4.2.1 分析塔的倾斜情况: 通过 4 个年份分别 4 次对古塔每个塔层的多点观测数据 (每个 年份对每一层的观测点都固定) ,,先将附件中的数据进行算术平均处理,然后把算术 平均处理后的数据进行一次多项式拟合后得出出 4 个年份的拟合方程组, 接着把塔底中 心点 O(XOK, YOK,0)的二维平面与垂直于 O 点的轴线 LK 构建空间直角坐标系,然后在塔 的中心轴线上随机取一个点为,然后再过 M 点做垂直于轴线 LK 的垂线于点 WK(XWK,YWK, ZWK) , 那么 WK 与 M 构成一个平面,通过勾股定理并且进行数角转换得出倾斜度 t。 4.2.2 分析塔的弯曲情况:从题目给定的附件中的数据进行 2 次多项式拟合,然后在 2 次多项式中选取不同高度的低(15m)、中(30m)、高(45m)三个层次来比较古塔 4 个年份古塔的弯曲程度,在判别弯曲程度的基础上我们采用曲率来分析古塔的弯曲程 度,曲率表示拟合二次曲线偏离拟合中心直线的程度,曲率越大,表示曲线的弯曲程度 越大 4.2.3 对于古塔的扭曲程度:先通过对于挠率,我们通过把算术平均处理的数据首先进行 三次多项式拟合得出 4 个年份的方程组, 然后分别对三次多项式 4 个年份的方程组进行 1 阶,2 阶和 3 阶求道,得多个求导公式,然后 分别对 Z=高(45m)/中(30m)/低(15m) 三个层次进行挠率求解代入 i=|
y'' z ' ' '? y'' ' z ' ' | 当中得出各个年份的三 y'' ^2 ? z' ' ^2 + ( y' z' '- y'' z' )^2

个层次挠率,然后再对各个年份的挠率进行算术平均,然后比较各个年份的平均挠率, 分析出古塔扭曲情况。 4.3 对于问题三: 在问题二的模型基础上,进行了数据分析。分别对古塔四个年份的三维平面的倾 斜角进行比较,估计出倾斜的趋势,古塔的各年曲率,估计出 4 个年份内古塔的弯曲走 势,最后对古塔的各年挠率进行比较,估计出塔 4 个年份内的扭曲状况。

五、 模型的建立与求解
5.1 各层的中心坐标 通过题目给出的附件 1 可知 1986 年 7 月、1996 年 8 月、2009 年 3 月和 2011 年 3 月这 4 次对该塔观测的数据。然而如何的求出各层的中心坐标,找出这个通式的方法, 我们建立模型一。并对每一年各个塔层的观测点进行算术平均处理(注:赛题附表当中 1986 年和 1996 年第 13 层只有 7 个观测点, 故在计算平均数时, 只除以七, 而不是八) , 第 14 层指的是塔尖部分。
3

以 1986 年第一层为例: 565.454 ? 562.058561 .39 ? 563.782 ? 567.941571 .25 ? 5571.938 ? 569.5 =566.66475 8 类似地可以用同样的方法得到其它数据;这些数据汇总后,可得到如下表 1。 表 1: 1986 年数据 1996 年数据 坐标 坐标 层 层 x(单位:m) y (单位: m) z (单位: m) x(单位:m) y(单位:m) Z (单位: m) 1 566.66475 522.7105 1.787375 1 566.66495 522.7102 1.783 2 566.719625 522.668375 7.32025 2 566.720525 522.6674375 7.314625 3 566.7735 522.62725 12.75525 3 566.7751 522.6256375 12.75075 4 566.816125 522.594375 17.07825 4 566.8183125 522.592175 17.07513 5 566.862125 522.559125 21.7205 5 566.8649125 522.556275 21.716 6 566.908375 522.524375 26.235125 6 566.911775 522.5209625 26.2295 7 566.94675 522.508125 29.836875 7 566.95055 522.504225 29.83225 8 566.98425 522.492375 33.350875 8 566.9884375 522.488075 33.34538 9 567.02175 522.476375 36.854875 9 567.026525 522.471375 36.84825 10 567.056875 522.462375 40.172125 10 567.0619875 522.457175 40.16763 11 567.1045 522.423 44.440875 11 567.1102375 522.4172625 44.43538 12 567.15175 522.383625 48.711875 12 567.157775 522.3775 48.70738 13 567.085 522.740286 52.8342857 13 567.0911571 522.7340429 52.83 14 567.24725 522.24375 55.12325 14 567.25435 522.23665 55.11975 2009 年数据 2011 年数据 坐标 坐标 层 层 x(单位:m) y (单位: m) z (单位: m) x(单位:m) y(单位:m) z (单位: m) 1 566.7268125 522.70146 1.7645 1 566.72695 522.7013625 1.76325 2 566.764 522.669275 7.309 2 566.7642 522.669025 7.2905 3 566.800075 522.6384 12.73225 3 566.8004375 522.6387442 12.72688 4 566.829313 522.613188 17.06975 4 566.8297125 522.612725 17.052 5 566.86035 522.586638 21.709375 5 566.8609625 522.586025 21.70388 6 566.947138 522.5342 26.211 6 566.9478375 522.53345 26.2045 7 566.979238 522.512325 29.824625 7 566.9800375 522.511525 29.817 8 567.0305 522.4797 33.339875 8 567.0313125 522.4788125 33.33663 9 567.081575 522.4466 36.84375 9 567.0825375 522.4456875 36.82225 10 567.136975 522.393688 40.161125 10 567.1380625 522.3925875 40.14413 11 567.179863 522.354663 44.432625 11 567.1809625 522.3535375 44.42488 12 567.222525 522.316013 48.69975 12 567.223825 522.3147125 48.68388 13 567.271213 522.271525 52.818375 13 567.2725125 522.2700875 52.81313 14 567.336 522.2148 55.091 14 567.3375 522.2135 55.087 再用拟合数据的方法可以拟合得出各年份的拟合中心位置的直线 P1、P2、P3 和 P4 (详细 matlab 编码见附录 4) 直线 P1:

4

?Y 1 = -0.5934204 72662X 1 + 858.971114 09774 ? ? Z1 ? 97.3298424 36X 1 ? - 55150.8630 515990

直线 P2:
?Y2 ? -0.5975321 41319X 2 ? 861.300620 492342 ? ? Z 2 ? 96.1669058 154X 2 - 54491.9082 965752

直线 P3:
?Y3 ? -0.7718499 50362 X 3 ? 960.128944 594815 ? ? Z 3 ? 84.3387697 304 X 3 - 47790.5081 085657

直线 P4:
?Y4 ? - 0.77247188 1683X 4 ? 960.481428 041587 ? ? Z 4 ? 84.1573793 412X 4 ? - 47687.7341 371109

其中把算术平均处理得出的 Z1 代进入通过 求出的拟合方程 Z1 =90.8972421008*X1+51506.6709048934 得出 X1 后,再把拟合出的 X1 的值代入 Y1=-0.669305985852 X1 +901.989584064071 得出的 Y1 值,如将第 1 层的中 55150.8630 515990 ? 1.787375 心点的 Z1 坐标带进②式求出 X1= =566.65714282, 即第 1 层 97.3298424 36 中心点的 X 坐标为 566.65714282 ;再把第 1 层中心点的 X 坐标带进①式求出 Y 拟合 =-0.593420472662*566.65714282 +858.97111409774 =522.705164563963,即第 1 层中 心点的 Y 坐标为 522.705164563963, 类似地可以用同样的方法得到其它数据;这些数 据汇总后,可得到表 2 如下。 表 2: 1986 年中心点 1996 年中心点 坐标 坐标 层 层 x(单位:m) y (单位: m) z (单位: m) x(单位:m) y(单位:m) Z (单位: m) 1 566.657142 522.70516 1.787375 1 566.6574 493.1002 1.783 2 566.71399 522.67143 7.32025 2 566.7149 493.0039 7.314625 3 566.76983 522.63829 12.75525 3 566.7715 492.9093 12.75075 4 566.8142 522.6119 17.07825 4 566.8164 492.834 17.07513 5 566.8619 522.5836 21.7205 5 566.8647 492.7533 21.716 6 566.90833 522.55611 26.235125 6 566.9116 492.6747 26.2295 7 566.94533 522.53415 29.836875 7 566.9491 492.612 29.83225 8 566.98144 522.51272 33.350875 8 566.9856 492.5509 33.34538 9 567.01744 522.49136 36.854875 9 567.0221 492.4899 36.84825 10 567.05152 522.47113 40.172125 10 567.0566 492.4322 40.16763 11 567.09538 522.44511 44.440875 11 567.101 492.3579 44.43538 12 567.13926 522.41907 48.711875 12 567.1454 492.2836 48.70738 13 567.18162 522.39393 52.834286 13 567.1882 492.2118 52.83 14 567.20513 522.37998 55.12325 14 567.2121 492.172 55.11975 2009 年中心点 2011 年中心点 层 坐标 层 坐标
5

x(单位:m) y (单位: m) z (单位: m) x(单位:m) y(单位:m) z (单位: m) 1 566.6705 522.7444 1.7645 1 566.6704 522.7445 1.76325 2 566.7361 522.6937 7.309 2 566.7361 522.6937 7.2905 3 566.8005 522.644 12.73225 3 566.8007 522.6438 12.72688 4 566.8518 522.6044 17.06975 4 566.8521 522.6041 17.052 5 566.9068 522.5619 21.70938 5 566.9074 522.5614 21.70388 6 566.9603 522.5206 26.211 6 566.9608 522.5201 26.2045 7 567.0031 522.4877 29.82463 7 567.0038 522.487 29.817 8 567.0447 522.4555 33.33988 8 567.0456 522.4547 33.33663 9 567.0862 522.4235 36.84375 9 567.087 522.4227 36.82225 10 567.1256 522.3931 40.16113 10 567.1265 522.3922 40.14413 11 567.1762 522.354 44.43263 11 567.1773 522.3529 44.42488 12 567.2269 522.3149 48.69975 12 567.228 522.3138 48.68388 13 567.2757 522.2772 52.81838 13 567.277 522.2759 52.81313 14 567.3029 522.2562 55.091 14 567.304 522.255 55.087 即 1986 年、1996、2009、2011 各层中心坐标分别为表 3 如下。 表 3: 1986 年各个塔层中心坐标 1996 年各个塔层中心坐标 塔 层 数 塔层数 塔层中心坐标 1 M1(566.6571 ,522.7052 ,1.7874) M1(566.6574 ,493.1002,1.7830) 2 M2(566.7140 ,522.6714 ,7.3203) M2(566.7149 ,493.0039,7.3146) 3 4 5 6 7 8 9 10 11 12 13 14 塔 层 数 1 2 M3(566.7698 ,522.6383 ,12.7553) M4(566.8142 ,522.6119 ,17.0783) M5(566.8619 ,522.5836 ,21.7205) M6(566.9083 ,522.5561 ,26.2351) M7(566.9453 ,522.5341 ,29.8369) M8(566.9814 ,522.5127 ,33.3509) M9(567.0174 ,522.4914 ,36.8549) M10(567.0515 ,522.4711 ,40.1721) M11(567.0954 ,522.4451 ,44.4409) M12(567.1393 ,522.4191 ,48.7119) M13(567.1816 ,522.3939 ,52.8343) M14(567.2051 ,522.3800 ,55.1233) 2009 年各个塔层中心坐标 M3(566.7715 ,492.9093,12.7508) M4(566.8164 ,492.8340,17.0751) M5(566.8647 ,492.7533 ,21.7160) M6(566.9116 ,492.6747,26.2295) M7(566.9491 ,492.6120 ,29.8323) M8(566.9856,492.5509 ,33.3454) M9(567.0221 ,492.4899 ,36.8483) M10(567.0566 ,492.4322,40.1676) M11(567.1010 ,492.3579,44.4354) M12(567.1454,492.2836 ,48.7074) M13(567.1882 ,492.2118 ,52.8300) M14(567.2121 ,492.1720 ,55.1198) 2011 年各个塔层中心坐标

塔层数 M1(566.6705 ,522.7444 ,1.7645) M2(566.7361 ,522.6937 ,7.3090)
6

塔层中心坐标 M1(566.6704 ,522.7445 ,1.7633) M2(566.7361 ,522.6937 ,7.2905)

3 4 5 6 7 8 9 10 11 12 13 14

M3(566.8005 ,522.6440 ,12.7323) M4(566.8518 ,522.6044 ,17.0698) M5(566.9068 ,522.5619 ,21.7094) M6(566.9603 ,522.5206 ,26.2110) M7(567.0031 ,522.4877 ,29.8246) M8(567.0447 ,522.4555 ,33.3399) M9(567.0862 ,522.4235 ,36.8438) M10(567.1256 ,522.3931 ,40.1611) M11(567.1762 ,522.3540 ,44.4326) M12(567.2269 ,522.3149 ,48.6998) M13(567.2757 ,522.2772 ,52.8184) M14(567.3029 ,522.2562 ,55.0910)

M3(566.8007 ,522.6438 ,12.7269) M4(566.8521 ,522.6041 ,17.0520) M5(566.9074 ,522.5614 ,21.7039) M6(566.9608 ,522.5201 ,26.2045) M7(567.0038 ,522.4870 ,29.8170) M8(567.0456 ,522.4547 ,33.3366) M9(567.0870 ,522.4227 ,36.8223) M10(567.1265 ,522.3922 ,40.1441) M11(567.1773 ,522.3529 ,44.4249) M12(567.2280 ,522.3138 ,48.6839) M13(567.2770 ,522.2759 ,52.8131) M14(567.3040 ,522.2550 ,55.0870)

为了更加精确, 减少误差故对表一的 4 个年份的观测数据再进行一次算术平均处理 (例如: 坐标 x 在第一层的算法是上表一 4 个年份在第 1 层 x 坐标相加后再除以 4 即是 (566.66475 ? 566.66495 ? 566.726812 5 ? 566.72695 ) =566.6959 ),类似计算可以算出以 4 下表 4 的数据。 表 4: 坐标 层 x (单位: m) y (单位: m) z (单位: m) 1 566.6959 522.7059 1.774531 2 566.7421 522.6685 7.308594 3 566.7873 522.6325 12.74128 4 566.8234 522.6031 17.06878 5 566.8621 522.572 21.71244 6 566.9288 522.5283 26.22003 7 566.9641 522.5091 29.82769 8 567.0086 522.4847 33.34319 9 567.0531 522.46 36.84228 10 567.0985 522.4265 40.16125 11 567.1439 522.3871 44.43344 12 567.189 522.348 48.70072 13 567.18 522.504 52.82395 14 567.2938 522.2272 55.10525 对数据进行拟合得出方程组

85852X ? 901.989584 064071 ?Y ? -0.6693059 ? 008X - 51506.6709 048934 ? Z ? 90.8972421

(详细 matlab 编码见附录 1)

其中把两次算术平均处理得出的 Z 两次算术平均代进入通过 求出的拟合方程 Z 两次算术平均=90.8972421008*X 拟合+ 51506.6709048934 得出 X 拟合后, 再把
7

拟合出的 X 拟合的值代入 Y 拟合=-0.669305985852 X 拟合 +901.989584064071 得出的 Y 拟合 值,如将第 1 层的中心点的 Z 两次算术平均坐标带进②式求出 X 拟合 51506.6709 048934 ? 1.774531 = =566.6668,即第 1 层中心点的 X 坐标为 566.6668 ;再 90.8972421 008 把 第 1 层 中 心 点 的 X 坐 标 带 进 ① 式 求 出 Y 拟 合 =-0.669305985852*566.6668 +901.989584064071=522.7161,即第 1 层中心点的 Y 坐标为 522.7161, 类似地可以用 同样的方法得到其它数据;这些数据汇总后,可得到如下表 5。 表 5 如下: 四年塔层拟合后的数据

层数 1 2 3 4 5 6 7 8 9 10 11 12 13 14

4 年平均拟 4 年平均拟 合后 Z 带进 合后 X 带进 两次算术平 中心线方程 中心线方程 均处理的 得到的 X 拟合 得到的 Y 拟合 Z 两次算术平均(单 (单位:m) (单位:m) 位:m) 566.6668 522.7161 1.774531 566.7276 522.6754 7.308594 566.7874 522.6354 12.74128 566.835 522.6035 17.06878 566.8861 522.5693 21.71244 566.9357 522.5361 26.22003 566.9754 522.5096 29.82769 567.0141 522.4837 33.34319 567.0526 522.4579 36.84228 567.0891 522.4335 40.16125 567.1361 522.402 44.43344 567.183 522.3706 48.70072 567.2284 522.3402 52.82395 567.2535 522.3234 55.10525

两间点的距离 D(每层层 高平方之差再开方)(单 位:m) D1 5.534547401 D2 5.433163519 D3 4.327879181 D4 4.644063133 D5 4.507988711 D6 3.607972357 D7 3.515808032 D8 3.499400345 D9 3.319259562 D10 4.272561834 D11 4.267655154 D12 4.123588961 (算层高时塔尖不考虑)

由上表可知塔与塔间的距离,为了能够更好的得出层高且让误差减少,所以求出 13 个 塔层间两两塔层的观测点间的距离,如此就可得出 D1,D2,??,D12。12 组塔层距离, 随后对这 12 组数据进行算术平均处理的得出层高 H=(D1+D2+D3+D4+D5+D6+D7+D8+D9+D10+D11+D12) ? 12=4.254491m 。因为塔层间的距 离相等 所以第一层的中心高度为 2.127245m, 第二层就是第一层的中心高度加上层高即 4.254491m +2.127245m =6.381736m,如此类推,可得出第 N 层的的中心点的 Z 坐标为
H ? H ? ( N ? 1) 2 。因此所有层的中心点的 Z 坐标如下表 6.

表 6: 层数 1 2 塔层中心点的 Z 坐标(单位:m) 2.127245 6.381735
8

3 4 5 6 7 8 9 10 11 12 13 14

10.636225 14.890715 19.145205 23.399695 27.654185 31.908675 36.163165 40.417655 44.672145 48.926635 53.181125 55.10525

把各个塔层的 Z 代入 4 个年份拟合出来的方程组
?Y拟合 ? -0.6693059 85852 X 拟合 ? 901.989584 064071 ① ,再得各个塔层中心的 X 和 Y ? ② ? Z 拟合 ? 90.8972421 008 X 拟合 - 51506.6709 048934

并 如 将 第 1 层 的 中 心 点 的 Z 坐 标 带 进 ② 式 求 出 X 拟 合 51506.6709 048934 ? 2.127245 = =566.6706377 , 即 第 1 层 中 心 点 的 X 坐 标 为 90.8972421 008 566.6706377 ; 再 把 第 1 层 中 心 点 的 X 坐 标 带 进 ① 式 求 出 Y 拟 合 =-0.669305985852*566.6706377 +901.989584064071=522.7135342, 即第 1 层中心点的 Y 坐标为 522.7135342,类似可得表 7 以下各组数据。 表 7 如下: 各层中心坐标 塔层中心 X 塔层中心 Y 塔层中心 Z 层 (单位:m) (单位:m) (单位:m) 1 566.6706377 522.7135342 2.127245341 2 566.7174432 522.682207 6.381736024 3 566.7642487 522.6508798 10.63622671 4 566.8110542 522.6195526 14.89071739 5 566.8578597 522.5882254 19.14520807 6 566.9046652 522.5568982 23.39969875 7 566.9514707 522.525571 27.65418944 8 566.9982762 522.4942438 31.90868012 9 567.0450817 522.4629166 36.1631708 10 567.0918872 522.4315894 40.41766148 11 567.1386927 522.4002622 44.67215217 12 567.1854982 522.368935 48.92664285 13 567.2323037 522.3376078 53.18113353 即各层的中心坐标分别为下表 8 表 8:
9

各个塔层中心坐标点 塔层 数 1 2 3 4 5 6 7 8 9 10 11 12 13 塔层中心坐标点 M1(566.6706377,522.7135342,2.127245341) M2(566.7174432,522.682207,6.381736024) M3(566.7642487,522.6508798,10.63622671) M4(566.8110542,522.6195526,14.89071739) M5(566.8578597,522.5882254,19.14520807) M6(566.9046652,522.5568982,23.39969875) M7(566.9514707,522.525571,27.65418944) M8(567.0450817,522.4629166,36.1631708) M9(567.0450817,522.4629166,36.1631708) M10(567.0918872,522.4315894,40.41766148) M11(567.1386927,522.4002622,44.67215217) M12(567.1854982,522.368935,48.92664285) M13(567.2323037,522.3376078,53.18113353)

14 M14(567.2535,522.3234,55.10525) 通过上图可以直观的看出各个塔层的水平方向的中心坐标没有发生较大的偏差, 故说明塔层越高偏移相对性大且呈现出倾斜的趋势 得出通用方法为: 求出层高后,第一层的中心位置的 Z 轴坐标为 H/2,然后将 H/2 代入拟合中心轴上,即 可得出第一层的中心位置的 X 轴与 Y 轴的坐标;第二层的中心的 Z 轴坐标为 H/2+H,然 后将 H/2 代入拟合中心线上,即可得出第二层的中心位置的 X 轴与 Y 轴的坐标;以此类 H H 推,易知,第 N 层的中心位置的 Z 轴坐标为 ? ( N ? 1) ? H ,同样把 ? ( N ? 1) ? H 代入 2 2 拟合中心轴上,得到第 N 层的中心位置 X 轴与 Y 轴的坐标,从而可以得出第 N 层的中心 位置的所有坐标。 5.1.2 各层中心坐标的优化模型二: 在模型一的求解过程中,发现运用加权平均的处理方法可以减少求出塔中心坐标 的误差,使模型更加贴近真实情况。因此建立模型二。 为了让模型更贴近现实,故对年份进行一个加权平均处理,求出一个组数据,而 这组数据将会是误差更小的。 根据古塔的塔龄增加,承受自重、气温、风力等各种作 用的时间加长,遭受地震、飓风的影响较多,古塔的倾斜、弯曲、扭曲等程度都会增大 考虑,故在加权平均处理上随着年份的增加,所占权重越小;我们把权重假设为 1986 年占 40%,1996 占 30%,2009 年占 20%,2011 年占 10%。因此我们利用表一的数据通过 以下公式(加权平均处理)把每一层的数据代进去,求得出四个年份中每一层的平均数 值。 公式如下: X 加权=( X1986 ? 40% ? X1996 ? 30% ? X 2006 ? 20% ? X 2011 ?10% ) ?4 ;

10

Y 加权=( Y1986 ? 40% ? Y1996 ? 30% ? Y2006 ? 20% ? Y2011 ?10% ) ?4 ; Z 加权=( Z1986 ? 40% ? Z1996 ? 30% ? Z2006 ? 20% ? Z2011 ?10% ) ?4 ; 例如 X 加权 566.66475 * 40% ? 566.66495 * 30% ? 566.726812 5 * 20% ? 566.72695 *10% = =566.6834; 4 522.7105 * 40% ? 522.7102 * 30% ? 522.70146 * 20% ? 522.701362 5 *10% Y 加权= =522.7077 4 1.787375 * 40% ? 1.783 * 30% ? 1.7645 * 20% ? 1.76325 *10% Z 加权= =1.779075 4 类似可以得出表 9 中 14 组数据如下: 表 9: 加权平均处理数据 坐标 层 X 加权(单位:m) Y加权(单位:m) Z加权(单位:m) 1 566.6834 522.7077 1.779075 2 566.7332 522.6683 7.313338 3 566.782 522.6301 12.74646 4 566.8208 522.5993 17.07299 5 566.8625 522.5665 21.71526 6 566.9211 522.5262 26.22555 7 566.9577 522.5081 29.83105 8 566.9995 522.4872 33.3456 9 567.0412 522.4659 36.8474 10 567.0825 522.4401 40.16578 11 567.1289 522.4007 44.43598 12 567.1749 522.3614 48.7053 13 567.1428 522.5976 52.8277 14 567.2762 522.2328 55.11213 得出加权处理后的各层数据在对其进行数据的拟合处理得出方程组 把各个塔层的 Z 加权代入 4 个年份拟合出来的方程组
?Y加权 ? -0.6349034 02438X 加权 ? 882.486100 074819 ? ? Z 加权 ? 93.4500334 137X 加权 - 52952.9634 028605 ① (详细 MATLAB 编码见附录 2 ) ②

,再得各个塔层中心的 X 加权和 Y 加权,如将第 1 层的中心点的Z加权坐标带进②式求出Z加权 52952.9634 028605 ? 1.779075 = =566.6637,即第 1 层中心点的 X 加权坐标为 566.6637 ; 93.4500334 137 再 把 第 1 层 中 心 点 的 X 加 权 坐 标 带 进 ① 式 求 出 Y 加 权 =-0.634903402438*566.6637 +882.486100074819=522.7094,即第 1 层中心点的 Y 加权坐标为 522.7094,类似可得表 10 以下各组数据。
11

表 10 四年塔层拟合后的数据

层数 1 2 3 4 5 6 7 8 9 10 11 12 13 14 塔 尖 567.2344 522.347 55.11213 (算层高时塔尖不考虑) 由上表可知塔与塔间的距离,为了能够更好的得出层高且让误差减少,所以求出 13 个 塔层间两两塔层的观测点间的距离,如此就可得出 D1,D2,??,D12。12 组塔层距离, 随后对这 12 组数据进行算术平均处理的得出层高 H=(D1+D2+D3+D4+D5+D6+D7+D8+D9+D10+D11+D12) ? 12= 4.254394m。因为塔层间的距 离相等 所以第一层的中心高度为 2.127197m, 第二层就是第一层的中心高度加上层高即 4.254394m +2.127197m =6.381591m,如此类推,可得出第 N 层的的中心点的 Z 坐标为
H ? H ? ( N ? 1) 2 。因此所有层的中心点的 Z 坐标如下表 11。

4 年平均拟 4 年平均拟 合后 Z 带进 合后 X 带进 两次算术平 中心线方程 中心线方程 均处理的 得到的 X 加权 得到的 Y 加权 Z 两次加权平均(单 (单位:m) (单位:m) 位:m) 566.6637 522.7094 1.779075 566.7229 522.6718 7.313338 566.7811 522.6349 12.74646 566.8274 522.6055 17.07299 566.877 522.5739 21.71526 566.9253 522.5433 26.22555 566.9639 522.5188 29.83105 567.0015 522.4949 33.3456 567.039 522.4711 36.8474 567.0745 522.4486 40.16578 567.1202 522.4196 44.43598 567.1659 522.3906 48.7053 567.21 522.3626 52.8277

两间点的距离 D(每层层 高平方之差再开方)(单 位:m) D1 5.534707 D2 5.433561 D3 4.326873 D4 4.642648 D5 4.51065 D6 3.60579 D7 3.514832 D8 3.502081 D9 3.318642 D10 4.270543 D11 4.269668 D12 4.122733

表 11: 层数 1 2 3 4 5 6 7 8 9 塔层中心点的 Z 加权坐标(单位:m) 2.127197 6.381591 10.63598 14.89038 19.14477 23.39917 27.65356 31.90795 36.16235
12

10 40.41674 11 44.67114 12 48.92553 13 53.17992 14 55.11213 把各个塔层中心点的 Z 加权代入 4 个年份拟合出来的方程组
?Y加权 ? -0.6349034 02438X 加权 ? 882.486100 074819 ? ? Z 加权 ? 93.4500334 137X 加权 - 52952.9634 028605 ① ②

,再得各个塔层中心的 X 加权和 Y 加权,如将第 1 层的中心点的Z加权坐标带进②式求出Z加权 52952.9634 028605 ? 2.127197 = =566.6674,即第 1 层中心点的 X 加权坐标为 566.6674 ; 93.4500334 137 再 把 第 1 层 中 心 点 的 X 加 权 坐 标 带 进 ① 式 求 出 Y 加 权 =-0.634903402438*566.6674 +882.486100074819=522.707,即第 1 层中心点的 Y 加权坐标为 522.707,类似可得表 12 以下各组数据。 表 12 如下: 各层中心坐标 塔层中心 X 塔层中心 Y 塔层中心 Z 层 (单位:m) (单位:m) (单位:m) 1 566.6674 522.707 2.127197 2 566.713 522.6781 6.381639 3 566.7585 522.6492 10.63603 4 566.804 522.6203 14.89043 5 566.8495 522.5914 19.14482 6 566.8951 522.5625 23.39922 7 566.9406 522.5336 27.65361 8 566.9861 522.5047 31.908 9 567.0316 522.4758 36.1624 10 567.0772 522.4469 40.41679 11 567.1227 522.418 44.67118 12 567.1682 522.3891 48.92558 13 567.2137 522.3602 53.17997

即各层的中心坐标分别为下表 13。 表 13: 各个塔层中心坐标 塔层 数 塔层中心坐标 1 M1(566.6674, 522.707,2.127197) 2 M2(566.7129599,522.6781136,6.381639306) 3 4 M3(566.7584858,522.6492091,10.63603327) M4(566.8040117,522.6203045,14.89042723)
13

5 6 7 8 9 10 11 12 13 14

M5(566.8495375,522.5914,19.1448212) M6(566.8950634,522.5624955,23.39921516) M7(566.9405893,522.533591,27.65360913) M8(566.9861152,522.5046864,31.90800309) M9(567.031641,522.4757819,36.16239706) M10(567.0771669,522.4468774,40.41679102) M11(567.1226928,522.4179728,44.67118498) M12(567.1682186,522.3890683,48.92557895) M13(567.2137445,522.3601638,53.17997291) M14(567.2344,522.347,55.11213)

由表 1 和表 8 数据通过 matlab 画出的三维图形(matlab 编程见附录 3)

上图可以直观的描绘古塔的立体图形且其中涵盖着空间中轴线,不但能够直观的看 出各个塔层段而且更能通过塔层中心坐标的连线段比较出各个塔的发展趋势 5.2 古塔的倾斜 对于问题 2 古塔倾斜的分析建立模型三,我们可以根据 1986 年、1996 年、2009 年 和 2011 年所测的数据先用算术平均数求每一层每个点的平均值, 例如取 1986 年第一层 观测的数据 X= (565.4540+562.0580+561.3900+563.7820+567.9410+571.2550+571.9380+569.5000) ? 8=566.66475 ; Y= (528.0120+525.5440+521.4470+518.1080+517.4070+519.8570+523.9530+527.3560) ? 8=522.7105 ;
14

Z=(1.7920+1.8180+1.7830+1.7690+1.7720+1.7700+1.7940+1.8010) ? 8 =1.787375 如此类推可以得到 1986 年、1996 年、2009 年和 2011 年的每层的相加平均数如下表 14 表 14: 1986 年数据 1996 年数据 坐标 坐标 层 x (单位: m) y (单位: m) z (单位: m) 层 x (单位: m) y (单位: m) z (单位: m) 1 566.6648 522.7105 1.787375 1 566.665 522.7102 1.783 2 566.7196 522.6684 7.32025 2 566.7205 522.6674 7.314625 3 566.7735 522.6273 12.75525 3 566.7751 522.6256 12.75075 4 566.8161 522.5944 17.07825 4 566.8183 522.5922 17.07513 5 566.8621 522.5591 21.7205 5 566.8649 522.5563 21.716 6 566.9084 522.5244 26.23513 6 566.9118 522.521 26.2295 7 566.9468 522.5081 29.83688 7 566.9506 522.5042 29.83225 8 566.9843 522.4924 33.35088 8 566.9884 522.4881 33.34538 9 567.0218 522.4764 36.85488 9 567.0265 522.4714 36.84825 10 567.0569 522.4624 40.17213 10 567.062 522.4572 40.16763 11 567.1045 522.423 44.44088 11 567.1102 522.4173 44.43538 12 567.1518 522.3836 48.71188 12 567.1578 522.3775 48.70738 13 567.085 522.7403 52.83429 13 567.0912 522.734 52.83 14 567.2473 522.2438 55.12325 14 567.2544 522.2367 55.11975 2009 年数据 2011 年数据 坐标 坐标 层 x (单位: m) y (单位: m) z (单位: m) 层 x (单位: m) y (单位: m) z (单位: m) 1 566.7268 522.7015 1.7645 1 566.727 522.7014 1.76325 2 566.764 522.6693 7.309 2 566.7642 522.669 7.2905 3 566.8001 522.6384 12.73225 3 566.8004 522.6387 12.72688 4 566.8293 522.6132 17.06975 4 566.8297 522.6127 17.052 5 566.8604 522.5866 21.70938 5 566.861 522.586 21.70388 6 566.9471 522.5342 26.211 6 566.9478 522.5335 26.2045 7 566.9792 522.5123 29.82463 7 566.98 522.5115 29.817 8 567.0305 522.4797 33.33988 8 567.0313 522.4788 33.33663 9 567.0816 522.4466 36.84375 9 567.0825 522.4457 36.82225 10 567.137 522.3937 40.16113 10 567.1381 522.3926 40.14413 11 567.1799 522.3547 44.43263 11 567.181 522.3535 44.42488 12 567.2225 522.316 48.69975 12 567.2238 522.3147 48.68388 13 567.2712 522.2715 52.81838 13 567.2725 522.2701 52.81313 14 567.336 522.2148 55.091 14 567.3375 522.2135 55.087 再用拟合数据的方法可以拟合得出各年份的拟合中心位置的直线 P1、P2、P3 和 P4(详 细 matlab 编码见附录 4) 直线 P1:

15

?Y 1 = -0.5934204 72662X 1 + 858.971114 09774 ? ? Z1 ? 97.3298424 36X 1 ? - 55150.8630 515990

直线 P2:
?Y2 ? -0.5975321 41319X 2 ? 861.300620 492342 ? ? Z 2 ? 96.1669058 154X 2 - 54491.9082 965752

直线 P3:
?Y3 ? -0.7718499 50362 X 3 ? 960.128944 594815 ? ? Z 3 ? 84.3387697 304 X 3 - 47790.5081 085657

直线 P4:
?Y4 ? - 0.77247188 1683X 4 ? 960.481428 041587 ? ? Z 4 ? 84.1573793 412X 4 ? - 47687.7341 371109

而这些直线分别经过塔与地面的切面圆心,以地面为参考系即是 XY 平面,此时塔底 切面圆心的 Z1、Z2、Z3、Z4 的坐标都为 0,将每年 Z 的坐标的值代入拟合塔层中心位 置直线得出塔底中心坐标的 X 轴和 Y 轴坐标,即塔底中心 O 的坐标为(XOK, YOK,0) ,我 们在该圆心点分别建立一条垂直于地面(即 X,Y 平面)的轴线 LK,在塔层中心位置直线 上任意取一个点 Mk(XMk,YMk,ZMk),过 M 点做垂直于轴线 LK 的垂线于点 WK,那么 WK 的坐标为 (XWK,YWK,ZWK) ,所以应有则 XWK = XOK 、YOk =YWK 、ZWK= ZMk ,即 WK=(XOK,YOK,ZMk) O,M,W ,这三个空间坐标点构成直角三角形, 所以角 MOW 就是拟合塔层中心位置直线 PK 与轴线 LK 所成的夹角 t。 OW= (( X W K ? X OK )^ 2 ? (YW K ? YOK )^ 2 ? ( ZW K ? 0)^ 2) = ZWK OM= (( X MK ? X OK )^ 2 ? (YMK ? YOK )^ 2 ? ( Z MK ? 0)^ 2) WM= (OW ^2 ? OM ^2)
WM OM 例如 1986 年的中心 O 坐标(XO1,YO1,0)求解即

tan t =

Y1=-0.593420472662X1+858.97111409774



Z1=97.329842436X1 + -55150.8630515990 ② 把 ZO1=0 带进②式得 XO1=566.6388 。把 XO1 代进①式得 YO1 = 522.7160,所以中心 O 坐 标(566.6388, 522.7160,0) 假设我们任意取一个点 M 的 Z 坐标都为 8,则 M1=(XM1,YM1,8) 把 Z1 带进②式可以得出 XM1=566.72097345549, 把 XM1 代进 ①式得 YM1 =522.6673。 即 M1 为 (566.72097, 522.6673, 8) ,W1=( 566.6388,522.716,8) 所以 OW= (( X W K ? X OK )^ 2 ? (YW K ? YOK )^ 2 ? ( ZW K ? 0)^ 2) = Z MK =8

16

OM= ((566.72097 ? 566.6388 )^2 ? ( 522.6673 ? 522.7160 )^2 ? (8 ? 0)^2) = 8.00057 WM=

(OW ^2 ? OM ^2) = (8.00057 ^2 ? 8^2) =0.0955776

WM 0.0955776 = =0.0119471 OM 8 类似可求得各年份的数据如下表 15 表 15 OM 的长 OW 的长 度 度 1986 8.000571 8 1996 8.000587 8 2009 8.000897 8 2011 8.000902 8 利用上表数据得到图 2 如下 图2

tan t =

WM 长度 0.095578 0.096908 0.119824 0.120119

tan t 0.011947 0.012114 0.014978 0.015015

斜率变化折线图 0.02

斜率

0.01 0 1986 1996 2009 2011 年份

tan t

5.3 古塔的弯曲 5.3.1 古塔的弯曲模型四 对于问题 2 古塔弯曲的分析建立模型四, 要分析塔的弯曲情况, 从古塔的曲率入手; 首先我们用 matlab 软件分别对 1986 年 7 月、1996 年 8 月、2009 年 3 月和 2011 年 3 月 对该塔进行了 4 次观测的数据进行拟合,拟合成二次多项式,分别为: (详细 Matlab 软 件程序见附录 5)
2 ?Y ? -0.574605251X 376X1 - 183835.161 349422 1 ? 650.946992 1986 年: ? 1 2 X1 ? 39129.9111 344X1 - 11119834.7 987881 ? Z1 ? - 34.4235380 2 ?Y2 ? -0.551656532X 572X2 - 176458.707 031079 2 ? 624.925445 1996 年: ? 2 7 X 2 ? 37452.4927 91X2 - 10644069.3 14419 ? Z 2 ? -32.945010 2 ?Y3 ? - 0.21885203 61* X 3 ? 247.415868 5235X 251517 3 - 69403.7772 2009 年: ? 2 1 * X 3 ? 54619.0845 833X3 - 15508981.5 815841 ? Z 3 ? -48.088762

17

2 ?Y4 ? - 0.02161803 943X 4 ? 244.385828 3476X4 - 68544.6481 083914 2011 年: ? 2 Z 4 ? 47.7721482 X 4 ? 54259.9242 496X4 -15407125. 6586547 ?

对 1986 年的拟合二次多项式曲线
2 ?Y1 ? -0.574605251X 376X1 - 183835.161 349422 1 ? 650.946992 ? 2 X1 ? 39129.9111 344X1 - 11119834.7 987881 ? Z1 ? - 34.4235380

进行空间曲线的曲率求算: 设 t=X1,则得出为: X1=t,Y1=-0.574605251t2+650.946992376t-183835.161349422,Z1=-34.4235380t2+3 9129.9111344t-11119834.7987881 (1) ① 讨论其在 P 点 (x1,y1,15) 的曲率,将 z1=15 代入①式, 解得 x1= 566.791589879, y1= 522.62741530, 把曲线看作是质点在空间中的运动轨迹,则质点的速度向量及加速度向量可以由曲 线方程中的一阶及二阶倒数给出: (求导 见附录 6) 速度向量 v=(1,-0.4159,107.9675) (2) 加速度向量 a=(0,-1.1492,-68.8471) (3) 其加速度可以分解为两部分:一个是与速度方向平行的切向加速度,另一个是与速 度方向垂直的法向加速度。 由圆周运动中速度与加速度的关系就可以确定圆周运动的半 径,也就是曲率半径,从而求出曲率。圆周运动关系式为:|a|=|v|2/2 (4) 其中 r 是圆周运动的半径。 切向加速度大小就是加速度向量在速度向量上的投影值, 由向量的运算规则及(2) (3) 可以求出切向加速度大小为: (v*a)/|v| 则法向加速度大小为: (| a |^2 ? (v * a)^2/ | v |^2) (5)

由(4) (5)计算出曲率半径: r=|v|^2/

(| a |^2 ? (v * a)^2/ | v |^2)









r=|v|^3/ (| a |^2* | v |^2 ? (v * a)^2) (6) 曲率 k1= (| a |^2* | v |^2 ? (v * a)^2) /|v|^3 把(2)(3)代入(7)中,解得 k1=0.000444。 同理可解得 1996 年、2009 年、2011 年的曲率分别为 k2=0.000481,k3=0.000141, k4=0.000142411。 ② 讨论其在 P 点 (x2,y2,30) 的曲率,将 z2=30 代入①式, 解得 x2=566.9372888,y2= 522.627415. 同 理 用 模 一 的 方 法 求 得 各 年 份 在 塔 高 30 米 的 曲 率 分 别 为 : K1=0.000423,K2=0.000459,K3=0.000142,K4=0.000144076. (7)

18

③ 讨 论 其 在 P 点 ( x3,y3,45 ) 的 曲 率 , 将 z3=45 代 入 ① 式 , 解 得 x3=567.099723,y3=522.44472.同理用模一的方法求得各年份在塔高 30 米的曲率分别 为:K1’=0.000405,K2’=0.000441,K3’=0.000145,K4’=0.000147087. 综上列出各年份的不同高度的曲率表 16 如下。 表 16: 年份 1986 1996 15 米高 的曲率 0.000444 0.000481 30 米高 的曲率 0.000423 0.000459 45 米高 的曲率 0.000405 0.000441 平均曲 率 0.000424 0.00046 根据上表所示,利用 图 3:

2009

2011

0.000141 0.000142 0.000142 0.000144 0.000145 0.000147 0.000143 0.000145

excel 做出各年份的不同高度的曲率变化趋势图图 3 如下.

0.002 0.0015 平均曲率 45米高的曲率 30米高的曲率 15米高的曲率

曲率

0.001 0.0005 0 1986 1996 2009 年份 2011

5.3.2 古塔的弯曲优化模型五 先将 4 个年份的算术平均处理数据 (即表 1) 通过 matlab (matlab 编码详见附录 8) 进行三次多项式数据拟合。为了研究出塔低、中、高三个不同高度的弯曲程度,我们选 取低(15m) 、中(30m) 、高(45m)三个层次来比较古塔 4 个年份古塔的扭曲程度,接 着分别把 Z=15,30,45 代入三次多项式拟合的三次项式如下: 1986 年:
3 2 ? Y1 ? - 8.14449X 6X1 - 853278.199 20X1 ? 1484098983 .65290 1 ? 13852.1760 ? 3 2 2 X1 - 199363050. 2574X 0.5327 1 ? 351670.899 1 ? 3767310478 ?Z1 ? - 206.7791X

19

1996 年:
3 2 ? Y2 ? 7.384537 X 2 ? 13343.5189 0 X 2 ? 7564951.90 784X2 ? 1429620872 .64414 ? 3 2 X 2 ? 338098.805 2 X 2 - 191670297. 0802X 2 ? 3621967052 2.5616 ?Z 2 ? - 198.7975

2009 年:
3 2 ?Y3 ? - 0.916012 X 3 ? 1558.00345 6 X 3 - 883313.361 250 X 3 ? 166932782. 340235 ? 3 2 X 3 - 32442.0551 6X3 ? 183983949. 2686X 768172 3 - 3.47800673 ? Z 3 ? 190.6846

2011 年:
3 2 ?Y4 ? - 0.901336 X 4 ? 1533.04368 8 X 4 - 869163.116 374X4 ? 164258753. 638376 曲 ? 3 2 X 4 - 321752.494 1 X 4 ? 182471151. 3122X4 - 3449414678 1.8903 ? Z 4 ? 189.1161

率计算公式: k=(

( y' z' ' - y'' z ' )^2 + z' ' ^2 + y'' ^2) )^(1/2) (曲率计算公式的得出见附录 7) (y'^2 ? z' ^2)^3

1986 年在 15m 处的曲率 k=(
( y' z' ' - y'' z ' )^2 + z' ' ^2 + y'' ^2) )^(1/2)= 8.71916E-06 (y'^2 ? z' ^2)^3

其他年份的不同高度的曲率类似可得如下表 17。 表 17: 15 米高的曲 30 米高的曲 率 率 1986 8.71916E-06 8.76396E-06 1996 1.73018E-12 1.72753E-12 2009 7.92269E-07 7.88483E-07 2011 7.69521E-06 7.65689E-06 根据上表所示,利 excel 做出曲率变化趋势图图 4 如下 图 4:
曲率变化趋势图 0.00003

45 米高的曲 率 8.80913E-06 1.72488E-12 7.84726E-07 7.61895E-06

不同塔高的挠率

0.000025 0.00002 0.000015 0.00001 0.000005 0 1986 1996 年份 2009 2011

45米高的曲率 30米高的曲率 15米高的曲率

由上表及趋势图分析从 1986 年至 1996 年曲率大幅降低, 1996 年至 2009 年以及 2009 年至 2011 年都呈增长趋势,可能是随着观测年份的增加,古塔的年龄增长,故弯曲的
20

程度也会随之增大。充分的证明出塔随着年份的增加,塔身的离 1986 年的基准越大, 且倾斜倾斜的幅度使先升高后减小,因此假如使塔继续保持这种势态的话塔就不会倒 塌。 5.4 古塔的扭曲 对于问题 2 古塔扭曲的分析,要分析塔的扭曲情况,从古塔的挠率入手。 (挠率:它 的绝对值度量了曲线上邻近两点的次法向量之间的夹角对弧长的变化率) 挠率 i=|
y'' z ' ' '? y'' ' z ' ' | y'' ^2 ? z' ' ^2 + ( y' z' '- y'' z' )^2

(8) (挠率的具体计算公式的见附录 7)

首先我们利用 matlab 分别对 1986 年 7 月、1996 年 8 月、2009 年 3 月和 2011 年 3 月对该塔进行了 4 次观测的数据进行拟合,拟合成三次多项式,分别为: (Matlabf 软件 程序见附录 7) 1986 年:
? Y1 ? - 8.14449X 3 ? 13852.1760 6X 2 - 853278.199 20X ? 1484098983 .65290 ? 3 2 ?Z1 ? - 206.7791X ? 351670.899 2 X - 199363050. 2574X ? 3767310478 0.5327

1996 年:
? Y2 ? 7.384537 X 3 ? 13343.5189 0 X 2 ? 7564951.90 784X ? 1429620872 .64414 ? 3 2 ?Z 2 ? - 198.7975 X ? 338098.805 2 X - 191670297. 0802 X ? 3621967052 2.5616

2009 年:
?Y3 ? - 0.916012 X 3 ? 1558.00345 6 X 2 - 883313.361 250 X ? 166932782. 340235 ? 3 2 ? Z 3 ? 190.6846 X - 32442.0551 6X ? 183983949. 2686X - 3.47800673 768172

2011 年:
?Y4 ? - 0.901336 X 3 ? 1533.04368 8 X 2 - 869163.116 374X ? 164258753. 638376 ? 3 2 ? Z 4 ? 189.1161 X - 321752.494 1 X ? 182471151. 3122X - 3449414678 1.8903

为了研究出塔低、中、高三个不同高度的弯曲程度,我们选取低(15m) 、中(30m) 、 高(45m)三个层次来比较古塔 4 个年份古塔的扭曲程度,接着分别把 Z=15,30,45 代 入三次多项式拟合的方程组; Ⅰ、对 1986 年的拟合三次多项式曲线,讨论其在 15m 处的挠率,将 z1=15 代入
? Y1 ? - 8.14449X 3 ? 13852.1760 6X 2 - 853278.199 20X ? 1484098983 .65290 中,解 ? 3 2 Z ? 206.7791X ? 351670.899 2 X 199363050. 2574X ? 3767310478 0.5327 ? 1

得 x=565.412646。 对拟合三次项式进行各阶求导(Matlab 软件程序见附录 9) ;将 x=565.412646 代入 求导后的各式中,得到各阶求导后的值列出表 18 如下. 表 18: Z1 Y1 的一阶导 Y1 的二阶导 Y1 的三阶导 Z1 的一阶导 Z1 的二阶导 Z1 的三阶导 15m -58.8647 74.36626 -48.8669 -1240.27 1848.69 -1240.67 30m -58.7624 74.29901 -48.8669 -1237.73 1846.982 -1240.67
21

45m -58.66 74.23164 -48.8669 -1235.18 1845.272 -1240.67 Z2 Y2 的一阶导 Y2 的二阶导 Y2 的三阶导 Z2 的一阶导 Z2 的二阶导 Z2 的三阶导 15m 29747407 51748.27 44.30722 -855.922 1529.564 -1192.79 30m 29747508 51748.36 44.30722 -852.938 1527.235 -1192.79 45m 29747609 51748.45 44.30722 -849.949 1524.899 -1192.79 Z3 Y3 的一阶导 Y3 的二阶导 Y3 的三阶导 Z3 的一阶导 Z3 的二阶导 Z3 的三阶导 15m -5.51544 7.213747 -5.49607 1287.117 -1689.12 1144.108 30m -5.52602 7.221808 -5.49607 1289.596 -1690.8 1144.108 45m -5.5366 7.229852 -5.49607 1292.072 -1692.47 1144.108 Z4 Y4 的一阶导 Y4 的二阶导 Y4 的三阶导 Z4 的一阶导 Z4 的二阶导 Z4 的三次导 15m -3.20905 5.516029 -5.40802 848.2398 5.516029 1134.697 30m -3.22167 5.528385 -5.40802 851.3123 5.528385 1134.697 45m -3.23427 5.540695 -5.40802 854.379 5.540695 1134.697 将 表 中 的 数 值 代 入 ( 8 ) 中 , 得 出 在 i=|
y'' z ' ' '? y'' ' z ' ' |=6.90788E-06 y'' ^2 ? z' ' ^2 + ( y' z' '- y'' z' )^2

15m

处 的 挠 率

Ⅱ、其他年份的不同高度的挠率类似可得如下表: 各年不同高度的挠率表 19 如下。 表 19: 年份 1986 1996 15 米高的挠率 30 米高的挠率 6.90788E-06 6.92213E-06 2.9789E-14 2.988E-14

2009 0.000361 0.0003602

2011 0.000285102 0.000282414

45 米高的挠率 6.93644E-06 2.9971E-14 0.0003595 0.000279767 根据上表所示,利用 excel 做出挠率变化趋势图图 5 如下: 由下表及趋势图图 5,分析从 1986 年至 1996 年挠率变化不大, 可能是塔龄还不能达 到扭曲度;1996 年至 2009 年增长较大,说明古塔在这个年间所受多方因素影响,造成 扭曲度大;2009 年至 2011 年又呈下降趋势,则可能是古塔弯曲与扭曲的双重影响。 5.5 古塔倾斜、弯曲和扭曲的分析 对于问题三: Ⅰ、倾斜程度的分析: 各年斜率如下表 20。 表 20: 年份 1986 1996 2009 2011 斜率 0.011947 0.012114 0.014978 0.015015
22

图 5:

挠率变化趋势图 0.0012

不同塔高的挠率

0.001 0.0008 0.0006 0.0004 0.0002 0 1986 1996 2009 年份 2011

45米高的曲率 30米高的曲率 15米高的曲率

各年斜率折线图图 6 如下. 图 6:

各年斜率趋势图 0.015015 0.014978

斜率

0.012114 0.011947

1986

1996 年份

2009

2011

由上表和趋势图可知塔的倾斜(斜率)随着的观测年份的增加,斜率也随之增加, 很明确的表明了塔的倾斜角度也在逐年增大。 Ⅱ、弯曲程度的分析: 各年不同高度的曲率如下表 21。 表 21 年份 1986 1996 2009 2011 15 米高的曲率 0.000444 0.000481 0.000141 0.000142411 30 米高的曲率 0.000423 0.000459 0.000142 0.000144076 45 米高的曲率 0.000405 0.000441 0.000145 0.000147087
23

塔的曲率由 1986 年到 1996 年增大,说明塔的弯曲程度变大;1996 年到 2009 年的 曲率大幅下降,说明塔的弯曲呈不规则状态,以及塔的扭曲影响;2009 年到 2011 年塔 的曲率小幅增加,可推断塔有回旋的趋势。 Ⅲ、扭曲程度的分析如下图 7。 图 7:

挠率变化趋势图 0.0012

不同塔高的挠率

0.001 0.0008 0.0006 0.0004 0.0002 0 1986 1996 2009 年份 2011

45米高的曲率 30米高的曲率 15米高的曲率

由挠率变化趋势图可知 1986 年和 1996 年的古塔高中低三个层段挠率变化大致相 同, 这充分的证明 1986 年到 1996 年间古塔高中低三个层段的扭曲程度由大致相同到升 高再降低,所以塔在 1986 年至 1996 年间大体上没有发生扭曲,然而从 1996 年至 2009 年间塔由于各种因素使塔发生了扭曲, 随后 2009 年至 2011 年间扭曲程度又相对性减小。 Ⅳ、对塔三个不同高度(15m 30m 45m )的挠率进行算术平均,得出下表,对古 塔进行进一步的扭曲度分析如下表 22。 表 22: 平均挠率变化表 1986 1996 2009 2011 6.92E-06 2.99E-14 3.60E-04 2.82E-04 根据平均挠率变化表利用 excel 画出平均挠率的变化趋势图如下图 8. 由下图可知 1986 年-2011 年间古塔整体扭曲程度呈现出变大的趋势,然而从折线 的 数 据 来 分 析 , 知 道 古 塔 这 三 个 年 分 段 间 的 挠 率 是 先 从 0.00000692215 减 少 到 0.00000000000002988 再 重 新 上 升 为 0.000360233333333333 再 下 降 为 0.000282427666666667, 这可充分说明了由于长时间承受自重、 气温、 风力等各种作用, 偶然还要受地震、飓风的影响,古塔会产生各种变形致使 1996 年到 2009 年间扭曲的跨 度比较大。

24

图 8:

各年平均挠率趋势图 3.60E-04 2.82E-04

挠率

2.99E-14 6.92E-06

系列2 系列1

1986

1996 年份

2009

2011

六、 模型评价
建模过程中运用算术平均以及加权平均,使多项数据简化,方便于模型求解。 我们建立了五个模型,并对其中的两个模型进行了优化,共建立了七个模型。在数 据处理上,运用 Matlab 的线性拟合、二次多项式拟合、三次多项式拟合等方法来求得 空间直线方程。并给出了斜率公式、曲率公式、挠率公式来计算古塔的倾斜、弯曲和扭 曲等变形情况,还能过画图、列表等方法对倾斜、弯曲、扭曲这三种情况进行了数据分 析。 不足的地方是,由于时间关系,未能对其它的模型进行优化处理。

七、 参考文献
[1]史峰,MATLAB 函数速查手册,北京:中国铁道出版社,2011。 [2]薛毅,数学建模基础,北京:北京工业大学出版社,2005。 [3]张德丰,MATLAB 数值分析与应用,北京:国防工业出版社,2009。 [4]叶其孝,大学生数学建模竞赛辅导材料[M],长沙:湖南教育出版社,2002 。 [5]张序,测量学,南京:东南大学出版社,2007

八、 附录
附录 1: 不同年份对古塔的观测数据拟合 MATLAB 编码如下: X 和 Y 的线性拟合 format long X=[566.6959 566.7421 566.7873 566.8234 566.8621 566.9288 566.9641 567.0086 567.0531 567.0985 567.1439 567.1890 567.1800 567.2938]; Y=[522.7059 522.6685 522.6325 522.6031 522.5720 522.5282 522.5091 522.4847 522.4600 522.4265 522.3871 522.3480 522.5040 522.2272];
25

A=polyfit(X,Y,1) B=polyval(A,X); plot(X,Y,'K+',X,B,'r') A = 1.0e+002 * -0.00669305985852 9.01989584064071 X 和 Z 的线性拟合 X=[566.6959 566.7421 566.7873 566.8234 566.8621 566.9288 566.9641 567.0086 567.0531 567.0985 567.1439 567.1890 567.1800 567.2938]; Z=[1.7745 7.3086 12.7413 17.0688 21.7124 26.2200 29.8277 33.3432 36.8423 40.1613 44.4334 48.7007 52.8239 55.1052]; C=polyfit(X,Z,1) D=polyval(C,X); plot(X,Z,'K+',X,D,'r') C = 1.0e+004 * 0.00908972421008 -5.15066709048934 附录 2: 模型优化中模型二的拟合 MTLAB 编码如下: X 和 Y 的线性拟合 format long X=[566.6834 566.7332 566.7820 566.8208 566.8625 566.9211 566.9577 566.9995 567.0412 567.0825 567.1289 567.1749 567.1428 567.2762]; Y=[522.7077 522.6683 522.6301 522.5993 522.5665 522.5262 522.5081 522.4872 522.4659 522.4401 522.4007 522.3614 522.5976 522.2328]; A=polyfit(X,Y,1) B=polyval(A,X); plot(X,Y,'K+',X,B,'r') A = 1.0e+002 * -0.00634903402438 8.82486100074819 X 和 Z 的线性拟合 format long X=[566.6834 566.7332 566.7820 566.8208 566.8625 566.9211 566.9577 566.9995 567.0412 567.0825 567.1289 567.1749 567.1428 567.2762]; Z=[1.7791 7.3133 12.7465 17.0730 21.7153 26.2255 29.8311 33.3456 36.8474 40.1658 44.4360 48.7053 52.8277 55.1121]; C=polyfit(X,Z,1) D2=polyval(C,X); plot(X,Z,'K+',X,D,'r') C = 1.0e+004 * 0.00934500334137 -5.29529634028605

26

附录 3: a=[565.454 528.012 1.792;562.058 525.544 1.818;561.39 521.447 1.783;563.782 518.108 1.769;567.941 517.407 1.772;571.255 519.857 1.77;571.938 523.953 1.794;569.5 527.356 1.801;565.454 528.012 1.792;565.48 527.764 7.326;562.238 525.364 7.351;561.663 521.42 7.314;564.001 518.226 7.301;567.995 517.563 7.306;571.165 519.961 7.304;571.801 523.908 7.324;569.414 527.141 7.336;565.48 527.764 7.326;565.506 527.52 12.761;562.415 525.188 12.786;561.931 521.394 12.749;564.216 518.343 12.736;568.048 517.716 12.741;571.076 520.063 12.74;571.666 523.864 12.758;569.33 526.93 12.771;565.506 527.52 12.761;565.526 527.327 17.084;562.555 525.047 17.109;562.144 521.373 17.072;564.387 518.435 17.059;568.091 517.838 17.064;571.005 520.144 17.063;571.558 523.829 17.081;569.263 526.762 17.094;565.526 527.327 17.084;565.548 527.119 21.726;562.706 524.896 21.751;562.373 521.351 21.714;564.571 518.534 21.701;568.136 517.969 21.705;570.929 520.232 21.708;571.443 523.791 21.723;569.191 526.581 21.736;565.548 527.119 21.726;565.57 526.915 26.267;562.854 524.748 26.309;562.6 521.329 26.308;564.752 518.632 26.264;568.18 518.095 26.189;570.857 520.315 26.136;571.333 523.755 26.164;569.121 526.406 26.244;565.57 526.915 26.267;565.671 526.652 29.869;563.132 524.585 29.911;562.883 521.356 29.91;564.949 518.846 29.866;568.172 518.346 29.791;570.679 520.441 29.737;571.094 523.672 29.765;568.994 526.167 29.846;565.671 526.652 29.869;565.77 526.397 33.383;563.403 524.427 33.425;563.158 521.382 33.424;565.141 519.055 33.38;568.164 518.59 33.305;570.506 520.564 33.251;570.862 523.591 33.279;568.87 525.933 33.36;565.77 526.397 33.383;565.868 526.141 36.887;563.674 524.268 36.929;563.433 521.408 36.928;565.333 519.263 36.884;568.156 518.834 36.809;570.333 520.686 36.755;570.63 523.51 36.783;568.747 525.701 36.864;565.868 526.141 36.887;565.961 525.9 40.201;563.927 524.12 40.214;563.693 521.433 40.244;565.516 519.462 40.223;568.148 519.068 40.171;570.171 520.801 40.038;570.408 523.433 40.129;568.631 525.482 40.157;565.961 525.9 40.201;566.078 525.628 44.472;564.193 523.95 44.485;563.958 521.463 44.505;565.649 519.607 44.486;568.094 519.242 44.442;570.013 520.885 44.309;570.236 523.35 44.4;568.615 525.259 44.428;566.078 525.628 44.472;566.195 525.355 48.743;564.459 523.78 48.756;564.224 521.492 48.776;565.782 519.753 48.757;568.039 519.415 48.713;569.854 520.969 48.58;570.063 523.268 48.671;568.598 525.037 48.699;566.195 525.355 48.743;566.308 525.092 52.866;564.716 523.616 52.878;564.481
27

521.521 52.897;565.91 519.893 52.88;569.701 521.05 52.703;569.897 523.188 52.794;568.582 524.822 52.822;566.308 525.092 52.866;567.255 522.238 55.128;567.235 522.242 55.108;567.247 522.251 55.128;567.252 522.244 55.129;567.255 522.238 55.128] b=[566.6706377 522.7135342 2.127245341;566.7174432 522.682207 6.381736024;566.7642487 522.6508798 10.63622671;566.8110542 522.6195526 14.89071739;566.8578597 522.5882254 19.14520807;566.9046652 522.5568982 23.39969875;566.9514707 522.525571 27.65418944;566.9982762 522.4942438 31.90868012;567.0450817 522.4629166 36.1631708;567.0918872 522.4315894 40.41766148;567.1386927 522.4002622 44.67215217;567.1854982 522.368935 48.92664285;567.2323037 522.3376078 53.18113353;567.2534718 522.3234399 55.10525] b1=b(:,1); b2=b(:,2); b3=b(:,3); x=a(:,1); y=a(:,2); z=a(:,3); plot3(x,y,z,'r.-'); hold on; plot3(b1,b2,b3,'bo-'); 附录 4: 1986 年拟合成 1 次多项式的 Matlab 编程 X 和 Y 的线性拟合 format long X1=[566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467 566.9842 567.0218 567.0569 567.1045 567.1517 567.0850 567.2473]; Y1=[522.7105 522.6684 522.6273 522.5944 522.5591 522.5244 522.5081 522.4924 522.4764 522.4624 522.4230 522.3836 522.7403 522.2437]; A1=polyfit(X1,Y1,1) B1=polyval(A1,X1); plot(X1,Y1,'K+',X1,B1,'r') A1 = 1.0e+002 * -0.00593420472662 8.58971114097744 X 和 Z 的线性拟合 format long X1=[566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467 566.9842 567.0218 567.0569 567.1045 567.1517 567.0850 567.2473]; Z1=[1.7874 7.3202 12.7553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119 52.8343 55.1232];
28

C1=polyfit(X1,Z1,1) D1=polyval(C1,X1); plot(X1,Z1,'K+',X1,D1,'r') C1 = 1.0e+004 * 0.00973298424365 -5.51508630515990 1996 年拟合成 1 次多项式的 Matlab 编程 X 和 Y 的线性拟合 format long X2=[566.6649 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.0620 567.1102 567.1578 567.0912 567.2544]; Y2=[522.7102 522.6674 522.6256 522.5922 522.5563 522.5210 522.5042 522.4881 522.4714 522.4572 522.4173 522.3775 522.7340 522.2367]; A2=polyfit(X2,Y2,1) B2=polyval(A2,X2); plot(X2,Y2,'K+',X2,B2,'r') A2 = 1.0e+002 * -0.00597532141319 8.61300620492342 X 和 Z 的线性拟合 format long X2=[566.6649 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.0620 567.1102 567.1578 567.0912 567.2544]; Z2=[1.7830 7.3146 12.7508 17.0751 21.7160 26.2295 29.8322 33.3454 36.8483 40.1676 44.4354 48.7074 52.8300 55.1198]; C2=polyfit(X2,Z2,1) D2=polyval(C2,X2); plot(X2,Z2,'K+',X2,D2,'r') C2 = 1.0e+004 * 0.00961669058154 -5.44919082965752 2009 年拟合成 1 次多项式的 Matlab 编程 X 和 Y 的线性拟合 X 和 Y 的线性拟合 format long X3=[566.7268 566.7640 566.8001 566.8293 566.8604 566.9471 566.9792 567.0305 567.0816 567.1370 567.1799 567.2225 567.2712 567.3360]; Y3=[522.7015 522.6693 522.6384 522.6132 522.5866 522.5342 522.5123 522.4797 522.4466 522.3937 522.3547 522.3160 522.2715 522.2148]; A3=polyfit(X3,Y3,1) B3=polyval(A3,X3); plot(X3,Y3,'K+',X3,B3,'r') A3 = 1.0e+002 *
29

-0.00771849950362 9.60128944594815 X 和 Z 的线性拟合 format long X3=[566.7268 566.7640 566.8001 566.8293 566.8604 566.9471 566.9792 567.0305 567.0816 567.1370 567.1799 567.2225 567.2712 567.3360]; Z3=[1.7645 7.3090 12.7323 17.0697 21.7094 26.2110 29.8246 33.3399 36.8438 40.1611 44.4326 48.6998 52.8184 55.0910]; C3=polyfit(X3,Z3,1) D3=polyval(C3,X3); plot(X3,Z3,'K+',X3,D3,'r') C3 = 1.0e+004 * 0.00843387697304 -4.77905081085657 2011 年拟合成 1 次多项式的 Matlab 编程 X 和 Y 的线性拟合 X 和 Y 的线性拟合 format long X4=[566.7269 566.7642 566.8004 566.8297 566.8610 566.9478 566.9800 567.0313 567.0825 567.1381 567.1810 567.2238 567.2725 567.3375]; Y4=[522.7014 522.6690 522.6387 522.6127 522.5860 522.5335 522.5115 522.4788 522.4457 522.3926 522.3535 522.3147 522.2701 522.2135]; A4=polyfit(X4,Y4,1) B4=polyval(A4,X4); plot(X4,Y4,'K+',X4,B4,'r') A4 = 1.0e+002 * -0.00772471881683 9.60481428041587 X 和 Z 的线性拟合 format long X4=[566.7269 566.7642 566.8004 566.8297 566.8610 566.9478 566.9800 567.0313 567.0825 567.1381 567.1810 567.2238 567.2725 567.3375]; Z4=[1.7633 7.2905 12.7269 17.0520 21.7039 26.2045 29.8170 33.3366 36.8222 40.1441 44.4249 48.6839 52.8131 55.0870]; C4=polyfit(X4,Z4,1) D4=polyval(C4,X4); plot(X4,Z4,'K+',X4,D4,'r') C4 = 1.0e+004 * 0.00841573793412 -4.76877341371109 附录 5: 1986 年拟合成二次多项式的 Matlab 编程: X 和 Y 的二次多项式拟合 format long X1=[566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467
30

566.9842 567.0218 567.0569 567.1045 567.1517 567.0850 567.2473]; Y1=[522.7105 522.6684 522.6273 522.5944 522.5591 522.5244 522.5081 522.4924 522.4764 522.4624 522.4230 522.3836 522.7403 522.2437]; A1=polyfit(X1,Y1,2) B1=polyval(A1,X1); plot(X1,Y1,'K+',X1,B1,'r') A1 = 1.0e+005 * -0.00000574605251 0.00650946992376 -1.83835161349422 X 和 Z 的二次多项式拟合 format long X1=[566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467 566.9842 567.0218 567.0569 567.1045 567.1517 567.0850 567.2473]; Z1=[1.7874 7.3202 12.7553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119 52.8343 55.1232]; C1=polyfit(X1,Z1,2) D1=polyval(C1,X1); plot(X1,Z1,'K+',X1,D1,'r') C1 = 1.0e+007 * -0.00000344235380 0.00391299111344 -1.11198347987881 1996 年拟合成二次多项式的 Matlab 编程: X 和 Y 的二次多项式拟合 format long X2=[566.6649 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.0620 567.1102 567.1578 567.0912 567.2544]; Y2=[522.7102 522.6674 522.6256 522.5922 522.5563 522.5210 522.5042 522.4881 522.4714 522.4572 522.4173 522.3775 522.7340 522.2367]; A2=polyfit(X2,Y2,2) B2=polyval(A2,X2); plot(X2,Y2,'K+',X2,B2,'r') A2 = 1.0e+005 * -0.00000551656532 0.00624925445572 -1.76458707031079 X 和 Z 的二次多项式拟合 format long X2=[566.6649 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.0620 567.1102 567.1578 567.0912 567.2544]; Z2=[1.7830 7.3146 12.7508 17.0751 21.7160 26.2295 29.8322 33.3454 36.8483 40.1676 44.4354 48.7074 52.8300 55.1198]; C2=polyfit(X2,Z2,2) D2=polyval(C2,X2); plot(X2,Z2,'K+',X2,D2,'r') C2 =
31

1.0e+007 * -0.00000329450107 0.00374524927910 -1.06440693144190 2009 年拟合成二次多项式的 Matlab 编程: X 和 Y 的二次多项式拟合 format long X3=[566.7268 566.7640 566.8001 566.8293 566.8604 566.9471 566.9792 567.0305 567.0816 567.1370 567.1799 567.2225 567.2712 567.3360]; Y3=[522.7015 522.6693 522.6384 522.6132 522.5866 522.5342 522.5123 522.4797 522.4466 522.3937 522.3547 522.3160 522.2715 522.2148]; A3=polyfit(X3,Y3,2) B3=polyval(A3,X3); plot(X3,Y3,'K+',X3,B3,'r') A3 = 1.0e+004 * -0.00002188520361 0.02474158685235 -6.94037772251517 X 和 Z 的二次多项式拟合 format long X3=[566.7268 566.7640 566.8001 566.8293 566.8604 566.9471 566.9792 567.0305 567.0816 567.1370 567.1799 567.2225 567.2712 567.3360]; Z3=[1.7645 7.3090 12.7323 17.0697 21.7094 26.2110 29.8246 33.3399 36.8438 40.1611 44.4326 48.6998 52.8184 55.0910]; C3=polyfit(X3,Z3,2) D3=polyval(C3,X3); plot(X3,Z3,'K+',X3,D3,'r') C3 = 1.0e+007 * -0.00000480887621 0.00546190845833 -1.55089815815841 2011 年拟合成二次多项式的 Matlab 编程: X 和 Y 的二次多项式拟合 format long X4=[566.7269 566.7642 566.8004 566.8297 566.8610 566.9478 566.9800 567.0313 567.0825 567.1381 567.1810 567.2238 567.2725 567.3375]; Y4=[522.7014 522.6690 522.6387 522.6127 522.5860 522.5335 522.5115 522.4788 522.4457 522.3926 522.3535 522.3147 522.2701 522.2135]; A4=polyfit(X4,Y4,2) B4=polyval(A4,X4); plot(X4,Y4,'K+',X4,B4,'r') A4 = 1.0e+004 * -0.00002161803943 0.02443858283476 -6.85446481083914 X 和 Z 的二次多项式拟合 format long X4=[566.7269 566.7642 566.8004 566.8297 566.8610 566.9478 566.9800 567.0313 567.0825 567.1381 567.1810 567.2238 567.2725 567.3375];
32

Z4=[1.7633 7.2905 12.7269 17.0520 21.7039 26.2045 29.8170 33.3366 36.8222 40.1441 44.4249 48.6839 52.8131 55.0870]; C4=polyfit(X4,Z4,2) D4=polyval(C4,X4); plot(X4,Z4,'K+',X4,D4,'r') C4 = 1.0e+007 * -0.00000477721482 0.00542599242496 -1.54071256586547 附录 6: 模型四求曲率的优化 一阶和二阶求导 MATLAB 编码如下: 曲线方程 (1)一阶求导: syms t; diff(-0.574605251*t^2+650.946992376*t-183835.161349422,t,1) diff(-34.4235380*t^2+39129.9111344*t-11119834.7987881,t,1) 曲线方程 (1)二阶求导: syms t; diff(-0.574605251*t^2+650.946992376*t-183835.161349422,t,2) diff(-34.4235380*t^2+39129.9111344*t-11119834.7987881,t,2) 附录 7
? x?t ? 设 t=x,则得出曲线方程为: ? y ? y (t ) ,即曲线为(t,y(t),z(t)) ? z ? z (t ) ?
?? 设? =曲线的一阶导数=(1,y’(t),z’(t))=(x 1 , y 1 , z 1 ) a ? ?? =曲线的一阶导数=(0,y’’(t),z’’(t))=(x 2 ,y 2 , z 2 ) b ? ?? =曲线的一阶导数=(0,y’’’(t),z’’’(t))=(x 3 , y 3 , z 3 ) c

所以有 x 1 =1,x 2 =0,x 3 =0;
?? ?? 所以( ? *? )^2=( y 1 z 2 - y 2 z 1 )^2+( x 2 z 1 - x 1 z 2 )^2+( x 1 y 2 - x 2 y 1 )^2; a b ?? ?? ?? ?? ?? ?? 因为 ( ? ,? ,? )=( ? *? )* ? =(y 1 z 2 - y 2 z 1 ,x 2 z 1 - x 1 z 2 , a b c a b c

x 1 y 2 - x 2 y 1 )*(x 3 , y 3 , z 3 )=( y 1 z 2 - y 2 z 1 )* x 3 +( x 2 z 1 - x 1 z 2 )* y 3 +( x 1 y 2 - x 2 y 1 )* z3; 所以 挠率
33

i=

(? ?? ,? ?? ,? ?? ) ( y1 z 2 - y 2 z 1) * x 3 + ( x 2 z1 - x 1z 2 ) * y 3 + ( x1 y 2 - x 2 y1 ) * z 3 a b c = ; ( y1 z 2 - y 2 z 1)^2 ? ( x 2 z1 - x 1z 2 )^2 + ( x1 y 2 - x 2 y1 )^2 (? ?? *? ?? )^2 a b
y 2 z3 ? y 3 z 2 y 2 ^2 ? z 2 ^2 + ( y1 z 2 - y 2 z1 )^2
y'' z ' ' '? y'' ' z ' ' | y'' ^2 ? z' ' ^2 + ( y' z' '- y'' z' )^2

所以 挠率 i=

即 挠率 i=|

关于曲率计算公式的计算: k=

|? ?? *? ?? | a b |? ?? |^3 a

=

(( y1 z 2 - y 2 z 1)^2 + ( x 2 z1 - x 1z 2 )^2 + ( x1 y 2 - x 2 y1 )^2) ^ (1/2) (x1 ^2 ? y1 ^2 + z1 ^2)^(3/2)
把 x1 =1, x 2 =0 代入上述方程得

曲率 k=(

( y1 z 2 - y 2 z 1)^2 + z 2 ^2 + y 2 ^2) )^(1/2) (y1 ^2 ? z1 ^2)^3

=(

( y' z' ' - y'' z ' )^2 + z' ' ^2 + y'' ^2) )^(1/2) (y'^2 ? z' ^2)^3

附录 8: 曲率模型 5 优化 matlab 编码 三次多项式拟合的 MATLAB 编码如下: 1986 三次多项式拟合编码 X 和 Y 的三次多项式拟合 format long X1=[566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467 566.9842 567.0218 567.0569 567.1045 567.1517 567.0850 567.2473]; Y1=[522.7105 522.6684 522.6273 522.5944 522.5591 522.5244 522.5081 522.4924 522.4764 522.4624 522.4230 522.3836 522.7403 522.2437]; A1=polyfit(X1,Y1,3) B1=polyval(A1,X1); plot(X1,Y1,'K+',X1,B1,'r') A1 = 1.0e+009 * -0.00000000814449 0.00001385217606 -0.00785327819920 1.48409898365290
34

X 和 Z 的三次多项式拟合 format long X1=[566.6648 566.7196 566.7735 566.8161 566.8621 566.9084 566.9467 566.9842 567.0218 567.0569 567.1045 567.1517 567.0850 567.2473]; Z1=[1.7874 7.3202 12.7553 17.0783 21.7205 26.2351 29.8369 33.3509 36.8549 40.1721 44.4409 48.7119 52.8343 55.1232]; C1=polyfit(X1,Z1,3) D1=polyval(C1,X1); plot(X1,Z1,'K+',X1,D1,'r') C1 = 1.0e+010 * -0.0000000206779 0.0000351670899 -0.01993630502574 3.76731047805327 1996 三次多项式拟合编码 X 和 Y 的三次多项式拟合 format long X2=[566.6649 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.0620 567.1102 567.1578 567.0912 567.2544]; Y2=[522.7102 522.6674 522.6256 522.5922 522.5563 522.5210 522.5042 522.4881 522.4714 522.4572 522.4173 522.3775 522.7340 522.2367]; A2=polyfit(X2,Y2,3) B2=polyval(A2,X2); plot(X2,Y2,'K+',X2,B2,'r') A2 = 1.0e+009 * -0.00000000784537 0.0000133435189 -0.00756495190784 1.42962087264414 X 和 Z 的三次多项式拟合 format long X2=[566.6649 566.7205 566.7751 566.8183 566.8649 566.9118 566.9506 566.9884 567.0265 567.0620 567.1102 567.1578 567.0912 567.2544]; Z2=[1.7830 7.3146 12.7508 17.0751 21.7160 26.2295 29.8322 33.3454 36.8483 40.1676 44.4354 48.7074 52.8300 55.1198]; C2=polyfit(X2,Z2,3) D2=polyval(C2,X2); plot(X2,Z2,'K+',X2,D2,'r') C2 = 1.0e+010 * -0.00000001987975 0.00003380988052 -0.01916702970802 3.62196705225616 2009 三次多项式拟合编码 X 和 Y 的三次多项式拟合 format long X3=[566.7268 566.7640 566.8001 566.8293 566.8604 566.9471 566.9792 567.0305 567.0816 567.1370 567.1799 567.2225 567.2712 567.3360]; Y3=[522.7015 522.6693 522.6384 522.6132 522.5866 522.5342 522.5123 522.4797 522.4466 522.3937 522.3547 522.3160 522.2715 522.2148];
35

A3=polyfit(X3,Y3,3) B3=polyval(A3,X3); plot(X3,Y3,'K+',X3,B3,'r') A3 = 1.0e+008 * -0.00000000916012 0.00001558003456 -0.00883313361250 1.66932782340235 X 和 Z 的三次多项式拟合 format long X3=[566.7268 566.7640 566.8001 566.8293 566.8604 566.9471 566.9792 567.0305 567.0816 567.1370 567.1799 567.2225 567.2712 567.3360]; Z3=[1.7645 7.3090 12.7323 17.0697 21.7094 26.2110 29.8246 33.3399 36.8438 40.1611 44.4326 48.6998 52.8184 55.0910]; C3=polyfit(X3,Z3,3) D3=polyval(C3,X3); plot(X3,Z3,'K+',X3,D3,'r') C3 = 1.0e+010 * 0.00000001906846 -0.00003244205516 0.01839839492686 -3.47800673768172 2011 三次多项式拟合编码 X 和 Y 的三次多项式拟合 format long X4=[566.7269 566.7642 566.8004 566.8297 566.8610 566.9478 566.9800 567.0313 567.0825 567.1381 567.1810 567.2238 567.2725 567.3375]; Y4=[522.7014 522.6690 522.6387 522.6127 522.5860 522.5335 522.5115 522.4788 522.4457 522.3926 522.3535 522.3147 522.2701 522.2135]; A4=polyfit(X4,Y4,3) B4=polyval(A4,X4); plot(X4,Y4,'K+',X4,B4,'r') A4 = 1.0e+008 * -0.00000000901336 0.00001533043688 -0.00869163116374 1.64258753638376 X 和 Z 的三次多项式拟合 format long X4=[566.7269 566.7642 566.8004 566.8297 566.8610 566.9478 566.9800 567.0313 567.0825 567.1381 567.1810 567.2238 567.2725 567.3375]; Z4=[1.7633 7.2905 12.7269 17.0520 21.7039 26.2045 29.8170 33.3366 36.8222 40.1441 44.4249 48.6839 52.8131 55.0870]; C4=polyfit(X4,Z4,3) D4=polyval(C4,X4); plot(X4,Z4,'K+',X4,D4,'r') C4 = 1.0e+010 * 0.00000001891161 -0.00003217524941 0.01824711513122 -3.44941467818903 附录 9
36

三次项式拟合的各年份各阶求导 matlab 编码如下: 1986 年: Y 的一阶导数: diff(-8.14449*X1^3+13852.17606*X1^2-7853278.19920*X1+1484098983.65290,X1, 1) Y 的二阶导数: diff( -8.14449*X1^3+13852.17606*X1^2-7853278.19920*X1+1484098983.65290,X1,2) Y 的三阶导数: diff( -8.14449*X1^3+13852.17606*X1^2-7853278.19920*X1+1484098983.65290,X1,3) Z 的一阶导数: diff(-206.7791*X1^3+351670.8992*X1^2-199363050.2574*X1+37673104780.5327,X1, 1) Z 的二阶导数: diff(-206.7791*X1^3+351670.8992*X1^2-199363050.2574*X1+37673104780.5327,X 1,2) Z 的三阶导数: diff(-206.7791*X1^3+351670.8992*X1^2-199363050.2574*X1+37673104780.5327,X1, 3) 1996 年: Y 的一阶导数: diff(7.384537*X2^3+13343.51890*X2^2+7564951.90784*X2+1429620872.64414,X2, 1) Y 的二阶导数: diff(7.384537*X2^3+13343.51890*X2^2+7564951.90784*X2+1429620872.64414,X2, 2) Y 的三阶导数: diff(7.384537*X2^3+13343.51890*X2^2+7564951.90784*X2+1429620872.64414,X2, 3) Z 的一阶导数: diff(-198.7975*X2^3+338098.8052*X2^2-191670297.0802*X2+36219670522.5616,X 2,1) Z 的二阶导数: diff(-198.7975*X2^3+338098.8052*X2^2-191670297.0802*X2+36219670522.5616,X 2,2) Z 的三阶导数: diff(-198.7975*X2^3+338098.8052*X2^2-191670297.0802*X2+36219670522.5616,X 2,3) 2009 年: Y 的一阶导数: diff( -0.916012*X3^3+1558.003456*X3^2-883313.361250*X3+166932782.340235,X 3,1) Y 的二阶导数:
37

diff( -0.916012*X3^3+1558.003456*X3^2-883313.361250*X3+166932782.340235,X3, 2) Y 的三阶导数: diff( -0.916012*X3^3+1558.003456*X3^2-883313.361250*X3+166932782.340235,X 3,3) Z 的一阶导数: diff(190.6846*X3^3-324420.5516*X3^2+183983949.2686*X3-34780067376.8172,X3 ,1) Z 的二阶导数: diff(190.6846*X3^3-324420.5516*X3^2+183983949.2686*X3-34780067376.8172,X3 ,2) Z 的三阶导数: diff(190.6846*X3^3-324420.5516*X3^2+183983949.2686*X3-34780067376.8172,X3,3 ) 2011 年: Y 的一阶导数: diff(-0.901336*X4^3+1533.043688*X4^2-869163.116374*X4+164258753.638376,X4 ,1) Y 的二阶导数: diff(-0.901336*X4^3+1533.043688*X4^2-869163.116374*X4+164258753.638376,X4,2 ) Y 的三阶导数: diff(-0.901336*X4^3+1533.043688*X4^2-869163.116374*X4+164258753.638376,X4 ,3) Z 的一阶导数: diff(189.1161*X4^3-321752.4941*X4^2+182471151.3122*X4-34494146781.8903,X4 ,1) Z 的二阶导数: diff(189.1161*X4^3-321752.4941*X4^2+182471151.3122*X4-34494146781.8903,X4 ,2) Z 的三阶导数: diff(189.1161*X4^3-321752.4941*X4^2+182471151.3122*X4-34494146781.8903,X4,3 )

38


相关文章:
古塔的变(胡、罗、钟)形
古塔的变(胡、罗、钟)形_学科竞赛_高中教育_教育专区。2013 数学建模 D题 广东省2等奖承 诺 书 我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生...
2014年大兴区一模语文卷及答案
征蓬出汉塞,归雁入天。 大漠孤烟直,长河落日圆...美不胜收的古典园林和宏伟壮观的坛庙、 帝陵、古塔...12.答案示例一: 从语气上看,钟子期评价都用感叹句...
山东省邹城市第一中学2016届高三语文4月模拟考试试题
古塔完全被树根抬举起来 ③终于觉察古塔是被树根簇拥...上的钟花女士,这笔善款将全部用于花女士的康复...岚罗、机场三条高速公路,形 成安全便捷、畅通高效...
盐城市2015届高三年级第一学期期中考试语文答案
译颇有学识,兼知律,善骑射,周武帝时, 译时丧...个尖角上,另一个尖角上有热那亚式古塔,古塔里住着...(5)无边落木萧萧下 (7)凌万顷之茫然 (2)声非...
填空
梁思成研究中国古塔,把塔分为楼阁式塔、密檐塔 空间...河南登封嵩岳寺塔 857 是我国现存最早的八边形塔...钟鼓楼、天王殿 北京故宫前三殿中的主殿奉天殿,...
考试卷
? 梁思成研究中国古塔,把塔分为楼阁式塔、密檐塔、...钟鼓楼、天王殿、大雄宝殿、法堂、藏经楼、僧人用...嵩岳寺塔平面和立面 (12 边形、15 层,高 40 米...
南京市2015~2016学年度第一学期期末学情调研测试卷语文
段横线处的词语,最恰当的一组是 在河南省开封市东北隅,有一座千年古塔,即有...这时情况又在变化。 光亮已经消逝, 那恐怖的迹象也荡然无存。雪下得正紧,大片...
2016届江苏省高考压轴预测卷 语文(word版)
另一 页 4第 个尖角上有热那亚式古塔,古塔里住着...(5)无边落木萧萧下 (6)风霜高洁 (7)凌万顷之...烘托守塔人的形 象。(2分)与下文环境描写相呼应...
人教版四年级上册语文全部课文
镇海的古塔、中山亭 和观潮台屹立在江边。远处,几...再近些,只见白浪翻滚,形 成一道两丈多高的白色...我一连看了 两个头,看得有些不耐烦了。 住宅...
图形和变换
把一个长方形作相似变换, 各条放大到原来的 3 ...钟示数,这时的实际时间应该是 18、如图,三角形 ...23:如图,小明欲测量一座古塔的高度,他站在该塔的...
更多相关标签:
古塔胡同 | 古塔胡同是什么 | 钟形曲线 | 钟形帽 | 钟形浮子式蒸汽疏水阀 | 钟形分布 | 钟形虫 | 杨殿钟和胡小芬图片 |