朱依民，田林亞，畢繼鑫，林 松

(1.河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 211100； 2.浙江華東測繪與工程安全技術(shù)有限公司，浙江 杭州 310014)

我國水能資源空間分布不均衡，具有明顯的空間聚集現(xiàn)象。水利水電工程多位于西南地區(qū)的四川、云南、貴州、西藏等地，這些地方具有山高坡陡、河谷狹窄、植被茂密、氣候復(fù)雜多變、測區(qū)地類地物復(fù)雜多樣等特點(diǎn)，致使傳統(tǒng)對地觀測技術(shù)無法及時(shí)滿足水電工程的勘測設(shè)計(jì)需求[1]。隨著無人機(jī)技術(shù)[2-3]的發(fā)展，基于無人機(jī)傾斜攝影測量技術(shù)的地形圖測繪得到了很好的應(yīng)用[4]，但該技術(shù)僅適用于地表較為裸漏的測區(qū)，而對于植被茂密的地區(qū)，植被下方無法通過影像匹配得到點(diǎn)云數(shù)據(jù)。機(jī)載激光雷達(dá)(light detection and ranging，LiDAR)測量技術(shù)作為一種主動式對地觀測技術(shù)，集激光測距技術(shù)、差分GPS技術(shù)、慣性測量技術(shù)于一體，具有受外界環(huán)境影響小、自動化程度高、采集數(shù)據(jù)時(shí)間短、測量精度高等優(yōu)點(diǎn)。與無人機(jī)傾斜影像立體匹配技術(shù)相比，機(jī)載LiDAR測量具有較強(qiáng)的穿透力，能快速采集包含大量地物三維坐標(biāo)信息的點(diǎn)云數(shù)據(jù)，特別是其獨(dú)特的掃描式測量方法使其在獲取高山區(qū)水電工程的豐富地表信息方面有著得天獨(dú)厚的優(yōu)勢，為快速、精準(zhǔn)、自動生成水電工程地表模型提供了可能[5]。

機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波是重建水電工程地表模型的重要環(huán)節(jié)，但高山區(qū)水電工程區(qū)域植被茂密、山區(qū)高度差大、河流寬度窄，致使基于數(shù)學(xué)形態(tài)學(xué)的濾波算法[6]、基于坡度變化的濾波算法[7]、基于漸進(jìn)三角網(wǎng)加密(progressive triangulated irregular network densification, PTD)的濾波算法[8-9]、基于曲面擬合的濾波算法[10]以及基于分割的濾波算法[11]在高山區(qū)機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波中的應(yīng)用效果不佳，后期仍需大量人工干預(yù)三角網(wǎng)建網(wǎng)成果。李樂林等[12]采用金字塔策略按分辨率從高到低的方式對點(diǎn)云數(shù)據(jù)進(jìn)行格網(wǎng)分層，通過局部統(tǒng)計(jì)分析的方法確定高差閾值，進(jìn)而利用薄板樣條插值法從上至下對不同層格網(wǎng)控制點(diǎn)的高程進(jìn)行插值，最終得到最下層格網(wǎng)的地面點(diǎn)云集合，但該濾波方法的精度依賴于頂層格網(wǎng)種子點(diǎn)的選取以及通過統(tǒng)計(jì)分析得到的高差閾值的合理性。楊銘等[13]通過改變八叉樹模型節(jié)點(diǎn)的尺寸，將復(fù)雜地形的點(diǎn)云分解成大量的山坡地形點(diǎn)云，然后對每個(gè)分解部分利用坡度進(jìn)行地面點(diǎn)濾波和地面點(diǎn)集合的合并，該方法中距離和坡度閾值的設(shè)置以及八叉樹的分割標(biāo)準(zhǔn)將直接影響最終的濾波效果。孫蒙等[14]利用統(tǒng)計(jì)學(xué)中的偏度與峰度思想實(shí)現(xiàn)點(diǎn)云濾波，首先對原始點(diǎn)云進(jìn)行基于偏度與峰度變化曲線的初次濾波，然后對初次濾波獲取的地面點(diǎn)進(jìn)行曲面擬合，最后對曲面擬合高差進(jìn)行統(tǒng)計(jì)和判斷，以此完成點(diǎn)云二次濾波獲得地面點(diǎn)，該方法需要持續(xù)將高程最大點(diǎn)作為非地面點(diǎn)進(jìn)行剔除直至偏度小于等于零為止，因此該方法僅適用于平坦區(qū)域的點(diǎn)云濾波，并不適用于起伏較大的山地點(diǎn)云濾波。鄭緝濤等[15]提出一種結(jié)合可變半徑滾動圓和B樣條擬合的點(diǎn)云濾波方法，該方法利用一個(gè)可變半徑圓環(huán)從面向地心側(cè)的方向滾過由點(diǎn)云構(gòu)成的掃描線，保留掃描線上被圓環(huán)滾過的點(diǎn)并作為地面點(diǎn)，對這些地面點(diǎn)進(jìn)行均勻差值后利用B樣條擬合方法進(jìn)行曲面擬合，最后對曲面擬合高差進(jìn)行閾值判斷完成點(diǎn)云二次濾波獲得地面點(diǎn)集合，但該方法在濾波過程中進(jìn)行的多次插值和采樣會丟失原始點(diǎn)云數(shù)據(jù)信息。針對上述問題，本文提出一種基于PTD和改進(jìn)曲面擬合的機(jī)載LiDAR點(diǎn)云濾波算法，算法思想是先對點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理，再使用PTD算法進(jìn)行初次濾波獲取部分地面點(diǎn)集合，然后基于格網(wǎng)劃分對每一個(gè)格網(wǎng)進(jìn)行曲面擬合實(shí)現(xiàn)二次濾波，最后通過點(diǎn)云的精細(xì)化處理得到精確地面點(diǎn)的集合，研究成果為機(jī)載LiDAR測量技術(shù)在高山區(qū)水電工程中的應(yīng)用提供技術(shù)支撐。

1 基于PTD和改進(jìn)曲面擬合的高山區(qū)點(diǎn)云濾波方法

1.1 點(diǎn)云數(shù)據(jù)的預(yù)處理

1.1.1點(diǎn)云去噪

機(jī)載LiDAR掃描系統(tǒng)在掃描時(shí)會受諸如多路徑效應(yīng)、飛鳥等影響而產(chǎn)生低位誤差和高位誤差等噪聲點(diǎn)[16]，絕大多數(shù)濾波算法都是將點(diǎn)云的最低點(diǎn)作為地面種子點(diǎn)進(jìn)行處理，為了避免噪聲點(diǎn)對后期點(diǎn)云數(shù)據(jù)濾波產(chǎn)生影響，應(yīng)先對原始點(diǎn)云數(shù)據(jù)進(jìn)行去噪處理。常用的去噪方法有局部平面擬合法、頻率域法、高程分布直方圖法等，由于噪聲點(diǎn)一般較少且在高程空間分布上比較孤立，本文結(jié)合高程分布直方圖和KD樹對點(diǎn)云數(shù)據(jù)進(jìn)行去噪。

1.1.2點(diǎn)云回波次數(shù)判斷

大多數(shù)商業(yè)LiDAR系統(tǒng)都可以從點(diǎn)云數(shù)據(jù)中獲取4次回波信息[17]。機(jī)載LiDAR測量系統(tǒng)發(fā)射的激光脈沖在高山區(qū)水電工程傳播路徑上會遇到植被、建筑物、地面等不同地物，進(jìn)而產(chǎn)生不同的反射信息，且激光脈沖具有一定的穿透力，穿過植被的樹葉或枝干會產(chǎn)生多次回波，但遇到建筑物或地面時(shí)只有一次回波信息[18]。

通過分析不同測量目標(biāo)點(diǎn)云的回波信息，可以發(fā)現(xiàn)單次回波的點(diǎn)云一般是地面點(diǎn)、建筑物點(diǎn)以及植被的冠層點(diǎn)和枝干點(diǎn)等；多次回波的首次回波點(diǎn)云通常為高大植被的冠層點(diǎn)，中間次回波點(diǎn)云主要為高大植被的枝葉點(diǎn)以及低矮植被點(diǎn)，末次回波點(diǎn)云一般是植被的中間枝葉點(diǎn)和地表反射的激光點(diǎn)。綜上，可以利用回波次數(shù)對點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理，在具有單次回波、多次回波的末次回波的點(diǎn)云集合中獲取地面點(diǎn)，減少非地面點(diǎn)的數(shù)量，可提高后期濾波操作的可靠性和算法的效率。

1.2 PTD算法

如圖1所示，P1、P2、P3為地面種子點(diǎn)，P4為待判斷點(diǎn)，ρ為P4到P1、P2、P3形成的平面的距離，β1、β2、β3分別為P1、P2、P3和P4的連線與P1、P2、P3形成的平面的夾角。為了方便數(shù)據(jù)索引以及后期計(jì)算，要先對經(jīng)過預(yù)處理的點(diǎn)云數(shù)據(jù)進(jìn)行格網(wǎng)劃分，格網(wǎng)邊長為區(qū)域內(nèi)建筑物最大尺寸，這可保證每一格網(wǎng)都含有地面點(diǎn)；其次，將每個(gè)格網(wǎng)內(nèi)的高程最低點(diǎn)作為地面種子點(diǎn)，并以種子點(diǎn)構(gòu)建初始PTD模型；遍歷區(qū)域內(nèi)所有待判斷點(diǎn)，對其進(jìn)行距離和角度的閾值判斷，將符合條件的納入至初始PTD模型中，并更新PTD模型，對剩余點(diǎn)進(jìn)行迭代判斷和PTD模型更新，PTD算法的具體流程見文獻(xiàn)[19]。在本文算法中，由于漸進(jìn)三角網(wǎng)加密是為下一步的曲面擬合算法提供一定數(shù)量的地面種子點(diǎn)，并且考慮到PTD算法的耗時(shí)較長，所以在此進(jìn)行兩次迭代計(jì)算，并把得到的地面點(diǎn)集合作為下一步改進(jìn)的曲面擬合算法的種子點(diǎn)。

圖1 PTD算法的原理

1.3 改進(jìn)的曲面擬合算法

曲面擬合的思想是將任意復(fù)雜空間曲面表示為一個(gè)二次曲面函數(shù)，函數(shù)值代表擬合得到的對應(yīng)平面點(diǎn)的高程值，再將擬合得到的高程值與實(shí)際的高程測量值相減便可得到擬合高程差。對于地面點(diǎn)而言，擬合高程差較小，而非地面點(diǎn)的擬合高程差較大，因此可根據(jù)這一特性設(shè)置一個(gè)合理的擬合高程差閾值來分離地面點(diǎn)和非地面點(diǎn)。

令xn、yn、zn為第n個(gè)地面點(diǎn)云的三維坐標(biāo)值，設(shè)地表二次曲面擬合模型的方程為

Z=a0+a1X+a2Y+a3X2+a4XY+a5Y2

(1)

其中X=(x1,x2,…,xn)TY=(y1,y2,…,yn)T

式中：a0、a1、a2、a3、a4、a5為方程的待求參數(shù)，其誤差方程矩陣形式為

(2)

l=(z1,z2,…,zn)T

當(dāng)n≥6時(shí)，可解算最佳未知參數(shù)進(jìn)而擬合該空間曲面，由于是同精度的獨(dú)立觀測，所以將由不同測量值的權(quán)構(gòu)成的權(quán)陣P設(shè)為單位陣，根據(jù)VTPV取最小值的原則，通過下式求解參數(shù)值：

(3)

將解得的參數(shù)值代入式(1)便可得到擬合曲面方程，將點(diǎn)云的平面坐標(biāo)值帶入擬合方程中可計(jì)算擬合高程值，再通過比較擬合高程差與閾值的大小實(shí)現(xiàn)點(diǎn)云的二次濾波。

圖2 基于地面種子點(diǎn)的格網(wǎng)劃分

圖3 單個(gè)格網(wǎng)的曲面擬合

1.4 地面點(diǎn)的精細(xì)化處理

經(jīng)過前面的點(diǎn)云去噪、回波判斷、漸進(jìn)三角網(wǎng)加密初次濾波以及改進(jìn)曲面擬合算法的二次濾波，原始的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)被分離成地面點(diǎn)集合和非地面點(diǎn)集合。此時(shí)地面點(diǎn)集合還包含少量的低矮植被點(diǎn)，究其原因是低矮植被點(diǎn)與地面點(diǎn)緊緊相連，進(jìn)行距離和坡度閾值判斷時(shí)很難完全將其去除。為了進(jìn)一步提高濾波的精度，本文利用點(diǎn)的幾何特征對這部分點(diǎn)云進(jìn)行進(jìn)一步剔除，對于經(jīng)過上述步驟處理得到的地面點(diǎn)集，遍歷計(jì)算每一點(diǎn)及其K鄰域內(nèi)點(diǎn)的高度方差σ2，該值在植被處較高，而在地面及平面建筑物屋頂處取值較小，因此可以設(shè)定一個(gè)合適的高度方差閾值進(jìn)一步剔除地面點(diǎn)集合中夾雜的低矮植被點(diǎn)。

2 實(shí)例驗(yàn)證

2.1 數(shù)據(jù)源

為了驗(yàn)證提出的濾波算法在高山區(qū)水電工程機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波中的可行性，選取DJI-M600搭載的HS-600機(jī)載LiDAR測量系統(tǒng)(圖4)的實(shí)測數(shù)據(jù)進(jìn)行試驗(yàn)，測區(qū)面積約323 072 m2，共掃描點(diǎn)云474 339個(gè)，點(diǎn)云密度1.47個(gè)/m2，測區(qū)地形起伏大、植被茂密，且存在部分建筑物(圖5中黃色方框區(qū)域)和水域(圖5中紅色方框區(qū)域)等地物類型，完全滿足試驗(yàn)的需要，原始點(diǎn)云數(shù)據(jù)如圖5所示。

圖4 DJI-M600搭載HS-600機(jī)載LiDAR測量系統(tǒng)

圖5 高山區(qū)水電工程原始點(diǎn)云數(shù)據(jù)

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

首先對原始點(diǎn)云數(shù)據(jù)進(jìn)行去噪和回波次數(shù)判斷處理，進(jìn)而基于PTD算法對經(jīng)過去噪和回波次數(shù)判斷的點(diǎn)云進(jìn)行初次濾波。在PTD算法中需要對點(diǎn)云進(jìn)行虛擬格網(wǎng)劃分，且格網(wǎng)的大小要大于等于測區(qū)內(nèi)最大建筑物的邊長才能保證將建筑物全部濾除，由于本試驗(yàn)實(shí)際測得測區(qū)內(nèi)建筑物的最大邊長為24.5 m，所以試驗(yàn)中將格網(wǎng)邊長設(shè)置為25 m。根據(jù)測區(qū)的實(shí)際情況，在試驗(yàn)中通過設(shè)置不同的閾值來觀察PTD的效果，試驗(yàn)發(fā)現(xiàn)，當(dāng)設(shè)置的三角網(wǎng)加密距離閾值ρ0=0.5 m，角度判斷閾值β0=10°時(shí)，采用PTD算法既能提取足夠數(shù)量的地面點(diǎn)集合，又能得到一個(gè)分布比較均勻的地面種子點(diǎn)集合，這有利于為接下來曲面擬合的二次濾波提供有效的地面種子點(diǎn)。利用初始種子點(diǎn)生成的三角網(wǎng)如圖6所示，經(jīng)過兩次迭代加密得到的地面點(diǎn)如圖7所示。

圖6 利用初始種子點(diǎn)生成的三角網(wǎng)

圖7 基于兩次三角網(wǎng)加密處理的地面點(diǎn)

圖8 基于改進(jìn)的曲面擬合算法濾波后的地面點(diǎn)

經(jīng)過上述步驟，原始的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)被分為地面點(diǎn)集合和非地面點(diǎn)集合，但此時(shí)地面點(diǎn)集合仍包含少量的低矮植被點(diǎn)，設(shè)置0.5 m的高度方差閾值剔除低矮植被點(diǎn)，實(shí)現(xiàn)高山區(qū)水電工程地面點(diǎn)的精細(xì)化處理，精細(xì)化處理結(jié)果如圖9所示，得到的數(shù)字高程模型如圖10所示。

圖9 基于精細(xì)化處理的地面點(diǎn)云

圖10 基于精細(xì)化地面點(diǎn)云的數(shù)字高程模型

目前，定量評價(jià)濾波算法誤差的指標(biāo)主要有兩個(gè)：①國際攝影測量與遙感學(xué)會(ISPRS)提出的濾波評價(jià)準(zhǔn)則，即點(diǎn)云誤分率；②濾波前后點(diǎn)云的均值與標(biāo)準(zhǔn)差的差值。本文采用常用的點(diǎn)云誤分率對提出的算法濾波效果進(jìn)行定量分析，第Ⅰ類誤差E1指將地面點(diǎn)分類為非地面點(diǎn)的誤差，第Ⅱ類誤差E2指將非地面點(diǎn)分類為地面點(diǎn)的誤差。E1和E2反映了算法的適應(yīng)性，總誤差Eall反映了算法的可行性，E1、E2、Eall的計(jì)算公式為

E1=b/(a+b)

(4)

E2=c/(c+d)

(5)

Eall=(b+c)/(a+b+c+d)

(6)

式中：a為濾波算法中正確分類的地面點(diǎn)數(shù)量；d為正確分類的非地面點(diǎn)數(shù)量；b為地面點(diǎn)誤分類成非地面點(diǎn)的數(shù)量；c為非地面點(diǎn)誤分類成地面點(diǎn)的數(shù)量。

為了進(jìn)行橫向?qū)Ρ确治?，采用傳統(tǒng)的形態(tài)學(xué)濾波算法、區(qū)域生長濾波算法、PTD算法和曲面擬合濾波算法分別對同一山地點(diǎn)云數(shù)據(jù)進(jìn)行濾波并進(jìn)行精度評價(jià)，計(jì)算結(jié)果如表1所示。

表1 各濾波算法的精度計(jì)算結(jié)果 %

由表1可知，與PTD算法、曲面擬合濾波算法、區(qū)域生長濾波算法和形態(tài)學(xué)濾波算法相比，本文提出的基于結(jié)合PTD和改進(jìn)曲面擬合的濾波算法的E1、E2和Eall的最大減小幅度分別達(dá)6.25%、1.82%和2.93%，表明本文提出的濾波算法更適用于高山區(qū)水電工程機(jī)載LiDAR點(diǎn)云數(shù)據(jù)的濾波處理。

3 結(jié) 語

本文提出了一種綜合點(diǎn)云回波特性、PTD算法、改進(jìn)的曲面擬合算法和地面點(diǎn)精細(xì)化處理的高山區(qū)深植被點(diǎn)云濾波方法，采用DJI-M600搭載HS-600的機(jī)載LiDAR測量系統(tǒng)的實(shí)測數(shù)據(jù)進(jìn)行試驗(yàn)，利用ISPRS提出的點(diǎn)云誤分率對濾波的精度進(jìn)行評定。結(jié)果表明，與傳統(tǒng)的PTD算法、曲面擬合濾波算法、區(qū)域生長濾波算法和形態(tài)學(xué)濾波算法相比，基于PTD和改進(jìn)的曲面擬合的濾波算法的第Ⅰ類誤差、第Ⅱ類誤差和總誤差的最大減小幅度分別為6.25%、1.82%和2.93%，驗(yàn)證了基于PTD和改進(jìn)曲面擬合的濾波算法在高山區(qū)水電工程機(jī)載LiDAR點(diǎn)云濾波中的適用性。