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

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

0 引言

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

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

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

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

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

本研究利用的LiDAR數(shù)據(jù)由ALS50-II系統(tǒng)(詳細(xì)參數(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],該類算法的應(yīng)用最為廣泛，常用的有不規(guī)則三角網(wǎng)濾波算法、數(shù)學(xué)形態(tài)法和基于坡度的濾波方法，倘若地形高差起伏較大(如陡崖、地裂縫等)，該類算法容易將坡面點分離為地物點； ②基于激光腳點回波強度信息的濾波算法[8],可以在指定回波強度閾值范圍內(nèi)剔除或者分離出感興趣的激光數(shù)據(jù),但是只是閾值范圍內(nèi)激光離散點云數(shù)據(jù)的分離，需結(jié)合其他算法進(jìn)一步的分類，以得到不同目標(biāo)對象。

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

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

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

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

2.2 影像的正射鑲嵌

航空數(shù)字影像經(jīng)過正射糾正得到的正射鑲嵌影像，已經(jīng)應(yīng)用到變形分析、土壤濕度分析以及裂隙構(gòu)造分析等領(lǐng)域[9]。實驗選取了17張影像進(jìn)行正射鑲嵌，初始外方位元素由POS系統(tǒng)獲取，差分GPS獲得像機航攝中心的瞬間位置(X，Y，Z)，慣導(dǎo)系統(tǒng)獲得像片姿態(tài)(φ，ω，κ)。正射鑲嵌的流程為： 首先，利用多項式進(jìn)行像機畸變校正[10]； 然后，加入控制點信息，空中三角測量獲取加密點，光束法區(qū)域網(wǎng)平差解算每張影像的外方位元素； 再由空三測量計算得到的地面離散點生成DEM，進(jìn)行單片正射糾正； 采用OrthoVista軟件進(jìn)行影像的自動勻色，調(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ù)，并挖掘這些參數(shù)對地裂縫信息提取及分析的能力。

3.1 山體陰影圖

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

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

3.2 坡度坡向圖

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

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

,

(1)

,

(2)

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

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

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

(4)

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

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

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

圖2 正射鑲嵌影像

首先，對激光點云數(shù)據(jù)進(jìn)行濾波處理，設(shè)定3個閾值參數(shù)分別為：Terrainangle=88°，Interationdistance=2 m ，Interationangle= 8°，進(jìn)行不規(guī)則三角網(wǎng)濾波算法獲取初始地面點； 然后，設(shè)定0.6 m的高程閾值參數(shù)，提取低矮植被加入到地面點集中； 最后，分離灰度反射強度為22～30、高程為0.6～1.2 m的植被點到地面點集，以保留陡峭的坡面點。采用目前應(yīng)用最廣泛的商業(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ù)字?jǐn)z影測量手段提取的地面點較少，構(gòu)建的DEM(圖5(b))比較粗糙，難以刻畫地表真實形態(tài)，尤其是地裂縫等地表微特征。利用TerraSoild軟件濾波分離地面點時，地裂縫的坡面點被誤分為植被點，因此生成的DEM(圖5(c))雖然在視覺效果上更能突顯地裂縫的位置，但是卻影響后續(xù)地裂縫提取的精度。圖5(d)為采用本文提出的濾波方法分離地面點派生的DEM，能夠更細(xì)致地刻畫地表形態(tài)，有效地突顯地裂縫等地表微特征。

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

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

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

基于生成的DEM，設(shè)定太陽方位角315°及太陽高度角45°，得到的定性參數(shù)山體陰影圖(圖6(a))，其表達(dá)地表形態(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)定性，陽坡的侵蝕強度往往遠(yuǎn)高于陰坡。

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

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

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

Fig.8CombinationofgroundfissurecontourandTIN(left)andDOM(right)

提取地裂縫的準(zhǔn)確位置之后，可通過點云的剖面信息瀏覽，來分析地裂縫的微細(xì)節(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處坡面都為坡度陡峭、雨水侵蝕較為嚴(yán)重的陽坡，坡面的土壤有塌陷，地裂縫有繼續(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)為了提取地裂縫的矢量線，得到地裂縫的長度及準(zhǔn)確位置信息，本文追蹤TIN剖面最小曲率等值線，獲取地裂縫矢量線，再與正射鑲嵌影像相疊合，定性分析確定其位置，結(jié)果比較精確。并通過對地裂縫剖面信息與輔助識別參數(shù)的綜合分析，得出了研究區(qū)地裂縫有繼續(xù)發(fā)育可能的結(jié)論。

參考文獻(xiàn)(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] 張小紅.機載激光雷達(dá)測量技術(shù)理論與方法[M].武漢:武漢大學(xué)出版社,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].測繪科學(xué),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.