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

(1. 湖北大學(xué)資源環(huán)境學(xué)院，湖北 武漢 430062； 2. 區(qū)域開發(fā)與環(huán)境響應(yīng)湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430062； 3. 武漢大學(xué)遙感信息工程學(xué)院，湖北 武漢 430079； 4. 浙江省第二測(cè)繪院，浙江 杭州 310012)

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

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

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

1 算法設(shè)計(jì)

1.1 地形開度算子

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

(1)

式中，φ表示均值天球球體的緯度；λ表示經(jīng)度。

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

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

(2)

Ωi=90-γi

(3)

(4)

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

圖2 天頂角計(jì)算示意圖

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

(5)

ΨiL=90+δi

(6)

(7)

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

圖3 天底角計(jì)算示意圖

1.2 基于地形開度算子的斷裂線自動(dòng)提取

大量的文獻(xiàn)指出[12-14]，地形開度算子具有兩大特點(diǎn)：①開度算子表現(xiàn)地表的主要變化，即地形結(jié)構(gòu)主要變化得到加強(qiáng)，而緩慢的連續(xù)形變或小范圍的突變則得到抑制，因此地形開度算子渲染后的地表具有極強(qiáng)的立體感和縱深感，也使得地形特征和周邊地形的區(qū)別更加明顯。②開度算子計(jì)算的地表亮度僅僅與搜索半徑大小有關(guān)，不會(huì)因?yàn)楣庹諚l件的不同和地形起伏而遺漏地表細(xì)節(jié)?；谏鲜鎏攸c(diǎn)，本文在獲取開度算子的基礎(chǔ)上，實(shí)現(xiàn)斷裂線的自動(dòng)提取。算法具體步驟如下：

步驟1：利用不規(guī)則三角網(wǎng)(Tin)加密算法實(shí)現(xiàn)點(diǎn)云濾波并提取DEM。

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

I=(Op-On)/2

(8)

圖4 正負(fù)開度和地表的關(guān)系

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

步驟4：利用邊緣檢測(cè)LOG算子對(duì)原始DEM進(jìn)行邊緣提取?？砂l(fā)現(xiàn)地形斷裂區(qū)域，地表起伏大的山坡、小的突起處都呈高亮狀態(tài)，但是地形平坦區(qū)域的亮度明顯低于其他區(qū)域，具有非常明顯的區(qū)分度?；谠撎攸c(diǎn)，通過亮度分割的方式提取地形平緩區(qū)域，對(duì)其進(jìn)行亮度反轉(zhuǎn)，然后利用形態(tài)學(xué)開運(yùn)算得到較完整的平緩區(qū)域。

步驟5：通過疊加求并運(yùn)算，利用步驟4提取的地形平坦區(qū)域進(jìn)一步消除步驟3提取的噪聲，提取結(jié)果即為斷裂線種子點(diǎn)。

步驟6：利用文獻(xiàn)[8]的方法，結(jié)合原始點(diǎn)云采用基于平面對(duì)的分段構(gòu)造法和基于三角網(wǎng)的直接構(gòu)造法提取完整的斷裂線矢量。

2 試驗(yàn)和精度評(píng)價(jià)

2.1 試驗(yàn)數(shù)據(jù)和試驗(yàn)工具

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

2.2 試驗(yàn)過程和結(jié)果分析

(1) 在LiDAR點(diǎn)云處理平臺(tái)LiDAR_Suite上實(shí)現(xiàn)點(diǎn)云精細(xì)分類，然后基于線性三角網(wǎng)插值算法獲取1 m分辨率的DEM，如圖5(b)所示。

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

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

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

(5) 基于步驟3和步驟4的提取結(jié)果進(jìn)行疊加分析，消除噪聲，得到完整的斷裂線種子點(diǎn)。在此基礎(chǔ)上，結(jié)合原始點(diǎn)云數(shù)據(jù)采用基于平面對(duì)的分段構(gòu)造法和基于三角網(wǎng)的直接構(gòu)造法提取斷裂線矢量。

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

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

3 總結(jié)和展望

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

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

參考文獻(xiàn)：

[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] 鄭輯濤，張濤.基于可變半徑圓環(huán)和B樣條擬合的機(jī)載LiDAR點(diǎn)云濾波[J].測(cè)繪學(xué)報(bào)，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] 易俐娜，劉鵬飛，喬小軍,等.結(jié)合遙感影像和DEM的線性體特征提取[J].測(cè)繪通報(bào)，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] 徐景中，萬(wàn)幼川，張圣望.基于機(jī)載激光雷達(dá)點(diǎn)云的斷裂線自動(dòng)提取方法[J].計(jì)算機(jī)應(yīng)用，2008，28(5):1214-1216.

[9] 彭檢貴，馬洪超，王宗躍,等.機(jī)載LiDAR點(diǎn)云的雙閾值自動(dòng)提取斷裂線方法[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，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，F(xiàn)ORNACIAI A.Visualization and Comparison of DEM-derived Parameters Application to Volcanic Areas[J].Geomorphology，2017，290: 69-84.