塔形建筑物輪廓點Z坐標點云傾斜監測方法

屈乾龍，朱慶偉，艾衛濤，溫 波，馬肇祥

(1.西安科技大學 測繪科學與技術學院，陜西 西安 710054；2.河南城建學院 測繪與城市空間信息學院，河南 平頂山 467000)

0 引 言

隨著時間的推移，古代高塔如今都是“十塔九斜”，塔體受基礎和結構本身的影響及諸多外界因素而產生彎曲變形，綜合表現為整個建筑體傾斜變形[1-3]。傾斜度是衡量施工技術質量和后期建筑物安全維護的重要指標，如何快速準確測量塔類建筑物傾斜度一直都是變形觀測領域的研究熱點。三維激光掃描技術以其無需直接接觸即可快速獲取物體表面高精度三維點云數據的特點被廣泛應用于塔類建筑物的變形監測[4-7]。國內對塔形建筑物的變形監測研究從未停止過，梁華等基于重慶某鐵塔點云數據，對鐵塔結構底部特征進行圓柱擬合獲取鐵塔中心坐標，再分別提取關鍵部位點云，利用特征擬合獲取幾何參數，計算鐵塔上、中、下端中心到鐵塔軸線偏差距離[8]。辛星等基于三維激光掃描技術對塔形構筑物進行傾斜監測，利用Imageware軟件對點云數據提取構筑物中軸線，計算出構筑物的傾斜值。最后與全站儀前方交會精度實驗對比驗證三維激光掃描技術在塔型構筑物傾斜監測中應用的可靠性[9]。王莉等以銅川市延昌塔為研究對象，采用三維激光掃描技術獲取塔體點云數據，利用Geomagic軟件進行建模，然后利用MATLAB軟件對點云數據擬合計算，得到塔體中軸線，繼而分析判斷塔體的傾斜變化狀況[10]。楊永林等以西安萬壽寺塔為例，基于三維激光掃描技術對塔體進行變形監測，提出最小二乘塔體傾斜度計算方法。該方法通過提取每層塔體的輪廓點云，根據最小二乘原理擬合出每層塔體的圓心，對每層圓心進行線性擬合提取塔體中軸線，繼而計算出塔體的傾斜量[11]。

以上學者在采用三維激光掃描技術做傾斜監測時，無論是橫斷面提取中心軸還是提取輪廓點擬合中心軸，均采用擬合中心軸的方法，數據處理方案較為單一，且均是將塔體看作一個剛體計算整體傾斜量，不能逐層計算。筆者提出一種點云輪廓點Z坐標傾斜監測數據處理方案，與提取中心軸法共同對西安大雁塔點云數據進行傾斜量計算，最后將計算的傾斜量進行對比。計算結果顯示，點云數據處理方法對塔體建筑物傾斜監測具有可行性和工程實際意義。

1 研究對象概括

大雁塔位于中國陜西省西安市南郊大慈恩寺內，大雁塔是磚仿木結構的四方形樓閣式磚塔，由塔基、塔身、塔剎組成，現通高為64.517 m。塔基高4.200 m，南北約48.700 m，東西45.700 m；塔體呈方錐形，平面呈正方形，底邊長為25.500 m，塔身高59.900 m，塔剎高4.870 m[12]。據記載，陜西測繪部門自1985年開始對其進行測量，其時大雁塔的傾斜速度處于加快過程中，到1996年左右塔傾斜達到最大程度，傾斜度達到1.011 m。如今，具有1 300多年歷史的大雁塔整體屬于動態平衡。大雁塔歷史悠久且傾斜明顯，故以大雁塔為例驗證文中方法的可行性，對大雁塔的安全運營及維修加固具有一定的參考意義。

2 點云數據獲取與處理

2.1 數據來源

三維激光掃描數據采集的誤差來源大致分為3種：與目標物反射面有關的誤差、儀器誤差、外界環境條件。王俊杰等通過以Trimble GX三維激光掃描儀為例研究數據采集誤差對數據精度的影響[13]。最終實驗得知測距誤差對點位精度影響微乎其微；由于掃描角的誤差，點位誤差隨測量距離的增大而增大，可通過點位精度估計公式來計算出掃描最大點位誤差，根據實際工作要求的點位精度來選取適當的掃描距離。文中數據是于2016年所測大雁塔點云數據。該數據是由Scan Station C10掃描儀采集的點云數據。設站情況如圖1所示，站點分布在塔的四周。近距離站點共8站，主要是獲取塔的中下部點云數據；遠距離的站點共4站，主要是獲取塔的頂部點云數據[12]。數據采集時，X軸在橫向掃描面內，Y軸在橫向掃描面內與X軸垂直，Z軸與橫向掃描面垂直。文中所涉及的所有Z值均是指Z軸坐標值。

圖1 設站分布

2.2 數據預處理

在獲得原始點云數據后，第一步要對多站掃描數據進行拼接配準處理。由于三維激光掃描過程中難免會掃描到塔體周圍各種地物而產生巨量的噪聲點，為了防止拼接過程中因為海量點云數據而導致軟件崩潰以及提高各站點云數據拼接速度，在原始數據導入Leica Cyclone軟件后，對每站除待測塔體和標靶球區域外大部分噪聲點云，直接在Cyclone軟件中進行區域裁剪，快速除去大量噪聲點云，這樣后續點云數據拼接和融合速度大大提高[14-16]。然后逐站標記同名點和標靶球，通過標靶球和同名點混合拼接方法進行數據拼接。為了減少個別標靶球或同名點的誤差傳播，采用每2站數據拼接一次，然后再把拼接好的結果進行拼接[17]，嚴格把控各站拼接精度，對于拼接誤差較大的測站，將該站標靶球和同名點刪除后重新標記同名點再進行拼接，直至拼接誤差在0.006 m左右。最終將各站基于儀器內部獨立坐標系統的點云數據轉化到同一坐標系中，形成一個整體。

拼接成一個整體后，仍存在大量的重復點云，所以需要通過設置平均采樣間隔或其它參數設置來進行點云數據融合處理，去除大量重復冗余點云[18-20]。

3 傾斜監測方法

3.1 已有三維激光掃描傾斜監測方法

文獻[9]提出基于中軸線上節點坐標偏移的方法，其思路是提取中軸線上節點的坐標[21-22]?；玖鞒倘鐖D2所示。

圖2 切片中心軸傾斜監測基本流程

文獻[10]利用Geomagic軟件進行建模，構建OBJ三維格網模型；再通過軟件提取每層塔體特征點坐標，計算出該層中心點坐標，利用 MATLAB軟件擬合計算塔的中心軸線，最終計算分析塔體傾斜變化狀況，基本原理如圖3所示。

圖3 擬合中心軸計算原理

文獻[11]提出基于最小二乘的塔體中軸線提取方法，其思路是首先提取出每層塔體外輪廓點云的圓心坐標，然后基于最小二乘原理對所有層圓心坐標進行線性擬合獲取塔體中軸線的方程，最后計算2期塔體的傾斜度，實現塔體的傾斜監測。

以上方法無論是對塔體進行區域分割還是基于塔體塔檐提取特征輪廓點，最后都是基于點云數據提取塔體中心軸的方法來確定塔體傾斜度，文中所提出的方法與該類方法有所不同，不用提取中心點即可計算出塔體傾斜度。

3.2 輪廓點Z坐標傾斜監測方法

筆者提出基于輪廓點云Z坐標提取計算法進行傾斜監測?；舅悸肥菙M合提取塔體固定輪廓特征點坐標，比如每一層的塔檐。如圖4所示，提取出每層塔檐外輪廓點云的坐標，A，B，C，D是塔檐的4個角，我們只需在提取出的外輪廓點云坐標中篩選出其Z坐標值最大點D和Z坐標值最小點B，此時LDO為D、B這2點Z坐標之差ΔZ，LDB可根據2點間距離公式求出，最后利用正弦定理公式(1)和反三角函數公式(2)即可算出該2點連線LDB與水平面夾角∠OBD，該角在數值上等于塔體傾斜角。

圖4 基于輪廓點Z坐標傾斜監測法原理

(1)

又∠DOB=90°，故公式(1)可化成

(2)

如圖5所示，這種是塔體傾斜特殊情況，塔體塔檐連線嚴格向另一側傾斜，這時LAD，LBC線上各點的Z坐標總體相同，此時可直接選取同側塔角點A，B，或者點D，C為Z坐標最大值和最小值點計算∠OEF。

圖5 塔體傾斜特殊情況

本方法同樣適用于圓形或六邊形等規則塔體建筑，原理相同，在此不再贅述。

4 實驗數據計算對比分析

4.1 擬合提取中心軸線傾斜監測法

擬合提取中心軸線傾斜監測法具體操作在此不再詳細介紹，基本操作流程如圖6所示。用Leica Cyclone軟件將點云數據預處理后，將拼接后的點云數據導入Imageware軟件中進行區域等距離切割點云，并將切割后的點云輪廓數據提取出來[23-25]，為了提高精度，準確計算出塔體各中心坐標，因此首先將輪廓點云擬合成等間距圓，如圖7所示，共提取7條輪廓線。

圖6 中心軸法實驗流程

圖7 切割輪廓點提取擬合圓

再計算出各層中心點圓心坐標并將計算得出的坐標和輪廓點坐標數據導入軟件具體效果如圖8所示。

圖8 頂點及各層中心點效果

計算出各層中心點坐標導入MATLAB軟件顯示后明顯發現各中心點與塔頂中心點并不在同一直線上，這是因為塔體建筑傾斜形變并不是線性形變的，因此中心軸法需要將塔體作為一個剛體整體分析。將中心點坐標和塔頂點坐標通過MATLAB編寫程序，基于最小二乘擬合原理提取塔體中軸線，結果如圖9所示。

圖9 塔體中心圓點擬合直線(單位：m)

最終得到擬合后頂點和各層中心點坐標，見表1。

表1 擬合后頂點和各層中心點坐標

最終通過擬合中心坐標計算，該年大雁塔塔體大約整體傾斜0.960 m。

4.2 基于輪廓點Z坐標傾斜監測法

流程如圖10所示。在點云拼接、去噪后，為了保證輪廓點提取精度，采用區域分割分層提取檐輪點云塔廓線，區域分割時，只需保證塔檐整體部分在分割區域內即可，如圖11所示。

圖10 輪廓點Z坐標提取法實驗流程

圖11 區域分割

因為在數據融合時為增加后期數據處理速度，故設置采樣間隔進行點云抽稀，為防止頂點和塔檐四角點云誤刪而導致Z坐標值最大最小點提取不準確，故在Cloudcompare軟件中選取各層塔檐四角和頂點坐標加入輪廓線提取坐標群，如圖12所示。

圖12 選取的塔角四點及塔頂

根據該方法進行輪廓點提取，最終各層塔檐共提取點云數量見表2。

表2 各層塔檐輪廓提取點云數量

提取各層輪廓點點云坐標后，導入表格從大到小排序，將Z坐標值最大最小點坐標單獨導出，根據公式(1)、(2)，求得各層傾斜角度以及基于該層傾斜度乘以塔的通高所求的塔整體傾斜量xi，詳細數據見表3。

表3 各層傾斜度及基于該層傾斜度的整體傾斜量

在求塔體整體傾斜量時，對表3中基于各層傾斜度整體傾斜量直接取算術平均值不能滿足測量中的精度要求，為了提高塔體整體傾斜計算精度，引入權要素，對表3中整體傾斜量做加權處理。對于塔形建筑物而言，底層較頂層相對傾斜量較小，且對整體傾斜影響最大，由于激光掃描點精度相同，且輪廓點提取點云密度直接影響能否提取出真正Z坐標值最大最小點，每一個輪廓點都相當于同精度觀測一次，因此用各層塔檐輪廓提取點數作為Ni次同精度觀測值參與定權計算。7次相同精度觀測塔體求解塔體傾斜量的過程，故可求同精度傾斜值算術平均值的權，具體公式見下式

(3)

(4)

(5)

式中Pi為第i層傾斜量的權；Ni為第i層輪廓點提取點云數，數據見表2；σ點為每一個點云的中誤差，每層測量精度相同，其中誤差為σ0，σi為第i層傾斜量中誤差。

將公式(4)(5)導入公式(3)中得

(6)

式中C為單位權輪廓點的點數，個。

取第一層輪廓點提取數做單位權觀測值，即C=2054。則根據公式(6)，可計算出每層傾斜量的權Pi。

(7)

由公式(7)可得傾斜量加權平均值即塔體整體傾斜量，結果為0.962 m。

4.3 數據分析與對比

通過提取的輪廓點最大值最小值點坐標數據對比后發現，文中方法所提取的Z坐標值最大、最小值點均為塔檐對角點，可根據文中方法原理可大致判斷出塔體傾斜方向為東南向西北方向傾斜?；谖闹蟹椒?，在能判斷塔體整體大致傾斜走向情況下，可直接從塔檐各角坐標中提取Z坐標值最大、最小點計算，省去提取輪廓點坐標這一步驟，減少計算量，可大大提高數據處理效率。為了驗證該規律的可靠性，在明確塔體大致傾斜方向后，隨機提取某層塔柱沿該傾斜方向對角點來進行驗證。最終，對塔體第3層和第4層塔柱的兩對角點進行提取坐標計算塔體的傾斜量以此來驗證該規律的可靠性，具體數據見表4。同樣，和提取中心軸法相比文中方法原理易懂，數據計算處理更為簡便。

根據表5各層傾斜度數據可以發現，每層傾斜角度并不相同，塔體底層傾斜量較小，隨著塔體升高各層傾斜量不一，這也是提取中心軸法不便將塔體分層計算傾斜量原因，塔各層的變形并不是線性的，采用不同層數據進行計算整體傾斜量會有不同的結果。因此現有基于點云數據提取中心軸法進行塔形建筑物傾斜監測分析時，均是將塔體假設為一個剛體進行分析，利用塔體中心點數據擬合出塔體的中軸線數據再計算塔體傾斜量。而文中提出的方法是通過分層計算塔體傾斜角度，再基于該傾斜角度進行整體偏移量的加權計算，既可計算塔整體傾斜量，亦可計算各層傾斜量。通過表4可以看出，無論是從塔檐提取輪廓點還是從塔柱提取輪廓點，都可以滿足文中輪廓點Z坐標傾斜監測方法。結合表5看出，2種方法計算塔體傾斜量盡管存在一些差異，所得出的大雁塔整體傾斜量相差0.002 m，約0.26%，但是在測量誤差范圍之內。總體來說，通過與擬合中心軸法整體傾斜量計算對比表明，文中所提的基于輪廓點Z坐標的塔形建筑物傾斜監測方法切實可行。

表4 數據特征驗證

表5 兩方法傾斜量計算結果對比

5 結 論

1)對大雁塔三維激光掃描數據進行研究計算，驗證文中所提出的基于輪廓點Z坐標傾斜監測法的有效性，該方法可以對基于三維激光掃描塔形建筑物傾斜監測提供一種新的數據處理方案。

2)采用基于塔體各層傾斜角度加權計算塔形建筑物整體傾斜量的方法，因此可以明確獲得塔體各層傾斜度，為塔體某層具體修繕維護提供數據參考。

3)在明確塔體基本傾斜方向后，可以直接提取該傾斜方向各層塔檐對角點云坐標數據，作為Z坐標最大最小值點參與傾斜量計算，減少點云計算量，可大大提高數據處理效率。

4)該方法適用性廣泛，圓形、六邊形等同類項目，可基于文中方法進行傾斜監測，對無中心軸不規則建筑物傾斜監測項目也具有一定的參考價值。