古塔变形分析模型

李 静，张 媛，顾 忠，贺 欣，王 鹏

(1.郑州铁路职业技术学院，河南 郑州 450052；2. 郑州市金水区教育体育局, 河南 郑州 450008)

古塔变形分析模型

李 静1，张 媛1，顾 忠1，贺 欣1，王 鹏2

(1.郑州铁路职业技术学院，河南 郑州 450052；2. 郑州市金水区教育体育局, 河南 郑州 450008)

根据2013年全国大学生数学建模比赛C题中某古塔4年的观测数据，给出了确定古塔各层中心位置的通用方法，建立了最优化模型，用Lingo软件求得4次测量的古塔的各层中心坐标。以斜率、曲率、投影、均方差、拟合等知识为基础，对倾斜度、弯曲度、扭曲度三个指标进行定义，结合使用Excel和Matlab软件对古塔变形情况进行量化分析，最后根据得到的数据对古塔的变形趋势进行预测。

倾斜度；弯曲度；扭曲度；均方差；优化模型；曲率；拟合

由于古塔长时间受自身或外界因素影响，会发生各种变形，包括倾斜、弯曲、扭曲。现对古塔4年中个别月份变形情况统计数据进行分析。

1 问题一的模型建立与求解

对于2013年全国大学生数学建模比赛C题(以下简称“C题”)中的问题一，从物理结构的稳定角度考虑，想象成对8个点结绳拉紧，绳长越短绳子越紧，求塔每层中心位置可以转化为求某一点到塔每层各个点距离之和最短。我们建立了最优化模型，运用Lingo软件求解出4年古塔每层的中心坐标，整理绘制成Excel表格并用Matlab软件制图。

1.1 绘制古塔三维立体图形

根据C题附件1给出1986年7月古塔x，y，z坐标，运用Matlab软件对数据进行预处理并且制出1986年古塔三维立体图(如图1)，以方便对问题进行思考与研究。

图1 古塔的基本形状

1.2 求古塔每层中心坐标

由于古塔常年受自身重量以及自然环境影响，发生倾斜、弯曲、扭曲等形变，致使原来在同一平面的8个顶点不在同一平面上。

我们想象成用一根绳子从中心位置处分别和8个点结绳，若要结构最稳定，也就是使结绳最紧，即所结绳长之和最短。利用这个思想，可以把求塔每层中心目标转化成求某一点到塔层各个端点距离之和最短，若求各层中心坐标可转化成某点到该层8个顶点距离之和最小，求这点坐标分量。

建立最优化模型一：

这是个非线性规划模型，运用Lingo软件编程求解，并将4年的各塔层中心目标数据用Excel软件进行统计制表(如表1)。

表1 各塔层中心坐标 单位：m

续表

层数 2009年3月xyz 2011年3月xyz1566．7089522．70261．7645566．7091522．70251．76332566．7490522．66667．3091566．7492522．66647．29063566．7885522．631412．7324566．7889522．631812．72704566．8201522．603217．0697566．8206522．602717．05215566．8542522．573021．7095566．8548522．572421．70406566．9384522．516826．2114566．9391522．516126．20507566．9693522．495429．8250566．9701522．494629．81738567．0192522．462133．3403567．0200522．461233．33709567．0688522．428636．8441567．0697522．427636．822610567．1248522．376540．1612567．1258522．375240．144211567．1681522．338944．4328567．1692522．337744．425112567．2112522．301548．7000567．2126522．300248．684013567．2598522．258152．8188567．2611522．256752．8135塔尖567．3360522．214855．0910567．3375522．213555．0870

2 问题二的模型建立与求解

对于C题问题二，我们首先定义了倾斜度、弯曲度、扭曲度3个指标，并给出计算这些指标的公式。接下来分别通过3个指标对古塔变形情况进行分析。

考虑倾斜变形时，我们细化出总偏移量、南北偏移量、东西偏移量、总倾斜度、南北倾斜度、东西倾斜度等指标，利用各层中心坐标，刻画出古塔4年倾斜程度，得出在2009年倾斜度和偏移量最大。

考虑弯曲变形时，借用均方差的思想，利用各层中心点坐标数据，将古塔14层分为13段，求出古塔每段斜率，用各段斜率与上下中心连线的斜率的差的平方和开方这一偏离程度作为弯曲度，偏离程度越大，弯曲度越大，得出在2009年弯曲度最大。

考虑扭曲变形时，把每年古塔各角在每层的观测点投影到平面xOy上，拟合出8条曲线，利用这8条曲线的平均曲率的总和来刻画古塔每年的扭曲度，得出在1996年扭曲度最大。

2.1 古塔倾斜情况分析

结合问题一得出的各塔层中心坐标，建立空间直角坐标系(如图2)，这里塔尖中心点P投影到塔底中心点O所在平面z=zO上，记为点Q，QM ⊥CM，QN ⊥CN，根据相应的空间几何知识可得出总偏移量、南北偏移量、东西偏移量、总偏移度、南北偏移度、东西偏移度。

图2 塔顶塔底中心坐标示意图

将4年的塔底和塔顶中心坐标代入以下各倾斜指标:

利用MATLAB软件对以上表达式进行编程求解，得出倾斜指标具体数据，同时根据中心坐标的移动得出古塔有南偏、东偏趋势，将变形数据汇总整理得到表2。

表2 古塔倾斜指标汇总

年份 总偏移量m南偏移量m东偏移量m总倾斜度°南倾斜度°东倾斜度°高度m19860．75060．46810．58680．80620．50280．630355．128119960．76080．47480．59440．81710．51000．638455．124020090．79450．48780．62710．85360．52410．673755．091020110．79620．48900．62840．85550．52540．675255．0870

为了清晰观察4年倾斜指标的变化，根据表2数据，运用Excel软件分别绘制古塔方向偏移量折线图(如图3)和古塔方向倾斜度折线图(如图4)。

图3 古塔4年偏移量走势图

图4 古塔4年倾斜度走势图

由表2整理数据结合图3、图4分析出：古塔在1986年、1996年、2009年发生倾斜形变的大小，随着古塔建立年代的久远，偏移量、倾斜度变化呈逐步上升趋势，在不采取对古塔变形预防或补救的情况下，由于古塔材质变形或损坏并且长时间承受自重、气温、风力等因素影响，古塔倾斜量只会逐年增加，满足理论依据；总偏移量和倾斜角度变化最大为1996—2009年，其余年份偏移变化差值均匀增加，另外，1996—2009年塔的高度减少也较大。

古塔每年承受自重、气温、风力等影响是自然的(忽略地震、飓风等自然灾害)，每年古塔发生形变量应大致呈现均匀增加趋势，而1996—2009年差值较大可能因为2008年四川地震，此次大型自然灾害使其发生倾斜程度较大，符合理论推断。

2.2 古塔弯曲情况分析

具体分析古塔弯曲情况，即将古塔弯曲程度进行量化，用具体指标反映古塔4年中的弯曲形变大小，自定义变量弯曲度β。

这里以塔的各层中心点坐标作为参考数据，将古塔分为13段，塔尖—13层、13层—12层、…、3层—2层、2层—1层，如图5。

图5 中心坐标曲线各段斜率图

结合图6分析得出：1986年第11层、12层、13层斜率有下降趋势，意味着古塔11层以上有弯曲，1996年第10层以上斜率减小，即古塔10层以上发生弯曲，而1986年、1996年前10层古塔斜率基本趋于稳定，说明古塔前10层没有发生弯曲；对于2009年、2011年斜率走势图基本吻合,在第5层、第7层、第9层斜率减小，表明古塔在这3层有不同程度的弯曲变形，其他各点没有太大波动，10层以上继续保持弯曲。

我们借用均方差思想，用各段斜率k(i)与k的偏差的平方和开方刻画弯曲度

图6 平均斜率走势图

用Maltab软件计算斜率，运用Excel软件对相应数据进行计算，同样方法处理每年每塔层中心坐标得出所需的弯曲度(如表3)。

表3 4年各塔层弯曲度

根据表3数据明显可以看出：随着古塔建立时间的增加，弯曲度也明显呈上升趋势，与古塔倾斜偏移量变化趋势大致一致。1996—2009年之间年弯曲度差值较大，2009—2011年时间较短，所以弯曲度没有发生太大变化；整体走势还是趋于稳定型增长，符合天气以及随着年代的久远古塔自身材质的改变与损坏规律。

2.3 古塔扭曲情况分析

针对古塔扭曲情况的分析，忽略弯曲等因素的影响。

把古塔各角每年的观测点投影到平面xOy上，拟合出8条曲线，利用这8条曲线的平均曲率的和来刻画古塔每年的扭曲度。

把拟合出的曲线y=ax3+bx2+cx+d的13个点(或12个点)代入曲率公式

计算出每条曲线的平均曲率，则每一年古塔的扭曲度定义为

γ越大，则扭曲程度越大；γ越小，扭曲程度越小。

把1986年7月13层的1号观察点，投影在xOy平面上，运用Matlab软件进行拟合，得出三次方程式

y=-8.3x3+1.4×104x2-7.9×106x+1.5×109， 然后将拟合方程代入曲率表达式，得出1号观察点的平均曲率为0.073 073 846×10-10；同时可以求出其他7个观察点的平均曲率，将这8个曲率相加即为扭曲度

γ=387.222 249 4×10-10。

同理将其他年份数据进行处理求得相应的扭曲度(如表4)。

表4 古塔4年扭曲度量化表

注：由于所得数据偏小，我们目的是分析古塔4年扭曲度变化，故将表中所得数据乘以 10-10。

我们得到了古塔4年的总扭曲度分别为：

387.222 249 4×10-10，68 503.476 02×10-10，93 806.860 49×10-10，93 821.802 53×10-10。

由表4分析得：扭曲度越来越大，古塔扭曲情况逐年严重，同样与倾斜、弯曲呈现的变化趋势相一致，而且从低层到高层扭曲度逐渐增加。

还可以看出：1986—1996年总体扭曲度变化最为明显，1996年、2009年、2011年这些年间扭曲度变化相应较小；在1986年和1996年测量中，4号观测点和7号观测点各自对应的两条棱的扭曲程度最显著，其中，4号观测点的棱的扭曲度变化幅度高于7号观测点；在2009年和2011年测量中，2号观测点和4号观测点各自对应的2条棱的扭曲程度最为显著，而7号观测点对应的棱在这两年的扭曲度反而有所下降。

[1]黄强. 古塔变形监测的探讨[J]. 测绘与空间地理信息，2013,36(6)：217-220.

[2]韩中庚. 数学建模实用教程[M]. 北京：高等教育出版社，2012.

[3](美)Mark M.Meerschaert.数学建模方法与分析[M].北京：机械工业出版社，2015.

[4]杨为民. 建筑物倾斜的计算方法的探讨[J]. 施工技术与测量技术，2010,30(5)：175-176.

[责任编辑：赵 伟]

The Deformation Analysis Model of the Ancient Pagoda

LI Jing1, ZHANG Yuan1,GU Zhong1,HE Xin1, WANG Peng2

(1. Zhengzhou Railway Vocational and Technical College, Zhengzhou 450052,China;2. Zhengzhou Jinshui District Education Sports Bureau,Zhengzhou 450008,China)

According to the 2013 national college students' mathematical modeling competition problem C a pagoda four years of observation data, gives the general method of the ancient towers layers center position, the optimization model is established, using Lingo software, obtained four times while hiking in the measurement of each layer center coordinates; And by slope, curvature, projection, the mean square error, fitting, such as knowledge, slope, bending and torsion degrees "define" the three indexes and combined with using Excel and Matlab software, the quantitative analysis of the deformation of the tower, and according to the data to predict the deformation trend of tower.

gradient; bending; distortion; mean square error; optimization model;curvature; fitting

2016 - 03 - 07

李静(1976— )，女，河北枣强人，郑州铁路职业技术学院公共教学部讲师,研究方向为数学教学。

O29

A

1008-6811(2016)-04-0015-06