基于空间直线位置关系的古塔变形研究

摘  要： 针对古塔变形问题，给出了一种基于极限思想的古塔各层中心的计算方法;利用空间直线的位置关系分析了古塔的各种变形，并利用MATLAB数值实验对模型进行了计算从而得到了数值结果。

关键词： 多边形的中心;极限思想;空间直线的位置关系;古塔变形;MATLAB

中图分类号： O29    文献标识码： A    DOI：10.3969/j.issn.1003-6970.2020.10.062

本文著录格式：游晋峰. 基于空间直线位置关系的古塔变形研究[J]. 软件，2020，41（10）：246249

【Abstract】： In this paper， the deformation problem of ancient pagoda is studied. A calculation method of the center of each layer of ancient pagoda based on limit thought is given; the various deformations of ancient pagoda is analyzed by using the position relationship between spatial lines， and the model is calculated by using MATLAB numerical experiment， and the numerical results are obtained.

【Key words】： The center of polygon; Limit thought; The position relationship between spatial lines; The deformations of ancient pagoda; MATLAB

0  引言

古塔长时间承受自重、气温、风力等作用和地震、飓风的影响，会产生倾斜、弯曲、扭曲等各种变形。2013年全国大学生数学建模竞赛C题“古塔的变形”[1]，给出了1986年7月、1996年8月、2009年3月和2011年3月4次观测数据，要求给出确定古塔各层中心位置的通用方法，并分析古塔的变形情况。

陈俊等人[2-4]利用投影多边形模型来确定古塔各层的中心，采用最小二乘法建立了线性模型，并借助三维高次曲线方程建立古塔的曲率模型和挠率模型。但这种用中轴变形曲线刻画物体的整体扭曲是有缺陷的[5]。蔡志杰等人[5]将古塔看作连续体，利用位移函数来刻画古塔的变形。其中所用知识均涉及了较专业的高等数学知识。而数学建模应遵循一个原则：尽量采用简单的数学工具[6-7]。

本文首先分析了题目所给数据，发现数据观测起止点不同。于是，对数据进行了处理，同时绘制了古塔的三维网面图和俯视图，便于进一步观察古塔。其次，给出了一种基于极限思想的古塔各层中心的计算方法。然后，利用空间直线的位置关系分析了古塔各种变形，利用Matlab软件完成了对模型的求解，得到了数值结果。

1  数据分析

通过分析题目所给的观测数据发现，1986年和1996年观测点的起止一致，2009年和2011年观测点的起止一致。但前两次观测点的起止和后两次的不同。

如圖1，对1986年第一层、1996年第二层、2009年第三层、2011年第四层的数据作图可以看到，1986年和1996年的起止一致，2009年和2011年的起止一致。但前两次观测和后两次观测的起止点不同。

为方便后续的计算和变形分析，需要对数据进行处理，使四次观测数据的起止点一致。

1.1  数据准备

对2009年和2011年观测点的起止顺序进行修改，即，将每层的第7个观测点、第8个观测点的数据依次移动到第1个观测点的数据之前。如图2可以看到，数据修改后，4次观测数据的起止保持一致。

1.2  古塔观测图

利用题目观测数据，绘制了古塔的三维网面图和俯视图，便于进一步观察古塔。由于1996年、2009年和2011年三次观测数据的作图结果和1986年的类似，故此处只给出1986年的古塔三维网面图（图3）和俯视图（图4）。

从1986年古塔的俯视图（如图4）可以看出，古塔各层的结构近似于正八边形，塔尖几乎投影到了一点。

2  基于极限思想的古塔各层中心的计算方法

2.1  确定古塔在xoy平面上各层的中心

对于某次测量数据的第层（）在x轴和y轴的投影，求该多边形的各边的中点，并依次连线形成新的多边形，再求新的多边形的各边的中点，……，如此循环多边形会变小，直到该多边形近似趋于一点，停止循环迭代。此时，得到的点为该多边形的中心。

算法如下：

%赋初值

x_old=0;

y_old=0;

x_new=x（i，：）;

y_new=y（i，：）;

%迭代循环，直至两次坐标差值比小

while （abs（x_old-x_new）>0.0000005） & （abs（y_ old-y_new）>0.0000005）

x_old=x_new;

y_old=y_new;

for j=1：size（x_new，2）-1

x_new（j）=（x_new（j）+x_new（j+1））/2;

y_new（j）=（y_new（j）+y_new（j+1））/2;

end

x_new（size（x_new，2））=x_new（1）;

y_new（size（y_new，2））=y_new（1）;

end

需要注意的是，由于顶层数据严重缺失，只考虑第1～13层。

本文中，考虑允许的观测误差为[8]，故设置当两次迭代的值相差小于时，停止迭代。如图5，对1986年的第一层的中心坐标利用上述算法进行计算，可以看到，多边形逐渐趋于一点。

2.2  确定古塔在z轴上各层的中心

2.1只考虑了古塔在xoy平面的中点，并未考虑z轴方向上的中心坐标。

利用求某次测量的第层（）的z轴方向上的坐标均值，来确定古塔第层（）的z轴方向上的中心坐标。

2.3  确定的古塔各层中心坐标

确定古塔各层中心坐标的最终方案为：

（1）利用无穷逼近的思想，确定古塔第层（）x轴和y轴方向上的中心坐标;

（2）利用z轴方向上的坐标均值，确定古塔第层（）的z轴方向上的中心坐标。

（3）利用四次测量结果的均值确定1986年和1996年塔尖的中心坐标。

（4）直接利用题目所给的一次测量结果确定2009年和2011年塔尖的中心坐标。

最终确定的四次测量的古塔各层中心坐标，见表1。

3  分析古塔的变形情况

针对古塔变形的情况，从古塔的墙体是否有裂缝或挤压、古塔整体的倾斜量和倾斜角、古塔各层的沉降[9]、古塔的扭曲四个角度展开研究。

3.1  古塔墙体的挤压与裂缝

检测墙体线段的长度，与前一年数据比较，根据发生的变化判定是否发生挤压或裂缝使墙体变形。

各层各线段的长度与前一年数据比较，结果如下。

通过分析图6～8发现，1996年到2009年测量数据的波动较大，即出现墙体有裂缝（图中正数）和压缩变形（图中负数）的现象。其他两个的比较数量级都在之间。

3.2  倾斜量和倾斜角

利用2.3节得到的古塔各层中心坐标计算倾斜量和倾斜角。

古塔中心坐标的曲线如图9。

3.2.1  整体的倾斜量和倾斜角

利用第1层的中心坐标和塔尖的中心坐标连线形成的直线的倾斜量和倾斜角，估计整体的倾斜量和倾斜角。

整体倾斜量的计算公式如下：

整体倾斜角的计算公式如下：

利用Matlab软件计算得到：

3.2.2  弯曲，即各层的倾斜量和倾斜角

对于第层（）的倾斜量和倾斜角，利用第层的中心坐标和第层的中心坐标连线形成的直线进行计算。

第层（）倾斜量的计算公式如下：

第层（）倾斜角的计算公式如下：

通過计算得到：

第1层～第13层的倾斜量依次为：0.0493、0.0472、0.0392、0.0411、0.1015、0.0390、0.0608、0.0610、0.0768、0.0580、0.0578、0.0660、0.0862。

第1层～第13层的倾斜角依次为：0.5113、0.4977、0.5188、0.5062、1.2924、0.6178、0.9901、1.0028、1.3249、0.7764、0.7780、0.9163、2.1702。

各层的倾斜量和倾斜角折线图，如图10和图11。

从图10、11可以看到，第5、9、13层的倾斜量和倾斜角相对比较明显，因此要对第5、9、13层加强防护。

3.3  古塔各层的沉降

利用2.3节得到的古塔各层中心的z轴方向的坐标值、、、的差值

分析古塔的各层沉降。

例如：比较1986年和1996年z轴方向的差值，计算公式如下：

1986年和1996年之间、1996年和2009年之间、2009年和2011年之间各层的沉降为见表2（负数表示下沉）。

观察表2发现，古塔各层均有不同程度的下沉，下沉范围在0.0013 m～0.0288 m之间。

3.4  古塔的扭曲变形

利用题目所给的观测数据，通过作图观察四次测量的第层（）的扭曲变形程度，见图12。

通过观察图像发现，第1、2、3层四年之间几乎没发生扭曲，第4、5、7～13层四年之间发生了轻微的的扭曲，第6层1996年和2009年之间的扭曲较严重。

下面以第六层1996年和2009年的数据为例，计算该层此两年间的扭曲变形情况，如图13。

对于第条线段的扭曲度，其计算公式如下：

计算得到，对应的八条线段的扭曲度依次为：1.1823、1.6083、0.3730、0.3603、0.8891、1.8980、2.1883、0.1329。

4  结语

数学建模讲求尽量采用简单的模型，不盲目追求所谓的“高大上”的方法和“全能新”的模型[10]。本文利用初等方法解决古塔变形问题，提出了利用取多边形各边中点的方法无限逼近古塔各层中心，从而得到古塔各层的中心坐标;利用空间直线的位置关系，分析了古塔可能的各种变形情况。对于古塔变形，不仅分析了题目要求的倾斜、弯曲、扭曲等变形，而且分析了题目古塔墙体是否有裂缝、古塔各层的沉降。该数据结果与陈俊的数据结果虽略有差距，但对于初接触数学建模的学生较容易理解。

参考文献

[1]2013年高教社杯全国大学生数学建模竞赛赛题. [DB/OL]. （2013-09-13）[2020-07-18]. http：//www.mcm.edu.cn/problem/ 2013/2013.html.

[2]陈俊. 古塔的形变研究[J]. 数学建模及其应用， 2014年2月第3卷第1期： 67-76.

[3]刘世杰. 基于多层中轴点拟合的古塔变形检测[J]. 同济大学学报（自然科学版）， 2018年3月， 第46卷第3期： 401-405.

[4]侯学慧. 古塔的变形问题研究[J]. 山西师范大学（自然科学版）， 2019年3月， 第33卷第1期： 47-53.

[5]蔡志杰， 谭永基. 用测量数据分析古塔的变形[J]. 数学建模及其应用， 2013年11月， 第2卷第5-6期： 14-19.

[6]姜启源， 谢金星， 叶俊. 数学模型（第四版）[M]. 北京： 高等教育出版社， 2011. 1， 第14页.

[7]徐仁旭， 孔亚仙. 数学实验与建模[M]. 长沙： 湖南师范大学出版社， 2011. 6.

[8]百度文库. 建筑施工测量技术要求及允许偏差. [DB/OL]. （2010-11-15）[2020-07-18]. http：//wenku.baidu.com/view/ 7c5c5f1aa8114431b90dd854.html.

[9]搜狐新闻. 地震加速西安大雁塔变形没有造成破坏性影响. [DB/OL]. （2008-05-26）[2020-07-18]. http：//news.sohu.com/ 20080526/n257078702.shtml.

[10]韩中庚. 数学建模竞赛论文的写作方法[J]. 数学建模及其应用， 2017年6月， 第6卷第2期： 42-48.