Kriging點云濾波改進算法及監測試驗研究

胡達,黎永索,張可能,梁曉東,梁小強,吳有平

（1. 湖南城市學院 城市地下基礎設施結構安全與防災湖南省工程研究中心；土木工程學院，湖南 益陽 413000；2. 中南大學 a.有色金屬成礦預測與地質環境監測教育部重點實驗室；b. 地球科學與信息物理學院，長沙 410083；3. 湖南聯智科技股份有限公司，長沙 410200；4. 中國電建集團中南勘測設計研究院有限公司，長沙 410014）

對隧道圍巖的穩定性進行監控及預報是確保隧道現場施工安全、提高施工信息化水平、優化圍巖支護設計參數必不可少的一項重要工作。目前，隧道監控量測的主要方法仍然是采用精密水準儀、全站儀、智能收斂計等儀器設備進行測量，這些測量工作首先需要在隧道洞壁上及時布設監控測點（三角掛鉤、反光片等）。這種傳統監控量測方法實際操作起來非常麻煩，效率極低，且誤差不容易控制，監測過程中極易造成對圍巖穩定性評價的錯誤判斷，使隧道施工安全性降低。

近年來，激光掃描技術逐漸興起，該技術具有高精度、高效率、高分辨率、全自動數字化采集、數據量豐富等特點[1]。將三維激光掃描技術應用于隧道監控量測是未來發展的一個重要研究方向，學者們在該技術的探索和研究上取得了大量的成果。王濤等[2]基于點云數據預處理問題，采用小波分析格網化濾波方法，并結合南京地鐵隧道進行了應用，取得了較好的濾波效果。許度等[3]提出了一種基于三維激光掃描技術的隧道變形非接觸原位測試技術，建立了基于掃描信息的圍巖結構面表征方法，并在中國錦屏地下深埋硐室群工程進行了應用。吳勇等[4]采用三維激光掃描技術進行運營期隧道檢測，為隧道的病害檢測及修復提供了一定參考。杜黎明等[5]提出了一種迭代橢圓擬合的方法，對隧道變形點云數據的提取進行了研究。張帆[6]提出了一種“聚合”算法，并應用于隧道管片結構變形可視化分析。徐衛東等[7]以北京地鐵7號線為依托，運用三維激光掃描技術對隧道開挖的現狀進行了測量。虞偉家[8]采用新型移動式三維激光測量技術，對盾構隧道進行了點云數據采集，并在蘇州地鐵進行了隧道監測試驗。黃帆等[9]對盾構隧道管片變形進行了隧道監測試驗，點云數據分析結果表明，三維激光掃描技術優勢明顯，對地鐵施工具有一定的指導意義。上述研究成果表明，目前精度最高的三維激光掃描儀在很多領域的應用上已經可以替代全站儀，能夠基本滿足隧道監控量測的精度要求。

大量研究表明[10-11]，點云濾波算法成為三維激光掃描精度是否滿足應用要求的關鍵。因此，對點云濾波進行針對性的研究十分必要。結合近年來各種濾波算法不同的理論背景，可以分成基于坡度、曲面擬合、分割、不規則三角網、形態學及機器學習的濾波算法6大類。基于坡度的濾波算法原理簡單、易實現，但過度依賴閾值的設定，在地形起伏較大區域不適合海量數據的實時處理，濾波效果不佳[12-13]。基于曲面擬合的濾波算法有賴于插值方法的選擇，所采用的多層級迭代方式會受到每一層級濾波結果的影響，容易出現誤差傳遞與累積[14-16]。基于分割的濾波算法其濾波效果過分依賴聚類分割的結果，而點云分割方法的選擇也對濾波結果有很大影響[17]。基于不規則三角網的濾波算法需要占用大量內存，并且對低位噪聲敏感，極易誤判低地勢物點[18]。基于形態學的濾波算法原理簡單、實現效率高，但此類算法在地形起伏較大區域的穩健性有待提高，如何提高其整體精度將是后續研究的重點[19]。基于機器學習的濾波算法需要大量的訓練樣本，其樣本必須覆蓋所有可能的地形特征，對計算機資源要求非常高，取得較好的濾波效果難度較大[20]。

綜上所述，以上各種點云濾波算法均各有局限性，如何改進和完善濾波算法，提高算法精度成為監測試驗能否成功的關鍵。由于Kriging插值算法在快速自動生成DEM方面精度高，整體精度評價方法具有可操作性[21]。筆者結合Kriging插值算法，通過現場試驗段監測，對點云數據進行提取分析，將試驗得出的變形數據與傳統方法的量測數據進行對比，以驗證其算法精度，為推進三維激光掃描技術在隧道監控量測中的應用提供理論依據。

1 Kriging濾波改進算法

Kriging插值法是一種運用結構分析與變異函數相關理論，在有限空間內針對區域化變量進行最優和無偏估值計算的預測方法，是研究空間變異和進行空間插值的一種線性無偏最優估計算法。Kriging插值法不僅考慮待估點位置與已知數據位置的相互關系，還考慮變量的空間相關性。假設區域化變量不是相互獨立的，具有一定的隨機性和結構特性，而且滿足二階平穩和內蘊平穩。同時，Kriging方法在數據網格化的過程中考慮了描述對象的空間相關性質，使插值結果更科學、更接近于實際情況，能給出插值的誤差（Kriging方差），使插值的可靠程度一目了然[22]。考慮時間因素，將經典的Kriging空間插值進行拓展[23-24]，建立時空Kriging插值模型，并將其應用于點云格網化濾波分析，提高三維激光掃描數據獲取的精度。

1.1 空間插值

Kriging插值法的基本數學模型[25]為

式中：Z′(X0)為預測點的估計值；Z(Xi)為預測點鄰域內參與預測的參考點的值；λi為Kriging權系數，在無偏性和最小方差性條件下，其值依賴變異函數的計算結果來確定。

式（1）是n個數值的線性組合，Kriging方法的原則就是保證估計量無偏，且在估計方差σ2E最小的前提下，求出n個權值系數λi。

在無偏性條件下，為了使估計方差最小，用拉格朗日乘子法求條件極值，見式（3）。

F是n個權 系 數 和μ的(n+1)元函 數，求出F對λi和μ的偏導數，最后得到用半變異函數值表示的Kriging方程組。

式 中：yi,j=y(xi,xj)=y(xi-xj)，為 半 變 異 函數值。

半變異函數或者半變差函數是從空間統計學中的方差概念演化而來的，將區域變量Z(x)在點x和（x+h）的值Z(x)和Z(x+h)之差的方差的1/2定義為區域變量Z(x)在x軸上的半變差函數，記為γ(h)，即

試驗半變異函數計算公式為

式中：N(h)為被增量h分隔的試驗數據對的數目。

為了估計區域化變量的未知值，需要將半變量試驗函數擬合到相應的理論半變異函數模型中。典型Kriging算法提供的半變異函數模型包括高斯模型、線性模型、球面模型、阻尼正弦模型、指數模型等，其中，原點附近的球形變異函數模型是線性的，指數變異函數模型是直線，高斯變異函數模型是拋物線。考慮到本研究是將其應用到隧道變形監測的濾波算法中，選用球形模型，這樣既考慮了儲層參數的隨機性，又考慮了儲層參數的相關性，在滿足插值方差最小的條件下，給出最佳線性無偏插值和方差。具體模型表達式為

式中：c0為塊金值；c0+c為基臺值；c為偏基臺值；a為變程；h為滯后距。

濾波算法流程如圖1所示。

圖1 Kriging算法流程圖Fig. 1 Kriging algorithm flow chart

1.2 時空插值

空間Kriging只能估計某一時間上未知區域的變形量，如果要估計任意時刻任意位置的變形量，需要進行時空Kriging插值。在隧道掘進過程中，監測點的布置、數據的采集以及位移值的變化不是空間和時間上的單一變化，而是時間與空間的統一。選用普通Kriging方法實現數據的時空插值[26-27]，計

算公式為

式中：Z*(s0,t0)為時空點(s0,t0)處的估計值；λi為臨近觀測值Z(si,ti)的加權系數，引入拉格朗日系數μ進行推導可得

式（9）中的加權系數λ的值可通過式（1）得到，繼而代入式（8），可得研究區域內任意點的插值估計值。

2 現場監測試驗方案

采用德國Z+F IMAGER 5010X三維激光掃描儀對長沙市城市固體廢棄物處理場管理生活區隧道工程進行監控量測，并將所量測的數據與傳統監控量測方法的數據進行對比研究。德國Z+F IMAGER 5010X三維激光掃描儀數據獲取速率超過100萬點每秒，50 m處精度最高可達0.8 mm，測程為0.3～187.3 m，是一臺高性能、高精度的旗艦掃描儀，非常適合應用于高精度的工業測量和隧道工程測量。該儀器精度信息如表1所示。

表1 Z+F IMAGER 5010X三維激光掃描儀精度信息Table 1 Precision information of Z+F IMAGER 5010X 3D laser scanner

2.1 工程概況

所采用的試驗數據來源于長沙市城市固體廢棄物處理場管理生活區隧道工程監控量測項目現場。隧道場地位于長沙市望城區橋驛鎮黑麋峰固體廢棄物處理場西側一帶的山間空地，為望湘巖基的南部地區，以侵蝕構造地貌為主，屬低山丘陵地貌，主要表現為低山—丘陵與沖溝形態，山脈脈絡清晰，山頂一般較圓滑，個別較尖，山脊大多數開闊舒緩，局部地段狹窄，形成窄陡山脊。隧道起始樁號為K0+195～K0+285，長為90 m。根據勘察結果，該場地內地形起伏大，山坡較為陡峭，場地內分布的地層主要有第四系殘坡積層，下伏基巖為燕山晚期花崗巖。

2.2 試驗方案

地面三維激光掃描外業實施計劃的制定需要根據測量任務、要求以及現場條件決定，主要包括坐標系、掃描儀及其配準靶標和掃描站的選擇等。根據實施計劃，提前布設好靶標，然后連接相關設備，設置掃描參數（如掃描范圍、掃描距離、掃描間隔等），在不同的掃描站上架設掃描儀進行掃描。

2.2.1 控制點的布設與測量 在隧道工程監測試驗中，由于現場條件限制，控制點只能布設在一側，如圖2所示。因此，在隧道洞口布設了3個自制的平面靶標，如圖3所示。作為控制坐標轉換點，要求這些點的位置相對固定，并且不易被施工破壞。布設好控制轉換點后，用全站儀測得靶標中心在施工控制坐標系下的坐標，用于后續轉換。對于站間轉換公共點，選擇不受入射角影響、反射率高、方便布設的標準球形靶標，如圖4所示。將3～4個球形靶標布設在兩相鄰測站之間，并在兩測站上均對其進行精掃。

圖2 平面控制點布設圖Fig. 2 Layout of plane control points

圖3 自制平面靶標Fig. 3 Homemade flat target

2.2.2 數據采集 外業測量是實際獲取數據的過程，掃描距離、物體表面材質、控制網、靶標測量精度、光斑大小、掃描點間距、點云拼接精度、全反射物質和外界環境等因素均對點云數據的質量產生直接影響。與其他光學電子儀器一樣，外界環境溫度、氣壓、空氣質量等因素對激光回波信號的影響顯著。因此，鑒于隧道內施工環境異常復雜，為避免數據采集質量受到影響，選取隧道內空氣可見度較高的時段進行試驗測試。

現場掃描試驗前，先在隧道內確定好儀器和靶標的位置，確保在進行掃描時不會被打斷或遮擋。布設好靶標球，如圖4所示。架好掃描儀，設置好掃描參數，然后進行掃描作業。掃描作業時靶標球需要保持固定狀態，避免工人觸摸或誤動，從而致使球的位置發生變化，影響后續點云配準精度。試驗掃描目標為隧道襯砌外表面的目標點，所以不必進行全景粗掃，而是先直接進行全景目標掃描，再在此基礎上對靶標球進行精掃（高精度高質量掃描）。試驗時可按圖5進行設站和靶標布設（站1～站3為測量站點），具體要求如下：兩站之間間距約25～30 m，保證兩站之間點云重合數不少于30%；靶標球或靶標板擺設與兩站中間至兩邊站的距離比較均勻；靶標盡量上下左右起伏較大，保證空間位置形狀良好；兩站之間至少有4個靶標，便于后期拼站的精度。圖6為三維激光掃描儀在隧道現場的掃描示意圖。

圖4 標準球形靶標Fig. 4 Standard spherical target

圖5 掃描測站及球形靶標位置分布圖Fig. 5 Location distribution of scanning station and spherical target

圖6 隧道掃描示意圖Fig. 6 Schematic diagram of tunnel scanning

2.3 數據處理

2.3.1 數據處理基本流程 對于三維激光掃描儀獲得的點云數據，并非每一個都能使用，過多的數據點會降低計算機運行的效率、增加存儲空間。隧道點云整體圖如圖7所示。為避免上述問題，需要刪除部分數據點，即對點云數據進行精簡處理，點云數據處理的主要流程為：

圖7 隧道點云整體圖Fig. 7 Overall drawing of tunnel point cloud

1）點云編輯。裁減掉點云粗差以及與目標無關的點云。

2）點云配準與拼接。利用公共點坐標將不同測站測量的點云數據轉換到同一坐標系中，實現點云的拼接，形成一個整體。目前，點云配準的方法主要有3種：人工靶標配準、點云自身配準和控制坐標系配準。配準后的點云圖如圖8所示。

圖8 配準后的兩個測站點云圖Fig. 8 Cloud map of two stations after registration

3）建立拓撲關系。點云通常是孤立的點，每個點只與其一定范圍內的周圍點相關。空間點云拓撲關系的建立主要有八叉樹法、網格法和k-d樹法等。

4）點云數據的精簡。點云數據精簡算法大致可以分為4類：包圍盒法、隨機采樣法、曲率采樣法和均勻網格法。

5）點云數據的去噪與光順。由于被測對象表面的粗糙度、波紋等缺陷和測量系統本身的影響，真實數據中混有噪聲點。可根據點云質量和后續建模要求靈活選擇合適的濾波算法。

6）孔洞的修補。激光掃描的過程中會因為各種原因（如局部遮擋）造成漏測，從而形成點云孔洞。修補算法可以分為2種：一種是在三角網格表面重建時進行三角網格面修補；另一種是先對散射點云數據孔洞進行修補，然后對三角網格面進行重建。

7）點云數據的分割。數據分割形成的不同曲面類型的子區域具有特征單一、凹凸一致的特點。對每一子區域進行單獨重構有利于曲面擬合時減小誤差和保持點云性質。目前，點云分割的算法主要有基于邊的算法、基于面的算法和基于聚類的算法。

8）點云數據的三維建模。在三維模型的建構過程中，曲面重構是最關鍵也最復雜的一步。目前，主要有兩類曲面重構方案：以三角網格面為基礎的自由曲面重構方案和以樣條曲面為基礎的自由曲面重構方案。

2.3.2 點云數據處理 點云濾波是點云數據處理中一個非常重要的環節，準確的濾波結果將有利于提高點云后處理結果的精度。根據普通Kriging插值的原理，采用Matlab程序實現Kriging濾波算法，并利用程序對點云數據進行處理。為減少濾波算法的工作量，選取最為接近垂直于靶標球表面的3個測點進行計算，取c0=0.01、c=0.2、a=20。取n=3，根據式（9），有μ=0.015 9、σ2E=0.054 7。

1）將外業掃描獲取的點云數據導入到數據處理程序中，進行點云數據預處理操作（濾波、去噪），為后續操作做準備。

2）數據預處理完成后，在配準流程（Registering）中，分別用平面靶標擬合工具和球形靶標擬合工具對靶標進行擬合，提取出靶標中心在掃描儀坐標系中的三維坐標。

3）靶標中心坐標提取完成后，再進行點云配準，將全站儀測得的控制靶標的坐標保存在TXT文件中，并將點云數據轉換到控制坐標系中。如果布設的公共點較多，配準時應該選擇質量較好的點進行點云配準。最后將配準后的點云數據坐標信息導出到TXT文件中，完成數據處理流程。

3 數據提取及分析

3.1 監測信息提取

根據隧道表面的點云數據坐標，可以提取出任意樁號處的隧道斷面數據，從而進行斷面擬合，提取出中心坐標，然后擬合出中線。在這些數據上可以對超欠挖情況、拱頂沉降情況、收斂狀態和軸線偏差等進行分析，對施工進行指導，及時發現不穩定及危險情況并上報，減少災害發生及其不良后果的產生。

3.1.1 斷面提取 在隧道監控量測中，所有的測量工作都是在斷面的基礎上完成的，所以，三維激光掃描數據應用在隧道中的第一步就是斷面提取，否則其他的監控量測項目將無法完成。可以根據隧道的設計資料來提取斷面。由斷面的定義可知，隧道中軸線上某點處的斷面與該點處的切線垂直，即該點處的切線向量就是該斷面的法向量。設該點在中軸線上的坐標為(x0,y0,z0)，該點處的切線向量為n=(A,B,C)，則該點處的斷面方程為

中軸線上某點的三維坐標可以根據隧道段的線路設計參數確定，該點處的切線向量n=(1,k,i/cos(atank))，k為xoy平面中的切線斜率，i為坡度。

不管掃描儀的分辨率有多高，測量的點與點之間總會有一定的間隔，所以，在實際提取斷面時，該斷面并不是嚴格的平面，而是有一定厚度的。在實踐應用中，一般提取1～2 cm厚度平面的數據作為斷面數據。圖9～圖11為利用該方法提取的同一樁號處的隧道斷面。圖9為在施工坐標系中提取的斷面；圖10為施工坐標系中旋轉到正面后的斷面，圖11為在獨立坐標系中的斷面，原點在軸線中心。

圖9 三維空間坐標系中提取的斷面Fig. 9 Extracted section in 3D coordinate system

圖10 三維空間坐標系中提取的斷面(橫截面)Fig. 10 Extracted section in 3D coordinate system(cross section)

圖11 平面坐標系中提取的斷面Fig. 11 Extracted section in plane coordinate system

3.1.2 變形信息提取 依據行業規范《公路隧道施工技術規范》（JTG/T 3660—2020），在隧道監測中，周邊位移和拱頂沉降是必測項目，并且每5～50 m需要測量一個斷面。三維激光掃描儀掃描的數據是連續、全面的，所以用掃描數據進行拱頂沉降和收斂變化的監測可以在任意斷面處進行。而傳統監控量測方法對隧道周邊的收斂量測一般采用隧道凈空變化測定計（簡稱收斂計）進行。拱頂沉降和收斂變化都是在斷面上選擇固定的點，對其進行連續的監測，收斂量測測點與拱頂下沉測點布置在同一斷面。傳統方法需要使用預埋件，而三維激光掃描方法不需要任何預埋件。根據隧道情況，拱頂沉降監測一般選取1～3個點，收斂變化一般選取2～3對點。

在斷面中提取監測點數據時，一般以監測點為中心的面狀數據為基礎，然后對區域內的點進行距離加權平均求取監測數據。設選取的監測點設計坐標為(x0,y0,z0)，選取區域內的點坐標為(xi,yi,zi),i=1,2,3…n，則監測點的測量坐標為

獲取了各個測點的監測數據后，就可以以第一期數據為基礎，作出時間變化曲線，以反映各斷面的收斂與沉降情況，還可以作出空間變化曲線，以反映隧道的整體變化情況。

3.2 數據分析

在隧道坐標文件的基礎上進行數據挖掘與分析。根據隧道表面的點云數據坐標，可以提取出任意樁號處的隧道斷面數據，從而進行斷面擬合，提取出中心坐標，然后擬合出中線。選擇該隧道工程具有代表性的3個斷面進行對比分析，分別采用反距離倒權法（IDW）、不規則三角網法（TIN）與克里金法（Kriging）進行對比研究。3個斷面在監測周期內的最終累計變形對比見表2、表3，圖12～圖17為所提取的K0+202、K0+207、K0+212斷面的位移收斂示意圖。

圖12 K0+202斷面(2016-09-30)Fig. 12 Section of K0+202 (2016-09-30)

圖17 K0+212斷面(2016-11-12)Fig. 17 Section of K0+212(2016-11-12)

表2 3種濾波方法對比Table 2 Comparison of three filtering methods

表3 試驗結果對比Table 3 Comparison of test results

圖13 K0+202斷面(2016-10-29)Fig. 13 Section of K0+202(2016-10-29)

圖14 K0+207斷面(2016-10-07)Fig. 14 Section of K0+207(2016-10-07)

圖15 K0+207斷面(2016-11-06)Fig. 15 Section of K0+207(2016-11-06)

圖16 K0+212斷面(2016-10-13)Fig. 16 Section of K0+212(2016-10-13)

從表2中的對比結果來看，由于三維激光掃描獲取了海量且較為均勻的監測數據，從而對IDW和Kriging方法十分有利，所得結果較相近；而TIN法的濾波結果誤差相對較大。總體情況表明，相對于IDW法，Kriging改進算法的濾波更接近真實值，算法精度更高，且更適合數據較多的情況，同時也說明該方法是行之有效的。

由表3可知，三維激光掃描量測的累計變形總體相對較小，其中K0+207、K0+212拱頂沉降數據較為接近傳統量測結果，周邊收斂數據差異較大。結合現場情況并對設備精度進行分析可知，隧道內環境因素對量測精度影響極大。對于拱頂沉降，掃描儀只需提取Z軸方向的坐標點數據就可以得出沉降結果，因而比水平坐標的換算更為穩定；對于周邊收斂，掃描儀必需提取X、Y、Z軸3個方向的坐標點數據才能有效得出收斂變形結果，由此造成的誤差較大。除數據提取外，隧道內環境等因素也對監測結果造成較大的不確定性，進而造成誤差累積。試驗數據同時也說明，在隧道變形監測中，三維激光掃描技術的抗干擾能力和精度還有待進一步優化和提高。在特定的監測環境下，三維激光掃描技術所量測的隧道累計變形值在一定程度上能較好地反映隧道變形的真實情況。

4 結論

采用基于Kriging的濾波算法，對三維激光掃描點云數據進行處理，通過現場試驗與傳統點云濾波數據處理方法的對比，說明該方法能夠高效識別并提取隧道輪廓斷面的可視化數據，為三維激光掃描技術在隧道監測領域的廣泛應用提供解決方案。試驗結果表明：

1）將經典的Kriging空間插值進行擴展，并將其應用于點云格網化濾波分析，計算結果表明，該算法在快速自動生成DEM方面精度較高，提高了三維激光掃描數據獲取的精度。然而，如何進一步提高針對海量點云數據的濾波效率、提升濾波算法自動化程度以及控制誤差，將是點云濾波算法今后的研究重點。

2）通過現場監測試驗，對點云數據進行提取分析，將試驗得出的變形數據與傳統方法量測數據進行對比，結果表明，基于Kriging濾波算法的三維激光掃描技術能夠較為高效、準確地獲得隧道變形數據。

3）三維激光掃描技術在隧道變形監測中的抗干擾能力和精度還有待進一步優化和提高。在一定的監測環境條件下，可以對拱頂沉降、周邊收斂和軸線偏差等進行分析，并及時預警圍巖穩定性異常及危險情況，防止地質災害的發生，為隧道工程安全施工保駕護航。