利用機載LiDAR數據提取與分析地裂縫

肖春蕾， 郭兆成， 張宗貴， 李遷， 尚博譞， 吳芳

(1.中國國土資源航空物探遙感中心，北京 100083; 2.國土資源部航空地球物理與遙感地質重點實驗室，北京 100083)

0 引言

地裂縫是內、外地質營力或人類活動綜合作用形成的具有一定長度、寬度和深度的地面裂隙[1]。地裂縫不僅對各類工程建筑、交通設施及土地資源造成破壞，而且會導致一系列的生態環境問題。地裂縫的形成機理主要有構造成因(如美國亞利桑那州南部的地震活動導致了已有斷裂破裂面的重新復活)、地下水(或地下礦產)開采過量成因及構造與地下水(或地下礦產)開采復合成因等3種[2]。

近年來，遙感手段逐步被應用于地裂縫調查及其變化速率、微地貌變化的研究。利用差分干涉合成孔徑雷達(InSAR)能夠獲取地表詳細的位移信息，但植被變化會造成信號的去相關，使其對植被覆蓋區活動地裂縫的識別與監測效果不佳[3]。同樣，在植被高覆蓋區使用高分辨率的光學圖像提取地裂縫的分布非常困難[4]。航空數字影像雖然可以提供地表的紋理信息，但是通過數字攝影測量方法獲取的數字地形模型(digital terrain model，DTM)往往不能消除植被的影響，難以獲取真實的地裂縫微地貌信息[5]。機載激光掃描系統(airborne laser scanning，ALS)不受陰影及太陽高度角限制，能快速獲取大面積、精細的地表三維數據。激光脈沖能夠穿透植被，在植被覆蓋區獲取的高精度DTM，明顯優于航攝圖像和衛星圖像。同時，LiDAR數據能夠最大程度地減小地形切割帶來的陰影影響[6]，通過其衍生的山體陰影圖能夠得到詳細的地表參數，這使得在植被覆蓋比較茂密的情況下，提取地裂縫等線性特征成為可能。

本文基于機載LiDAR激光點云數據，構建能刻畫微地貌特征的高精度數字高程模型(digital elevation model,DEM)，設定定性及定量的識別參數，挖掘其在地裂縫提取及分析方面的能力，依據所提取的參數和微特征信息分析地裂縫的穩定性。

1 研究區概況及數據源

研究區位于湖南冷水江市，屬于湖南省地質災害比較典型的區域，地面塌陷、地裂縫現象比較常見。浪石灘位于冷水江市西南2 km處，是與新化縣接壤的資江右畔丘陵區域。該區地裂縫、地表塌陷始于1987年，由于地下采煤造成坡腳和坡體深部采空，以及相應的礦區排水疏干，誘發出現地裂縫和地表塌陷等災害。近年來，暴雨、洪水等因素加快了地表的變形發展，其地裂陷變形機制以沉陷拉裂為主。準確調查和評估地裂縫的分布及其發展狀態是避讓和防止次生災害的前提，但由于該區域植被茂密，利用地面調查及其他光學影像解譯等方法難以全面地獲取區域地裂縫與地表塌陷的空間分布狀況。

本研究利用的LiDAR數據由ALS50-II系統(詳細參數見表1)所獲取，點云密度6.8點/ m2，總數為6 039 290。該系統配備RCD中幅面像機，同時獲取了測區的影像數據(空間分辨率10 cm)。

表1 獲取實驗數據采用的機載激光掃描參數

2 數據處理

2.1 激光點云濾波

激光點云濾波的目的是剔除地物點數據，得到裸露的地面點，以便于DEM的建立。機載LiDAR數據的濾波算法可以分為2大類： ①基于高程突變的濾波算法[7],該類算法的應用最為廣泛，常用的有不規則三角網濾波算法、數學形態法和基于坡度的濾波方法，倘若地形高差起伏較大(如陡崖、地裂縫等)，該類算法容易將坡面點分離為地物點； ②基于激光腳點回波強度信息的濾波算法[8],可以在指定回波強度閾值范圍內剔除或者分離出感興趣的激光數據,但是只是閾值范圍內激光離散點云數據的分離，需結合其他算法進一步的分類，以得到不同目標對象。

本文的目的是提取植被覆蓋茂密區域地裂縫的信息，所以采用基于高程突變和激光回波強度信息相結合的濾波算法來保留地裂縫的坡面點，且分離低矮的植被點加入到地面點集中，以保證植被茂密地區地表形態的完整性。具體實現過程如下：

1)基于不規則三角網的濾波。利用不規則三角網濾波算法提取地面點，首先獲取一些認為是地面點的低點構建一個初始的稀疏不規則三角格網，設定3個閾值：①最大地形坡度(terrain angle)，表達地形的起伏情況； ②最大內插距離值(interation distance)，定義內插點到三角網的距離； ③夾角(interation angle)，表達內插點和其最鄰近地面點連線與三角網之間的夾角。對非地面點進行判斷，每次將滿足最大內插角度，到三角面的距離小于給定閾值的點納入三角網中，同時高程閾值的選擇要隨著迭代次數的增加而適當減小，重復多次直到不再有新點加入為止。

2)高程濾波。在不規則三角網濾波得到非地面點的基礎上，通過設定高程閾值，獲取一定低高度的低矮植被點作為地面點來進行DEM的構建。將數字地面模型(digital surface model,DSM)離散點云數據減去地面點數據可得到nDSM(normalized DSM)。在nDSM中，地表地物(如建筑物、車、樹等)視為在同一水平面上，可以通過選取一定的高度閾值從nDSM中提取所需的低矮植被點。

3)利用激光點云回波強度信息分離坡面點。分析由不規則三角網及高程濾波得到的非地面點的回波強度信息，設定合理的回波強度閾值，提取陡峭坡度的激光腳點。

2.2 影像的正射鑲嵌

航空數字影像經過正射糾正得到的正射鑲嵌影像，已經應用到變形分析、土壤濕度分析以及裂隙構造分析等領域[9]。實驗選取了17張影像進行正射鑲嵌，初始外方位元素由POS系統獲取，差分GPS獲得像機航攝中心的瞬間位置(X，Y，Z)，慣導系統獲得像片姿態(φ，ω，κ)。正射鑲嵌的流程為： 首先，利用多項式進行像機畸變校正[10]； 然后，加入控制點信息，空中三角測量獲取加密點，光束法區域網平差解算每張影像的外方位元素； 再由空三測量計算得到的地面離散點生成DEM，進行單片正射糾正； 采用OrthoVista軟件進行影像的自動勻色，調節匹配相鄰影像的顏色和亮度，計算得到每張片子的輻射校正參數[11]； 最后，利用OrthoVista軟件的縫合線和自動羽化功能，得到無縫的正射鑲嵌圖。

3 基于LiDAR的地裂縫識別參數獲取

通過濾波算法分離出的地面點數據，可構建不規則三角網(triangulated irregular network,TIN)，反距離加權內插可以得到DEM。本文以上述方法構建的TIN和DEM數據為基礎，設定定性及定量識別參數，并挖掘這些參數對地裂縫信息提取及分析的能力。

3.1 山體陰影圖

山體陰影圖利用假想光源對地表進行照射，產生地形表面的陰影圖，對實際地形特征進行“逼真”的立體模擬，增強地面的起伏感。它能夠凸顯很多微地貌特征，尤其是一些線性特征，比如地裂縫、線性構造、斷層崖等。本文采用山體陰影圖作為地裂縫識別的定性參數。

山體陰影圖的計算有3個重要參數： ①太陽方位角，即假想光源發出的光線在水平面上的投影與正北方向的夾角，確定其值需分析試驗區地裂縫的主體走向及河流山脈等的走向； ②太陽高度角，即假想光源發出的光線與水平面的夾角，一般選45°； ③地表灰度值，取值范圍為0～255，或者依據高程變化情況選擇最為合適的值。

3.2 坡度坡向圖

地面上某點的坡度表示了地表在該點的傾斜程度。根據坡度起伏變化可以確定地裂縫、泥石流或嚴重的土壤侵蝕區。坡向(aspect)定義為坡面法線在水平面上的投影與正北方向的夾角，可確定地裂縫的走向。本文將坡度(slope)和坡向圖作為地裂縫識別的定量參數。

地面上某點的坡度(S)和坡向(A)為地形曲率在東西(Y)、南北(X)方向上高程變化率的函數，即

,

(1)

,

(2)

式中：fx為X方向高程變化率；fy為Y方向高程變化率。通常，fx和fy的求解，是在3×3的移動窗口中(圖1)，通過數值微分或局部曲面擬合方法進行。

圖1 差分DEM示意圖

坡度運用三階反距離平方權差分算法計算，即

(3)

式中：d為格網分辨率；Zi(i=1,2,…，9)為中心點“5”周圍i格網點的高程。

3.3 線性提取

目前，利用LiDAR點云數據自動或半自動提取地表斷裂線的研究比較少。Brügelmann利用LiDAR點云高程值重采樣為距離圖像，以圖像處理的方法計算像素的二階導數來提取斷裂點的候選點，并進一步擬合成斷裂線[12]，這種方法不適應于有高差起伏的山區。Briese通過局部平面擬合相交的方法來提取斷裂線[13]，該方法直接對離散LiDAR點云進行處理，較好地保持了斷裂線的精度，但是需要人工給出斷裂線的起始位置和方向，并且進行面片擬合相交逐步提取，計算量較大。本文采用基于曲率等值線的方法提取地裂縫矢量線，利用離散的三維地面點云數據構建TIN格網，追蹤剖面最小曲率等值線，該方法即能較好地保持地裂縫的精度，又可減少相對于離散點云的計算量。

地形表面曲率是局部地形曲面在各個界面方向上的形狀、凹凸變化的反映，顯示出地形曲面在不同方向上的結構和形態特征。Wood提出用縱向曲率(longitudinal curvature)和斷面曲率(cross section curvature)來進行地形特征的識別和提取[14]。地形剖面曲率的實質為地面坡度的變化率，即

(4)

4 實驗結果與分析

4.1 正射鑲嵌及點云濾波精度

研究區位于資江右畔的丘陵區域，由正射鑲嵌后的影像(圖2)可以看出，研究區植被覆蓋茂密，建筑物、農田較多。在研究區布設了6個控制點來評定正射影像的精度，X，Y，Z這3個方向的均方根誤差RMSE分別為0.282 m，0.259 m和0.883 m。

圖2 正射鑲嵌影像

首先，對激光點云數據進行濾波處理，設定3個閾值參數分別為：Terrainangle=88°，Interationdistance=2 m ，Interationangle= 8°，進行不規則三角網濾波算法獲取初始地面點； 然后，設定0.6 m的高程閾值參數，提取低矮植被加入到地面點集中； 最后，分離灰度反射強度為22～30、高程為0.6～1.2 m的植被點到地面點集，以保留陡峭的坡面點。采用目前應用最廣泛的商業軟件Terrasoild提取的地面點僅占點云總數的13.1%(表2)，且有地裂縫的坡面點被誤分為植被點的現象(圖3)。

表2 不同方式獲取的地面點對比

(a) TerraSoild軟件提取的地面點和植被點(b) 激光點云的空間分布

為了更好地識別與提取地裂縫，本文在基于高程突變濾波的基礎上，分析裸地反射強度信息，保留了地裂縫的坡面點(圖4)，分離的地面點數目占點云總數的43.2%(表2)。

(a) 本文方法提取的地面點和植被點(b) 激光點云的空間分布

利用濾波分離出的地面點構建DEM，為后續地裂縫的提取與識別提供數據基礎。原始激光點云數據內插生成的DSM見圖5(a)。由于地表地物(植被、房屋等)的影響，很難辨識地裂縫的位置。在植被覆蓋的區域，利用數字攝影測量手段提取的地面點較少，構建的DEM(圖5(b))比較粗糙，難以刻畫地表真實形態，尤其是地裂縫等地表微特征。利用TerraSoild軟件濾波分離地面點時，地裂縫的坡面點被誤分為植被點，因此生成的DEM(圖5(c))雖然在視覺效果上更能突顯地裂縫的位置，但是卻影響后續地裂縫提取的精度。圖5(d)為采用本文提出的濾波方法分離地面點派生的DEM，能夠更細致地刻畫地表形態，有效地突顯地裂縫等地表微特征。

(a) 原始激光點云派生出的DSM(b) 傳統數字攝影測量方法提取的DEM

(c) Terrasoild軟件濾波提取地面點所派生出的DEM(d) 本文方法提取地面點所派生的DEM

4.2 識別參數及地裂縫微特征分析

基于生成的DEM，設定太陽方位角315°及太陽高度角45°，得到的定性參數山體陰影圖(圖6(a))，其表達地表形態的立體能力比較強。對于植被覆蓋茂密的區域，較正射鑲嵌影像(圖6(b))而言，山體陰影圖(圖6(c))更能清晰地反映地裂縫的起伏及破碎特征。

(a) 實驗區山體陰影圖

(b) 局部有分布地裂縫處的正射影像(c) 局部有分布地裂縫處的山體陰影圖((a)圖黑色方框范圍內)((a)圖黑色方框范圍內)

圖7為研究區的坡度和坡向圖。地裂縫分布的區域為坡度比較大(約在30°～50°之間)的地區，分布規律性較強，地裂縫處的坡度最小為40°，最大為65°，大部分位于50°～60°之間。坡向表示斜坡體表面的任一高程值的變化量的最大變化方向，由坡向圖可以看出，地裂縫的坡向為NW方向，位于南坡區域，即陽坡，而坡向決定斜坡的穩定性，陽坡的侵蝕強度往往遠高于陰坡。

對于地裂縫線性提取，首先用獲取的地面點數據構建TIN，追蹤曲率的最小值Kmin，得到地裂縫矢量線(圖8)，長度分別為(從左到右)24.8 m，43.2 m和81.8 m，總長149.8 m。將矢量線與正射影像圖相疊合，如圖8(b)所示，可以看出地裂縫提取的位置比較準確。

(a) 地裂縫矢量線與TIN格網疊合圖(b) 地裂縫矢量線與DOM疊合圖

圖8地裂縫矢量線與TIN格網及DOM疊合圖

Fig.8CombinationofgroundfissurecontourandTIN(left)andDOM(right)

提取地裂縫的準確位置之后，可通過點云的剖面信息瀏覽，來分析地裂縫的微細節，如深度、寬度及幾何形態等(圖9)。

(a) 地裂縫A處影像(左)及其點云剖面(右)(b) 地裂縫B處影像(左)及其點云剖面(右)

(c) 地裂縫C處影像(左)及其點云剖面(右)(d) 地裂縫D處影像(左)及其點云剖面(右)

地裂縫的形態為規則的V字形(圖9(a)(c)右)； 由圖9(b)(c)可以發現，這2處地裂縫的坡面有變緩的現象，而根據地裂縫提取輔助識別參數判別，2處坡面都為坡度陡峭、雨水侵蝕較為嚴重的陽坡，坡面的土壤有塌陷，地裂縫有繼續發育的可能。

5 結論

本文以機載LiDAR(ALS)數據為基礎，開展激光點云數據對地裂縫識別提取能力的研究。主要結論如下：

1)本文依次基于不規則三角網濾波、高程濾波及回波信息強度濾波提取地面點，很好地保留了地裂縫陡峭邊的坡面點，保證了地表微特征的完整性。

2)提取了地裂縫輔助識別參數，定性參數山體陰影圖能突顯地裂縫的地貌特征，以確定地裂縫分布區域，由定量識別參數坡度坡向圖可確定地裂縫發生在坡度為30°～50°的區域，坡向為NW方向，分布規律性較強。

3)為了提取地裂縫的矢量線，得到地裂縫的長度及準確位置信息，本文追蹤TIN剖面最小曲率等值線，獲取地裂縫矢量線，再與正射鑲嵌影像相疊合，定性分析確定其位置，結果比較精確。并通過對地裂縫剖面信息與輔助識別參數的綜合分析，得出了研究區地裂縫有繼續發育可能的結論。

參考文獻(References)：

[1] Jachens R C,Holzer T L.Differential compaction mechanism for earth fissures near Casa Grande,Arizona[J].Geological Society of America Bulletin,1982,93(10):998-1012.

[2] Holzer T L,Davis S N,Lofgren B E.Faulting caused by groundwater extraction in south-central Arizona[J].Journal of Geological Research,1979,84(B2):603-612.

[3] Belardinelli M E,Sandri L,Baldi P.The major event of the 1997 Umbria-Marche(Italy)sequence:What could we learn from DInSAR and GPS data[J].Geophysical Journal International,2003,153(1):242-252.

[4] Niebergall S,Loew A,Mauser W.Object-orientated analysis of very high resolution QuickBird data for mega city research in Delphi/India[C]//Proceedings of the Urban Remote Sensing Joint Event,Paris.ISBN:1-4244-0712-5.2007:8-13,IEEE07EX1577.

[5] Baltsavias E P.A comparison between photogrammetry and laser scanning[J].ISPRS J Photogrammetry and Remote Sensing,1999,54(2):83-94.

[6] Carter W E,Shrestha R L,Slatton K C.Geodetic laser scanning[J].Phys Today,2007:60(12):41-47.

[7] 張小紅.機載激光雷達測量技術理論與方法[M].武漢:武漢大學出版社,2007:93-116。

Zhang X H.The Theory and Technical Method of Airborne Laser Radar Measurement[M].Wuhan:Wuhan University Process,2007:93-116.

[8] 路興昌,張雪霞.基于回波強度和采樣點距離的點云濾波研究[J].測繪科學,2009,34(6):196-197.

Lu X C,Zhang X X.Study on points cloud filtering based on reflectance intensity and range[J].Science of Surveying and Mapping,2009,34(6):196-197.

[9] Niethammer U,Rothmund S,Joswig M.UAV-based remote sensing of the slow moving landslide Super-Sauze[C]//Malet J P,Rema tre A,Boogard T.Proceedings of the International Conference on Landslide Processes:From Geomorphologic Mapping to Dynamic Modelling.Strasbourg:CERG Press,2009:69-74.

[10]Niethammer U,Rothmund S,James M R,et al.UAV-based remote sensing of landslides[C]//International Archives of Photogrammetry,Remote Sensing and Spatial Information Sciences,Vol.XXXVIII,Part 5 Commission V Symposium,Newcastle upon Tyne,UK.2010:496-501.

[11]OrthoVista.Official OrthoVista software homepage[EB/OL].http://www.orthovista.com.(accessed 1 August 2010).

[12]Brügelmann R.Automatic breakline detection from airborne laser range data[J].International Archives of Photogrammetry and Remote Sensing,2000,33(B3):109-115.

[13]Briese C.Three-dimensional modelling of breaklines from airborne laser scanner data[J].International Archives of Photogrammetry Remote Sensing,2004,35(B3):109-I102.

[14]Wood J.The geomorphological characterization of digital elevation models[D].Leicester,UK:University of Leicester,1996.