古塔觀測(cè)數(shù)據(jù)的計(jì)算和變形分析

付小娟 吳洪坤

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

1 引言

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

2 建立模型

2.1 問(wèn)題分析

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

圖1 第一層坐標(biāo)圖.藍(lán):1986，綠:1996，紅:2009，黑:2011

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

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

圖2 1986 年測(cè)試第一層坐標(biāo)點(diǎn)擬合的空間平面圖

2.2 模型建立

2.2.1 塔心坐標(biāo)模型的建立與求解

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

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

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

表1 四次測(cè)試第一層八個(gè)測(cè)試點(diǎn)之間的距離

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

圖3 1986 年測(cè)試各層測(cè)試點(diǎn)空間圖

圖4 平面多邊形重心坐標(biāo)的計(jì)算

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

表2 四次測(cè)試各層的塔心坐標(biāo)

2.2.2 塔的傾斜度模型

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

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

表3 國(guó)外主流BIM 服務(wù)器對(duì)比

2.2.3 塔的彎曲度模型

把同一年的13 個(gè)塔心坐標(biāo)擬合出一條空間曲線，共可以擬合出四條曲線如圖5 所示。塔的彎曲度可以利用各層塔心坐標(biāo)連線的曲率來(lái)刻畫，因此，先把各層塔心坐標(biāo)分別投影到XOZ 面(圖6)和YOZ 面上(圖7)，利用投影點(diǎn)擬合成一個(gè)圓(9)，因?yàn)橛卸鄠€(gè)塔心，所以把圓方程寫成矩陣方程(10)

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

表4 古塔彎曲的曲率年份

2.2.4 塔的扭曲度

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

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

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

圖7 塔心在YOZ 面投影曲線

3 結(jié)論

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

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

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

［3］周博，謝東來(lái)等.MATLAB 科學(xué)計(jì)算［M］.北京，機(jī)械工業(yè)出版社.

［4］郭幼操.從四邊形重心到多邊形的重心［J］.浙仁農(nóng)村技術(shù)師專學(xué)報(bào)，1996(1-2):30-33.

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

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