常金勇，崔艷星，李海增

（長治學院 數(shù)學系，山西 長治 046011）

1 問題重述

某測繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對某千年古塔進行了四次觀測。所得數(shù)據(jù)有助于對該塔的傾斜、彎曲、扭曲等變形情況加以研究。問題可以歸結(jié)為：

(1)通過找出確定古塔各層中心位置的通用方法，給出古塔的各層中心坐標列表；

(2)再分析該塔傾斜、彎曲、扭曲等變形情況；

(3)最后分析該塔的變形趨勢。

2 問題的分析

對于問題（1），我們可以從古塔的每層上隨機取n個點，繪出塔的輪廓曲線，然后利用MATLAB插值方法[1]37-381求出該曲線的中心，最終找出計算各層中心位置的通用方法，進而計算出古塔各層及塔尖的中心坐標。

對于問題（2），我們從建筑測量學來定義傾斜、彎曲、扭曲等概念[2]67-78。對于古塔這種底面積小、高度不高的建筑，我們用測水平角法[2]來度量它的傾斜。

撓度是指結(jié)構(gòu)構(gòu)件的軸線或中面由于彎曲引起垂直于軸線或中面方向的線位移。求解方法是通過求出中心位移量、以及中心位移量在平均扭曲方向垂直面上的投影，進一步繪出撓度曲線圖來直觀給出古塔的彎曲與扭曲程度。

對于問題（3）我們可以根據(jù)四次各層的中心坐標用MATLAB[3]101-176繪出插值圖，根據(jù)圖形來分析該塔的彎曲情況。

3 模型的建立與求解

3.1 問題1 的解決方案及模型

在每層取出n個點，分別計算這n個點X,Y,Z坐標的平均值。最后我們得到了求各層中心坐標的通用方法：

我們在實際操作中運用Maple程序[4]24-36計算了各層的中心坐標，以及古塔各層的中心坐標變化。

3.2 問題2 的解決方案及模型

根據(jù)建筑測量學知識，古塔主體幾何中心就是其底面的中心。要求觀測塔尖相對于底面中心的平面位移ΔDi，塔尖的中心與第一層中心的距離hi，由ΔDi與hi的比值來計算出tanθi，運用反三角函數(shù)知識來求出每次測量的古塔傾斜弧度δi，再用MATLAB程序[5]求出傾斜角度θi，就可以確定古塔的傾斜情況。（見下面公式）

古塔的傾斜角度結(jié)果如下表：

?

結(jié)論1：在前十年間該古塔傾斜度變化了0.0103度，平均每年會傾0.00103度。在接下來的十三年里古塔的傾斜度又變化了0.0256度，則平均每年古塔會傾斜0.00197度。在最后的兩年里，古塔傾斜了0.0020度，平均每年會傾斜0.0010度。所以古塔開始的一段時間內(nèi)傾斜的速度相對較慢。在后來的一段時間里，古塔的傾斜速度開始變快，最近幾年內(nèi)古塔的傾斜速度又變的相對緩慢。所以文物部門需適時對古塔進行觀測了解各種變形量，以制定必要的加固措施。

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 古塔第四次測量的撓度曲線

結(jié)論2：撓度曲線在截面位置坐標dertaki處的斜率，或撓度Hi對坐標dertaki的一階導數(shù)，等于該截面的轉(zhuǎn)角。從圖中可以看出隨著層數(shù)的增加撓度逐漸變大。說明塔身隨著層數(shù)的增加導致扭曲程度越來越劇烈。而從曲率曲線可以看出同樣的規(guī)律。

3.4 問題3 的分析結(jié)果：

通過層與層的中心位移量插值可知，每年層與層之間的彎曲程度相近，但是在1986年，塔尖比其它各層的彎曲程度大，1996年時，塔尖對第十三層的彎曲程度相對較大。2009年，第六層對第五層的彎曲程度相對較大，而2011年第六層對第五層的彎曲程度比較大。從曲率曲線圖可看出1986年和1 996年均是在高度為50米處的曲率大，則說明塔尖的彎曲程度較大。2009年和2011年均是在高度為26米處的曲率大，表明在第六層處的彎曲程度大。

根據(jù)撓度得知曲線上某點的切線的傾斜角就是以該點為高度的平面的扭曲度，由撓度曲線圖和古塔中心位移量在平均扭曲方向垂直面上的投影表分析得出：1986年，在接近50米處即古塔的12、13層的扭曲度較大；1996年，同樣在50米左右處古塔的扭曲度較大；到了2009年時，撓度曲線變的平緩，則古塔扭曲的角度差不多；2011年的撓度曲線較2009年的撓度曲線沒有太大差異。根據(jù)各層中心相對于第一層中心的扭曲方向表，發(fā)現(xiàn)θi的值均為負值，說明古塔是沿順時針方向扭曲的。

根據(jù)各層中心相對于第一層中心的位移量表，我們可以看出每年的Δk 在逐漸地增大，說明古塔由下而上彎曲的程度越來越大。這四年里第8層的位移量相對于第7層的位移量和第12層的位移量相對于第11層的位移量差異不大，說明第8層對于第7層和第12層對于第11層的彎曲度很小。

總體上，隨著層數(shù)的增加，塔身逐漸偏離z 軸。而且從四次模擬的各層中心我們可以發(fā)現(xiàn)：中心隨著年份的增加在加速偏移z 軸。所以文物部門需要對該古塔實施搶救性保護措施，以配合我國重點古建筑的保護工作。

［1］任玉杰.數(shù)值分析及其M ATLAB 實現(xiàn)［M］.北京：高等教育出版社,2007.37-381.

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

［3］龔純, 王正林.M ATLAB 語言常用算法程序集［M］.北京：電子工業(yè)出版社,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.