一種局部稀疏地面點(diǎn)云與已有DEM的融合方法

雷麗珍，林 超

(廣東省國土資源技術(shù)中心，廣東 廣州 510075)

利用機(jī)載LiDAR點(diǎn)云數(shù)據(jù)快速完成高精度數(shù)字高程模型(DEM)的大規(guī)模生產(chǎn)，已成為近幾年自然資源管理、城鎮(zhèn)規(guī)劃與建設(shè)等諸多領(lǐng)域的一個(gè)研究和應(yīng)用熱點(diǎn)[1-5]。機(jī)載LiDAR點(diǎn)云數(shù)據(jù)質(zhì)量好壞對后期DEM生產(chǎn)具有重要影響。由于設(shè)備特點(diǎn)、測區(qū)地形地貌、植被覆蓋、成本和天氣等諸多因素影響[6]，實(shí)際獲取的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)存在地面點(diǎn)局部過于稀疏的問題，導(dǎo)致采用三角網(wǎng)內(nèi)插算法構(gòu)建DEM時(shí)，出現(xiàn)較嚴(yán)重的“三角面片化”問題，在一定程度上降低了成果精度[7]，影響DEM的制圖與使用，因此需要對這些LiDAR點(diǎn)數(shù)據(jù)進(jìn)行數(shù)據(jù)質(zhì)量改善。而進(jìn)行局部區(qū)域的重飛或補(bǔ)飛，要消耗較大的成本、時(shí)間，且后期獲取的數(shù)據(jù)與前期獲取的數(shù)據(jù)之間可能存在一定的誤差，數(shù)據(jù)一致性與整合也會(huì)存在一定的問題[8]。針對以上難題，本文提出一種基于高斯核函數(shù)加權(quán)迭代插值算法的局部稀疏地面點(diǎn)云與已有DEM融合方法，以實(shí)現(xiàn)LiDAR點(diǎn)云生成DEM時(shí)局部地面點(diǎn)稀疏區(qū)域的地形特征保持及高程精度的改善。

1 研究方法

1.1 基本思路

LiDAR地面點(diǎn)云與已有DEM之間存在一定的高程差異，且高程差異為非剛性的偏差，不同位置的差異值不一致。本文提出一種內(nèi)插算法，在LiDAR點(diǎn)云稀疏區(qū)域，以濾波后少量正確的地面點(diǎn)作為控制點(diǎn)，利用高斯核函數(shù)加權(quán)迭代插值方法改正低精度的DEM數(shù)據(jù)，對填補(bǔ)在LiDAR地面點(diǎn)較稀疏區(qū)域的DEM點(diǎn)云的高程值進(jìn)行局部改善，使其高程數(shù)值與其鄰近的LiDAR地面點(diǎn)的高程值保持一致，同時(shí)盡量保持原始DEM的地形形態(tài)特征[9]。

1.2 融合插值原理

地表形態(tài)通常具有連續(xù)性，在一定的空間鄰域內(nèi)具有一定的相似性。地理學(xué)第一定律所描述的“越鄰近的事物越相關(guān)”，為空間數(shù)據(jù)插值等方法應(yīng)用提供了基礎(chǔ)的理論依據(jù)。空間插值方法是將不同粒度、不同尺度的空間數(shù)據(jù)進(jìn)行有效融合的一種常用方法，它使用一組采樣點(diǎn)的觀測值的加權(quán)組合來確定待觀測點(diǎn)的值。本文采用加權(quán)插值方法實(shí)現(xiàn)稀疏地面點(diǎn)云與已有DEM的數(shù)據(jù)融合，將稀疏地面點(diǎn)云作為高精度的控制點(diǎn)，通過插值原理估計(jì)控制點(diǎn)一定鄰域范圍內(nèi)DEM高程的局部改正值，同時(shí)在對DEM高程進(jìn)行糾正時(shí)盡量保持DEM的整體地形形態(tài)特征。其原理如圖1所示。

如圖1所示，估計(jì)p點(diǎn)處的高程改正值時(shí)，首先搜索其一定距離范圍內(nèi)的LiDAR地面點(diǎn)；然后，依據(jù)這些地面點(diǎn)云的真實(shí)高程與其估計(jì)高程值間的高程差異(即dH)，通過高斯核函數(shù)加權(quán)來估計(jì)p點(diǎn)位置的高程改正。計(jì)算公式如下

(1)

式中，p表示待估計(jì)高程改正值的DEM點(diǎn)；q表示LiDAR地面點(diǎn)云G中的一個(gè)點(diǎn)；d(q,p)表示q與p之間的歐氏距離。上述公式值中，權(quán)值wq的確定是影響高程改正的重要因素。本文采用待估點(diǎn)與控制點(diǎn)間距離的高斯核函數(shù)作為控制點(diǎn)的加權(quán)權(quán)重，即

(2)

式中，h為高斯核函數(shù)的帶寬。

由于稀疏地面點(diǎn)云數(shù)據(jù)的空間分布密度不均，采用全局高斯核函數(shù)帶寬參數(shù)很難獲得比較理想的插值結(jié)果，為此，本文提出基于Delaunay三角網(wǎng)的高斯核函數(shù)帶寬h的自適應(yīng)確定方法[10]。首先采用稀疏地面點(diǎn)云構(gòu)建Delaunay三角網(wǎng)，然后統(tǒng)計(jì)與每個(gè)地面點(diǎn)(如pi)相連的三角網(wǎng)邊的平均長度(記為mi)和Delaunay三角網(wǎng)中所有邊長的標(biāo)準(zhǔn)方差(記為s)，則對于點(diǎn)pi的高斯核加權(quán)核函數(shù)的帶寬h計(jì)算為

h(pi)=mi+3s

(3)

進(jìn)而，按照以下步驟迭代對已有DEM的高程進(jìn)行局部改正：

(1) 根據(jù)已有的DEM數(shù)據(jù)，構(gòu)建一定距離間隔(如2 m×2 m)的規(guī)則格網(wǎng)，通過雙線性插值算法獲得每個(gè)格網(wǎng)角點(diǎn)的高程值，將這些格網(wǎng)的角點(diǎn)保存為帶高程的DEM點(diǎn)并輸出。

(2) 針對每個(gè)DEM點(diǎn)，通過最鄰近查詢，搜索DEM點(diǎn)附近的LiDAR地面點(diǎn)，計(jì)算這些DEM點(diǎn)的平均高程作為LiDAR點(diǎn)的估計(jì)高程值，記為H1，LiDAR點(diǎn)原始高程值記為H0，計(jì)算兩者高程差異

dH=H0-H1

(4)

(3) 統(tǒng)計(jì)所有LiDAR點(diǎn)的高程差異dH的平均值、最大差異值及大部分點(diǎn)的差異值。

(4) 利用高斯核函數(shù)加權(quán)插值方法，根據(jù)LiDAR點(diǎn)云平面位置、dH數(shù)值和DEM點(diǎn)的平面位置計(jì)算出DEM點(diǎn)位置處的高程改正項(xiàng)correctH，將原始DEM高程數(shù)值加上改正后的高差作為改正后DEM格網(wǎng)點(diǎn)的高程值[11]。

(5) 迭代執(zhí)行以上步驟，直到DEM點(diǎn)的高程值收斂(或改變很微小)時(shí)停止。

1.3 技術(shù)流程

本文方法在LAS點(diǎn)云格式下進(jìn)行運(yùn)算，技術(shù)流程如圖2所示。選取需處理的稀疏地面點(diǎn)云區(qū)域的數(shù)據(jù)，與轉(zhuǎn)換為LAS格式的DEM點(diǎn)進(jìn)行融合處理生成插值點(diǎn)，與LiDAR點(diǎn)云地面點(diǎn)疊加內(nèi)插生成DEM數(shù)據(jù)。

1.4 算法實(shí)現(xiàn)

基于以上原理開發(fā)設(shè)計(jì)了點(diǎn)云數(shù)據(jù)融合工具軟件，界面設(shè)計(jì)如圖3所示。

點(diǎn)云數(shù)據(jù)融合工具主要包括以下功能：① 數(shù)據(jù)輸入：選擇需要進(jìn)行處理的LiDAR點(diǎn)云數(shù)據(jù)文件和DEM點(diǎn)云數(shù)據(jù)文件；② 參數(shù)設(shè)置：設(shè)置點(diǎn)搜索的鄰域半徑，LiDAR點(diǎn)按類別讀取及高程改正的平滑參數(shù)等；③ 結(jié)果輸出與保存；④ 運(yùn)行狀態(tài)提示，反饋中間計(jì)算的運(yùn)行情況。

2 試驗(yàn)結(jié)果

試驗(yàn)采樣了某測區(qū)內(nèi)3.9 km×1.5 km范圍內(nèi)的點(diǎn)云數(shù)據(jù)作為原始待融合的LiDAR點(diǎn)云數(shù)據(jù)。該區(qū)域LiDAR數(shù)據(jù)中一些區(qū)域地面點(diǎn)較密，也有一些區(qū)域地面點(diǎn)非常稀少，這些區(qū)域并非平地，過于稀少的LiDAR地面點(diǎn)會(huì)影響這些區(qū)域的DEM構(gòu)建。為此，采樣同一區(qū)域已有的較低精度DEM點(diǎn)云數(shù)據(jù)進(jìn)行LiDAR地面點(diǎn)的填補(bǔ)與改善。

測試結(jié)果的局部區(qū)域顯示結(jié)果如圖4所示。

從測試結(jié)果可以發(fā)現(xiàn)，該算法可以有效地對原始稀疏LiDAR地面點(diǎn)進(jìn)行填補(bǔ)和改善，在對稀疏LiDAR點(diǎn)進(jìn)行插值的同時(shí)，也盡可能地保持了原始地形走勢和總體的地貌形態(tài)特征。根據(jù)插值后的點(diǎn)云生產(chǎn)的DEM如圖5所示。

通過比較融合前后地面點(diǎn)生產(chǎn)DEM可以發(fā)現(xiàn)，基于原始地面點(diǎn)生產(chǎn)的DEM在地面點(diǎn)稀少的局部區(qū)域存在較嚴(yán)重的“面片化”現(xiàn)象，地形過渡不自然、不流暢；而經(jīng)過插值后的點(diǎn)云在這些區(qū)域的地形更自然，貼近自然的地貌形態(tài)特征。

3 精度及可靠性驗(yàn)證

為了綜合驗(yàn)證本文方法插值的精度及可靠性，選取了某區(qū)域作為驗(yàn)證區(qū)(如圖6(b)所示)。該區(qū)域范圍為700 m×700 m，數(shù)據(jù)獲取時(shí)設(shè)備穿透性高，分類后地面點(diǎn)分布密集，點(diǎn)云數(shù)量為4 223 559，統(tǒng)計(jì)密度為5.2點(diǎn)/平方米，地面點(diǎn)數(shù)量為215 440，該測區(qū)實(shí)際檢驗(yàn)的點(diǎn)云高程中誤差為0.2 m，通過人工剔除120 133個(gè)(占55.8%)地面點(diǎn)模擬實(shí)際中可能出現(xiàn)的稀疏點(diǎn)云狀態(tài)(如圖6(c)、圖6(d)所示)，并與該區(qū)域已有低精度DEM(如圖6(a)所示)數(shù)據(jù)進(jìn)行本文融合算法處理，并生成2 m格網(wǎng)DEM成果(如圖6(e)、圖6(f)所示)。

由圖6可以看出，已有低精度DEM與基于點(diǎn)云生成的高精度DEM在山體地形走勢上基本一致，但高精度DEM在細(xì)節(jié)上更為豐富和準(zhǔn)確，模擬稀疏點(diǎn)云剔除部分關(guān)鍵地面點(diǎn)后局部山脊及山谷明顯失真，通過高斯核函數(shù)加權(quán)內(nèi)插后，失真情況在視覺上有明顯改善。從剔除的地面點(diǎn)中均勻選取1881個(gè)模型關(guān)鍵點(diǎn)作為本次驗(yàn)證的檢查點(diǎn)，對以上數(shù)據(jù)進(jìn)行誤差統(tǒng)計(jì)。檢查點(diǎn)如圖6(g)所示，主要分布在模擬稀疏點(diǎn)云區(qū)域，且數(shù)量較多，對試驗(yàn)中各個(gè)DEM數(shù)據(jù)進(jìn)行誤差計(jì)算分析。統(tǒng)計(jì)后中誤差情況見表1，對各區(qū)間殘差分布的百分比對比情況分析如圖7所示。

表1 中誤差統(tǒng)計(jì) m

利用高斯核函數(shù)加權(quán)內(nèi)插融合后的DEM的中誤差為1.796 m，在精度上相對于融合前的稀疏地面點(diǎn)情況有一定改善，且誤差呈正態(tài)分布，與稀疏地面點(diǎn)DEM誤差分布曲線相比，內(nèi)插后粗差有明顯減小，DEM在弱精度區(qū)域的可靠性顯著提升，鑒于對點(diǎn)云細(xì)微地物完全刪除后難以通過內(nèi)插恢復(fù)，誤差在±2 m區(qū)間內(nèi)的變化差異不大。

4 結(jié) 語

本文針對機(jī)載LiDAR點(diǎn)云數(shù)據(jù)中地面點(diǎn)稀少問題，提出了基于已有低精度DEM數(shù)據(jù)與LiDAR稀疏地面點(diǎn)進(jìn)行融合的技術(shù)方法，并通過試驗(yàn)對融合效果、精度及可靠性進(jìn)行了分析驗(yàn)證。結(jié)果表明，使用融合后的點(diǎn)云數(shù)據(jù)生產(chǎn)的DEM產(chǎn)品，可以減少稀疏地面點(diǎn)云生成的DEM中存在的“三角面片化”問題，并提高了DEM的精度[12-13]，尤其是在存在較大誤差的區(qū)域，對實(shí)際生產(chǎn)中出現(xiàn)的地面點(diǎn)稀疏區(qū)域進(jìn)行局部融合，能夠改善DEM的整體質(zhì)量。