付小娟 吳洪坤

(1.廣州民航職業(yè)技術學院 人文社科學院，廣東 510403;2.廣州民航職業(yè)技術學院 飛機維修工程學院，廣東 510403)

1 引言

某古塔已經有上千年歷史，是我國重點保護文物。但由于古塔長時間承受自重、氣溫、風力等各種作用，偶然還要受地震、颶風的影響，古塔會產生各種變形，諸如傾斜、彎曲、扭曲等。為保護古塔，文物部門適時對古塔進行觀測，每次都得到一組數(shù)據(jù)，以制定必要的保護措施?，F(xiàn)根據(jù)2013 年全國大學生數(shù)學建模競賽C 題附錄1 的數(shù)據(jù)，提出確定古塔各層中心位置的通用方法，并列表給出各次測量的古塔各層中心坐標。分析該塔傾斜、彎曲、扭曲等變形情況。

2 建立模型

2.1 問題分析

古塔在傾斜、彎曲、扭曲等變形中，每一層的層面結構、塔心坐標等都發(fā)生了變化，為方便起見，首先只計算各次測量的第一層塔心，把每次測量得到的第一層的測試點作為研究對象，并且按照題目給出的順序編號為1，2，3.....8 號點，為了解八個點的相對位置，利用matlab 畫出四次測試第一層八個點的空間圖形，如圖1 所示，發(fā)現(xiàn)測試點不在一個水平面上，組成一個類似于八邊形的結構，經計算八條連線的長度如表1 所示，每條連線是不等長的，所以測試點并不是均勻取定的，但是前兩次測試的數(shù)據(jù)比較接近，各條連線對應幾乎相等，后兩次測試的數(shù)據(jù)也比較接近，各條連線也對應幾乎相等，由此推斷，前兩次測試點是一一對應的，而后兩次測試點也是一一對應的。

圖1 第一層坐標圖.藍:1986，綠:1996，紅:2009，黑:2011

為了組成一個平面多邊形，設想把測試點投影到某個平面上，用投影八邊形的重心坐標表示塔心坐標，由于各次測試點并不在一個水平面上，所以采用空間傾斜平面作為投影面，首先利用各層八個測試點擬合一個空間平面，如圖2 表示1986 年測試第一層八個測試點擬合的空間平面，然后把測試點投影到這個空間平面上，這樣就可以得到一個空間平面上的八邊形［1］，利用matlab 求得八邊形的重心，就得到了各層的塔心坐標。

根據(jù)問題一求出的四組塔心的坐標，算出頂層與底層的傾斜位移，進而求出傾斜量。再用各次測量每一層的塔心坐標擬合出四條空間曲線，用曲線的曲率表示塔的彎曲情況。塔的扭曲度可以利用同一個測試點相對于塔心發(fā)生的扭轉來描述。通過前面變形數(shù)據(jù)的分析，可以得到塔在以后的變形趨勢。根據(jù)塔的變形趨勢，我們提出了幾點建議。

圖2 1986 年測試第一層坐標點擬合的空間平面圖

2.2 模型建立

2.2.1 塔心坐標模型的建立與求解

圖1 所示第一層的八個測試點不在一個水平面上，把1986 年測得的原始數(shù)據(jù)繪成圖3，發(fā)現(xiàn)形如塔狀，其他各層測試點也不在一個水平面上;所以采用先根據(jù)測試點擬合一個空間平面Π，Π 的方程常用公式(1)表示，當有n 個測試點時，要擬合這個平面，可以表示成矩陣［2］方程(2)的形式，利用matlab［3］求解出各層測試點擬合的平面方程，把每次測的各測試點投影到這個擬合平面上，設空間中任意一個點是B1(x1，y1，z1)，其在平面Π 上的投影點設為B2(x2，y2，z2)，則投影方程如公式(3)所示，算出k 值，帶入(3)式得到各個測試點的投影點坐標，連成投影八邊形，求其重心坐標。

多邊形的重心計算［4］，設有圖4 所示的n 邊形A1A2…Ai…An，各頂點坐標分別是Ai(xi，yi，zi)(i=1，2...n)，連接AnA2，AnA3，…AnAi，…，AnAn-2，得到n-2 個三角形，設各個三角形的重心坐標分別是Gi()(i=1，2...n -2)，其中每個坐標的計算如公式(4)所示。

設第i 個三角形的面積是σi，先利用兩點之間距離公式，求出第i 個三角形的三邊長ai，bi，ci，再套入公式(5)求出σi。設多邊形的重心坐標為O ()，計算公式如公式(6)所示。

表1 四次測試第一層八個測試點之間的距離

由公式(1)到(6)，利用matlab 算出第一、七、十三和塔尖的塔心坐標，如表2 所示。

圖3 1986 年測試各層測試點空間圖

圖4 平面多邊形重心坐標的計算

表2 中數(shù)據(jù)與測試的原始數(shù)據(jù)比較吻合，非常接近各層測試值，1986 年和1996 年的測試數(shù)據(jù)及計算數(shù)據(jù)都更接近，說明前兩次的測試點是一一對應的，而且這10 年間，古塔變形情況比較小，維護較好;而2009 年和2011 年的數(shù)據(jù)更加接近，這說明后兩次的測試點也是一一對應的，說明這期間的維護比較頻繁，古塔變形較小;但是從1996 年到2009 年這15 年間，發(fā)現(xiàn)古塔的變形比較大，同一層的塔心橫坐標x 的變化范圍是(0 -0.4961)m，縱坐標y 的變化范圍是(0 -0.6551)m，豎坐標的變化范圍是(0 -0.0146)m。同一層塔心的橫坐標值越來越大，說明古塔的傾斜朝著x 軸的正方向;縱坐標的變化規(guī)律比較復雜，這和古塔的扭曲變形嚴密相關;同一層塔心的豎坐標z 的值隨著時間的推移變得越來越小，這表明了古塔的實際彎曲變形。

表2 四次測試各層的塔心坐標

2.2.2 塔的傾斜度模型

塔的傾斜是由于基礎立柱頂面高低不平而引起塔中心偏離鉛垂線位置的現(xiàn)象，采用頂層塔心相對于底層塔心的水平位移變化來刻畫古塔的傾斜［2］［5］。設古塔第一層的塔心坐標為(x1，y1，z1)，頂層的塔心坐標為 (x14，y14，z14)，塔底到塔頂塔心的高度為H，則傾斜位移δ 如公式(7)所示，傾斜度β 如公式(8)。計算結果如表3所示。

該古塔的傾斜度隨著時間的推移越來越大，查詢數(shù)據(jù)資料可得當傾斜度小于0.004 是正常允許的，但是表四中數(shù)據(jù)都大于0.004，所以該古塔需要盡快進行修繕，在傾斜方向加固立柱。

表3 國外主流BIM 服務器對比

2.2.3 塔的彎曲度模型

把同一年的13 個塔心坐標擬合出一條空間曲線，共可以擬合出四條曲線如圖5 所示。塔的彎曲度可以利用各層塔心坐標連線的曲率來刻畫，因此，先把各層塔心坐標分別投影到XOZ 面(圖6)和YOZ 面上(圖7)，利用投影點擬合成一個圓(9)，因為有多個塔心，所以把圓方程寫成矩陣方程(10)

利用matlab 計算出a，b，c，由此找到半徑R如公式(11)，曲率ρ=1/R。由此可以計算出曲線在X=0 和Y=0 兩個投影面上的綜合曲率(如表4)。發(fā)現(xiàn)1986 年和1996 年彎曲度基本相同，2009 年和2011 年也相同，塔朝著x 軸正方向的的彎曲度從1996 年到2011 年有很大遞減，而朝著y 軸負方向的彎曲度卻有所遞增，說明在這期間對塔的彎曲變形采取了一定保護措施，使之逐漸回歸直立的狀態(tài)。

表4 古塔彎曲的曲率年份

2.2.4 塔的扭曲度

平面的扭曲變形定義為“所考慮平面的兩個對邊旋轉角度之差隨其距離的變化率”來描述【6】。采用同一個測試點與塔心連線發(fā)生的扭轉角度來計算塔的扭曲度，根據(jù)前邊證明1986年和1996 年的測試點是一一對應的，設1986年第一層塔心坐標為O，第i 次測量第一層的第j 個測試點為Aij，連接，求出兩條線之間的夾角，然后求另外七對測試點與O 點連線的夾角，……，，算出八個夾角的平均值，表示古塔第一層從1986 年到1996 年所發(fā)生的扭曲度，同理可求其他扭曲度，篇幅所限，本文不做計算。

圖5 各層塔心擬合空間曲線

圖6 塔心在XOZ 面投影曲線圖

圖7 塔心在YOZ 面投影曲線

3 結論

該古塔的傾斜變形比較嚴重，而且逐年增長，建議今后加強傾斜檢測與防護;彎曲變形比較小，說明在彎曲防護方面做得比較好。古塔監(jiān)測內容的實施以及頻率按各塔的健康程度和委托方的要求而定，監(jiān)測周期根據(jù)現(xiàn)場實際情況可適當調整，建議今后每隔3 -4 年進行一次古塔平面和高程的監(jiān)測，如發(fā)現(xiàn)周邊建筑或者地鐵等施工，地質狀況發(fā)生劇烈變化，可進行加密測量，以便監(jiān)測出古塔隨時間的變化產生變形量，從而保證對古塔的安全維護提供最新最可靠的數(shù)據(jù)和更好地進行古塔保護。

［1］黃強.古塔變形監(jiān)測的探討［J］.測繪與空間地理信息，2013，36(6):217 -220.

［2］胡志曉.古塔傾斜觀測和數(shù)據(jù)分析［J］.江蘇建筑，2011(6):34-36.

［3］周博，謝東來等.MATLAB 科學計算［M］.北京，機械工業(yè)出版社.

［4］郭幼操.從四邊形重心到多邊形的重心［J］.浙仁農村技術師專學報，1996(1-2):30-33.

［5］余周佑.安慶振風塔傾斜測量與數(shù)據(jù)分析［J］.施工技術研究與應用，2003(5):34.

［6］韓煊，Jamie R Standing，李寧.地鐵施工引起的建筑物扭曲變形分析［J］.土木工程學報，2010，43(1):82-88.