常金勇，崔艳星，李海增

（长治学院 数学系，山西 长治 046011）

1 问题重述

某测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对某千年古塔进行了四次观测。所得数据有助于对该塔的倾斜、弯曲、扭曲等变形情况加以研究。问题可以归结为：

(1)通过找出确定古塔各层中心位置的通用方法，给出古塔的各层中心坐标列表；

(2)再分析该塔倾斜、弯曲、扭曲等变形情况；

(3)最后分析该塔的变形趋势。

2 问题的分析

对于问题（1），我们可以从古塔的每层上随机取n个点，绘出塔的轮廓曲线，然后利用MATLAB插值方法[1]37-381求出该曲线的中心，最终找出计算各层中心位置的通用方法，进而计算出古塔各层及塔尖的中心坐标。

对于问题（2），我们从建筑测量学来定义倾斜、弯曲、扭曲等概念[2]67-78。对于古塔这种底面积小、高度不高的建筑，我们用测水平角法[2]来度量它的倾斜。

挠度是指结构构件的轴线或中面由于弯曲引起垂直于轴线或中面方向的线位移。求解方法是通过求出中心位移量、以及中心位移量在平均扭曲方向垂直面上的投影，进一步绘出挠度曲线图来直观给出古塔的弯曲与扭曲程度。

对于问题（3）我们可以根据四次各层的中心坐标用MATLAB[3]101-176绘出插值图，根据图形来分析该塔的弯曲情况。

3 模型的建立与求解

3.1 问题1 的解决方案及模型

在每层取出n个点，分别计算这n个点X,Y,Z坐标的平均值。最后我们得到了求各层中心坐标的通用方法：

我们在实际操作中运用Maple程序[4]24-36计算了各层的中心坐标，以及古塔各层的中心坐标变化。

3.2 问题2 的解决方案及模型

根据建筑测量学知识，古塔主体几何中心就是其底面的中心。要求观测塔尖相对于底面中心的平面位移ΔDi，塔尖的中心与第一层中心的距离hi，由ΔDi与hi的比值来计算出tanθi，运用反三角函数知识来求出每次测量的古塔倾斜弧度δi，再用MATLAB程序[5]求出倾斜角度θi，就可以确定古塔的倾斜情况。（见下面公式）

古塔的倾斜角度结果如下表：

?

结论1：在前十年间该古塔倾斜度变化了0.0103度，平均每年会倾0.00103度。在接下来的十三年里古塔的倾斜度又变化了0.0256度，则平均每年古塔会倾斜0.00197度。在最后的两年里，古塔倾斜了0.0020度，平均每年会倾斜0.0010度。所以古塔开始的一段时间内倾斜的速度相对较慢。在后来的一段时间里，古塔的倾斜速度开始变快，最近几年内古塔的倾斜速度又变的相对缓慢。所以文物部门需适时对古塔进行观测了解各种变形量，以制定必要的加固措施。

3.3 对古塔扭曲与弯曲的研究：

图1 古塔第一次测量的曲率曲线图

图2 古塔第二次测量的曲率曲线图

我们通过挠度和曲度观测得到了如下算法步骤：

步骤1：按下式计算各层中心相对于第一层中心的位移量：

步骤2：按下式用Maple[5]45-49计算各层中心相对于第一层中心的扭曲方向：

步骤3：取各层中心扭曲方向的平均值，作为该古塔主体的扭曲方向：

步骤4：计算各层中心位移量在平均扭曲方向垂直面上的投影：

步骤5：绘制主体曲率曲线图与挠度曲线图(单位为米)。以第一层中心为原点，以观测点的高度为纵坐标，以中心点的位移量为横坐标，用MATLAB软件的三次样条插值的方法[1]，运用文献[6]中方法绘出主体曲率曲线图（见图1-图4）。

图3 古塔第三次测量的曲率曲线图

图4 古塔第四次测量的曲率曲线图

以第一层中心为原点，以观测点的高度为纵坐标，以观测点的平均扭曲方向上的位移量为横坐标，用MATLAB的三次样条插值的方法[1]，运用文献[6]中方法绘出主体挠度曲线图（见图5-图8）。

图5 古塔第一次测量的挠度曲线

图6 古塔第二次测量的挠度曲线

图7 古塔第三次测量的挠度曲线

图8 古塔第四次测量的挠度曲线

结论2：挠度曲线在截面位置坐标dertaki处的斜率，或挠度Hi对坐标dertaki的一阶导数，等于该截面的转角。从图中可以看出随着层数的增加挠度逐渐变大。说明塔身随着层数的增加导致扭曲程度越来越剧烈。而从曲率曲线可以看出同样的规律。

3.4 问题3 的分析结果：

通过层与层的中心位移量插值可知，每年层与层之间的弯曲程度相近，但是在1986年，塔尖比其它各层的弯曲程度大，1996年时，塔尖对第十三层的弯曲程度相对较大。2009年，第六层对第五层的弯曲程度相对较大，而2011年第六层对第五层的弯曲程度比较大。从曲率曲线图可看出1986年和1 996年均是在高度为50米处的曲率大，则说明塔尖的弯曲程度较大。2009年和2011年均是在高度为26米处的曲率大，表明在第六层处的弯曲程度大。

根据挠度得知曲线上某点的切线的倾斜角就是以该点为高度的平面的扭曲度，由挠度曲线图和古塔中心位移量在平均扭曲方向垂直面上的投影表分析得出：1986年，在接近50米处即古塔的12、13层的扭曲度较大；1996年，同样在50米左右处古塔的扭曲度较大；到了2009年时，挠度曲线变的平缓，则古塔扭曲的角度差不多；2011年的挠度曲线较2009年的挠度曲线没有太大差异。根据各层中心相对于第一层中心的扭曲方向表，发现θi的值均为负值，说明古塔是沿顺时针方向扭曲的。

根据各层中心相对于第一层中心的位移量表，我们可以看出每年的Δk 在逐渐地增大，说明古塔由下而上弯曲的程度越来越大。这四年里第8层的位移量相对于第7层的位移量和第12层的位移量相对于第11层的位移量差异不大，说明第8层对于第7层和第12层对于第11层的弯曲度很小。

总体上，随着层数的增加，塔身逐渐偏离z 轴。而且从四次模拟的各层中心我们可以发现：中心随着年份的增加在加速偏移z 轴。所以文物部门需要对该古塔实施抢救性保护措施，以配合我国重点古建筑的保护工作。

［1］任玉杰.数值分析及其M ATLAB 实现［M］.北京：高等教育出版社,2007.37-381.

［2］程大锦.建筑：形式、空间和秩序(第二版)［M］.天津：天津大学出版社,2005.67-78.

［3］龚纯, 王正林.M ATLAB 语言常用算法程序集［M］.北京：电子工业出版社,2011.101-176.

［4］何青,王丽芬.M aple 教程［M］.北京：科学出版社,2006.24-36.

［5］张磊,毕靖,郭莲英.M ATLAB 实用教程［M］.北京：人民邮电出版社,2008.45-49.

［6］Yanxing Cui, Chuanlong Wang, Ruiping Wen. On the convergence of generalized parallel multisplitting iterative methods for semidefinite linear systems［J］. Numerical Algebra, Control and Optimization,2012,2(4)：863-873.