基于Matlab的古塔變形分析

王雪萍， 楊慶生

（烏魯木齊職業(yè)大學(xué) 基礎(chǔ)教育部，新疆 烏魯木齊 830002）

0 引 言

古代建筑物中的古塔造型在我國許多地區(qū)隨處可見，一個時期一個地方的古塔建造體現(xiàn)一定的歷史文化背景；一些古塔歷經(jīng)幾個朝代，任憑風(fēng)吹雨打，巍然矗立至今，體現(xiàn)著我國古代建筑設(shè)計師和建造師的聰明智慧。研究古塔的結(jié)構(gòu)，分析古塔的變形對保護(hù)我國古建筑，研究我國歷史文化和建筑設(shè)計有著重要的價值和意義。

本問題是2013年全國大學(xué)生數(shù)學(xué)建模競賽乙組C題［1］，內(nèi)容是：由于長時間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會產(chǎn)生各種變形，諸如傾斜、彎曲、扭曲等。為保護(hù)古塔，文物部門需適時對古塔進(jìn)行觀測，了解各種變形量，以制定必要的保護(hù)措施。某古塔已有上千年歷史，是我國重點(diǎn)保護(hù)文物。管理部門委托測繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對該塔進(jìn)行了4次觀測。試給出確定古塔各層中心位置的通用方法，古塔各層的中心坐標(biāo)；分析該塔傾斜、彎曲、扭曲等變形情況；分析該塔的變形趨勢。

1 古塔形狀的確定和缺失數(shù)據(jù)的預(yù)測

1.1 對所給的古塔測量數(shù)據(jù)進(jìn)行分析，確定古塔形狀

假設(shè)每次測量觀測方法相同，用同一儀器和設(shè)備，在基本相同的環(huán)境和條件下觀測，所測量的數(shù)據(jù)都是真實(shí)、精確、可靠的。對附件中提供的某古塔1986年、1996年、2009年、2011年4次每層8個觀測點(diǎn)的測量數(shù)據(jù)進(jìn)行觀察，發(fā)現(xiàn)該古塔共有13層，每層有8個觀測點(diǎn)，塔尖有4個觀測點(diǎn)，且只有1986年和1996年的觀測數(shù)據(jù)，2009年和2011年數(shù)據(jù)缺失，故不考慮塔尖的變形。對給與的觀測數(shù)據(jù)的橫坐標(biāo)和縱坐標(biāo)應(yīng)用excel繪圖功能做散點(diǎn)圖，發(fā)現(xiàn)每層8個測量點(diǎn)近似呈八邊形分布，1986年和1996年的每層8個測量點(diǎn)觀測順序和2009年、2011年測量點(diǎn)觀測順序不一致，如圖1和圖2所示。

圖1 1986年第1層8個觀測點(diǎn)坐標(biāo)平面圖

圖2 2009年第1層8個觀測點(diǎn)坐標(biāo)平面圖

為了方便分析，我們從橫坐標(biāo)值最小的測量點(diǎn)開始編號，即按照2009年和2011年的測量點(diǎn)順序調(diào)整1986年和1996年每層8個觀測點(diǎn)的順序，即1986年和1996年各層的觀測點(diǎn)順序由12345678改為78123456。即古塔每層近似呈八邊形，且上小下大。1986年古塔空間立體圖如圖3所示。

圖3 1986年古塔空間立體圖

1.2 對1986年和1996年第13層觀測點(diǎn)3的缺失數(shù)據(jù)預(yù)測

由于所提供的1986年和1996年第13層觀測點(diǎn)3的坐標(biāo)數(shù)據(jù)缺失，為了便于分析，需要對缺失的數(shù)據(jù)進(jìn)行預(yù)測。1986年、1996年、2009年、2011年每層觀測點(diǎn)3的空間坐標(biāo)用Matlab7.0繪出的圖形如圖4所示。

圖4 1986年、1996年、2009年、2011年各層第3個觀測點(diǎn)空間圖形

可以從圖形直觀看出，古塔在第3個觀測點(diǎn)上，1986年和1996年、2009年和2011年測量的數(shù)據(jù)空間圖形幾乎一致，變化不大，但是，1996年到2009年這14年間古塔在第3個觀測點(diǎn)處變化比較大，在10層到13層上這4次變化趨勢是一致的，呈線性趨勢。

用xi，yi表示第i層的橫坐標(biāo)和縱坐標(biāo)，Δxi＝xi－xi－1，Δyi＝y(tǒng)i－yi－1分別表示相鄰兩層橫坐標(biāo)和縱坐標(biāo)的改變量［2］，則

對1986年和1996年第13層第3個觀測點(diǎn)的缺失的橫坐標(biāo)和縱坐標(biāo)用式（1）預(yù)測，結(jié)果見表1。

表1 1986年、1996年第13層觀測點(diǎn)3缺失數(shù)據(jù)預(yù)測 m

由表1可以看出，各層的絕對改變量Δxi和Δyi很小，相差1mm，這樣可以近似用第12層的改變量來推測第13層的坐標(biāo)。

1986年：

1996年：對于觀測點(diǎn)3的z3值，如果把各層看做在同一平面上的話，運(yùn)用均值法，即用該塔第13層其余7個觀測點(diǎn)的豎坐標(biāo)的平均值來作為z3的預(yù)測值，則

以此類推：得出1986年第13層測量點(diǎn)3坐標(biāo)（567.984，519.588，52.834）、1996年第13層測量點(diǎn)3坐標(biāo)（568.1，519.581 6，52.83）。

另外一種方法可以應(yīng)用Matlab軟件做線性擬合。

1986年第13層觀測點(diǎn)3的x坐標(biāo)預(yù)測，從第10層開始取x坐標(biāo)值，Matlab命令：

＞＞tr＝［10 11 12；568.148 568.094 568.039］；

＞＞t＝tr（1，：）＇；

＞＞x＝tr（2，：）＇；

＞＞p＝polyfit（t，x，1）

則橫坐標(biāo)x的擬合函數(shù)x＝568.693－0.054 5t，把t＝13代入，x＝567.984 67。

縱坐標(biāo)y擬合函數(shù)y＝517.333＋0.173 5t，把t＝13代入，y＝519.588 667。

兩種方法預(yù)測的x和y坐標(biāo)誤差分別是0.67mm和0.667mm，非常小，故這兩種方法用來預(yù)測缺失數(shù)據(jù)都可行。

2 古塔變形的分析

2.1 各層中心坐標(biāo)的確定

古塔建造初期應(yīng)該是一上小下大的正八邊形空間立體，除了第1層外，其它各層是一個截面為八邊形的正棱臺，經(jīng)過上千年，古塔發(fā)生了變形，不再是一個規(guī)則的正八邊形，由于所提供的數(shù)據(jù)是每層8個點(diǎn)的測量坐標(biāo)，為了簡化問題，方便問題的解決，我們把每層看做一個八邊形正棱柱。正八邊形內(nèi)接于圓，內(nèi)接圓的圓心即可作為所在層的中心，可以用正八邊形擬合圓的方程，則圓心坐標(biāo)是各層中心的橫坐標(biāo)和縱坐標(biāo)。第k層擬合的圓方程：

根據(jù)所給各觀測點(diǎn)的橫坐標(biāo)和縱坐標(biāo)用Matlab軟件擬合圓的方程。1986年第1層的中心橫坐標(biāo)和縱坐標(biāo)的Matlab命令［3－4］：

則可得擬合圓心坐標(biāo)（566.665，522.709）。假設(shè)測量點(diǎn)位于各層的頂部，則各層中心應(yīng)該位于相鄰兩層的中間位置，第1層中心的豎坐標(biāo)取層高的平均值的二分之一，第k層中心豎坐標(biāo)記為zk0，其中

則可得1986年第1層的中心坐標(biāo)為（566.665，522.709，0.894），同 理，可 以 計 算 出1986年第2層的中心坐標(biāo)（566.723，522.671，4.554）。

另外，計算各層中心的簡單方法可以用各層8個測量點(diǎn)的橫坐標(biāo)和縱坐標(biāo)的平均值來近似表示各層的中心坐標(biāo)，即

則可得1986年第1層的中心坐標(biāo)為（566.665，522.711，0.894），同理，可計算出1986年第2層中心坐標(biāo)（566.720，522.668，4.554），兩種方法計算的中心坐標(biāo)偏差均小于0.005m，因此兩種方法計算中心坐標(biāo)都是可行的。

2.2 古塔的傾斜、彎曲、扭曲分析

2.2.1 古塔的傾斜分析

如果古塔沒有傾斜，各層的中心坐標(biāo)應(yīng)該在古塔中心軸上，即其投影在xoy平面上的同一點(diǎn)處，根據(jù)建筑施工誤差標(biāo)準(zhǔn)，軸線位置偏移＜10mm是允許的［5－8］。由1986年各層的中心坐標(biāo)可以看出，古塔的中心坐標(biāo)不在同一點(diǎn)處，做1986年各層中心在xoy面的投影坐標(biāo)圖，如圖5所示。

中心坐標(biāo)在x軸上的極差0.542m，在y軸上的極差0.348m，且橫坐標(biāo)隨著層數(shù)的增加遞增，縱坐標(biāo)隨著層數(shù)的增加遞減，即向x軸正向傾斜。記第1層與第13層中心的投影距離為l1，13，塔身與xoy平面的夾角α，則

即1986年塔的傾斜角為89.31°。同理，可計算出1996年塔的傾斜角為89.303°；2009年塔的傾斜角為89.251°；2011年塔的傾斜角為89.246°?？梢钥闯?，1986年和1996年塔的傾斜變化不大，1996年到2009年間塔的傾斜變化較大，2009年和2011年間塔的傾斜變化不大。

圖5 1986年古塔各層中心坐標(biāo)在xoy平面的投影

對各層8個觀測點(diǎn)的高計算層高差、平均層高、層高均方差，見表2。

表2 1986年各層8個觀測點(diǎn)層高數(shù)據(jù)統(tǒng)計

對于1986年的古塔各層8個觀測點(diǎn)的z坐標(biāo)進(jìn)行分析，得出各層高的變化趨勢，第1層到第5層極差較小，傾斜程度不大，第6層到第10層極差較大，傾斜較大，第11層和第12層次之，第13層傾斜最大。

其它3年的分析方法相同。

2.2.2 古塔的彎曲與扭曲分析古塔建造初期應(yīng)該是一個正八棱臺，每層8個觀測點(diǎn)在xoy平面的投影應(yīng)該是均勻分布的直線，如果彎曲，則在xoy平面的投影不再是均勻分布的直線，如果沒有發(fā)生扭曲，則在xoy平面各層的觀測點(diǎn)的投影的距離就是非均勻分布。如果發(fā)生扭曲，則投影就不在同一條直線上。1986年古塔13層8個觀測點(diǎn)空間曲線和每層8個觀測點(diǎn)在xoy平面的投影分別如圖6和圖7所示。

圖6 1986年古塔各層8個觀測點(diǎn)空間線圖

圖7 1986年古塔各層8個觀測點(diǎn)在xoy平面投影點(diǎn)

從1986年各層8個點(diǎn)的空間線圖可以直觀看出，1～6層變化不大，從第7層開始發(fā)生彎曲，第7層正好處于該塔的中心位置，可能受重力的作用開始發(fā)生彎曲。

1986年每層8個點(diǎn)在xoy平面上的投影不在一條直線上，可以看出發(fā)生了扭曲，觀測點(diǎn)7，8從第7層開始向x軸正向發(fā)生扭曲，觀測點(diǎn)6到第10層扭曲較大，觀測點(diǎn)5扭曲不明顯，觀測點(diǎn)4和觀測點(diǎn)3、觀測點(diǎn)2都從第7層開始扭曲，其中觀測點(diǎn)3扭曲較嚴(yán)重，觀測點(diǎn)1扭曲不明顯。

2.2.3 古塔變形趨勢分析

根據(jù)測量的坐標(biāo)數(shù)據(jù)，用Matlab7.0做出1986年、1996年、2009年、2011年觀測點(diǎn)1的空間曲線如圖8所示。

圖8 1986年、1996年、2009年、2011年觀測點(diǎn)1的空間曲線

從圖形觀察可以看出，1986年和1996年這兩次測量，古塔的變形不大，2009年和2011年兩次測量古塔變形也不大，但在1996年和2009年期間，也許發(fā)生了地震、颶風(fēng)等原因，古塔發(fā)生了很大的變形。從1986年和1996年這兩次測量可以看出，從第7層開始發(fā)生了彎曲，其它形變不明顯，而2009年，從第1層開始到第5層觀測點(diǎn)1發(fā)生了位移和整體傾斜，向著y軸正向傾斜，從第6層開始到7，8，9層發(fā)生了彎曲和扭曲，從第10層開始再一次發(fā)生了彎曲和扭曲現(xiàn)象，之后的2011年測量，幾乎變化不大。結(jié)合圖4的第3個觀測點(diǎn)的空間圖形，可以看出變形情況基本一致。

3 結(jié) 語

應(yīng)用Matlab7.0對某古塔4次測量數(shù)據(jù)進(jìn)行分析，給出了每層中心坐標(biāo)的2種算法及塔與xoy平面的傾斜角，結(jié)合觀測點(diǎn)在平面的投影和空間曲線圖，得出1986年和1996年兩次測量，古塔變化不大，1996年到2009年期間，古塔發(fā)生了嚴(yán)重的變形，2011年相對于2009年古塔的變化不大，說明1996年至2009年期間有可能發(fā)生了地震、颶風(fēng)等自然災(zāi)害。在沒有地震、颶風(fēng)等發(fā)生的情況下，由于受重力等影響，古塔易發(fā)生傾斜，在其高度的1／2處易發(fā)生較明顯的彎曲，在2／3處也易發(fā)生輕微彎曲；由于受風(fēng)力和風(fēng)向影響，古塔在其高度1／2之上處容易發(fā)生扭曲，高度的1／2處之下一般比較穩(wěn)固，對于文物保護(hù)部門來說，古塔高度的1／2處及以上樓層要多加關(guān)注，及時采取保護(hù)措施。

［1］ 全國大學(xué)生數(shù)學(xué)建模競賽網(wǎng).2013年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽賽題［EB／OL］.［2014－12－22］.http：／／www.mcm.edu.cn／.

［2］ 李麗，王振領(lǐng).MATLAB工程計算及應(yīng)用［M］.北京：人民郵電出版社，2003.

［3］ 同濟(jì)大學(xué)數(shù)學(xué)教研室.高等數(shù)學(xué)［M］.6版.北京：高等教育出版社，2007.

［4］ 蘇金明，阮沈永.MATLAB實(shí)用教程［M］.北京：電子工業(yè)出版社，2005.

［5］ 卜璞，向亭譯.用整體方差分析法對變形監(jiān)測數(shù)據(jù)進(jìn)行處理［J］.科學(xué)技術(shù)與工程，2012，12（2）：407－409.

［6］ 丁寧，孫英俊.高層建筑物變形監(jiān)測數(shù)據(jù)處理與分析［J］.測繪科學(xué)，2011，36（5）：93－95.

［7］ 黃強(qiáng).古塔變形檢測的探討［J］.測繪與空間地理信息，2013，36（6）：217－220.

［8］ 張方程，李曉天.同心圓亞像素中心定位方法［J］.長春工業(yè)大學(xué)學(xué)報：自然科學(xué)版，2014，35（5）：500－505.