許 穎，岳東杰，孟德友

（1.河南財(cái)經(jīng)政法大學(xué)資源與環(huán)境學(xué)院，鄭州 450046；2.河海大學(xué)地球科學(xué)與工程學(xué)院，南京 211100）

1 研究背景

侵蝕溝作為一種常見的土壤侵蝕類型，侵蝕量大，速度快，在土壤侵蝕中占據(jù)著重要的位置［1］。隨著對地觀測技術(shù)和點(diǎn)云數(shù)據(jù)處理能力的發(fā)展，機(jī)載激光雷達(dá)（Light Detection and Ranging，LiDAR）的出現(xiàn)，為獲取高時(shí)空分辨率的地球空間信息提供了全新的技術(shù)手段［2－3］。但是，點(diǎn)云數(shù)據(jù)量、特征差異和自然地物的復(fù)雜性嚴(yán)重制約著機(jī)載LiDAR技術(shù)數(shù)據(jù)后處理的發(fā)展。不像建筑物、樹木、道路、橋梁、電力線［4］等有著明顯地表特征的地物［5］，侵蝕溝嵌入在地表中，很難形成明顯的表面特征，并且它們的變化范圍、形成過程和外觀不規(guī)則性更是增加了檢測的難度［6］。

在LiDAR技術(shù)出現(xiàn)的初期，Betts等［7］借助由Li-DAR數(shù)據(jù)構(gòu)建的數(shù)字高程模型（DEM）提取局部溪谷的形態(tài)信息。Perroy等［8］應(yīng)用機(jī)載和地面激光雷達(dá)相結(jié)合的方法進(jìn)行溝蝕評估，由于地面激光雷達(dá)掃描侵蝕溝底部和側(cè)部有一定的限制，因此機(jī)載激光雷達(dá)技術(shù)更適合大范圍侵蝕溝監(jiān)測。Evans等［9］提出從高精度DEM中提取侵蝕溝，并在侵蝕溝邊緣內(nèi)插出其深度。Mason等［10］針對點(diǎn)云利用邊界檢測算子得到高梯度值，此方法在復(fù)雜地形中有一定的局限性。Rutzinger等［11］基于最大曲率分割算法提出對象驅(qū)動(dòng)方法進(jìn)行興趣點(diǎn)檢測，算法適于深且寬的輪廓特征，對紋理特征復(fù)雜的溝壑地區(qū)效果不佳。

以上基于DEM對提取類似侵蝕溝的微地貌特征還有一定的局限性，基于點(diǎn)云本身的操作會(huì)更符合地物地形特征?；诖?，本文提出一種基于機(jī)載LiDAR點(diǎn)云的不同尺度下的表面特征差異侵蝕溝點(diǎn)云分割方法?？紤]到侵蝕溝點(diǎn)云信息位于表面特征復(fù)雜山區(qū)，改進(jìn)基于點(diǎn)云的邊緣分割算法，基于法向量和曲率的變化關(guān)系來進(jìn)行侵蝕溝點(diǎn)云分割。

2 特征差異分割方法

特征差異分割方法的思想是基于侵蝕溝處于地表發(fā)生突變的地方，選取不同半徑，得到不同表面特征參數(shù)，設(shè)定一定閾值分割得到侵蝕溝點(diǎn)云。曲面的幾何信息是指能夠反映曲面局部特性的數(shù)據(jù)量［12］，其主要包括曲面的法向量、高斯曲率、平均曲率和曲面變化度等［13］。本文將法向量和曲率作為點(diǎn)云的表面特征參數(shù)。

2.1 點(diǎn)云特征分析

侵蝕溝的地貌形態(tài)一般表現(xiàn)為縱剖面與所在坡面基本一致，橫剖面一般為U或V字形，因此通過切溝兩邊的溝沿線及U字形切溝溝底的兩邊或V字形切溝的溝底線可很好地反映切溝的三維形態(tài)［14］。由于溝壁斗力，溝谷深壑，因而侵蝕溝在水力、重力的雙動(dòng)力混合作用下，往往具有隨機(jī)性特征?？偨Y(jié)來看，侵蝕溝的點(diǎn)云特征主要表現(xiàn)以下幾個(gè)方面：①從侵蝕溝形態(tài)上來看，其溝沿線點(diǎn)云表現(xiàn)出不規(guī)則、不對稱的特征，溝壁點(diǎn)云密度分布不均勻，增加了侵蝕溝提取的難度，溝沿線和溝底點(diǎn)云高程差別明顯；②從與其他地物差別來看，比如與建筑物、樹木、道路、橋梁等地物表現(xiàn)出的點(diǎn)云特征進(jìn)行比較，侵蝕溝表現(xiàn)出的點(diǎn)云嵌入在地形特征中，很難形成明顯的表面特征，并且一般侵蝕溝發(fā)生在山區(qū)，地形起伏比較大，地物特征不規(guī)則；③侵蝕溝的動(dòng)態(tài)變化更是增大了檢測的難度。

2.2 法向量與曲率估計(jì)

法向量反映了曲面在某點(diǎn)處的幾何特征，原始的散亂點(diǎn)云本身沒有法向量，點(diǎn)云數(shù)據(jù)的法向量主要是通過點(diǎn)的鄰域來確定［15］。因此，點(diǎn)的法向量近似為鄰域表面的切平面法向量估計(jì)問題，對于點(diǎn)云中的每個(gè)點(diǎn)pi＝( xi，yi，zi)T，獲取由距離半徑確定的k個(gè)相鄰點(diǎn)，每個(gè)點(diǎn)pi的協(xié)方差矩陣可以表示為

式中：k代表點(diǎn)pi鄰域中所有點(diǎn)的個(gè)數(shù)；p－表示鄰域中所有點(diǎn)的質(zhì)心；λj是協(xié)方差矩陣的第j個(gè)特征值；v→—j是第j個(gè)特征向量。

用數(shù)學(xué)方法計(jì)算出的法向量解決不了它的符號(hào)問題，因此得到法向量的方向是任意的。為了解決這個(gè)問題，我們需要進(jìn)行法向量一致性調(diào)整，選取視點(diǎn)vp，使其滿足方程：

2.3 特征差異原理

特征差異處理為無序大量三維點(diǎn)云數(shù)據(jù)的處理提供了一種高效率、多尺度方法，以點(diǎn)集P中的每個(gè)點(diǎn)pi為圓心，選取 2個(gè)不同的半徑r1，r2（r1＞r2），由上述方法會(huì)分別得到2個(gè)不同的法向量n1( pi，r1)，n2( pi，r2)和曲率值，如圖1所示。表面特征差異定義為

2.4 基于特征差異的侵蝕溝提取

基于特征差異的侵蝕溝分割方法的主要工作是在特征差異的估算上。本文鄰域點(diǎn)的確定采用固定距離法，其中距離的選擇可以根據(jù)點(diǎn)云分辨率（點(diǎn)間距的大?。﹣盱`活確定，一般情況下，領(lǐng)域點(diǎn)數(shù)量設(shè)定在10～20之間比較好。

基于特征差異的侵蝕溝提取具體的算法如下：

圖1 半徑不同產(chǎn)生的特征差異Fig.1 Characteristic difference produced by different radii

（1）選取半徑r1得到鄰域點(diǎn)，采用主成分分析方法擬合得到點(diǎn)云法向量和曲率值，并對法向量進(jìn)行一致性調(diào)整。

（2）選取半徑r2得到鄰域點(diǎn)，采用主成分分析方法擬合得到點(diǎn)云法向量和曲率值，并對法向量進(jìn)行一致性調(diào)整（r1＞r2）。

（3）根據(jù)法向量和曲率的不同，計(jì)算每個(gè)點(diǎn)的特征差異值。

（4）根據(jù)侵蝕溝的點(diǎn)云特征，設(shè)定誤差閾值（0.25），提取差異特征大的值。

（5）基于距離（3 m）聚類的方法，剔除小塊特征差異大的點(diǎn)，進(jìn)而得到侵蝕溝點(diǎn)云。

因此，基于不同尺度下的表面特征差異侵蝕溝分割算法流程如圖2所示。

圖2 點(diǎn)云分割流程Fig.2 Flow chart of point cloud segmentation

3 試驗(yàn)及結(jié)果分析

3.1 數(shù) 據(jù)

為了測試本文算法的可行性，在VC＋＋環(huán)境下實(shí)現(xiàn)了該算法。試驗(yàn)數(shù)據(jù)為甘肅省嘉峪關(guān)市某區(qū)域，選取區(qū)域共358 977個(gè)點(diǎn)中的一部分，經(jīng)過噪聲處理和濾波后地面點(diǎn)數(shù)為11 038，點(diǎn)云平均密度為14個(gè)點(diǎn)／m2，點(diǎn)云數(shù)據(jù)的三維顯示如圖3所示。

圖3 點(diǎn)云三維顯示圖Fig.3 Three-dimensional display of point cloud

3.2 試驗(yàn)過程

本文根據(jù)固定距離法選擇鄰域點(diǎn)，點(diǎn)云密度14個(gè)點(diǎn)／m2，點(diǎn)云密度比較高，半徑分別取1 m和2 m選擇鄰域點(diǎn)。試驗(yàn)過程如下所述。

（1）對地面點(diǎn)云進(jìn)行歸一化處理，提高計(jì)算效率。

（2）對每一點(diǎn)云取鄰域半徑分別為1 m和2 m，計(jì)算法向量和曲率值，進(jìn)而得到特征差異值，設(shè)定誤差閾值為0.02，得到分割后結(jié)果，如圖4（a）所示。

（3）設(shè)定聚類距離為2 m，得到侵蝕溝點(diǎn)云，如圖 4（b）所示。

圖4 侵蝕溝點(diǎn)云提取過程Fig.4 Process of extracting point cloud of gully

對圖3截取的原始點(diǎn)云進(jìn)行極大曲率特征點(diǎn)提?。?4］，提取結(jié)果如圖5（a）所示，于是特征點(diǎn)與圖 4分割得到的侵蝕溝點(diǎn)云合并得到完整溝壑點(diǎn)云，如圖 5（b）所示。

基于提取的侵蝕溝點(diǎn)云，對其進(jìn)行定量分析，分別得到侵蝕溝的長度約為32 m，寬度約為2.5 m，如圖6所示。侵蝕溝定量數(shù)據(jù)的獲得，為以后侵蝕溝定量化評價(jià)和水土流失侵蝕量計(jì)算提供了一種新的方法。

圖5 完整溝壑點(diǎn)云Fig.5 Integrated point cloud of gully

圖6 侵蝕溝參數(shù)計(jì)算Fig.6 Calculation of parameters of gully

3.3 結(jié)果分析

視覺上來看，本文提取的侵蝕溝點(diǎn)云與原始點(diǎn)云形態(tài)相似，符合溝的形態(tài)特征，如圖5所示。另一方面，對原始三維點(diǎn)云和分割得到的溝壑點(diǎn)云疊加，如圖7（a）所示，可以看出所分割點(diǎn)云完全落在原始三維點(diǎn)云的溝壑位置上，符合溝壑點(diǎn)云特征；對原始點(diǎn)云建立DEM，通過侵蝕溝三維點(diǎn)云與DEM疊加，如圖7（b）所示，可以看出，侵蝕溝溝沿線與DEM的高程變化相吻合。以上試驗(yàn)結(jié)果表明，采用特征差異方法能從點(diǎn)云中準(zhǔn)確分割出比較完整的侵蝕溝點(diǎn)云信息。

圖7 點(diǎn)云分割結(jié)果驗(yàn)證Fig.7 Verification of point cloud segmentation results

4 結(jié) 論

本文提出的基于特征差異的機(jī)載LiDAR侵蝕溝點(diǎn)云分割方法，考慮到侵蝕溝點(diǎn)云特征以及所處地形特性的特殊性，從僅僅考慮點(diǎn)云法向量和曲率單因素微分信息擴(kuò)展到鄰域的多尺度分析，即從它們變化鄰域的尺寸可以得到更多的曲面信息。從結(jié)果分析來看，多尺度的算法可以從散亂點(diǎn)云中分割精確的侵蝕溝點(diǎn)云，很好地保持了侵蝕溝形態(tài)特征。所提取的侵蝕溝點(diǎn)云能合理表達(dá)真實(shí)形態(tài)，為進(jìn)一步進(jìn)行侵蝕量評價(jià)提供基礎(chǔ)。

［1］ WU Y，CHENG H.Monitoring of Gully Erosion on the Loess Plateau of China Using a Global Positioning System［J］.Catena，2005，63（2／3）：154－166.

［2］ 鄭輯濤，張 濤.基于可變半徑圓環(huán)和 B樣條擬合的機(jī)載 LiDAR點(diǎn)云濾波［J］.測繪學(xué)報(bào)，2015，44（12）：1359－1366.

［3］ 張 煜，竇延娟，張曉東.機(jī)載激光雷達(dá)數(shù)據(jù)采集及數(shù)據(jù)處理［J］.長江科學(xué)院院報(bào)，2010，27（1）：13－16.

［4］ 韓文軍.機(jī)載LiDAR數(shù)據(jù)的電力塔建模方法研究［J］.長江科學(xué)院院報(bào)，2012，29（10）：122－126.

［5］ 張春亢，趙學(xué)勝，王洪斌.采用Morse理論的小尺度地形特征提取方法［J］.測繪科學(xué)技術(shù)學(xué)報(bào)，2015，（3）：266－270.

［6］ FILIN S，AVNI Y，BARUCH A，et al.Characterization of Land Degradation along the Receding Dead Sea Coastal Zone Using Airborne Laser Scanning［J］.Geomorphology，2014，206：403－420.

［7］ BETTSH D，DEROSE R C.Digital Elevation Models as a Tool for Monitoring and Measuring Gully Erosion［J］.International Journal of Applied Earth Observation and Geoinformation，1999，1（2）：91－101.

［8］ PERROY R L，BOOKHAGEN B，ASNER G P，et al.Comparison of Gully Erosion Estimates Using Airborne and Ground-based LiDAR on Santa Cruz Island，California［J］.Geomorphology，2010，118（3）：288－300.

［9］ EVANSM，LINDSAY J.High Resolution Quantification of Gully Erosion in Upland Peatlands at the Landscape Scale［J］.Earth Surface Processes and Landforms，2010，35（8）：876－886.

［10］MASON D C，SCOTT T R，WANG H J.Extraction of Tidal Channel Networks from Airborne Scanning Laser Altimetry［J］.ISPRS Journal of Photogrammetry＆Remote Sensing，2006，61（2）：67－83.

［11］RUTZINGER M，H?FLE B，PFEIFER N，et al.Objectbased Analysis of Airborne Laser Scanning Data for Natural Hazard Purposes Using Open Source Components［C／CD］∥ISPRS.International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences.Salzburg，Austria，July 4－5，2006：36（part 4／C42）.

［12］于海洋，余鵬磊，謝秋平，等.機(jī)載LiDAR數(shù)據(jù)建筑物頂面點(diǎn)云分割方法研究［J］.測繪通報(bào)，2014，（6）：20－23.

［13］劉 宇.基于微分信息的散亂點(diǎn)云拼合和分割［D］.武漢：華中科技大學(xué)，2008.

［14］BARUCH A，F(xiàn)ILIN S.Detection of Gullies in Roughly Textured Terrain Using Airborne Laser Scanning Data［J］.ISPRSJournal of Photogrammetry and Remote Sensing，2011，66（5）：564－578.

［15］TANG C K，MEDIONI G.Curvature-augmented Tensor Voting for Shape Inference from Noisy 3D Data［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2002，24（6）：858－864.