利用機載LiDAR數(shù)據(jù)提取與分析地裂縫

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

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

0 引言

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

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

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

1 研究區(qū)概況及數(shù)據(jù)源

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

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

表1 獲取實驗數(shù)據(jù)采用的機載激光掃描參數(shù)

2 數(shù)據(jù)處理

2.1 激光點云濾波

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

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

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

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

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

2.2 影像的正射鑲嵌

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

3 基于LiDAR的地裂縫識別參數(shù)獲取

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

3.1 山體陰影圖

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

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

3.2 坡度坡向圖

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

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

,

(1)

,

(2)

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

圖1 差分DEM示意圖

坡度運用三階反距離平方權(quán)差分算法計算，即

(3)

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

3.3 線性提取

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

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

(4)

4 實驗結(jié)果與分析

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

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

圖2 正射鑲嵌影像

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

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

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

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

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

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

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

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

4.2 識別參數(shù)及地裂縫微特征分析

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

(a) 實驗區(qū)山體陰影圖

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

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

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

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

圖8地裂縫矢量線與TIN格網(wǎng)及DOM疊合圖

Fig.8CombinationofgroundfissurecontourandTIN(left)andDOM(right)

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

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

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

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

5 結(jié)論

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

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

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

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

參考文獻(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] 張小紅.機載激光雷達測量技術(shù)理論與方法[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.