王杰棟 俞涵

（浙江省第二測(cè)繪院，浙江 杭州 310012）

當(dāng)前，測(cè)繪生產(chǎn)部門利用機(jī)載LiDAR 數(shù)據(jù)生產(chǎn)DEM 主要采用自動(dòng)濾波加后期人工編輯的方式[1]。基于單一LiDAR 點(diǎn)云數(shù)據(jù)進(jìn)行濾波的結(jié)果往往存在大量的錯(cuò)分和漏分點(diǎn)云現(xiàn)象，后期需要大量的人工編輯來(lái)提高DEM精度，大大增加了DEM的生產(chǎn)周期和生產(chǎn)成本。

考慮到許多數(shù)據(jù)管理部門已經(jīng)存有測(cè)區(qū)DEM 的實(shí)際情況，在獲得最新的點(diǎn)云數(shù)據(jù)后，可以利用這些數(shù)據(jù)直接檢測(cè)地形的變化區(qū)域（變化區(qū)域一般遠(yuǎn)小于整個(gè)測(cè)區(qū)范圍），實(shí)現(xiàn)DEM 的快速更新[2]。因此，本文提出一種基于變化檢測(cè)的LiDAR 點(diǎn)云DEM 更新方法。

1 算法原理

1.1 預(yù)處理

預(yù)處理主要包括LiDAR 點(diǎn)云數(shù)據(jù)與已有DEM 坐標(biāo)系統(tǒng)一和LiDAR 點(diǎn)云濾波。本文以已有DEM 所定義的坐標(biāo)系為參考坐標(biāo)系，將LiDAR 點(diǎn)云數(shù)據(jù)統(tǒng)一到參考坐標(biāo)系中。利用商用軟件Terrasolid 對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行初始濾波處理，獲得初始地面點(diǎn)及DEM。

1.2 候選變化區(qū)域獲取

針對(duì)預(yù)處理獲得的DEM 數(shù)據(jù)及已有的DEM 數(shù)據(jù)，利用公式（1）獲取高程差值模型（dDEM）。

公式（1）中，t1、t2表示前后兩個(gè)不同時(shí)相。通過(guò)對(duì)dDEM模型進(jìn)行判斷，可以定位測(cè)區(qū)內(nèi)的變化區(qū)域。初始的dDEM除了表示空間形態(tài)變化的大區(qū)域范圍外，還存在少量細(xì)小的圖斑。這些細(xì)小的圖斑是因?yàn)辄c(diǎn)云自動(dòng)濾波獲取的數(shù)據(jù)存在錯(cuò)分及漏分現(xiàn)象，導(dǎo)致生成的DEM 在部分區(qū)域存在偏差。為了有效去除干擾區(qū)域，確定真正發(fā)生變化的區(qū)域，本文采用高程差閾值處理和形態(tài)學(xué)濾波算法對(duì)dDEM模型進(jìn)行進(jìn)一步處理。

1.3 偽變化區(qū)域去除

1.3.1 高差閾值處理

通過(guò)設(shè)定兩期DEM 相對(duì)高差閾值，可以消除大部分偽變化區(qū)域。

公式（2）中，P 表示像素值，（c，r）代表像素點(diǎn)所在位置，T 為高程差閾值。

根據(jù)機(jī)載LiDAR 數(shù)據(jù)后處理技術(shù)規(guī)定，比例尺為1∶2000 的DEM 在平地和丘陵地的高程中誤差分別為0.4m 和0.5m。本文采用0.4m 作為高程差閾值。

1.3.2 多尺度形態(tài)學(xué)開(kāi)重建濾波處理

由于濾波結(jié)果存在少量突起或凹陷錯(cuò)分點(diǎn)、DEM插值誤差等問(wèn)題，進(jìn)行高差閾值去噪后的dDEM（0.4）中依舊有一部分偽變化區(qū)域。針對(duì)這些偽變化區(qū)域，本文通過(guò)將dDEM（0.4）轉(zhuǎn)化為二值圖像，采用多尺度形態(tài)學(xué)開(kāi)重建濾波算法進(jìn)行去除。

形態(tài)學(xué)開(kāi)重建方法是一種基于測(cè)地膨脹的形態(tài)學(xué)重建方法[3]。測(cè)地膨脹不需要直接選擇特定的結(jié)構(gòu)元素，而是基于標(biāo)記影像和掩膜影像進(jìn)行迭代計(jì)算。假定f 為標(biāo)記影像，p 為掩膜影像（Df＝Dp），B 表示結(jié)構(gòu)元素，則形態(tài)學(xué)膨脹可定義為：

基于上述形態(tài)學(xué)膨脹重建的定義，形態(tài)學(xué)開(kāi)重建可定義為：

公式（4）中O 表示形態(tài)學(xué)開(kāi)運(yùn)算。形態(tài)學(xué)開(kāi)重建是在開(kāi)運(yùn)算的基礎(chǔ)上增加了重建過(guò)程。

由于偽變化區(qū)域的形態(tài)大小各不相同，如果采用單一的結(jié)構(gòu)元素進(jìn)行濾波，往往不能有效去除不同尺度的偽變化區(qū)域噪聲。本文采用有限個(gè)不同大小的結(jié)構(gòu)元素參與影像濾波運(yùn)算，利用不同尺度結(jié)構(gòu)元素能去除對(duì)應(yīng)尺度噪聲的特點(diǎn)，構(gòu)建多尺度形態(tài)學(xué)濾波器。多尺度形態(tài)學(xué)濾波定義如下：

公式（5）中，nB 表示尺度為n 對(duì)應(yīng)的結(jié)構(gòu)元素。

1.4 區(qū)域?yàn)V波

利用已有的DEM數(shù)據(jù)輔助變化區(qū)域點(diǎn)云進(jìn)行濾波，主要采用優(yōu)化初始地面種子點(diǎn)的選取方式。其基本流程如下：

（1）對(duì)變化區(qū)域內(nèi)點(diǎn)云數(shù)據(jù)建立規(guī)則格網(wǎng)索引，格網(wǎng)尺度根據(jù)經(jīng)驗(yàn)設(shè)置為5m。

（2）設(shè)置高程差閾值 Hthreshold，判斷每個(gè)點(diǎn)與已有DEM的高程差值，將滿足閾值的點(diǎn)云記為待定地面點(diǎn)。

（3）遍歷格網(wǎng)內(nèi)的待定地面點(diǎn)，選擇高程值最小的待定地面點(diǎn)作為初始地面點(diǎn)。

（4）從已有DEM 中提取斷裂線區(qū)域，對(duì)落入斷裂線區(qū)域內(nèi)的格網(wǎng)進(jìn)行分裂處理，獲得小格網(wǎng)（小格網(wǎng)尺度為格網(wǎng)的一半），將小格網(wǎng)內(nèi)高程值最小的待定地面點(diǎn)加入到初始地面點(diǎn)集中。

（5）利用初始地面點(diǎn)數(shù)據(jù)，采用TIN 加密濾波算法獲得變化區(qū)域的地面點(diǎn)云及DEM。

在進(jìn)行區(qū)域?yàn)V波之后，還需對(duì)變化區(qū)域進(jìn)行少量人工交互編輯修正，從而獲得精確的DEM。

2 試驗(yàn)與分析

2.1 試驗(yàn)數(shù)據(jù)說(shuō)明

試驗(yàn)數(shù)據(jù)取自吉林省長(zhǎng)春市。前一期DEM 數(shù)據(jù)，如圖1（a）所示，由采集自2009 年5 月的LiDAR 點(diǎn)云數(shù)據(jù)生成，DEM 分辨率是1m；后一期點(diǎn)云數(shù)據(jù)，如圖1（b）所示，采集自2010 年10 月，數(shù)據(jù)范圍約為500×600 m2，平均點(diǎn)密度約為6.20 個(gè)/m2。2009 年，該地區(qū)很多地塊處在開(kāi)發(fā)初級(jí)階段，2010 年，隨著大量房地產(chǎn)項(xiàng)目的啟動(dòng)，地表變化劇烈，地面出現(xiàn)了大量土堆和較深的建筑物地基。圖1（c）是利用商業(yè)軟件Terrasolid 進(jìn)行自動(dòng)濾波處理后生成的概略DEM。

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

利用兩期DEM 數(shù)據(jù)相減，可以獲取該時(shí)間段的地形地貌變化情況，如圖2（a）所示；圖2（b）是對(duì)圖2（a）進(jìn)行0.4m 高程差閾值處理后的結(jié)果。

圖1 試驗(yàn)數(shù)據(jù)

圖2 dDEM

針對(duì)部分高程差值和真正變化區(qū)域類似的偽變化區(qū)域，采用多次度形態(tài)學(xué)開(kāi)重建進(jìn)行去除，如圖3 所示，其中圖3（a）是圖3（b）的二值化結(jié)果，圖3（b）是形態(tài)學(xué)運(yùn)算后的結(jié)果。

圖3 形態(tài)學(xué)運(yùn)算結(jié)果

選取圖3（b）中紅色區(qū)域作為變化區(qū)域，采用區(qū)域?yàn)V波算法獲得的變化區(qū)域DEM，如圖4（a），并將早期DEM 中的變化區(qū)域替換為新生成的區(qū)域DEM，從而實(shí)現(xiàn)DEM 更新，如圖4（b）。在地形變化區(qū)域采用基于DEM 輔助的濾波方法能有效提高濾波的準(zhǔn)確度，減少錯(cuò)分和漏分點(diǎn)云現(xiàn)象，降低后期的人工編輯量，且地形變化區(qū)域范圍遠(yuǎn)小于整個(gè)測(cè)區(qū)范圍，導(dǎo)致大部分測(cè)區(qū)可以使用原有的DEM 數(shù)據(jù)，提高了DEM 的生產(chǎn)效率，降低了生產(chǎn)成本。

圖4 DEM

3 結(jié)束語(yǔ)

本文提出了一種基于變化檢測(cè)的LiDAR 點(diǎn)云DEM更新方法，該方法充分利用已有的DEM，通過(guò)比較分析兩期DEM 數(shù)據(jù)提取地形變化區(qū)域，再對(duì)這些區(qū)域進(jìn)行濾波和人工編輯生成區(qū)域DEM，實(shí)現(xiàn)DEM的快速更新。試驗(yàn)表明，該方法在兩次量測(cè)期間測(cè)區(qū)地形沒(méi)有發(fā)生重大變化的情況下，利用點(diǎn)云數(shù)據(jù)和已有DEM 生成新的DEM，能有效提高生產(chǎn)及更新DEM 的速度， 降低生產(chǎn)成本。