基于橢圓柱模型的隧道斷面擬合算法研究

陳瀟，王鑫森，俞雪薇

(天津市勘察院，天津 300191)

1 引 言

隧道斷面測量是調線調坡設計和結構限界檢測中的一個關鍵環節，可用于檢驗施工過程中產生的施工誤差以及結構變形等不利因素，作為后續工程的參考依據，因此其測量的精度與數據的準確性對于整條地鐵線路的施工建設至關重要。傳統的全站儀測量方法采集數據速度慢、效率低，且僅能反映不連續的結構特征，難以作為長距離、大范圍的數據采樣手段。而地面三維激光掃描技術的出現為解決該問題帶來了革命性突破，它通過高速激光測距迅速獲取大量的目標對象表面數據點，可提供全面的點云信息，具有非接觸性、高精度、高密度等優點，克服了傳統測量技術的局限性和片面性[1，2，3，7]。

利用三維激光掃描儀在盾構隧道內部施測往往需面對較為復雜的現場環境，施工器具、金屬支架、臨時軌道均會對掃描激光準確捕捉隧道內壁點產生不利影響[4，5]。因此，如何通過濾波有效剔除多余的噪聲點，同時快速、準確地求解隧道斷面參數、提取隧道中心線都是工程中亟須解決的技術難題。雖然國外的一些廠商(如瑞士Amberg)針對此類問題都提供過專業的解決方案，但其硬件設備和軟件系統的價格十分昂貴，難以在隧道工程中普及應用。為此本文提出了一套可行的濾波及擬合算法，并已在工程實踐中進行驗證。

2 主要技術流程

受施工工藝影響，地鐵盾構隧道的斷面形狀一般皆非標準圓，而在數學模型中，正圓可視為橢圓的一個特例，故本文以通用的橢圓柱體作為地鐵區間隧道的數學模型[8]。本文的主要技術流程為：

(1)對掃描點云數據進行配準，并手動刪除明顯、易剔的雜點；

(2)利用天寶公司Realworks軟件對隧道點云進行橫斷面切片處理[6]，該步驟需提供隧道中心線的初始值，一般經兩至三次迭代切片便可批量獲取高質量的隧道橫斷面切片；

(3)建立實地測量坐標系與斷面坐標系之間的轉換關系，利用最小二乘原理擬合模型參數和隧道中心坐標，其計算結果作為初始值進入下一輪迭代，同時每次迭代之前需先完成一次點云濾波(濾波閾值逐次減小)，直至迭代收斂；

(4)記錄、保存各斷面的擬合參數，包括橢圓模型參數、旋轉參數、平移參數及擬合中心點三維坐標，以固定格式輸出成果。其流程圖如圖1所示：

圖1 總體技術流程圖

3 數學模型及算法關鍵點

3.1 橢圓柱模型建立

建立隧道模型的關鍵在于求解模型的幾何參數以及確定實地測量坐標系與模型坐標系(即斷面坐標系)之間的轉換關系。如圖2所示，其中O-XYZ為測量坐標系，O′-X′Y′Z′為斷面坐標系，O′點為斷面(即橢圓)中心，X′軸與橢圓柱中軸線重合，Y′、Z′軸在隧道斷面內。

圖2 測量坐標系與斷面坐標系的示意圖

以上兩種坐標系均為左手系，二者間的坐標轉換關系為：

(1)

式中，△x、△y、△z為坐標系平移參數，R為旋轉矩陣，包含9個旋轉參數：

隧道斷面在O′-X′Y′Z′系下的方程為：

y′2/a2+z′2/b2=1

(2)

寫成一般式為：

λy′2+μz′2=1

(3)

將式(1)代入式(3)，按泰勒公式展開后，可構建下列誤差方程：

v=a1dλ+a2dμ+a3dr21+a4dr22+a5dr23+a6dr31+a7dr32+a8dr33+a9△y+a10△z+l0

(4)

式中：

l0=λ(△y+r21x+r22y+r23z)2

+μ(△z+r31x+r32y+r33z)2-1

a1=r21x2+r22y2+r23z2+2r21r22xy+2r21r23xz+2r22r23yz

+2r21x△y+2r22y△y+2r23z△y+△y2

a2=r31x2+r32y2+r33z2+2r31r32xy+2r31r32xz+2r32r33yz

+2r31x△z+2r32y△z+2r33z△z+△z2

a3=2λr21x2+2λr22xy+2λr23xz+2λx△y

a4=2λr22y2+2λr21xy+2λr23yz+2λy△y

a5=2λr23z2+2λr21xy+2λr22yz+2λz△y

a6=2μr31x2+2μr32xy+2μr33xz+2μx△z

a7=2μr32y2+2μr31xy+2μr33yz+2μy△z

a8=2μr33z2+2μr31xz+2μr32xz+2μz△z

a9=2λ△y+2λr21x+2λr22y+2λr23z

a10=2μ△z+2μr31x+2μr32y+2μr33z

以上各系數均以未知參數的近似值計算。

由于橢圓方程(2)不含x項，式(4)中不含△x和r11、r12、r13。旋轉矩陣為正交矩陣，故式中的6個元素必須滿足如下限制條件：

(5)

將式(5)線性化展開后，與式(4)共同構建附有限制條件的間接平差[9]。解算該方程組，應至少選取10個隧道斷面點的數據。

3.2 坐標重心化

測量坐標系一般采用國家坐標系或地方坐標系，其平面坐標值遠大于高程值，亦遠大于斷面坐標系下的坐標值。為避免由此造成的數值運算精度損失，同時為增強兩種坐標間的可對比性，應對點云的測量坐標實施坐標重心化處理，即用各斷面點的三維坐標分別減去同一斷面內所有點的對應坐標分量均值，然后以重心化坐標代替測量坐標進行平差計算。

3.3 參數初始化

上述附有限制條件的間接平差方程包含10個未知參數，分別為2個橢圓模型參數λ、μ， 2個平移參數△y、△z， 6個旋轉參數r21～r33。另外為求取隧道中心線，X方向上的平移參數△x應一并求解，所以待求的參數初始值共計11個。

(1)橢圓模型參數

橢圓模型參數的初始化一般采用隧道的設計參數，即令橢圓的長軸a與短軸b皆等于設計半徑，在初次平差時視其為標準圓進行處理。

(2)旋轉參數

旋轉矩陣的9個參數實際是由3個獨立的旋轉角所確定，但因缺少公共點，無法利用式(1)反算旋轉參數(或旋轉角)的初始值。若以單位矩陣作為其初始值，即假定旋轉角的初始值皆為零，則試驗證明：當測量坐標系與斷面坐標系的x軸(或y軸)夾角小于30°時，迭代平差最終趨于收斂；否則，結果將無法收斂[8，10]。因此該方法不具備普適性。

為實現旋轉矩陣的初始化，本文需基于下列兩點假設：

①假設兩坐標系的Z軸平行。考慮在豎直方向上，Z軸與Z′軸的夾角往往較小，故暫且視其為平行軸，如此可將三軸旋轉問題簡化為單軸旋轉問題(繞Z軸旋轉)，僅剩求解1個旋轉角θ。

②假設斷面切片上所有點的x′坐標值均為0。截取隧道斷面時嚴格控制切片厚度(約 2 cm)，計算時僅考慮y′、z′分量變化，而令x′分量始終為0。

如此可利用式(1)中的第一個方程建立唯一的函數關系式：

x′=△x+cosθ·x+sinθ·y=0

(6)

值得注意的是，式中x、y均為重心化測量坐標。因重心化測量坐標系與斷面坐標系的原點在同一平面內，故△x=0，式(6)可變換為：

sinθ/cosθ=-x/y

(7)

正弦與余弦的比值易受誤差影響而產生較大偏差，對此應通過多點求解再取均值，然后代入下列限制條件方程求解：

sin2θ+cos2θ=1

(8)

該二次方程具有正、負兩根，取正或取負的差別僅取決于斷面坐標系X軸的不同指向(大里程或小里程)，不會影響擬合結果，因此計算時可任選其一，進而求出全部旋轉參數的初始值。

(3)平移參數

本文采用的濾波方法是基于斷面幾何中心的迭代篩選，平移參數初始值的偏差過大很可能導致大量有效點云被自動剔除，選擇合適的初始值對于后續的點云濾波具有十分重要的意義。

由上文可知平移參數△x始終為零，其余的△y、△z可通過式(4)和式(5)平差求解。考慮到首次平差時點云數據未經濾波，若噪聲點在局部較為集中(如鋪軌后的軌道及道床掃描點均視為噪聲點)，易引起隧道模型的長軸、短軸發生較大的形變，極大地影響擬合精度。因此，本文在初始化平移參數時，省略式(4)中dλ、dμ兩項改正，即使用固定的橢圓模型參數計算平移參數△y、△z，得到更準確的初始值。

3.4 迭代濾波

掃描點云中存在大量的噪聲點必然會嚴重影響模型參數和中心坐標的擬合精度，因此本文針對隧道噪聲點的分布特點，并結合隧道設計參數(半徑)，代入逐級縮小的閾值進行迭代濾波，直至收斂穩定為止。

根據平面橢圓的數學特性，可定義如下濾波函數：f(y′，z′)=(λy′2+μz′2)1/2。對于橢圓平面內的任意一點，若該點恰恰處在橢圓曲線上，則f(y′，z′)=1；當該點落在橢圓內部時，有0≤f(y′，z′)<1；而當該點位于橢圓外部時，有f(y′，z′)≥1。濾波閾值的設置應按從大到小的次序逐級遞減(如 30 cm、 20 cm、10 cm、5 cm)，具體數值視實際需求而定。由于濾波函數為無量綱表達式，其閾值也應改用無量綱的數值，即：(1±p/r)，其中p代表某一閾值，r為設計半徑值。

點云濾波的基本方法步驟為：

(1)利用既有的平移參數和旋轉參數將點云的重心化測量坐標轉換為斷面坐標；

(2)將各點一一代入濾波函數進行檢驗，剔除超限點，再以過濾后的點云數據重新平差，計算全部未知參數；

(3)重復上述步驟(1)和(2)，直到相鄰兩次迭代的改正數小于0.001，停止迭代；

(4)采用更小的閾值，重復上述步驟(1)～(3)，直到最小閾值條件下的收斂完畢。

4 工程試驗

本文選擇天津地鐵3號線某段隧道作為研究對象，該區間位于市內繁華地帶，采用盾構法施工，隧道斷面的設計內徑為 5.5 m。本次試驗使用瑞士Amberg的隧道狀態激光全息成像掃描系統進行點云數據采集，該系統集成最新的Profiler 6012型360°激光掃描儀與GRP 5000隧道檢測小車，其采樣頻率可達100萬點/秒，測距精度優于 ±1.0 mm(相對儀器中心 5 m半徑)。為驗證本文算法的擬合效果，掃描得到的原始點云采用兩種方式處理：一是用系統配套的專業限界及結構變形分析軟件，二是用本文闡述的方法及編制的C#軟件。

4.1 濾波分析

本次濾波采用的閾值依次為R±20 cm、R±10 cm、R±5 cm，其中R為設計半徑 2.75 m。以上閾值分別除以設計半徑可得到對應的比例閾值，即：(1±0.07)、(1±0.04)、(1±0.02)。圖3所示為某斷面作逐次濾波時先后得到效果示意圖：

圖3逐次濾波效果圖

從圖3中可以看出，隧道的原始點云中混雜著鋼支架、電纜、道床等大量噪聲點，利用本文的方法經三次濾波，已將主要的噪聲點濾除，剩余的點云數據(除道床區域外)皆均勻分布于隧道內壁，這說明本方法具有較好的去噪效果。

4.2 擬合分析

本次試驗掃描的隧道長度約400 m，分別采用TMS軟件與Realworks軟件基于相同的設計中心線、在相同里程處進行斷面切片，相鄰切片的間距約為 8 m(彎道及特征點處需要加密)。限于篇幅，僅選取部分數據作對比分析。

圖4 擬合中心對比

如圖4所示，圖中橫軸表示隧道斷面中心點號C1～C15，豎軸表示由TMS軟件或自編程序擬合的實際中心點與設計中心點的點位偏差△P。不難發現二者的擬合結果存在一定的差別：TMS軟件解算結果在數值上與采用的設計中心線更為吻合，其△P均值為 1.9 mm，最大值也僅為 5.8 mm；而自編軟件的對比結果多介于4 mm～8 mm，僅一點的偏差超過 8 mm(10.8 mm)。該差異產生的原因主要在于：兩種軟件在指定位置截取的切片會存在細小的差別，進而影響其中心點的擬合坐標。對比而言，作為專業的隧道點云處理軟件，TMS明顯更具有精度優勢及自動化優勢，但以本文的方法作為一種更經濟的替代方案，也基本滿足工程需求。

各斷面的橢圓度 表1

表1為通過本方法平差計算出的上述15個斷面的長、短半軸及橢圓度。各斷面的長半軸與短半軸的較差約 1 cm，其橢圓度值基本處于3.0～5.0的區間浮動，這充分說明較之標準圓柱體，本文選用橢圓柱體作為盾構隧道的模型是更符合實際的。另外由于長、短軸的差別較小，點云濾波環節若作為正圓模型進行處理，會令運算更為簡潔高效。

5 結 語

針對地鐵隧道的形態及其掃描點云的分布特點，本文提出了采用橢圓柱體進行建模的方法，并推導參數擬合及點云濾波的計算公式，編制相關的計算程序。經對比試驗證明：本方法的擬合結果與預期效果相符。在顧及成本的前提下，以本方法替代專業廠商系統軟件的部分關鍵功能具有可行性，可滿足工程需求。此外，本方法中的隧道斷面切片仍然依靠Realworks等處理軟件以半自動化的方式提取，因此探索基于點云的隧道斷面自動化截取將是下一步研究工作的重點。