蘇永濤 邢艷秋 焦義濤 邢萬里 尤號田

(東北林業(yè)大學(xué)，哈爾濱 150040)

機(jī)載激光雷達(dá)林區(qū)數(shù)字高程模型的建立1)

蘇永濤 邢艷秋 焦義濤 邢萬里 尤號田

(東北林業(yè)大學(xué)，哈爾濱 150040)

選擇地處草原和森林過渡地帶的上庫力農(nóng)場作為研究區(qū)(E120°36′50.48″～120°52′56.53″，N50°21′11.08″～50°24′32″)，由機(jī)載激光雷達(dá)Leica ALS60采集實(shí)驗(yàn)數(shù)據(jù)，對TerraSolid分類獲取的地形點(diǎn)建立數(shù)字高程模型(DEM)；利用IDL編譯一次樣條有限元內(nèi)插法對點(diǎn)云數(shù)據(jù)進(jìn)行分塊處理，分析生產(chǎn)DEM的精度。結(jié)果表明：1.0、1.5、2.0 m三種不同分辨率的DEM精度，分別為0.034、0.078、0.096 m。

機(jī)載激光雷達(dá)；點(diǎn)云數(shù)據(jù)；數(shù)字高程模型；一次樣條有限元法

The grassland and forest transition zone library force the Shang Ku Li farm as study area (E120°50.48′36″-120°52′56.53″, N50°21′11.08″-50°24′32″) were selected. The point cloud data was acquired by the airborne laser scanning (ALS) technique of Leica ALS60. The ground points were classified by the TerraSolid software, and the digital elevation model (DEM) was created by a spline finite element interpolation method which was compiled by IDL. The accuracy of DEM were 0.034, 0.078, and 0.096 m when the resolution were 1.0, 1.5 and 2.0 m, respectively.

機(jī)載激光雷達(dá)(ALS)能夠?qū)崟r地獲取三維信息，是一個以飛機(jī)為平臺的主動遙感系統(tǒng)[1]。ALS采集的點(diǎn)云數(shù)據(jù)具有精度高、密度大等優(yōu)勢，因此從點(diǎn)云數(shù)據(jù)中獲取地形點(diǎn)建立DEM，成為了一種快速獲取高精度數(shù)字高程模型(DEM)產(chǎn)品的有效途徑[2-4]。3D點(diǎn)云數(shù)據(jù)數(shù)據(jù)量較大，在處理、應(yīng)用時會有許多不便。在一個DEM項(xiàng)目中，并不需要點(diǎn)云數(shù)據(jù)所表達(dá)出的全部信息，只需要獲取相應(yīng)的數(shù)據(jù)點(diǎn)，以達(dá)到一定的可信度和精度[5-6]。但是，在實(shí)際的點(diǎn)云數(shù)據(jù)采集過程中，實(shí)地地形并不能同采樣密度很好的匹配，不可避免地出現(xiàn)過度采樣問題，增加了數(shù)據(jù)的存儲量和點(diǎn)云數(shù)據(jù)的處理時間[7]。所以，如何進(jìn)行正確有效的DEM數(shù)據(jù)提取，并獲取一種便于處理和操作點(diǎn)云數(shù)據(jù)的方法，對于建立DEM的地形數(shù)據(jù)集以及數(shù)學(xué)模型具有非常重要的實(shí)際意義。

DEM可分為不規(guī)則格網(wǎng)DEM和規(guī)則格網(wǎng)DEM[8]。不規(guī)則格網(wǎng)DEM的格網(wǎng)，可以是三角形、四邊形或是其它的多邊形格網(wǎng)；規(guī)則格網(wǎng)DEM的格網(wǎng)，一般是正方形。不規(guī)則格網(wǎng)DEM可很好的表示復(fù)雜地形，但數(shù)據(jù)量巨大、數(shù)據(jù)結(jié)構(gòu)復(fù)雜，使用和管理較復(fù)雜[9]；由于規(guī)則格網(wǎng)DEM具有數(shù)據(jù)量小、結(jié)構(gòu)簡單、便于存儲和管理的特點(diǎn)，成為了目前被廣泛使用的一種DEM[10]。目前無論是為了商業(yè)目的，還是科學(xué)研究，雖然對應(yīng)用激光雷達(dá)點(diǎn)云數(shù)據(jù)建立高精度DEM做了很多研究，但是，對林區(qū)地形的準(zhǔn)確模擬仍然還需要面對許多挑戰(zhàn)[11]。例如：如何準(zhǔn)確識別并剔除林下地表的草本植物，以獲取真實(shí)的地面點(diǎn)，建立真實(shí)反映林下地形的高精度DEM；如何針對復(fù)雜的林下地形，選擇合適的內(nèi)插方法建立DEM，以便后期的使用和管理等等。

本研究依據(jù)機(jī)載小光斑激光雷達(dá)點(diǎn)云數(shù)據(jù)，采用一次樣條有限元內(nèi)插法[12]，結(jié)合點(diǎn)云分塊處理的原理建立DEM。使用TerraSolid軟件對點(diǎn)云數(shù)據(jù)進(jìn)行分類，通過迭代處理，盡可能準(zhǔn)確地提取真實(shí)地面點(diǎn)；然后，對提取的地面點(diǎn)分塊，采用一次樣條有限元內(nèi)插法，分別將分塊后的地面點(diǎn)劃分成互連的格網(wǎng)單元，進(jìn)行內(nèi)插，建立規(guī)則的格網(wǎng)DEM。生成的DEM，數(shù)據(jù)結(jié)構(gòu)簡單、精度高，而且占據(jù)的存儲空間小。

1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古呼倫貝爾市西北部、額爾古納市東南部的上庫力農(nóng)場(E120°36′50.48″～120°52′56.53″，N50°21′11.08″～50°24′32″)，地處草原和森林過渡地帶，屬寒溫帶大陸季風(fēng)氣候，海拔600～700 m。研究區(qū)內(nèi)地形起伏落差小，但較多為丘陵，地形地貌比較復(fù)雜。山脈、丘陵的陰坡，廣泛分布著以白樺(Betulaplatyphylla)為主的天然次生林，混生樹種包括落葉松(Larixgmelinii)和樟子松(Pinussylvestris)等；樹下灌木層，主要由石棒繡線菊(Spiraeamedia)、筐柳(Salixlinearistipularis)等中層灌木類組成。

2 研究方法

2.1 數(shù)據(jù)獲取

實(shí)驗(yàn)數(shù)據(jù)是由機(jī)載激光雷達(dá)Leica ALS60在2012年8月26日采集的。LeicaALS60發(fā)射的激光脈沖為近紅外波段，波長為1 064 nm，最大脈沖頻率為200 kHz，最大掃描頻率為100 Hz；采用線性掃描方式，能夠記錄4次回波信息和3次強(qiáng)度信息，平均點(diǎn)云密度為8點(diǎn)/m2，光斑直徑為0.22 m。數(shù)據(jù)存儲格式采用標(biāo)準(zhǔn)las1.2格式。該數(shù)據(jù)格式可以存儲點(diǎn)云的坐標(biāo)值、高程值、回波強(qiáng)度、回波次數(shù)、自定義分類等信息。一同搭載的還有Leica RCD105相機(jī)，65 mm鏡頭，同步獲取高清照片1 203幅，圖像分辨率為20 cm，圖像重疊率均達(dá)到50%以上。

2.2 點(diǎn)云去噪

在對研究區(qū)進(jìn)行測量過程中，會產(chǎn)生各種測量誤差和隨機(jī)誤差，例如，測量儀器的震蕩、被測物表面粗糙不平、鏡面反射、遮擋物遮擋等各種因素，這些因素均會使測點(diǎn)云數(shù)據(jù)中不可避免的存在不合理的噪聲點(diǎn)。所以，在對點(diǎn)云數(shù)據(jù)進(jìn)行處理之前，首先要進(jìn)行去噪與平滑工作。

對于點(diǎn)云去噪的研究，國內(nèi)外已經(jīng)提出了許多行而有效的方法，這些處理方法主要與獲取點(diǎn)云的排列形式密切相關(guān)。點(diǎn)云排列形式，一般分為有序或部分有序排列和無序排列[13]。對于有序或部分有序排列的點(diǎn)云數(shù)據(jù)當(dāng)中的噪聲點(diǎn)的處理，通常按照掃描線采用平滑濾波[14]的方法進(jìn)行處理；而對于無序點(diǎn)云，由于點(diǎn)云自身點(diǎn)與點(diǎn)的拓?fù)潢P(guān)系尚未建立，所以，目前對于無序點(diǎn)云數(shù)據(jù)的去噪處理還沒有一個簡單、便捷的方法，在現(xiàn)階段，去噪聲點(diǎn)一般是自動與手工結(jié)合的方法。

本研究采用孤立點(diǎn)算法對點(diǎn)云去噪基本原理，首先選擇一點(diǎn)作為球心，然后在其一定閾值范圍內(nèi)的球形空間中判斷有沒有其它點(diǎn)云存在，有則當(dāng)前點(diǎn)云為非噪聲點(diǎn)，否則為噪聲點(diǎn)或體外孤點(diǎn)，最后人工剔除靠近地表的低噪聲點(diǎn)。

2.3 點(diǎn)云分類

本研究使用TerraSolid軟件對點(diǎn)云進(jìn)行分類，采用不規(guī)則三角網(wǎng)(TIN)對點(diǎn)云進(jìn)行迭代處理，分類出地形點(diǎn)和非地形點(diǎn)。首先，選擇最低點(diǎn)作為初始的種子點(diǎn)，構(gòu)建初始的三角網(wǎng)；然后，對點(diǎn)云進(jìn)行篩選，將符合條件的點(diǎn)作為種子點(diǎn)，同原先的種子點(diǎn)重新構(gòu)建新的三角網(wǎng)；再次，對點(diǎn)云進(jìn)行篩選，直到分類出所有的地形點(diǎn)和非地形點(diǎn)，迭代結(jié)束。

2.4 選擇擬合點(diǎn)和檢查點(diǎn)

對于擬合點(diǎn)和檢查點(diǎn)的選擇，本研究采用密度選擇法。密度選擇法，是選擇對擬合精度貢獻(xiàn)大的點(diǎn)云為擬合點(diǎn)，是以擬合誤差最小為目標(biāo)函數(shù)。先將測區(qū)分成合適長度(δ)的若干單元格，δ為點(diǎn)云密度分辨率；然后，確定單元格內(nèi)的點(diǎn)；最后，選擇擬合點(diǎn)和檢查點(diǎn)。首先，求得單元格中心點(diǎn)的坐標(biāo)；然后，求單元格內(nèi)每個點(diǎn)同中心點(diǎn)的距離，距離中心點(diǎn)最近的點(diǎn)選為擬合點(diǎn)，最遠(yuǎn)點(diǎn)作為檢查點(diǎn)。單元格內(nèi)可能有多個點(diǎn)，也可能一個點(diǎn)都沒有，每個單元格我們只選一個點(diǎn)或兩個點(diǎn)，一個作為擬合點(diǎn)，一個作為檢查點(diǎn)。若只有一個點(diǎn)，其作為擬合點(diǎn)。若單元格內(nèi)沒有點(diǎn)，選擇距單元格最近的鄰近點(diǎn)作為擬合點(diǎn)。

2.5 一次樣條有限元內(nèi)插法

地形的表面可以看作是由許多個局部的獨(dú)立單元組成，由相應(yīng)的邊界條件將這些獨(dú)立的單元相互連接，構(gòu)成一個連續(xù)的整體，即用這些相關(guān)的局部區(qū)域的解求得整個區(qū)域的解，這是有限元內(nèi)插法[12]構(gòu)建數(shù)字高程模型的基本內(nèi)涵。用有限元DEM內(nèi)插解算高程(見圖1)，圖1中dx、dy的計(jì)算方法見公式(1)。

圖1 已知點(diǎn)P與DEM格網(wǎng)

(1)

如公式(2)所示，Δx、Δy為格網(wǎng)邊長單位化的坐標(biāo)增量。

Δx=(xp-xi)/dx，0<Δx<1；

Δy=(yp-yi)/dy，0<Δy<1。

(2)

如圖1所示，點(diǎn)P的函數(shù)值Φ(x，y)可由其所在格網(wǎng)4個頂點(diǎn)的函數(shù)值Ci，j、Ci+1，j、Ci+1，j+1、Ci，j+1，按一次樣條函數(shù)表示。

Φ(x，y)=(1-Δx)(1-Δy)Ci，j+Δx(1-Δy)Ci+1，j+

(1-Δx)ΔyCi，j+1+ΔxΔyCi+1，j+1。

(3)

使用公式(3)可以根據(jù)已知高程的數(shù)據(jù)點(diǎn)建立DEM。若P點(diǎn)為已知高程的數(shù)據(jù)點(diǎn)，則可用其高程hp作為觀測值。以格網(wǎng)高程hi，j…作為待定的未知數(shù)，由式(3)列出誤差方程，如式(4)。

vp=(1-Δx)(1-Δy)hi，j+Δx(1-Δy)hi+1，j+(1-Δx)Δyhi，j+1+

ΔxΔyhi+1，j+1-hp。

(4)

式(4)稱作一類觀測值誤差方程式。為了保證地面的圓滑，可利用x和y方向上的二次差分條件，構(gòu)成第二類虛擬觀測值誤差方程式[15]，見公式(5)。

vx(i，j)=hi-1，j-2hi，j+hi+1，j-0；

vy(i，j)=hi，j-1-2hi，j+hi+1，j+1-0。

(5)

第二類誤差方程式的列法，需要分析多種情況：如圖2所示，首先由于格網(wǎng)的4個頂點(diǎn)，無論在x方向，還是y方向，都無法保證它的圓滑，因?yàn)?，頂點(diǎn)在x、y方向都是單向的，所以不存在第二類誤差方程式。其次，在格網(wǎng)的四周，x方向的上下兩條邊，只能保證在x方向保持圓滑的條件，在y方向，無法滿足圓滑的條件；同樣，在y方向也是如此，只能保證y方向的圓滑，而無法滿足x方向。最后，是中間(M-2)×(N-2)個格網(wǎng)點(diǎn)全部滿足第二類誤差方程式所有條件，由此可以列出方程。尤其注意兩類誤差方程式的結(jié)構(gòu)是不同的，其中第一類誤差方程式的個數(shù)等于數(shù)據(jù)點(diǎn)的個數(shù)(s)，其結(jié)構(gòu)與該數(shù)據(jù)點(diǎn)所在的DEM格網(wǎng)編號有關(guān)，即有多少個數(shù)據(jù)點(diǎn)，就有多少個第一類誤差方程。第二類誤差方程的個數(shù)等于m(n-2)+n(m-2)=2(mn-m-n)，它的結(jié)構(gòu)與DEM節(jié)點(diǎn)序號有關(guān)。每個格網(wǎng)中一共含有m×n個DEM格網(wǎng)點(diǎn)(見圖2)?，F(xiàn)把高程格網(wǎng)點(diǎn)按照先列后行的順序排列成一維的序列：1、2、3、…、n、n+1、n+2、…、2n、…、(m-1)n+1、…、mn。

圖2 m行n列DEM格網(wǎng)點(diǎn)排列順序

根據(jù)間接平差的數(shù)學(xué)模型V=Bx-l(B為系數(shù)矩陣，x為參數(shù)矩陣，l為常數(shù)矩陣)，對于第一類誤差方程式，其間接平差的常數(shù)項(xiàng)為輸入數(shù)據(jù)點(diǎn)的z坐標(biāo)值；對于第二類誤差方程式，其常數(shù)項(xiàng)為0。所以，可以將兩類方程合并，一起列誤差方程式。

2.6 DEM的精度評定

一次樣條有限元內(nèi)插法的精度評定，是從原始點(diǎn)云數(shù)據(jù)中選取一部分用做檢查點(diǎn)，其余部分最接近小格網(wǎng)單元中心的點(diǎn)用來作為擬合點(diǎn)。通過擬合點(diǎn)擬合出一個曲面后，將檢查點(diǎn)代入曲面中，計(jì)算檢查點(diǎn)z值與擬合值的差值。由公式(6)計(jì)算數(shù)字高程模型中誤差。

(6)

式中：Δh為檢查點(diǎn)擬合值與檢查點(diǎn)z值的差值；n為檢查點(diǎn)的個數(shù)；Mh為擬合精度。

3 結(jié)果與分析

機(jī)載激光雷達(dá)在掃描過程中，不可避免會出現(xiàn)噪聲點(diǎn)，高空中的噪聲點(diǎn)對建立高精度的DEM影響較小，但是地面以下的噪聲點(diǎn)，對建立高精度的DEM影響較大，所以必須盡可能地將其剔除，以減少其對DEM精度的影響。圖3、圖4分別為去噪前和去噪后的點(diǎn)云數(shù)據(jù)，圖3中林分點(diǎn)云上方和下方位置離散的點(diǎn)云即為噪聲點(diǎn)或體外孤點(diǎn)。經(jīng)過觀察發(fā)現(xiàn)，軟件可以很好地剔除較明顯的噪聲點(diǎn)，但是還有部分地面以下的噪聲點(diǎn)不能剔除。這部分噪聲點(diǎn)的產(chǎn)生是由于激光在掃描過程中，部分到達(dá)地面的激光經(jīng)過反射打到別的物體上，經(jīng)過再次反射被激光接收器接收，經(jīng)過計(jì)算其距離大于實(shí)際距離，解算出的高程便低于地面。這部分點(diǎn)如果不能夠很好的剔除，在建立的DEM中會出現(xiàn)坑洼的區(qū)域，對DEM精度影響較大。對其的剔除，現(xiàn)在還沒有很好的方法，只能通過人工進(jìn)行刪除，一般會比較費(fèi)時費(fèi)力。

圖3 原始點(diǎn)云

圖4 去噪后的點(diǎn)云

由于本研究目的是建立DEM，所以在進(jìn)行分類過程中只將點(diǎn)云數(shù)據(jù)分成了兩類，一類是地面點(diǎn)，另一類是非地面點(diǎn)。采用TIN方法，對原始點(diǎn)云數(shù)據(jù)分類，提出地形點(diǎn)(見圖5)。

在實(shí)際生產(chǎn)中，并不是DEM產(chǎn)品的精度越高越好，DEM精度越高在生產(chǎn)的過程中耗時和耗財(cái)越多。因此，在選擇DEM時，一般根據(jù)自己的需要，建立不同精度的DEM。本研究依據(jù)提取的地形點(diǎn)，采用一次樣條有限元內(nèi)插法，分別建立了3種分辨率的DEM產(chǎn)品，分別為1.0 m(見圖6)、1.5 m(見圖7)、2.0 m(見圖8)。經(jīng)對比觀察發(fā)現(xiàn)，1.0 m分辨率的DEM更能表現(xiàn)實(shí)際地形，分辨率越小，其表示的實(shí)際地形越詳細(xì)。在實(shí)驗(yàn)過程中，建立1.0 m分辨率DEM的時間，遠(yuǎn)長于建立1.5、2.0 m分辨率的DEM花費(fèi)的時間，其精度也高于1.5、2.0 m分辨率的DEM(見表1)。

圖5 提取的地形點(diǎn)

圖6 1.0 m分辨率DEM

圖7 1.5 m分辨率DEM

圖8 2.0 m分辨率DEM

DEM分辨率/m擬合點(diǎn)數(shù)檢查點(diǎn)數(shù)擬合精度/m1.0220416110.0341.5111710180.0782.06226100.096

一般對DEM進(jìn)行精度評定時，需要使用GPS、全站儀或其他的測量儀器，獲取部分地面控制點(diǎn)，以檢驗(yàn)DEM的精度；但是，由于林區(qū)特殊的環(huán)境，無論是使用GPS，還是全站儀，都不能獲取精度較高的控制點(diǎn)。因?yàn)?，在林區(qū)由于林木冠層的遮擋，GPS接收的信號較弱，會出現(xiàn)單點(diǎn)解或無解；因此，測量誤差較大，如果用這樣測得的控制點(diǎn)用來檢驗(yàn)建立的DEM，會造成較大的點(diǎn)位誤差和高程中誤差。如果使用全站儀或其它儀器測控制點(diǎn)，受林區(qū)地形和環(huán)境的影響，不但工作量大，而且費(fèi)時費(fèi)力，精度也不一定高。

本研究從提取的地面點(diǎn)中選取部分點(diǎn)云作為檢查點(diǎn)，以檢核建立的DEM精度，高程中誤差如表1所示。從表1中可看出，不同分辨DEM的擬合點(diǎn)數(shù)和檢查點(diǎn)數(shù)是不同的，這是因?yàn)樵谕粔K區(qū)域，分辨率的大小代表了每個單元格的大小，分辨率越高，單元格越小，其數(shù)量越多，選擇的擬合點(diǎn)和檢查點(diǎn)越多，在計(jì)算時耗的時間越長。分辨率1.0、1.5、2.0 m的擬合精度，分別0.034、0.078、0.096 m，可見，在進(jìn)行建立DEM時，分辨率越低，擬合精度越低，這主要因?yàn)槊總€小格網(wǎng)單元越大，其表示的地面越不精確，因此，在進(jìn)行驗(yàn)證時產(chǎn)生的誤差越大。

在建立DEM的過程中，DEM的精度往往受地形的影響。本研究區(qū)屬于丘陵地區(qū)，起伏不大，但林下地形較復(fù)雜，既有平坦的地方，同時也有地勢較陡的地方。激光雷達(dá)在掃描過程中，由于林木的遮擋不能測到一些地形特征點(diǎn)；或在選擇擬合點(diǎn)的過程中，也會忽略一些地形的特征點(diǎn)，如山脊、山谷等。這常常會影響建立的DEM的精度，即生成較大的高程中誤差，而且建立的DEM并不能較好地模擬林下實(shí)際地形狀況。解決這個問題，一是需要提高激光雷達(dá)掃描點(diǎn)的密度，二是需要優(yōu)化選擇擬合點(diǎn)的方法，盡可能多的選擇地形特征點(diǎn)作為擬合點(diǎn)，或是在點(diǎn)云分類過程中直接分出地形特征點(diǎn)作為擬合點(diǎn)，進(jìn)行內(nèi)插建立DEM。

以上3種分辨的的擬合精度，即高程中誤差都在厘米級。因此，實(shí)際應(yīng)用中，作為林區(qū)的DEM參與求解樹高和冠層高度，精度完全是可以的。所以，機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)使用一次樣條有限元內(nèi)插法建立DEM，在林區(qū)應(yīng)用中模擬林下地形以及求解樹木結(jié)構(gòu)參數(shù)是非常實(shí)用的。一次樣條有限元內(nèi)插法建立的DEM屬于規(guī)則格網(wǎng)DEM，在地形較復(fù)雜地區(qū)不可避免的會忽略一些地形的結(jié)構(gòu)及細(xì)部，這會造成模擬的地形不精確，需要附加些地形特征點(diǎn)，如山脊線、山谷線上的點(diǎn)參與內(nèi)插計(jì)算建立DEM。但如何能夠準(zhǔn)確地提取出林下地形的特征點(diǎn)，并參與內(nèi)插計(jì)算建立DEM還需要做較多的研究。

4 結(jié)論

結(jié)果表明：使用一次樣條有限元內(nèi)插法生產(chǎn)林區(qū)DEM可以達(dá)到較高精度，對于1.0 m分辨率的DEM，精度可達(dá)到0.034 m。因此，使用一次樣條有限元內(nèi)插法建立DEM，可以很好地模擬林下地形，進(jìn)而提高估測的森林參數(shù)精度。在使用一次樣條有限元內(nèi)插法建立DEM的過程中，對點(diǎn)云數(shù)據(jù)分塊處理建立DEM是非常重要的。由于受計(jì)算機(jī)內(nèi)存的限制，如果分的區(qū)域太大，不僅建立DEM的速度會慢，還會因?yàn)閮?nèi)存不足而導(dǎo)致建立DEM失敗。一次樣條有限元內(nèi)插法建立的是規(guī)則格網(wǎng)的DEM，數(shù)據(jù)結(jié)構(gòu)簡單、便于存儲，但輸出的DEM中沒有保留原始的點(diǎn)，對于復(fù)雜的地形，尤其是落差比較大的地形，其精度并不能保持很好。所以，對于較復(fù)雜的地形，還需要考慮其它的擬合算法。

[2] SITHOLE G， VOSSELMAN G. Experimental comparison of filter algorithms for bare-Earth extraction from airborne laser scanning point clouds[J]. ISPRS Journal of Photogrammetry and Remote Sensing，2004，59(1/2)：85-101.

[3] GOBAKKEN T， BOLLANDS?S O M， NAESSET E. Comparing biophysical forest characteristics estimated from photogrammetric matching of aerial images and airborne laser scanning data[J]. Scandinavian Journal of Forest Research，2014，30(1)：73-86.

[5] LIU X Y. Airborne LiDAR for DEM generation： some critical issues[J]. Progress in Physical Geography，2008，32(1)：31-49.

[6] NAESSET E. Vertical height errors in digital terrain models derived from airborne laser scanner data in a boreal-alpine ecotone in Norway[J]. Remote Sensing，2015，7(4)：4702-4725.

[7] HANSEN E H， GOBAKKEN T， NAESSET E. Effects of pulse density on digital terrain models and canopy metrics using airborne laser scanning in a tropical rainforest[J]. Remote Sensing，2015，7(7)：8453-8468.

[8] PETZOLD B， REISS P， ST?SSEL W. Laser scanning-Surveying and mapping agencies are using a new technique for the derivation of digital terrain models[J]. ISPRS Journal of Photogrammetry and remote Sensing，1999，54(2/3)：95-104.

[9] WACK R， WIMMER A. Digital terrain models from airborne laserscanner data-a grid based approach[J]. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences，2002，34(3/B)：293-296.

[10] RAZAK K A， STRAATSMA M W， VAN WESTEN C J， et al. Airborne laser scanning of forested landslides characterization： Terrain model quality and visualization[J]. Geomorphology，2011，126(1)：186-200.

[11] CLARK M L， CLARK D B， ROBERTS D A. Small-footprint lidar estimation of sub-canopy elevation and tree height in a tropical rain forest landscape[J]. Remote Sensing of Environment，2004，91(1)：68-89.

[12] 劉萬林，李天文.建立GPS水準(zhǔn)綜合模型的探討[J].西北大學(xué)學(xué)報(自然科學(xué)版)，2004，34(5)：611-614.

[13] 李元旺，黃文明，溫佩芝，等.空間超限鄰域點(diǎn)云去噪算法的研究與實(shí)現(xiàn)[J].計(jì)算機(jī)系統(tǒng)應(yīng)用，2010，19(3)：35-38，52.

[14] HU H， DING Y L， ZHU Q， et al. An adaptive surface filter for airborne laser scanning point clouds by means of regularization and bending energy[J]. ISPRS Journal of Photogrammetry and Remote Sensing，2014，92：98-111.

[15] 姚吉利，褚麗麗，于志路.GPS水準(zhǔn)面擬合方法研究[J].測繪工程，2005，14(4)：23-26.

EstablishmentofDEMinForestAreawithAirborneLaserRadar

//Su Yongtao, Xing Yanqiu, Jiao Yitao, Xing Wanli, You Haotian

(Northeast Forestry University, Harbin 150040, P. R. China)

Airborne laser scanning; Point cloud data; Digital elevation model; A spline finite element method

P231；P237

1)黑龍江省博士后基金項(xiàng)目(LBH-Z15194)；林業(yè)公益性行業(yè)科研專項(xiàng)經(jīng)費(fèi)(201504319)；國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目(2013CB733404)。

蘇永濤，男，1973年3月生，東北林業(yè)大學(xué)文法學(xué)院，副教授；黑龍江省林業(yè)科學(xué)院博士后科研工作站，博士后科學(xué)研究。E-mail：syt200810@163.com。

邢艷秋，東北林業(yè)大學(xué)森林作業(yè)與環(huán)境研究中心，教授。E-mail：yanqiuxing@nefu.edu.cn。

2017年6月16日。

責(zé)任編輯：張 玉。

//Journal of Northeast Forestry University，2017,45(10)：30-34.