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

陳瀟，王鑫森，俞雪薇

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

1 引 言

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

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

2 主要技術(shù)流程

受施工工藝影響，地鐵盾構(gòu)隧道的斷面形狀一般皆非標(biāo)準(zhǔn)圓，而在數(shù)學(xué)模型中，正圓可視為橢圓的一個(gè)特例，故本文以通用的橢圓柱體作為地鐵區(qū)間隧道的數(shù)學(xué)模型[8]。本文的主要技術(shù)流程為：

(1)對(duì)掃描點(diǎn)云數(shù)據(jù)進(jìn)行配準(zhǔn)，并手動(dòng)刪除明顯、易剔的雜點(diǎn)；

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

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

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

圖1 總體技術(shù)流程圖

3 數(shù)學(xué)模型及算法關(guān)鍵點(diǎn)

3.1 橢圓柱模型建立

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

圖2 測(cè)量坐標(biāo)系與斷面坐標(biāo)系的示意圖

以上兩種坐標(biāo)系均為左手系，二者間的坐標(biāo)轉(zhuǎn)換關(guān)系為：

(1)

式中，△x、△y、△z為坐標(biāo)系平移參數(shù)，R為旋轉(zhuǎn)矩陣，包含9個(gè)旋轉(zhuǎn)參數(shù)：

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

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

(2)

寫成一般式為：

λy′2+μz′2=1

(3)

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

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

以上各系數(shù)均以未知參數(shù)的近似值計(jì)算。

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

(5)

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

3.2 坐標(biāo)重心化

測(cè)量坐標(biāo)系一般采用國(guó)家坐標(biāo)系或地方坐標(biāo)系，其平面坐標(biāo)值遠(yuǎn)大于高程值，亦遠(yuǎn)大于斷面坐標(biāo)系下的坐標(biāo)值。為避免由此造成的數(shù)值運(yùn)算精度損失，同時(shí)為增強(qiáng)兩種坐標(biāo)間的可對(duì)比性，應(yīng)對(duì)點(diǎn)云的測(cè)量坐標(biāo)實(shí)施坐標(biāo)重心化處理，即用各斷面點(diǎn)的三維坐標(biāo)分別減去同一斷面內(nèi)所有點(diǎn)的對(duì)應(yīng)坐標(biāo)分量均值，然后以重心化坐標(biāo)代替測(cè)量坐標(biāo)進(jìn)行平差計(jì)算。

3.3 參數(shù)初始化

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

(1)橢圓模型參數(shù)

橢圓模型參數(shù)的初始化一般采用隧道的設(shè)計(jì)參數(shù)，即令橢圓的長(zhǎng)軸a與短軸b皆等于設(shè)計(jì)半徑，在初次平差時(shí)視其為標(biāo)準(zhǔn)圓進(jìn)行處理。

(2)旋轉(zhuǎn)參數(shù)

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

為實(shí)現(xiàn)旋轉(zhuǎn)矩陣的初始化，本文需基于下列兩點(diǎn)假設(shè)：

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

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

如此可利用式(1)中的第一個(gè)方程建立唯一的函數(shù)關(guān)系式：

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

(6)

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

sinθ/cosθ=-x/y

(7)

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

sin2θ+cos2θ=1

(8)

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

(3)平移參數(shù)

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

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

3.4 迭代濾波

掃描點(diǎn)云中存在大量的噪聲點(diǎn)必然會(huì)嚴(yán)重影響模型參數(shù)和中心坐標(biāo)的擬合精度，因此本文針對(duì)隧道噪聲點(diǎn)的分布特點(diǎn)，并結(jié)合隧道設(shè)計(jì)參數(shù)(半徑)，代入逐級(jí)縮小的閾值進(jìn)行迭代濾波，直至收斂穩(wěn)定為止。

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

點(diǎn)云濾波的基本方法步驟為：

(1)利用既有的平移參數(shù)和旋轉(zhuǎn)參數(shù)將點(diǎn)云的重心化測(cè)量坐標(biāo)轉(zhuǎn)換為斷面坐標(biāo)；

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

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

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

4 工程試驗(yàn)

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

4.1 濾波分析

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

圖3逐次濾波效果圖

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

4.2 擬合分析

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

圖4 擬合中心對(duì)比

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

各斷面的橢圓度 表1

表1為通過本方法平差計(jì)算出的上述15個(gè)斷面的長(zhǎng)、短半軸及橢圓度。各斷面的長(zhǎng)半軸與短半軸的較差約 1 cm，其橢圓度值基本處于3.0～5.0的區(qū)間浮動(dòng)，這充分說明較之標(biāo)準(zhǔn)圓柱體，本文選用橢圓柱體作為盾構(gòu)隧道的模型是更符合實(shí)際的。另外由于長(zhǎng)、短軸的差別較小，點(diǎn)云濾波環(huán)節(jié)若作為正圓模型進(jìn)行處理，會(huì)令運(yùn)算更為簡(jiǎn)潔高效。

5 結(jié) 語(yǔ)

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