王效蓋 王 健 劉翔宇 曹 一

(山東科技大學(xué) 測繪與空間信息學(xué)院, 山東 青島 266590)

0 引言

機(jī)載激光雷達(dá)(light detection and ranging,LiDAR)技術(shù)是一種高效、快速、大區(qū)域采集三維空間點(diǎn)云的主動(dòng)遙感技術(shù),已廣泛應(yīng)用于數(shù)字模型構(gòu)建、實(shí)景三維建設(shè)、森林資源管理、地質(zhì)災(zāi)害監(jiān)測等領(lǐng)域[1]。其中,機(jī)載點(diǎn)云濾波是分離地面點(diǎn)和非地面點(diǎn)的關(guān)鍵步驟,也是LiDAR點(diǎn)云數(shù)據(jù)后續(xù)實(shí)際應(yīng)用的基礎(chǔ)。目前地面濾波算法主要分為坡度濾波、形態(tài)學(xué)濾波、分割濾波和插值濾波。

坡度濾波是利用點(diǎn)云之間坡度變化分離地面點(diǎn)的一種濾波方法,這類方法十分依賴坡度閾值的設(shè)定,其對連續(xù)地形有較好的分割結(jié)果,但在陡坡、樹木與地面連接處等坡度突變區(qū)域容易出現(xiàn)錯(cuò)誤濾波現(xiàn)象[2-3]。形態(tài)學(xué)濾波主要原理是LiDAR點(diǎn)云在形態(tài)學(xué)開閉運(yùn)算前后存在高程差異,高程變化小的為地面點(diǎn),反之為非地面點(diǎn),其精度依賴窗口大小的設(shè)置,需要使用者設(shè)定合理的窗口參數(shù),對非專業(yè)人員不友好[4]。分割濾波首先根據(jù)LiDAR點(diǎn)云某種屬性對其進(jìn)行分割操作,并根據(jù)分割屬性之間的差異分離地面點(diǎn)和非地面點(diǎn),在特殊地形中有較好的分割結(jié)果,但其精度受分割結(jié)果影響較大[5-6]。插值濾波首先需要選擇地面種子點(diǎn),并使用插值方法生成初始參考面,若待分類點(diǎn)與參考面之間滿足一定的關(guān)系則將其劃歸為地面點(diǎn),基于插值的濾波方法原理簡單、濾波結(jié)果可靠、魯棒性強(qiáng),但是其結(jié)果容易受種子點(diǎn)構(gòu)造初始參考面精度的影響[7-8]。Axelsson提出的漸進(jìn)加密三角網(wǎng)(triangulated irregular network,TIN)濾波方法(progressive TIN densification,PTD)是插值濾波中的經(jīng)典算法,因其較高的濾波精度和較快的濾波速度而被集成于多個(gè)開源或商業(yè)軟件中,如TerraSolid等[9]。實(shí)驗(yàn)發(fā)現(xiàn),PTD算法在對平坦或地形簡單的區(qū)域有較好的濾波效果,但在復(fù)雜林區(qū)中因獲取初始地面種子點(diǎn)數(shù)量變少且可靠性降低,導(dǎo)致在進(jìn)行迭代運(yùn)算時(shí)誤差累積增大,使其濾波效果不太理想。

為此,在PTD算法的基礎(chǔ)上,本文提出了一種適用于林區(qū)LiDAR點(diǎn)云改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波方法(improved progressive TIN densification,IPTD)。IPTD算法在地面種子點(diǎn)選取策略上做出一定改進(jìn),借助布料模擬算法和局部薄板樣條插值獲得大量分布均勻且可靠的地面種子點(diǎn),用以提高初始TIN的精度,然后對待分類點(diǎn)進(jìn)行迭代濾波,從而提高算法在林區(qū)環(huán)境中LiDAR點(diǎn)云的濾波精度。

1 漸進(jìn)加密三角網(wǎng)濾波

漸進(jìn)加密三角網(wǎng)濾波算法是Axelsson在2000年提出的一種經(jīng)典插值濾波方法并且得到廣泛的使用[9]。算法流程中選取局部窗口的最低點(diǎn)作為初始地面種子點(diǎn),以此構(gòu)建初始三角網(wǎng),通過計(jì)算待定點(diǎn)P到ABC平面的距離d及連線之間的夾角α、β、γ,判斷距離和角度是否滿足設(shè)定閾值,若滿足則將P點(diǎn)加入地面點(diǎn)并構(gòu)建新的三角網(wǎng),反復(fù)迭代至沒有新點(diǎn)加入,濾波完成。漸進(jìn)加密三角網(wǎng)濾波算法原理如圖1所示。

圖1 漸進(jìn)加密三角網(wǎng)濾波原理

該算法具體步驟為：

(1)選擇地面種子點(diǎn)。對研究區(qū)域進(jìn)行規(guī)則格網(wǎng)劃分,格網(wǎng)大小設(shè)置為區(qū)域最大建筑物大小,選取格網(wǎng)中最低點(diǎn)作為初始地面種子點(diǎn)。

(2)構(gòu)建初始三角網(wǎng)。根據(jù)(1)中選擇的地面種子點(diǎn)構(gòu)建一個(gè)稀疏的初始三角網(wǎng)。

(3)待定點(diǎn)分類。依次遍歷所有待定點(diǎn),分別計(jì)算待定點(diǎn)到三角網(wǎng)平面的距離和夾角,若距離和角度均滿足設(shè)定閾值,將該點(diǎn)劃分為地面點(diǎn)。

(4)更新三角網(wǎng)。使用新加入的地面點(diǎn)對三角網(wǎng)進(jìn)行更新,生成新的三角網(wǎng)。

(5)重復(fù)步驟(3)和(4),迭代至沒有新增的地面點(diǎn),算法結(jié)束。

該算法原理簡單,效果穩(wěn)定,對于城市和林區(qū)都有較強(qiáng)的適用性。其算法精度主要取決于初始三角網(wǎng)構(gòu)建的準(zhǔn)確性和迭代閾值的設(shè)定。在植被密集覆蓋且地形崎嶇復(fù)雜的林區(qū),采用局部最小值獲取的種子點(diǎn)可靠性降低,從而影響初始三角網(wǎng)的準(zhǔn)確性,最終使得濾波精度較低濾波效果不理想。

2 改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波方法

2.1 噪聲剔除

原始機(jī)載LiDAR數(shù)據(jù)中噪聲來源于掃描儀本身或受環(huán)境影響產(chǎn)生的多路徑效應(yīng)。這對使用布料模擬算法獲取潛在地面種子點(diǎn)時(shí)產(chǎn)生嚴(yán)重影響,導(dǎo)致獲取的潛在地面種子點(diǎn)結(jié)果不可靠。因此,在進(jìn)行后續(xù)處理前需要剔除原始點(diǎn)云中遠(yuǎn)離主體點(diǎn)云的離群噪聲點(diǎn)。基于空間位置分布的去噪算法(statistical outlier removal,SOR)能夠穩(wěn)健的檢測這些離群點(diǎn)。SOR濾波器對每個(gè)點(diǎn)進(jìn)行K鄰域統(tǒng)計(jì)分析,計(jì)算K個(gè)鄰近點(diǎn)平均距離分布的均值和標(biāo)準(zhǔn)差,將平均距離在給定閾值范圍之外的點(diǎn)視為離群點(diǎn)并剔除。如式(1)所示,本文使用12個(gè)最近鄰計(jì)算每個(gè)點(diǎn)的均值和中值距離。

Max(d)>Mean(d)+Std(d)

(1)

2.2 獲取潛在地面種子點(diǎn)

布料模擬濾波算法是由Zhang等提出的一種新型濾波方法,其基本思想：首先將機(jī)載LiDAR點(diǎn)云進(jìn)行倒置并在測量點(diǎn)上方覆蓋一塊柔性布料,柔性布料在重力的作用下逐漸黏附在測量點(diǎn)上,最終柔性布料模擬出一個(gè)近似的地形表面,如果測量點(diǎn)到柔性布料的距離滿足閾值要求則標(biāo)記為地面點(diǎn)[10]。布料模擬的優(yōu)點(diǎn)在于：①參數(shù)設(shè)置簡單,僅需較少的參數(shù)即可完成濾波;②地形適應(yīng)性好,可用于連續(xù)地形、陡坡和斷崖;③運(yùn)算效率高。PTD算法中只有少數(shù)點(diǎn)被選為地面點(diǎn),導(dǎo)致初始三角網(wǎng)精度較低。相比之下,布料模擬算法可以獲取大量分布均勻的潛在地面種子點(diǎn)。布料模擬算法過程如圖2所示。

圖2 布料模擬算法過程

2.3 獲取地面種子點(diǎn)

在應(yīng)用布料模擬算法濾波后,某些非地面點(diǎn)會(huì)錯(cuò)誤的包含在潛在地面點(diǎn)中,為了提高原始三角網(wǎng)構(gòu)建的精度,需要識別并剔除這些非地面點(diǎn),本文使用薄板樣條插值方法來檢驗(yàn)并消除非地面點(diǎn)。薄板樣條插值是一種空間離散點(diǎn)擬合精度高、魯棒性強(qiáng)的空間插值方法,在對復(fù)雜地形的模擬中表現(xiàn)出較好效果[11]。

對于3維空間中的測量點(diǎn)集(X,Y,Z),薄板樣條插值曲面可由式(2)表示。

(2)

(3)

其中,N是用于TPS插值的控制點(diǎn)數(shù);p(x,y)是趨勢函數(shù),a是它的系數(shù);wi是與控制點(diǎn)相關(guān)的權(quán)重;q(ri)是徑向基函數(shù)。徑向基函數(shù)定義為

q(r)=r2ln(r)

(4)

r為兩點(diǎn)間的歐氏距離。矩陣表達(dá)式如下：

(5)

實(shí)驗(yàn)發(fā)現(xiàn),非地面點(diǎn)存在局部高程突變的現(xiàn)象,本文認(rèn)為在鄰域點(diǎn)云中高程值異常的點(diǎn)為非地面點(diǎn),并構(gòu)造插值曲面從潛在種子點(diǎn)中篩選并去除非地面點(diǎn)。因?yàn)榍蠼獗“鍢訔l插值函數(shù)系數(shù)的時(shí)間和插值點(diǎn)的數(shù)量成正比,所以本文采用局部薄板樣條插值代替整體薄板樣條插值來提高效率。具體步驟為：①使用規(guī)則格網(wǎng)(20×20)組織潛在地面點(diǎn)數(shù)據(jù),并選擇局部最小值作為插值點(diǎn);②搜索待估計(jì)潛在種子點(diǎn)K鄰域的插值點(diǎn),構(gòu)造局部薄板樣條插值曲面,K設(shè)置為12,并將激光腳點(diǎn)坐標(biāo)代入求解插值高程;③計(jì)算真實(shí)高程與插值高程的殘差,若高程殘差大于設(shè)定的殘差閾值,則認(rèn)為該點(diǎn)為非地面點(diǎn)并從潛在種子點(diǎn)中剔除;④逐點(diǎn)篩選得到地面種子點(diǎn)。

2.4 改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波

傳統(tǒng)漸進(jìn)加密三角網(wǎng)濾波算法選擇種子點(diǎn)采用格網(wǎng)區(qū)域高程最低點(diǎn),在林區(qū)環(huán)境下種子點(diǎn)選取較少且可靠性低,在進(jìn)行迭代運(yùn)算時(shí)造成誤差累積。因此,本文在種子點(diǎn)選取策略上做出改進(jìn),使用布料模擬算法和局部薄板樣條插值獲得大量分布均勻和可靠的地面種子點(diǎn),解決傳統(tǒng)漸進(jìn)加密三角網(wǎng)濾波在林區(qū)環(huán)境下種子點(diǎn)稀疏的問題,提高林區(qū)下的適應(yīng)性。改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波步驟如下：

(1)種子點(diǎn)選取。首先使用布料模擬算法提取機(jī)載LiDAR點(diǎn)云中的潛在地面點(diǎn),然后通過局部薄板樣條插值剔除潛在地面點(diǎn)中的非地面點(diǎn),獲得分布均勻和可靠的地面種子點(diǎn)。

(2)初始三角網(wǎng)構(gòu)建。在機(jī)載LiDAR點(diǎn)云四周設(shè)置寬度為m的緩沖區(qū),并以固定間隔n構(gòu)造輔助點(diǎn),輔助點(diǎn)高程為最近種子點(diǎn)的高程[12]。通過種子點(diǎn)和輔助點(diǎn)構(gòu)建初始三角網(wǎng)。

(3)三角網(wǎng)濾波。計(jì)算待定點(diǎn)到鄰近三角形的距離和角度,若其距離和角度小于設(shè)定的距離閾值和角度閾值,則分為地面點(diǎn)。

(4)更新三角網(wǎng)。將步驟(3)地面點(diǎn)與已有種子點(diǎn)重新構(gòu)建三角網(wǎng)。

(5)迭代運(yùn)算。重復(fù)步驟(3)和(4),迭代運(yùn)算至沒有新的地面點(diǎn)加入。

3 實(shí)驗(yàn)及結(jié)果分析

采用國際攝影測量與遙感學(xué)會(huì)(International Society for Photogrammetry and Remote Sensing，ISPRS)標(biāo)準(zhǔn)數(shù)據(jù)和林區(qū)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)并對濾波結(jié)果進(jìn)行精度分析。本文采用的濾波精度評價(jià)指標(biāo)為Ⅰ類誤差(T.Ⅰ)、Ⅱ類誤差(T.Ⅱ)、總誤差(T.E)和Kappa系數(shù),其中三種誤差數(shù)值越小和Kappa系數(shù)數(shù)值越大表示濾波精度越高。計(jì)算方法如表1所示,其中a和b分別為正確分類和錯(cuò)誤分類的地面點(diǎn)數(shù);d和c分別是正確和錯(cuò)誤分類的非地面點(diǎn)數(shù);e是a、b、c和d的總和。

表1 地面濾波精度評價(jià)方法

3.1 ISPRS標(biāo)準(zhǔn)數(shù)據(jù)

采用ISPRS提供的標(biāo)準(zhǔn)數(shù)據(jù)集中山地部分(samp51-samp71)作為測試數(shù)據(jù),驗(yàn)證本文所提方法的有效性。每個(gè)樣本的參考數(shù)據(jù)都是人工濾波生成的,樣本點(diǎn)被分類為地面點(diǎn)和非地面點(diǎn)。IPTD對ISPRS標(biāo)準(zhǔn)數(shù)據(jù)集的濾波結(jié)果如表2所示。

表2 ISPRS標(biāo)準(zhǔn)數(shù)據(jù)濾波精度

從表2可以看出：IPTD算法濾波結(jié)果Ⅰ類誤差為1.23%,Ⅱ類誤差為12.29%,樣本平均Kappa系數(shù)和總誤差分別為84.96%和2.16%,說明本文算法有效控制了Ⅰ類誤差的產(chǎn)生。samp51、samp61和samp71獲得的Ⅰ類誤差均小于1%,其中samp51的Ⅰ類誤差最小,僅為0.09%,說明本文方法在緩坡濾波中能較較好地保留真實(shí)地面點(diǎn),減小其誤分類的可能性。Samp53Ⅱ類誤差最大,高達(dá)34.06%,其余樣本數(shù)據(jù)都控制在一定范圍內(nèi),本文算法在斷崖處的濾波還有待優(yōu)化。表中samp51的Kappa系數(shù)最高,高達(dá)96.04%,samp61的總誤差最小,僅為0.85%。samp52和samp53的總誤差都較大,samp53的Kappa系數(shù)最小,僅為56.06%。圖3和圖4展示了兩個(gè)樣本的數(shù)字地表模型(digital surface model，DSM)、濾波后的數(shù)字高程模型(digital elevation model，DEM)及濾波誤差的空間分布情況,其中藍(lán)色為Ⅰ類誤差,紅色為Ⅱ類誤差。由圖3可見,samp52的數(shù)據(jù)中存在斷崖和陡坡,沿河分布低矮的植被,這些高差突變區(qū)域?qū)﹂撝涤幸欢ㄌ魬?zhàn),判斷錯(cuò)誤而造成Ⅱ類誤差較大。由圖4可見,samp53的數(shù)據(jù)中分布有梯田等落差較大區(qū)域,其中誤差分布較多的地方地形變化較大,導(dǎo)致Ⅱ類誤差較大。

圖3 Samp52濾波誤差空間分布

圖4 Samp53濾波誤差空間分布

將本文方法與近年來提出的5種方法進(jìn)行對比[8,13-16],如表3所示。本文方法在6種方法中總誤差平均值最小,較原有濾波方法總誤差的平均值減小0.75%,且本文方法在所有樣本數(shù)據(jù)中有3個(gè)總誤差最小,說明本文方法在濾波精度上具有一定的優(yōu)勢。

表3 本文方法與其他5種方法的總誤差對比 單位：%

3.2 林區(qū)數(shù)據(jù)

ISPRS標(biāo)準(zhǔn)數(shù)據(jù)集點(diǎn)云密度較小且植被特征不夠明顯,為了進(jìn)一步驗(yàn)證本文方法在森林地區(qū)的適用性,本研究從OpenTopography網(wǎng)站(http://www.opentopography.org/)下載3組高密度林區(qū)機(jī)載LiDAR數(shù)據(jù)。3組林區(qū)數(shù)據(jù)具有不同地地形特征：數(shù)據(jù)Site1位于美國加利福尼亞州內(nèi)華達(dá)山脈南部,為坡度較小和植被覆蓋密集的緩坡地形;Site2位于美國猶他州中南部山區(qū),為覆蓋分散低矮植被的陡坡地形;Site3位于美國科羅拉多州,地形起伏明顯,覆蓋較多密集的植被。3組林區(qū)數(shù)據(jù)點(diǎn)云如圖5所示。

圖5 林區(qū)點(diǎn)云數(shù)據(jù)

采用本文方法和PTD算法對3組樣本數(shù)據(jù)進(jìn)行濾波處理并進(jìn)行對比分析,其結(jié)果如表4所示。根據(jù)表中數(shù)據(jù)可以看出,本文方法在Ⅰ類誤差和Ⅱ類誤差上都有大幅的降低,說明本文方法在林區(qū)濾波中能夠大幅度減少錯(cuò)誤分類現(xiàn)象,獲得較準(zhǔn)確的地面模型。從總誤差來看,本文方法相較于PTD算法在Site1區(qū)域降低了3.51%,在Site2區(qū)域降低8.68%,在Site3區(qū)域降低4.02%,說明本文方法在總誤差控制上有較好的表現(xiàn)。本文方法Kappa系數(shù)最大為94.57,最小為80.14,在Kappa系數(shù)有了大幅的提高,其中Site2區(qū)域提高僅14.50%。由圖6、圖7和圖8可以看出,本文方法生成的DEM與參考DEM接近。綜上所述,相較于傳統(tǒng)的PTD算法,本文方法可以應(yīng)用于林區(qū)點(diǎn)云濾波中。

表4 林區(qū)數(shù)據(jù)濾波精度對比

圖6 Site1參考數(shù)據(jù)、PTD和本文方法生成的DEM

圖7 Site2生成的DEM

圖8 Site3生成的DEM

4 結(jié)束語

為了提高漸進(jìn)加密三角網(wǎng)濾波算法在林區(qū)機(jī)載點(diǎn)云數(shù)據(jù)中的濾波精度,提出了一種改進(jìn)的漸進(jìn)加密三角網(wǎng)濾波方法。所提方法使用布料模擬濾波和局部薄板樣條插值方法代替局部最小值法獲取大量均勻可靠的地面種子點(diǎn),解決了因種子點(diǎn)稀疏導(dǎo)致初始三角網(wǎng)精度較低的問題。使用ISPRS標(biāo)準(zhǔn)數(shù)據(jù)和林區(qū)數(shù)據(jù)對方法進(jìn)行驗(yàn)證,并與其他濾波方法進(jìn)行對比分析,本文方法有較低的總誤差和較高的Kappa系數(shù),證明本文方法能有效應(yīng)用于林區(qū)點(diǎn)云濾波中。本文方法在陡坡等復(fù)雜地形中濾波誤差仍然較大,未來就提高不同地形下林區(qū)點(diǎn)云的濾波精度繼續(xù)展開工作。