基于机载全波形LiDAR数据的森林地上生物量估测算法研究

邢艳秋，姚松涛，李梦颖，谢杰，闫灿

(东北林业大学 森林作业与环境研究中心，哈尔滨 150040)

基于机载全波形LiDAR数据的森林地上生物量估测算法研究

邢艳秋，姚松涛，李梦颖，谢杰，闫灿

(东北林业大学 森林作业与环境研究中心，哈尔滨 150040)

为提高森林地上生物量(AGB)的估测精度，本研究以白桦林为研究对象，以机载全波形激光雷达(LiDAR)数据为研究数据，首先提出了机载全波形LiDAR数据读取与全波形特征信息提取的相关算法，然后结合具体算法的实现提取出每条全波形数据对应的各波形分量的能量信息，进而依据波形能量信息计算出每个样地的全波形激光穿透指数(FiLPI)，之后通过全波形激光穿透指数FiLPI建立其与对应样地实测森林AGB的统计回归模型，同时将森林AGB的估测值与森林AGB的实测值进行对比。结果表明全波形激光穿透指数FiLPI与森林AGB具有很好的相关性(R2为0.885，RMSE为0.095)，并且森林AGB的估测值与实测值之间误差的波动较小，提高了森林AGB的估测精度，弥补和提供了机载全波形LiDAR数据估测森林AGB的方法和思路。

机载激光雷达；全波形数据；波形数据读取；波形特征信息提取；白桦林；地上生物量

0 引言

在林业调查中，森林地上生物量(Aboveground Biomass，AGB)通常指的是单位面积内植被的绝干物质量[1-2]。传统森林AGB的测定方法耗时耗力耗材且更新较慢。激光雷达(Light Detection And Ranging，LiDAR)作为一种主动式遥感技术，广泛应用于森林经营管理与生态系统研究等领域[3]。按回波模式，LiDAR系统通常可分为离散回波和全波形两种，其中，离散回波LiDAR通过记录冠层的单个或多个离散回波信号以获得不同密度的点云数据；全波形LiDAR则通过记录返回脉冲的全部能量来得到亚米级森林垂直剖面。LiDAR遥感技术估测森林AGB是非破坏性的，有利于生态环境的保护。

大多研究学者基于离散回波LiDAR点云数据对森林AGB进行了大量的估测研究，Bortolot等[4]提出了一种适用于不同林型和不同区域的森林AGB估测的可视化算法，该算法以LiDAR点云数据提取的林木参数和外业实测数据为训练样本同AGB建立回归关系，进而估测其它区域的森林AGB，结果表明森林AGB的外业实测值同估测值的相关性在0.59～0.82之间，能够成功估测不同区域的森林AGB。刘清旺等[5]研究了离散回波LiDAR点云数据估测单株树AGB，由点云数据生成的冠层高度模型结合优化的单株树冠特征识别算法估测相关的结构参数，结果表明，由LIDAR点云数据得到的单株树估测参数与实测参数显著相关(R2=0.729)，可以估测单株树AGB。王涛等[6]研究表明，森林AGB同叶面积指数等参数相关，然而有关机载全波形LiDAR数据估测森林AGB的研究相对较少。

为研究机载全波形LiDAR数据在森林AGB估测研究中的应用潜能，本研究提出了全波形数据读取和波形特征信息提取算法，进而利用波形能量信息参数进行森林AGB估测研究，弥补机载全波形LiDAR数据估测森林AGB的方法和思路。

1 研究区与数据获取

1.1 研究区概况

本研究选择内蒙古自治区额尔古纳市东南部的依根农林交错区作为研究区域，地理位置坐标为120°36′50.48″E～120°52′56.53″E，50°21′11.08″N～ 50°24′32.00″N。植物资源以天然次生白桦林(Betulaplatyphylla)为主，混生树种包括落叶松(Larixgmelinii)、樟子松(Pinussylvestrisvar.mongolica)等，林下灌木层主要由石棒绣线菊(Spiraeamedia)、筐柳(Salixlinearistipularis)等组成。

1.2 外业样地数据

本研究于2012年8月至9月对研究区天然次生白桦林进行了相应的地面调查。综合考虑林地的森林类型及生长状况，一共布设了80个边长为10 m的方形样地，测量了树高、胸径、冠幅等单木参数。采用Vertex IV超声波测高测距仪测量了每木树高，30 m范围内其高度分辨率为0.1 m；采用胸径尺测量每木胸径(1.30 m树高处周长)，并利用皮尺测量东西和南北方向的单木冠幅。根据研究需求，借助Trimble GeoXT6000手持GPS定位仪对样地中心点进行定位，定位精度达50 cm。外业样地数据统计结果见表1。

表1 外业样地数据基本统计量

1.3 机载全波形LiDAR数据

本研究机载全波形LiDAR数据采集时间为2012年8月26日，天气状况良好，采用的飞机型号为运-5，相对飞行高度1300 m，搭载的LiDAR扫描系统为Leica ALS60，采用的大地坐标系统为WGS84，投影坐标系统为UTM。Leica ALS60系统最高扫描频率为200 kHz，最大视场角为75°。每个激光脉冲在地面上所形成的光斑直径大约为0.2～0.3 m，全波形数据是返回脉冲被接收望远镜接收后由采样器对其进行数字化采样而成，记录了激光光斑内目标物对应的完整波形数据[7-8]。

2 数据处理与研究方法

2.1 样地森林地上生物量计算

王正文等[9]对位于内蒙古自治区呼伦贝尔市额尔古纳右旗境内白桦林生物量进行了评估，并确定了白桦林AGB相对生长模型如公式(1)所示，模型相关系数达到0.985。由于本研究区与前者地理位置相邻且树种相同，因此采用公式(1)计算样地内单木AGB。

(1)

则样地AGB平均值为公式(2)所示：

(2)

式中：Wi表示第i株单木生物量，kg；Di表示第i株单木胸径，cm；Hi表示第i株单木树高，m；Wm表示样地平均地上生物量，kg；n表示样地内单木株数。

2.2 全波形数据读取

全波形文件的读取是分析处理波形数据的首要前提。本研究采集的全波形数据经过坐标转换等处理后，以LAS1.3格式存储在文件中，LAS1.3增加了对全波形数据的支持，以二进制形式存储，记录了接收到的完整波形数据，包括采样间隔和实时记录的方位等。通过访问一个激光点数据，能够以映射字段访问到与该点关联的返回脉冲波形，读取到对应的波形信息。计算机系统访问内存速度比访问硬盘速度快得多，在读取数据时如果少访问硬盘而更多地访问内存，那么读取速度将会快速提升。基于这种思路，本研究使用IDL编程语言借助内存映射方法[10]实现了研究区LeicaALS60LAS1.3全波形数据的快速读取。

2.3 分析处理与波形特征信息提取

本研究将原始全波形数据近似地看作不同目标物后向散射高斯波形的叠加，是噪声与一系列高斯分量之和在时间轴上的叠加，即原始全波形数据是等时间间隔采样的序列{tm，f(tm)}，其中{tm，m=1，…，n}表示采样时间，{f(tm)，m=1，…，n}表示采样时间tm对应的振幅值，符合混合高斯分布。使用公式(3)来表示由一系列简单高斯分量组成的原始波形数据。

(3)

2.3.1 噪声去除

每条完整全波形数据对应着19.2 m的沿脉冲距离，在实际工作环境中若不考虑异常噪声点，则可穿透目标物的分布范围通常要小于这个长度[11]。据此可选某个四分位数振幅值作为去噪阈值，从而剔除整个波形数据中的噪声振幅值。四分位数方法适合于具有脉冲特征信号数据的平滑去噪，四分位数位置的计算如公式(4)所示。

(4)

2.3.2 波形分量个数确定

全波形数据高斯分解的准确性和效率性与初始参数的选取有很大关系。本研究依据获得的数据以离散点数据的形式给出，提出了离散型函数求解拐点。根据拐点特性，拐点前后曲线的上升(或下降)趋势不一样，若点P为该曲线的拐点，则其肯定满足公式(5)。

(kAB-kBP)(kPC-kPD)<0。

(5)

式中：A、B表示点P前任意两点，C、D表示点P后任意两点，kAB，kBP，kPC，kCD分别表示直线AB、直线BP、直线PC、直线CD的斜率。该算法尤其适用于离散点数据横坐标为递增等间距的情况，采集的LeicaALS60全波形数据本身是递增等间隔采样的数据。所以在离散点数据递增等间隔采集时，若点(tm，f(tm))为波形的拐点，则满足的判断条件由公式(5)变换为了简化公式(6)。

G(m)=(f(tm)+f(tm-2)-f(tm-1))(f(tm+2)+f(tm)-2f(tm+1))<0。

(6)

式中：m=3，4，5，…，N-3，N表示全波形数据采样点个数，通过全波形数据的读取结果可知N为128。全波形数据可能存在多个拐点，当算法开始时，需要将所有的序列点{tm，f(tm)}依次代入公式(6)，计算出G(m)，当G(m)小于0时算法即可终止并存储记录m时的点(tm，f(tm))。由于高斯分布关于直线t=μ对称，基于此将高斯分布的两个拐点分别称为左拐点和右拐点，自然地，一个高斯分布由其一个左拐点(tleft=μ-σ，f(tleft))及其一个右拐点(tright=μ-σ，f(tright))唯一确定，进而确定该波形数据中的高斯分量个数[12-13]。

2.3.3 高斯分量初始特征参数估计

为解决实际波形数据中的左右拐点划分，本研究在去除背景噪声后将检测出的拐点分别与其右边第一个点连成直线，如果直线斜率大于零，则该拐点为左拐点；如果直线斜率小于零，则该拐点为右拐点。从而可将判断出的所有拐点分为左、右拐点。

小学时期正是学生成长的重要时期。在这个时期，学生的身心都在发展过程中，是对于学生进行思想培养、促进思维成长的重要阶段。在过去的教学当中，课堂上的主体是教师，学生知识被动的进行学习，学生更多是死记硬背，使得学生的思想固化，阻碍了学生多种思维的发展，以至于影响了学生未来的成就。因此在现今的小学教学当中，应该提出学生的主体性地位，培养学生多方面能力，促进学生的发展。

(7)

在确定了一个高斯分量的左拐点与其右拐点后，可以反推出其位置tm和半波宽σm的初始估计值，同时，将其左拐点和右拐点之间点的最大振幅值作为其峰值Am的初始估计值，它们的初始估计值如公式(7)所示。据此可以得到每个高斯分量的初始参数(Am、tm、σm)，从而完成该波形数据高斯分量初始特征参数估计步骤。

2.3.4 波形能量信息计算

高斯分量的能量信息是其波形拟合曲线与采样时间轴所形成的面积，因而其计算方法如公式(8)所示。

(8)

式中：Em表示第m个高斯分量的能量值；Am表示通过波形拟合计算出的第m个高斯分量的波峰振幅值；tm表示通过波形拟合计算出的第m个高斯分量的位置；σm表示通过波形拟合计算出的第m个高斯分量的标准差(半宽)。

2.4 全波形数据激光穿透指数

本研究采用统计模型借助全波形激光穿透指数LPIfi对森林地上生物量进行分析和验证。全波形数据激光穿透指数反映了激光脉冲的穿透能力，由以上内容得知，本研究第i块样地内的全波形激光穿透指数FiLPI恰好如公式(9)所示。

(9)

式中：i=1，2…80，N表示样地内个数；Eci和Egi分别表示第i块样地内的全波形激光脉冲中的冠层返回脉冲能量和地面返回脉冲能量，需要通过波形能量信息计算公式(8)计算得到。

从外业采集的80个样地数据中随机选取60个作为建模样本，余下20个作为检验样本分别进行森林地上生物量回归模型的建模和验证。建模精度用决定系数R2评价，其值越大，说明模型构建结果越好；模型预测精度用RMSE(Root Mean Square Error，RMSE)进行评价，其值越小，说明模型预测效果越好。

3 结果与分析

3.1 样地生物量计算

根据外业样地实测数据按照2.1中的方法得到样地AGB统计结果见表2。

表2 样地生物量统计结果

3.2 全波形数据读取结果与显示

本研究按照上述方法借助IDL8.3平台实现了研究区机载全波形LiDAR数据的读取，并将其以波形的方式显示出来，结果如图1所示(某一条全波形数据)。此外，不仅顺利读取了全波形数据，还同时获取了与全波形数据相关的参数，如采样的起始时间、结束时间、采样间隔、采样个数和坐标位置参数等，这些数据为后续的波形分析处理与波形特征信息提取提供基础参数信息。

图1 编译程序读取并显示的某条全波形数据Fig.1 A waveform record of this compiled program

由结果可知波形采样数为128，相较于传统离散回波LiDAR，全波形LiDAR不仅能记录目标物的三维坐标信息[14-15]，还能存储整个返回波形。由于读取全波形数据文件需要专业软件，本研究基于采集的LAS1.3机载全波形LiDAR数据，实现了全波形数据的快速读取及可视化，为后续的波形分析处理与波形特征信息提取提供基础。

3.3 全波形数据分析处理与波形特征信息提取结果

本研究根据后续工作的需要，通过上述阐述的方法对研究区样地内的Leica ALS60 LAS1.3格式的全波形数据进行了处理。共有80个边长为10 m的方形样地，记录了约 8.2 万条的全波形数据，其中大约99.93%的都能得到很好处理，其余的根据局部最大值法进行了处理。某条原始全波形数据的各步骤处理结果如图2～4所示。

图2 某条原始全波形数据的四分位数法噪声估计Fig.2 Noise estimation of the quartiles of the raw full-waveform data

图3 拐点计算结果与原始全波形数据的比较Fig.3 Comparison between result of inflection point calculation and the raw full-waveform data

图4 高斯分量波形拟合结果与原始全波形数据的比较Fig.4 Comparison between component waveform fitting and the raw full-waveform data

在全波形扫描器的情况下，波形分解还可提供脉冲宽度和脉冲强度作为属性，其中脉冲强度是返回信号的积分，因此该强度值物理上是脉冲能量。相比之下，传统离散回波LiDAR系统仅提供返回信号的三维坐标，并且在许多情况下也提供强度，然而一般情况下对该强度值是如何在这些系统中得到的知之甚少，通常它是返回脉冲的最大振幅值，其被错误地称为强度[14]。

通过波形数据的高斯分解，可以提取每个波形分量的振幅、宽度、位置等重要的波形特征信息，进而计算每个波形分量的能量信息，即图5中所示的起波点与止波点之间围城的面积，且其计算方法如公式(8)所示。

图5 Leica ALS60全波形数据及波形特征信息Fig.5 Leica ALS60 full-waveform data and waveform parameters

3.4 模型构建结果与精度分析

经统计分析得知，通过机载全波形LiDAR数据求出的全波形激光穿透指数FiLPI与实测AGB有较好的相关性(R2=0.885，RMSE=0.095)，模型构建结果如图6所示，建模结果比较理想且实测AGB与全波形激光穿透指数FiLPI具有明显的负相关关系：AGB实测值随着FiLPI的增大而减小，即当森林地上生物量越大或森林越茂密时激光脉冲的穿透能力反而越小，这与森林样地实际情况基本吻合。

图6 FiLPI与AGB实测值的统计回归关系Fig.6 Statistical Regression Relation between FiLPI and Measured AGB

通过验证样本进行预测，将计算出的AGB的估测值与AGB的实测值进行对比，可以直观地得到估测值与实测值之间的误差，如图7所示，可知AGB实测值与AGB估测值之间存在一定误差，但出现较大误差的情况极少。

通过计算得到验证样地的RMSE为21.12 kg，且误差波动较小，这说明该方法能够有效地进行森林AGB的估测。因此，基于本研究的机载全波形LiDAR数据的相应算法与研究方法可以准确地估测森林AGB，弥补了全波形数据估测森林AGB的方法和思路。

图7 实测AGB与预测AGB的预测误差Fig.7 The prediction error of the field-measured AGB against the predicted AGB

4 结论与讨论

本研究基于机载LiDAR全波形数据，通过全波形数据相关处理算法的实现，能够准确地提取每个波形分量以及相应的波形特征信息，在实际处理中表现出较好的适应性。此外在每条全波形数据处理过程中，存储记录了其中每个波形分量的波峰振幅值、位置和波宽等波形特征信息(参数)，尤其是首末次波形分量的，并计算存储对应的能量信息。进而提取出全波形激光穿透指数FiLPI估测参数，并通过统计模型估测了森林样地AGB，表明本文基于机载全波形LiDAR数据的相应算法研究可以准确地估测森林AGB，弥补了全波形数据估测森林AGB的方法和思路。

由于本研究的林型为阔叶林，树种为白桦，在林型和树种的多样性方面比较单一。因此，在下一步工作中，希望将本研究中的算法应用到不同类型的林区，以使机载全波形LiDAR数据森林AGB估测算法研究能够得到推广应用验证，以提高相关算法研究及统计模型的理论性和可推广性。

[1]尤号田，邢艳秋，王萌，等.小光斑激光雷达数据估测森林生物量研究进展[J].森林工程，2014，30(3)：39-42.

[2]李妍，徐新良，张超.中国乔木林碳储量变化研究[J].森林工程，2015，31(4)：50-55.

[3]李增元，刘清旺，庞勇.激光雷达森林参数反演研究进展[J].遥感学报，2016，20(5)：1138-1150.

[4]Bortolot Z J，Wynne R H.Estimating Forest Biomass Using Small Footprint LiDAR Data：An Individual Tree-Based Approach that Incorporates Training Data[J].Isprs Journal of Photogrammetry & Remote Sensing，2005，59(6)：342-360.

[5]刘清旺，李增元，陈尔学，等.机载LIDAR点云数据估测单株木生物量[J].高技术通讯，2010，20(7)：765-770.

[6]王涛，龚建华，张利辉，等.基于机载激光雷达点云数据提取林木参数方法研究[J].测绘科学，2010，35(6)：47-49.

[7]王素元，马洪超，王杰栋，等.基于分组LM算法的全波形LiDAR高斯分解[J].测绘与空间地理信息，2016(7)：144-147.

[8]邱赛，邢艳秋，田静，等.星载LiDAR与HJ-1A/HSI高光谱数据联合估测区域森林冠层高度[J].林业科学，2016，52(5)：142-149.[9]王正文，王德利，臧传来，等.大兴安岭次生林白桦对林下日阴菅及其它主要草本植物的影响[J].生态学报，2001，21(8)：1301-1307.

[10]杨群.在VB利用内存映射文件读取大数据量原始遥感影像图[J].计算机工程与应用，2003，39(6)：135-137.

[11]卢昊，庞勇，徐光彩，等.机载激光雷达全波形数据与系统点云差异的定量分析[J].武汉大学学报信息科学版，2015，40(5)：588-593.

[12]段乙好，张爱武，刘诏，等.一种用于机载LiDAR波形数据高斯分解的高斯拐点匹配法[C]//全国地理信息科学博士生学术论坛，2014-11-21，2014：189-197.

[13]Mallet C，Bretar F.Full-Waveform Topographic Lidar：State-of-the-Art[J].Isprs Journal of Photogrammetry & Remote Sensing，2009，64(1)：1-16.

[14]Reitberger J，Schnörr C，Krzystek P，et al.3D Segmentation of Single Trees Exploiting Full Waveform LIDAR data[J].Isprs Journal of Photogrammetry & Remote Sensing，2009，64(6)：561-574.

[15]Wagner W，Ullrich A，Ducic V，et al.Gaussian Decomposition and Calibration of a Novel Small-Footprint Full-Waveform Digitising Airborne Laser Scanner[J].Isprs Journal of Photogrammetry & Remote Sensing，2006，60(2)：100-112.

Estimation Algorithm of Forest Aboveground Biomass Based on Airborne Full-waveform LiDAR Data

Xing Yanqiu，Yao Songtao，Li Mengying，Xie Jie，Yan Can

(Center for Forest Operational Environment，Northeast Forestry University，Harbin 150040)

In order to improve the estimation accuracy of forest aboveground biomass(AGB)，this paper takes birch forest as the research object and airborne full-waveform light detection and ranging(LiDAR)data as the research data.First，the algorithms of LiDAR data reading and full-waveform information extraction were proposed.Then the waveform energy corresponding to each sample of full-waveform data was extracted and the full-waveform laser penetration indexFiLPIwas calculated.After that，the statistical regression model of AGB was established by using the indexFiLPI.Finally，the estimated value of forest AGB was compared with the measured value.Results showed that the indexFiLPIhas a good correlation with forest AGB(R2is 0.885，RMSEis 0.095)，and the error between the estimated value and the measured value of forest AGB was small，which improved the estimation accuracy of forest AGB and provided the method and idea of forest AGB estimated by airborne full-waveform LiDAR data.

Airborne LiDAR；full-waveform data；waveform data reading；waveform feature information extraction；birch forest；aboveground biomass

2016-04-18

国家林业局林业公益性行业科技专项基金(201504319)；国家重点基础研究发展计划项目(2013CB733404)

邢艳秋，博士，教授。研究方向：森工管理与林业信息工程。E-mail：yanqiuxing@nefu.edu.cn

邢艳秋，姚松涛，李梦颖，等.基于机载全波形LiDAR数据的森林地上生物量估测算法研究[J].森林工程，2017，33(4)：21-26.

S 758.5

A

1001-005X(2017)04-0021-06