机载激光雷达林区数字高程模型的建立1)

苏永涛 邢艳秋 焦义涛 邢万里 尤号田

(东北林业大学，哈尔滨 150040)

机载激光雷达林区数字高程模型的建立1)

苏永涛 邢艳秋 焦义涛 邢万里 尤号田

(东北林业大学，哈尔滨 150040)

选择地处草原和森林过渡地带的上库力农场作为研究区(E120°36′50.48″～120°52′56.53″，N50°21′11.08″～50°24′32″)，由机载激光雷达Leica ALS60采集实验数据，对TerraSolid分类获取的地形点建立数字高程模型(DEM)；利用IDL编译一次样条有限元内插法对点云数据进行分块处理，分析生产DEM的精度。结果表明：1.0、1.5、2.0 m三种不同分辨率的DEM精度，分别为0.034、0.078、0.096 m。

机载激光雷达；点云数据；数字高程模型；一次样条有限元法

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.

机载激光雷达(ALS)能够实时地获取三维信息，是一个以飞机为平台的主动遥感系统[1]。ALS采集的点云数据具有精度高、密度大等优势，因此从点云数据中获取地形点建立DEM，成为了一种快速获取高精度数字高程模型(DEM)产品的有效途径[2-4]。3D点云数据数据量较大，在处理、应用时会有许多不便。在一个DEM项目中，并不需要点云数据所表达出的全部信息，只需要获取相应的数据点，以达到一定的可信度和精度[5-6]。但是，在实际的点云数据采集过程中，实地地形并不能同采样密度很好的匹配，不可避免地出现过度采样问题，增加了数据的存储量和点云数据的处理时间[7]。所以，如何进行正确有效的DEM数据提取，并获取一种便于处理和操作点云数据的方法，对于建立DEM的地形数据集以及数学模型具有非常重要的实际意义。

DEM可分为不规则格网DEM和规则格网DEM[8]。不规则格网DEM的格网，可以是三角形、四边形或是其它的多边形格网；规则格网DEM的格网，一般是正方形。不规则格网DEM可很好的表示复杂地形，但数据量巨大、数据结构复杂，使用和管理较复杂[9]；由于规则格网DEM具有数据量小、结构简单、便于存储和管理的特点，成为了目前被广泛使用的一种DEM[10]。目前无论是为了商业目的，还是科学研究，虽然对应用激光雷达点云数据建立高精度DEM做了很多研究，但是，对林区地形的准确模拟仍然还需要面对许多挑战[11]。例如：如何准确识别并剔除林下地表的草本植物，以获取真实的地面点，建立真实反映林下地形的高精度DEM；如何针对复杂的林下地形，选择合适的内插方法建立DEM，以便后期的使用和管理等等。

本研究依据机载小光斑激光雷达点云数据，采用一次样条有限元内插法[12]，结合点云分块处理的原理建立DEM。使用TerraSolid软件对点云数据进行分类，通过迭代处理，尽可能准确地提取真实地面点；然后，对提取的地面点分块，采用一次样条有限元内插法，分别将分块后的地面点划分成互连的格网单元，进行内插，建立规则的格网DEM。生成的DEM，数据结构简单、精度高，而且占据的存储空间小。

1 研究区概况

研究区位于内蒙古呼伦贝尔市西北部、额尔古纳市东南部的上库力农场(E120°36′50.48″～120°52′56.53″，N50°21′11.08″～50°24′32″)，地处草原和森林过渡地带，属寒温带大陆季风气候，海拔600～700 m。研究区内地形起伏落差小，但较多为丘陵，地形地貌比较复杂。山脉、丘陵的阴坡，广泛分布着以白桦(Betulaplatyphylla)为主的天然次生林，混生树种包括落叶松(Larixgmelinii)和樟子松(Pinussylvestris)等；树下灌木层，主要由石棒绣线菊(Spiraeamedia)、筐柳(Salixlinearistipularis)等中层灌木类组成。

2 研究方法

2.1 数据获取

实验数据是由机载激光雷达Leica ALS60在2012年8月26日采集的。LeicaALS60发射的激光脉冲为近红外波段，波长为1 064 nm，最大脉冲频率为200 kHz，最大扫描频率为100 Hz；采用线性扫描方式，能够记录4次回波信息和3次强度信息，平均点云密度为8点/m2，光斑直径为0.22 m。数据存储格式采用标准las1.2格式。该数据格式可以存储点云的坐标值、高程值、回波强度、回波次数、自定义分类等信息。一同搭载的还有Leica RCD105相机，65 mm镜头，同步获取高清照片1 203幅，图像分辨率为20 cm，图像重叠率均达到50%以上。

2.2 点云去噪

在对研究区进行测量过程中，会产生各种测量误差和随机误差，例如，测量仪器的震荡、被测物表面粗糙不平、镜面反射、遮挡物遮挡等各种因素，这些因素均会使测点云数据中不可避免的存在不合理的噪声点。所以，在对点云数据进行处理之前，首先要进行去噪与平滑工作。

对于点云去噪的研究，国内外已经提出了许多行而有效的方法，这些处理方法主要与获取点云的排列形式密切相关。点云排列形式，一般分为有序或部分有序排列和无序排列[13]。对于有序或部分有序排列的点云数据当中的噪声点的处理，通常按照扫描线采用平滑滤波[14]的方法进行处理；而对于无序点云，由于点云自身点与点的拓扑关系尚未建立，所以，目前对于无序点云数据的去噪处理还没有一个简单、便捷的方法，在现阶段，去噪声点一般是自动与手工结合的方法。

本研究采用孤立点算法对点云去噪基本原理，首先选择一点作为球心，然后在其一定阈值范围内的球形空间中判断有没有其它点云存在，有则当前点云为非噪声点，否则为噪声点或体外孤点，最后人工剔除靠近地表的低噪声点。

2.3 点云分类

本研究使用TerraSolid软件对点云进行分类，采用不规则三角网(TIN)对点云进行迭代处理，分类出地形点和非地形点。首先，选择最低点作为初始的种子点，构建初始的三角网；然后，对点云进行筛选，将符合条件的点作为种子点，同原先的种子点重新构建新的三角网；再次，对点云进行筛选，直到分类出所有的地形点和非地形点，迭代结束。

2.4 选择拟合点和检查点

对于拟合点和检查点的选择，本研究采用密度选择法。密度选择法，是选择对拟合精度贡献大的点云为拟合点，是以拟合误差最小为目标函数。先将测区分成合适长度(δ)的若干单元格，δ为点云密度分辨率；然后，确定单元格内的点；最后，选择拟合点和检查点。首先，求得单元格中心点的坐标；然后，求单元格内每个点同中心点的距离，距离中心点最近的点选为拟合点，最远点作为检查点。单元格内可能有多个点，也可能一个点都没有，每个单元格我们只选一个点或两个点，一个作为拟合点，一个作为检查点。若只有一个点，其作为拟合点。若单元格内没有点，选择距单元格最近的邻近点作为拟合点。

2.5 一次样条有限元内插法

地形的表面可以看作是由许多个局部的独立单元组成，由相应的边界条件将这些独立的单元相互连接，构成一个连续的整体，即用这些相关的局部区域的解求得整个区域的解，这是有限元内插法[12]构建数字高程模型的基本内涵。用有限元DEM内插解算高程(见图1)，图1中dx、dy的计算方法见公式(1)。

图1 已知点P与DEM格网

(1)

如公式(2)所示，Δx、Δy为格网边长单位化的坐标增量。

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

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

(2)

如图1所示，点P的函数值Φ(x，y)可由其所在格网4个顶点的函数值Ci，j、Ci+1，j、Ci+1，j+1、Ci，j+1，按一次样条函数表示。

Φ(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)可以根据已知高程的数据点建立DEM。若P点为已知高程的数据点，则可用其高程hp作为观测值。以格网高程hi，j…作为待定的未知数，由式(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方向上的二次差分条件，构成第二类虚拟观测值误差方程式[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所示，首先由于格网的4个顶点，无论在x方向，还是y方向，都无法保证它的圆滑，因为，顶点在x、y方向都是单向的，所以不存在第二类误差方程式。其次，在格网的四周，x方向的上下两条边，只能保证在x方向保持圆滑的条件，在y方向，无法满足圆滑的条件；同样，在y方向也是如此，只能保证y方向的圆滑，而无法满足x方向。最后，是中间(M-2)×(N-2)个格网点全部满足第二类误差方程式所有条件，由此可以列出方程。尤其注意两类误差方程式的结构是不同的，其中第一类误差方程式的个数等于数据点的个数(s)，其结构与该数据点所在的DEM格网编号有关，即有多少个数据点，就有多少个第一类误差方程。第二类误差方程的个数等于m(n-2)+n(m-2)=2(mn-m-n)，它的结构与DEM节点序号有关。每个格网中一共含有m×n个DEM格网点(见图2)。现把高程格网点按照先列后行的顺序排列成一维的序列：1、2、3、…、n、n+1、n+2、…、2n、…、(m-1)n+1、…、mn。

图2 m行n列DEM格网点排列顺序

根据间接平差的数学模型V=Bx-l(B为系数矩阵，x为参数矩阵，l为常数矩阵)，对于第一类误差方程式，其间接平差的常数项为输入数据点的z坐标值；对于第二类误差方程式，其常数项为0。所以，可以将两类方程合并，一起列误差方程式。

2.6 DEM的精度评定

一次样条有限元内插法的精度评定，是从原始点云数据中选取一部分用做检查点，其余部分最接近小格网单元中心的点用来作为拟合点。通过拟合点拟合出一个曲面后，将检查点代入曲面中，计算检查点z值与拟合值的差值。由公式(6)计算数字高程模型中误差。

(6)

式中：Δh为检查点拟合值与检查点z值的差值；n为检查点的个数；Mh为拟合精度。

3 结果与分析

机载激光雷达在扫描过程中，不可避免会出现噪声点，高空中的噪声点对建立高精度的DEM影响较小，但是地面以下的噪声点，对建立高精度的DEM影响较大，所以必须尽可能地将其剔除，以减少其对DEM精度的影响。图3、图4分别为去噪前和去噪后的点云数据，图3中林分点云上方和下方位置离散的点云即为噪声点或体外孤点。经过观察发现，软件可以很好地剔除较明显的噪声点，但是还有部分地面以下的噪声点不能剔除。这部分噪声点的产生是由于激光在扫描过程中，部分到达地面的激光经过反射打到别的物体上，经过再次反射被激光接收器接收，经过计算其距离大于实际距离，解算出的高程便低于地面。这部分点如果不能够很好的剔除，在建立的DEM中会出现坑洼的区域，对DEM精度影响较大。对其的剔除，现在还没有很好的方法，只能通过人工进行删除，一般会比较费时费力。

图3 原始点云

图4 去噪后的点云

由于本研究目的是建立DEM，所以在进行分类过程中只将点云数据分成了两类，一类是地面点，另一类是非地面点。采用TIN方法，对原始点云数据分类，提出地形点(见图5)。

在实际生产中，并不是DEM产品的精度越高越好，DEM精度越高在生产的过程中耗时和耗财越多。因此，在选择DEM时，一般根据自己的需要，建立不同精度的DEM。本研究依据提取的地形点，采用一次样条有限元内插法，分别建立了3种分辨率的DEM产品，分别为1.0 m(见图6)、1.5 m(见图7)、2.0 m(见图8)。经对比观察发现，1.0 m分辨率的DEM更能表现实际地形，分辨率越小，其表示的实际地形越详细。在实验过程中，建立1.0 m分辨率DEM的时间，远长于建立1.5、2.0 m分辨率的DEM花费的时间，其精度也高于1.5、2.0 m分辨率的DEM(见表1)。

图5 提取的地形点

图6 1.0 m分辨率DEM

图7 1.5 m分辨率DEM

图8 2.0 m分辨率DEM

DEM分辨率/m拟合点数检查点数拟合精度/m1.0220416110.0341.5111710180.0782.06226100.096

一般对DEM进行精度评定时，需要使用GPS、全站仪或其他的测量仪器，获取部分地面控制点，以检验DEM的精度；但是，由于林区特殊的环境，无论是使用GPS，还是全站仪，都不能获取精度较高的控制点。因为，在林区由于林木冠层的遮挡，GPS接收的信号较弱，会出现单点解或无解；因此，测量误差较大，如果用这样测得的控制点用来检验建立的DEM，会造成较大的点位误差和高程中误差。如果使用全站仪或其它仪器测控制点，受林区地形和环境的影响，不但工作量大，而且费时费力，精度也不一定高。

本研究从提取的地面点中选取部分点云作为检查点，以检核建立的DEM精度，高程中误差如表1所示。从表1中可看出，不同分辨DEM的拟合点数和检查点数是不同的，这是因为在同一块区域，分辨率的大小代表了每个单元格的大小，分辨率越高，单元格越小，其数量越多，选择的拟合点和检查点越多，在计算时耗的时间越长。分辨率1.0、1.5、2.0 m的拟合精度，分别0.034、0.078、0.096 m，可见，在进行建立DEM时，分辨率越低，拟合精度越低，这主要因为每个小格网单元越大，其表示的地面越不精确，因此，在进行验证时产生的误差越大。

在建立DEM的过程中，DEM的精度往往受地形的影响。本研究区属于丘陵地区，起伏不大，但林下地形较复杂，既有平坦的地方，同时也有地势较陡的地方。激光雷达在扫描过程中，由于林木的遮挡不能测到一些地形特征点；或在选择拟合点的过程中，也会忽略一些地形的特征点，如山脊、山谷等。这常常会影响建立的DEM的精度，即生成较大的高程中误差，而且建立的DEM并不能较好地模拟林下实际地形状况。解决这个问题，一是需要提高激光雷达扫描点的密度，二是需要优化选择拟合点的方法，尽可能多的选择地形特征点作为拟合点，或是在点云分类过程中直接分出地形特征点作为拟合点，进行内插建立DEM。

以上3种分辨的的拟合精度，即高程中误差都在厘米级。因此，实际应用中，作为林区的DEM参与求解树高和冠层高度，精度完全是可以的。所以，机载激光雷达点云数据使用一次样条有限元内插法建立DEM，在林区应用中模拟林下地形以及求解树木结构参数是非常实用的。一次样条有限元内插法建立的DEM属于规则格网DEM，在地形较复杂地区不可避免的会忽略一些地形的结构及细部，这会造成模拟的地形不精确，需要附加些地形特征点，如山脊线、山谷线上的点参与内插计算建立DEM。但如何能够准确地提取出林下地形的特征点，并参与内插计算建立DEM还需要做较多的研究。

4 结论

结果表明：使用一次样条有限元内插法生产林区DEM可以达到较高精度，对于1.0 m分辨率的DEM，精度可达到0.034 m。因此，使用一次样条有限元内插法建立DEM，可以很好地模拟林下地形，进而提高估测的森林参数精度。在使用一次样条有限元内插法建立DEM的过程中，对点云数据分块处理建立DEM是非常重要的。由于受计算机内存的限制，如果分的区域太大，不仅建立DEM的速度会慢，还会因为内存不足而导致建立DEM失败。一次样条有限元内插法建立的是规则格网的DEM，数据结构简单、便于存储，但输出的DEM中没有保留原始的点，对于复杂的地形，尤其是落差比较大的地形，其精度并不能保持很好。所以，对于较复杂的地形，还需要考虑其它的拟合算法。

[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水准综合模型的探讨[J].西北大学学报(自然科学版)，2004，34(5)：611-614.

[13] 李元旺，黄文明，温佩芝，等.空间超限邻域点云去噪算法的研究与实现[J].计算机系统应用，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水准面拟合方法研究[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)黑龙江省博士后基金项目(LBH-Z15194)；林业公益性行业科研专项经费(201504319)；国家重点基础研究发展计划项目(2013CB733404)。

苏永涛，男，1973年3月生，东北林业大学文法学院，副教授；黑龙江省林业科学院博士后科研工作站，博士后科学研究。E-mail：syt200810@163.com。

邢艳秋，东北林业大学森林作业与环境研究中心，教授。E-mail：yanqiuxing@nefu.edu.cn。

2017年6月16日。

责任编辑：张 玉。

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