利用地形開度算子自動提取激光雷達點云斷裂線

張 良，張 帆，姜曉琦，周佳雯，王杰棟

(1. 湖北大學資源環境學院，湖北 武漢 430062； 2. 區域開發與環境響應湖北省重點實驗室，湖北 武漢 430062； 3. 武漢大學遙感信息工程學院，湖北 武漢 430079； 4. 浙江省第二測繪院，浙江 杭州 310012)

LiDAR系統是一種集激光、全球定位系統(GPS)和慣性導航系統(INS)3種技術于一身的新型傳感器，無需大量地面控制點，便能夠快速準確地獲取地表高密度、高精度的三維點云，同時具備全天時觀測能力、受天氣影響小等優點，是高精度DEM的優質數據源之一[1-2]。地形表面通常不都是光滑的或均勻變化的，在陡坎、梯田、沖溝、堤岸等處，常出現轉折或突變，被稱為地形斷裂線[3]。斷裂線是建立DEM非常重要的依據，對數字高程模型的精度有相當大的影響，尤其在地形變化大的山區或水系多的地區[4]。因此，為了有效保留地表的關鍵地形特征，為用戶提供精度更高的DEM，需要直接從LiDAR點云數據中提取地形斷裂線作為構建DEM的約束條件。

近10年來，已有學者對基于LiDAR點云的斷裂線提取進行了研究并取得了一定的成果。現有的斷裂線自動提取算法包括斷裂線候選點提取和斷裂線矢量跟蹤兩大步驟，其中斷裂線候選點的精度是斷裂線提取精度和完整性的前提和關鍵，也是本文重點探討的地方。文獻[5—6]利用二階導數等圖像邊緣算子處理LiDAR距離圖像獲取斷裂線候選點并擬合斷裂線。這類算法自動化程度高，但是自動提取的斷裂線比較破碎且錯提取率高。文獻[7]通過局部平面相交的方法直接從離散點云提取較精確的斷裂線候選點，該方法一般需要人工給出斷裂線的起始位置和方向作為斷裂線生長的初始值。文獻[8]嘗試利用三角網面片之間的法向差異提取候選斷裂點，再基于方向優先策略追蹤斷裂線以提高算法的自動化程度。文獻[9]首先利用雙閾值從距離圖像提取斷裂線概略位置，再結合原始點云數據用分段構造法獲取更精確的斷裂線。雙閾值法提取的斷裂線完整性好，但大、小閾值的設置需要多次試驗才能有效確定。

如前文所述，目前大部分斷裂線提取算法都是基于地表局部鄰域的變化開展，因此對地表突變和點云錯分敏感，精度和完整性無法滿足生產需求，因此目前有經驗的處理人員更信任結合可視化、光照分析等工具判斷地形的整體變化，進而人工繪制斷裂線。目前大部分商業軟件(如Qtmodeler)利用陰影(或與高程結合)渲染的方式表現三維地形細節。陰影渲染下，地表亮度與光源位置、地形起伏密切相關，地形起伏大的區域，部分地表細節因無法接受足夠光照而被遺漏，而部分區域則可能曝光過度[10-11]。因此很難基于陰影渲染圖提取完整的斷裂線。地形開度(topoghaphic openness)算子利用均質虛擬輻射源代替單方向點陣直射光源，從根本上突破了基于陰影渲染采用平行直射光源導致的陰影和曝光過度問題。近幾年來，作為一種新型的高分辨率地形可視化工具，在軍事、考古、地質[12-15]等需要對微地貌進行詳細解譯的領域得到廣泛應用。本文從人工解譯的角度出發，將地形開度作為一種定量描述地表整體變化的地形特征，提出基于地形開度(topographic openness)的斷裂線自動提取算法，以提高復雜地區高精度DEM的提取能力。

1 算法設計

1.1 地形開度算子

地形開度算子或SVF(sky view factor)算子是近幾年發展的利用均質輻射源代替單方向點陣直射光源的地形渲染算法。SVF算法從能量傳播的角度出發，假設光源為一個虛擬的半圓形天球，待渲染的地表位于該天球的中心并受天球中均勻的散射光所照射，在該假設下，可以認為地表上某個點的亮度取決于該點及其鄰域暴露在散射光源下的體積或表面積(如圖1(b)所示)。圖1(a)表示起伏地表的一個剖面截圖，其中A為谷底待計算亮度值的像素，扇形區域代表A能夠接收均質散射光輻射區域。因此，通過式(1)可估算出該點暴露在天球光線照射的體積Ω即該點的亮度。

(1)

式中，φ表示均值天球球體的緯度；λ表示經度。

圖1 基于輻射度的地形渲染原理

地形開度[10]則相當于SVF的簡化算子。地形開度包括正開度和負開度兩個部分，正開度用于表示山脊、背脊等隆起地貌，負開度表示地形的整體凹陷程度，著重表達山谷、溝壑等地形。正開度通常取多個方位的天頂角的平均值(一般取8個方向)。圖2表示方位角為90°和270°的天頂角剖面示意圖。首先在一定半徑內尋找點B，使得天頂角Ωi最小，即點A和點B之間的仰角γi最大。利用式(2)計算γi，再基于式(3)即可計算天頂角Ωi。通過同樣的原理計算其余方向天頂角值之后，可根據式(4)計算出該點的正開度值Op。

(2)

Ωi=90-γi

(3)

(4)

式中，γi為點A和點B之間的仰角；HA、HB分別為點A和點B的高程；dAB為點A和點B之間的水平距離；Op為正開度值；n表示所取方向的數量，一般n值取8。

圖2 天頂角計算示意圖

負開度的計算原理和正開度類似，通常取8個方位的天底角的平均值。首先在搜索區域內找到點C，使得天底角Ψi最小，即點A和點C之間的俯視角δi最大，然后利用式(5)、式(6)計算該方向的天底角Ψi(圖3為方位角90°和270°下天底角的示意圖)。計算其余方向的天底角之后，利用式(7)統計負開度值On。

(5)

ΨiL=90+δi

(6)

(7)

式中，δi為點A和點B之間的俯視角；HA、HC分別為點A和點C的高程；dAC為點A和點C之間的水平距離；On為負開度值；n為所取方向的數量，一般n值取8。

圖3 天底角計算示意圖

1.2 基于地形開度算子的斷裂線自動提取

大量的文獻指出[12-14]，地形開度算子具有兩大特點：①開度算子表現地表的主要變化，即地形結構主要變化得到加強，而緩慢的連續形變或小范圍的突變則得到抑制，因此地形開度算子渲染后的地表具有極強的立體感和縱深感，也使得地形特征和周邊地形的區別更加明顯。②開度算子計算的地表亮度僅僅與搜索半徑大小有關，不會因為光照條件的不同和地形起伏而遺漏地表細節。基于上述特點，本文在獲取開度算子的基礎上，實現斷裂線的自動提取。算法具體步驟如下：

步驟1：利用不規則三角網(Tin)加密算法實現點云濾波并提取DEM。

步驟2：基于DEM分別提取正負開度Op和On，并且通過正負開度相減得到IDEM(式(8))。如圖4(a)所示，若點在地表凹陷處，正開度小于負開度，IDEM取負值。如圖4(b)所示，點落在起伏較為平緩的表面，其正、負開度值相近，IDEM值接近零。如圖4(c)所示，若點落在斷裂附近或山脊等突起地表附近，正開度大于負開度，IDEM取正值。因此，根據IDEM可區分凸出地表、較平緩地表(含傾斜地表)和凹陷地表。

I=(Op-On)/2

(8)

圖4 正負開度和地表的關系

步驟3：通過邊緣檢測LOG算子對IDEM進行邊緣提取。對該邊緣提取結果進行二值化，利用面積閾值刪除連通區域小的邊緣像素，得到的高亮區域即為斷裂線種子點及少量噪聲。

步驟4：利用邊緣檢測LOG算子對原始DEM進行邊緣提取。可發現地形斷裂區域，地表起伏大的山坡、小的突起處都呈高亮狀態，但是地形平坦區域的亮度明顯低于其他區域，具有非常明顯的區分度。基于該特點，通過亮度分割的方式提取地形平緩區域，對其進行亮度反轉，然后利用形態學開運算得到較完整的平緩區域。

步驟5：通過疊加求并運算，利用步驟4提取的地形平坦區域進一步消除步驟3提取的噪聲，提取結果即為斷裂線種子點。

步驟6：利用文獻[8]的方法，結合原始點云采用基于平面對的分段構造法和基于三角網的直接構造法提取完整的斷裂線矢量。

2 試驗和精度評價

2.1 試驗數據和試驗工具

為了檢驗算法的有效性，筆者選取一塊中國陜西某地區的點云數據(如圖5(a)所示)進行試驗。試驗數據位于黃土高原，由于植被覆蓋較少，土質疏松，氣候干旱，導致水土流失嚴重，侵蝕切割強烈，斷裂和裂隙縱橫交錯。試驗區域范圍為0.6 km×0.6 km，平均點間距為0.6 m，激光腳點數量為1 418 228，試驗區內包含大量的不同種類的斷裂、村莊和低矮植被。利用LiDAR點云數據處理平臺LiDAR_Suite對試驗數據進行濾波和人工分類并基于LiDAR_Suite軟件進行二次開發，利用C++語言實現了本文算法。

2.2 試驗過程和結果分析

(1) 在LiDAR點云處理平臺LiDAR_Suite上實現點云精細分類，然后基于線性三角網插值算法獲取1 m分辨率的DEM，如圖5(b)所示。

(2) 基于DEM提取開度圖像IDEM。根據文獻[11]的建議，開度半徑設置為50像素。如圖5(c)所示，開度圖像中，突起地形(如斷裂線或道路橋梁邊緣)呈明顯的高亮狀態，山坡、山谷則亮度值非常低，具有非常明顯的區分度。

(3) 利用LOG算子對IDEM進行邊緣提取，其中高斯拉普拉斯卷積核的大小為25×25像素，高斯方差設為2.5，并且利用灰度直方圖自動進行二值分割，提取結果如圖5(d)所示，二值化后的結果基本上涵蓋了全部地形斷裂，以及少量由于地形突變和錯分類導致的噪聲。

(4) 利用LOG算子直接對DEM進行邊緣提取，LOG算子參數與步驟3相同，處理結果如圖5(e)所示，可見在黃土高原這種特殊地形條件下，直接基于LOG算子也能取得較好的效果，但是斷裂線完整度不如基于IDEM的提取結果，噪聲數量也明顯更大。另一方面，如圖5(e)所示，利用LOG算子處理后的DEM，地形平坦區域的亮度明顯低于其他區域，利用該特征對邊緣提取結果進行二值化處理，并且通過形態學算子提取較完整的平緩區域，結果如圖5(f)所示。

(5) 基于步驟3和步驟4的提取結果進行疊加分析，消除噪聲，得到完整的斷裂線種子點。在此基礎上，結合原始點云數據采用基于平面對的分段構造法和基于三角網的直接構造法提取斷裂線矢量。

圖5 基于開度的斷裂線提取

(6) 從兩個方面來評價本文方法自動提取斷裂線的效果：①將最終提取的斷裂線矢量和開度圖疊加，結果如圖5(g)所示，可見本文所提取的斷裂線具備非常好的完整性，同時與實際的斷裂邊緣吻合度也較好；另一方面，地形結構變化不夠明顯的區域，如圖5(g)左下角的橋梁邊緣沒有識別，主要原因在于為了更好地表達地表整體的變化，開度半徑設置偏大，因此小的結構變化不夠突出。②在LiDAR_Suite平臺通過人機交互的方式提取斷裂線，與本文算法進行對比，對比結果如圖5(h)所示，其中左側為本文算法提取的斷裂線，右側為人工提取的斷裂線，通過對比發現，兩種方式提取的斷裂線基本一致，通過本文提取的斷裂線細節更加豐富，保真度更高。

3 總結和展望

DEM的制作一直是測繪地理信息系統領域的重點之一，提取地形斷裂線并作為構建DEM的約束條件，可以有效減少地貌的失真，提供更高精度的DEM。本文從人工解譯的角度出發，引入地形開度作為一種定量描述地表整體變化的地形特征，提出了基于開度算子的斷裂線自動提取算法，突破了傳統方法僅僅基于地形局部變化考慮的局限，結合邊緣提取算子和形態學算子，能夠快速提取完整準確的斷裂線。算法簡潔有效，不需要人工干預，具備較好的推廣價值。

開度算子表示的地形細節的豐富程度與搜索半徑相關，為了提取更高精度的斷裂線，需要研究基于不同地形地貌自適應調整開度半徑的方法；同時，本文著重探討了斷裂線候選點的提取方法，如何結合開度算子，基于斷裂線候選點提取更精確的斷裂線矢量，還需要進一步研究。另外，目前對斷裂線提取精度的評價僅僅是與影像和手工提取的結果進行比較，還難以進行定量分析，這也是今后需要深入研究的問題。

參考文獻：

[1] AXELSSON P.Processing of Laser Scanner Data-Algorithms and Applications[J].Isprs Journal of Photogrammetry & Remote Sensing，1999，54(2-3):138-147.

[2] 鄭輯濤，張濤.基于可變半徑圓環和B樣條擬合的機載LiDAR點云濾波[J].測繪學報，2015，44(12):1359-1366.

[3] YANG B，HUANG R，DONG Z，et al.Two-step Adaptive Extraction Method for Ground Points and Breaklines from Lidar Point Clouds[J].Isprs Journal of Photogrammetry & Remote Sensing，2016，119:373-389.

[4] 易俐娜，劉鵬飛，喬小軍,等.結合遙感影像和DEM的線性體特征提取[J].測繪通報，2014(10):19-22.

[5] UGELMANN R B.Automatic Breakline Detection from Airborne Laser Range Data[J].International Archives of Photogrammetry and Remotes Sensing,2000,33(B3):109-115.

[6] TOSCANO G J，GOPALAM U，DEVARAJAN V.A Novel Method for Automation of 3D Hydro Break Line Generation from LIDAR Data Using Matlab[J].ISPRS-International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences，2013，XL-2/W2(1):99-104.

[7] BRIESE C.Three-dimensional Modelling of Breaklines from Airborne Laser Scanner Data[J].International Archives of Photogrammetry Remote Sensing & Spatial Information Sciences，2011，35:1097-1102.

[8] 徐景中，萬幼川，張圣望.基于機載激光雷達點云的斷裂線自動提取方法[J].計算機應用，2008，28(5):1214-1216.

[9] 彭檢貴，馬洪超，王宗躍,等.機載LiDAR點云的雙閾值自動提取斷裂線方法[J].測繪科學技術學報，2010，27(4):275-279.

[10] SHIRASAWA M，YOKOYAMA R.Visualizing Topography by Openness: A New Application of Image Processing to Digital Elevation Models[J].Photogrammetric Engineering & Remote Sensing，2002，68(3):257-266.

[11] NGUYEN H T，PEARCE J M.Incorporating Shading Losses in Solar Photovoltaic Potential Assessment at the Municipal Scale[J].Solar Energy，2012，86(5):1245-1260.

[12] DONEUS M.Openness as Visualization Technique for Interpretative Mapping of Airborne Lidar Derived Digital Terrain Models[J].Remote Sensing，2013，5(12):6427-6442.

[13] ELMAHDY S I.Topographic Openness Algorithm for Characterizing Geologic Fractures of Kuala Lumpur Limestone Bedrock Using DEM[J].Journal of Geomatics，2010，4(2):61-68.

[14] HODUL M，KNUDBY A，HO H.Estimation of Continuous Urban Sky View Factor from Landsat Data Using Shadow Detection[J].Remote Sensing，2016，8(7):568.

[15] FAVALLI M，FORNACIAI A.Visualization and Comparison of DEM-derived Parameters Application to Volcanic Areas[J].Geomorphology，2017，290: 69-84.