基于激光雷达数据的东北林区航空林分材积表编制

曾伟生,孙乡楠,王六如,王 威,蒲 莹

(国家林业和草原局调查规划设计院,北京 100714)

激光雷达是一项新兴的主动遥感技术,它通过发射激光来测定传感器与目标物之间的距离[1-3]。激光雷达分为对空和对地两大类,而对地激光雷达又可分为星载、机载、车载、定位4类,其中最为常用的是机载激光雷达系统[2,4]。机载激光雷达通过获取较大范围内森林的三维扫描数据,可用来估计蓄积量等主要森林参数[1,2,4]。由于激光雷达在估计林木高度和林分结构方面具有独特优势,自20世纪80年代中期[5]首次成功用于估计森林蓄积量以后,基于激光雷达开展森林蓄积量估计方面的研究,已经取得了很多成果[6-14],国外已经开始将其应用于森林资源清查[15-18]。我国开展激光雷达技术应用于林业方面的研究相对较晚[2-4],目前,利用机载激光雷达技术估计森林蓄积量仍基本处于研究阶段[10,11,14,19],尚未见到其应用于森林资源调查实践的报道。

森林蓄积量是森林资源调查的核心指标,已作为两项森林资源约束性指标之一纳入生态文明建设评价考核指标体系。激光雷达技术的出现,为森林资源调查领域的森林蓄积量估计开辟了一条新的途径。在国家机构改革的新形势下,如何充分发挥激光雷达技术的优势,为新时期森林资源调查监测体系的发展提供技术支撑,是当前面临的一项重要课题。本研究将利用东北林区10种森林类型的1 947个样地的地面实测蓄积量数据和机载激光雷达数据,通过利用多元线性回归、非线性回归及哑变量建模方法[20-22],建立基于相同激光雷达变量和不同森林类型参数的蓄积量模型,既为森林资源调查提供计量依据,也为规范森林蓄积量建模与评价方法提供科学参考。

1 数据与方法

1.1 数据资料

本研究所用数据来自中国国土勘测规划院的招投标项目“主要树种航空林分材积表编制”(编号:GXTC-A-19070081),包括地面样地调查数据和机载激光雷达数据两部分。

1)地面样地调查数据。该数据涉及东北林区的云冷杉林(Piceaspp.&Abiesspp.)、落叶松林(Larixspp.)、樟子松林(Pinussylvestrisvar.mongolica)、红松林(Pinuskoraiensis)、栎树林(Quercusspp.)、桦树林(Betulaspp.)、杨树林(Populusspp.)、榆树林(Ulmusspp.)、椴树林(Tiliaspp.)和水胡黄林(Fraxinusmandshurica,Juglansmandshurica&Phellodendronamurense)等10种主要森林类型,共计2 020个样地,其中每个类型不少于200个。调查范围覆盖内蒙古大兴安岭、黑龙江大兴安岭、伊春、松花江、牡丹江、长白山等6个林区的12个片区,调查时间为2019年9月至11月。样地为600m2圆形样地,中心点用RTK技术定位,定位精度在1m以内。除每株样木测量胸径以外,还测量了15株不同径阶的样木树高,以此为基础建立树高-胸径回归模型,推算每株样木的树高,并依据部颁二元立木材积表计算样木材积,从而得到样地的每公顷蓄积量,作为建模的目标变量(Y,m3/hm2)。

2)机载激光雷达数据。该数据范围与地面样地调查相同,获取时间为2019年9月至10月,数据点云密度大于4个点/m2。在对原始数据进行预处理的基础上,利用LiDAR360软件进行激光雷达点云数据的分类、平差,以及数字高程模型(DEM)、数字表面模型(DNS)和冠层高度模型(CHM)的处理,最后提取反映林分高度、密度和结构信息的98个变量作为建模的解释变量,其中高度变量56个(X01—X56),强度变量42个(X57—X98)。

在建模之前,需要对2 020个样地中的异常数据作剔除处理。由于涉及的解释变量个数达到98个之多,很难逐一绘制目标变量与解释变量之间的散点图来判定和剔除异常数据,本研究首先采用全部数据建模,再根据主要解释变量的残差图来判定和剔除异常数据,但剔除异常数据的比例原则上要控制在5%以内,最后参与建模的样地数为1 947个。表1为10种森林类型参与建模样地的每公顷蓄积量的特征值。

表1 建模样地每公顷蓄积量的特征值

1.2 建模方法

基于前述10种森林类型1 947个样地的机载激光雷达数据和地面实测数据,首先建立多元线性和非线性回归模型,并通过对比分析,研究确定合适的基础回归模型；然后利用哑变量建模方法,建立基于相同激光雷达变量和不同森林类型参数的蓄积量标准模型,并编制林分蓄积量表。

1.2.1多元线性回归

采用逐步回归方法确定有显著相关的激光雷达信息变量,拟合如下多元线性回归模型:

Y=a0+a1X1+a2X2+…+akXk+ε

(1)

式中:Y为目标变量(每公顷蓄积量,m3/hm2)；X1,X2,…,Xk为解释变量；a0,a1,…,ak为模型参数,参数的t值原则上应该大于2,否则视为无统计学意义(即与0无显著差异),从模型中剔除；ε为误差项,假定其服从均值为0的正态分布。

1.2.2非线性回归

在线性模型(1)基础上,再建立非线性回归模型:

Y=b0X1b1X2b2…Xkbk+ε

(2)

式中:b0,b1,…,bk为非线性模型参数,其他符号的意义同模型(1)；模型(2)的参数求解采用非线性回归估计方法。

通过对多元线性模型与非线性模型的相关评价指标进行对比分析,确定与森林蓄积量相关最紧密的激光雷达信息变量,并提出统一形式的模型作为基础回归模型。

1.2.3哑变量模型

在确定蓄积量基础模型后,采用哑变量建模方法[20-22],通过引入哑变量代表不同的森林类型,建立基于相同激光雷达变量和不同森林类型参数的蓄积量标准模型。以二元非线性模型为例,其表达式如下:

Y=(∑b0i×Si)×X1(∑b1i×Si)×X2(∑b2i×Si)+ε

(3)

式中:Si为反映不同森林类型的哑变量(i=1,2,…,10)；b0i,b1i,b2i为不同森林类型的参数；其他符号同前。模型(3)的参数求解方法同模型(2)。

1.2.4模型评价

用于模型评价的指标包括以下6项:确定系数R2、估计值的标准差SEE、总体相对误差TRE、平均系统误差ASE、平均预估误差MPE和平均百分标准误差MPSE[22-23]。其中MPE和MPSE的计算公式如下:

(4)

(5)

对建立的蓄积量回归模型,计算以上6项指标值,根据指标大小进行模型评价。从实用性角度考虑,一般要求模型的TRE和ASE均在±5%以内,MPE小于5%,MPSE小于20%。另外,残差图也是评价模型的重要参考依据。一个好的模型,残差应当呈随机分布。

2 结果与分析

利用10种森林类型1 947个样地的蓄积量数据和激光雷达数据,分别拟合多元线性模型(1)和非线性模型(2)。基于篇幅考虑,此处仅列出模型的评价指标如表2所示,每个模型的解释变量及相应的参数估计值从略。

表2 多元线性和非线性蓄积量模型的评价指标

从表2中的解释变量个数k可知,非线性模型相比线性模型要简单一些,多数森林类型仅用2个变量就能得到满意的结果。从表3的评价指标对比来看,总体相对误差TRE几乎无差异；平均预估误差MPE差异也不大,且都达到了预期的5%以内；确定系数R2和估计值的标准差SEE线性模型好于非线性模型；但平均系统误差ASE和平均百分标准误差MPSE则大部分非线性模型好于线性模型,尤其是云冷杉林、红松林、水胡黄林3种类型的ASE和MPSE,其差异更为突出。究其原因,主要是线性模型的估计值出现极小值甚至负值导致的。综合来看,非线性模型(2)的拟合效果要好于线性模型(1),特别是可以避免出现负值。因此,应该优先考虑采用非线性模型。

关于非线性模型(2),10种森林类型中有6种仅包含2个解释变量,其中作用最大的变量是反映激光雷达点云高度的变量(X24,X26,X35,X36或X44),第二个变量是反映点云强度的变量(X81,X82,X86,X90或X92)。如果将云冷杉林、落叶松林、红松林和水胡黄林的非线性模型从三元或四元模型简化为二元模型,其确定系数R2分别从0.771,0.731,0.846和0.466减少至0.761,0.672,0.821和0.403,平均预估误差MPE分别从3.65%,3.71%,2.66%和2.82%增加至3.72%,4.09%,2.85%和2.98%。从模型的实用性考虑,这样的变化幅度是可以接受的。因此,从模型的简约性和实用性出发,最终选定二元非线性模型作为基础回归模型。

再进一步分析二元模型中包含的点云高度变量和强度变量,一个为显著正相关变量,另一个为显著负相关变量。将全部10种森林类型的2 020个样地数据合在一起进行相关分析后发现,正相关最大的变量集中分布在3个区段:X35—X38(点云高度平均值、二次幂平均值、三次幂平均值、中位数),X23—X26(点云高度50%,60%,70%,75%分位数),X04—X07(点云累积高度20%,25%,30%,40%分位数)；负相关最大的变量也分布在3个区段:X90—X91(点云强度平均值、中位数),X79—X83(点云强度50%,60%,70%,75%,80%分位数),X46(点云高度偏斜度)(图1)。这些变量应该是需要重点关注的激光雷达变量,尤其是排在正负相关两端的变量更是核心变量。

注图中29个实心柱形对应的变量为10种森林类型的线性和非线性模型所涉及的29个解释变量

经建模比较,用正负相关最大的2个变量X35和X82(相关系数分别为0.709和-0.420)建立的二元模型,其确定系数反而低于用X35和X90(相关系数-0.407)组合的二元模型。由于X35和X90是最符合预期的激光雷达点云信息变量,本文把基于这2个变量的二元蓄积量模型定义为标准模型。剔除部分异常样地后,利用1 947个样地数据拟合基于X35和X90的二元模型,其确定系数达到0.675,仅低于表2中云冷杉林、落叶松林和红松林非线性蓄积量模型的确定系数,高于其他7种森林类型非线性模型的确定系数。与标准模型相对应的哑变量模型为:

Y=(∑b0i×Si)×X35(∑b1i×Si)×X90(∑b2i×Si)+ε

(6)

式中:X35为点云平均高度,X90为点云平均强度,其他符号同模型(3)。哑变量模型(6)的确定系数达到0.792,10种森林类型的参数估计值及模型评价指标如表3所示。

表3 哑变量模型(6)的参数估计值和模型评价指标

从表3的各项评价指标可知,采用标准模型估计10种森林类型的蓄积量,其拟合效果与表2中的非线性模型(2)差异不大。但是,表2中不同森林类型的非线性模型(2)不仅变量个数有差异,参数值也相差悬殊；而表3中的标准模型(6),不同森林类型都是基于相同的2个变量X35和X90(仅水胡黄林例外),每个变量对应的参数值大小也比较一致。所有模型的平均预估误差MPE均在5%以内,总体相对误差TRE和平均系统误差ASE都在±3%以内且接近于0,平均百分标准误差MPSE为20%左右,但模型的确定系数R2还是存在较大差异,仅红松林达到0.8以上,云冷杉林、栎树林、杨树林在0.7～0.8之间,落叶松林、桦树林、榆树林在0.6～0.7之间,樟子松林、椴树林和水胡黄林在0.6以下。尤其是水胡黄林,其标准模型的确定系数仅为0.111,改用另外2个变量X13和X62后,模型的确定系数也仅为0.440。究其原因,除了水胡黄林这个森林类型树种组成较为复杂、混交程度也大于其他类型外,还有一个重要原因就是样本的结构不够理想,这主要也是取样范围受到了限制,难以找到合适的样地。表1中的数据也表明水胡黄林的蓄积量变动系数最小,最小值和最大值之间的分布区间最窄。

另外,表3中10个森林类型蓄积量模型的残差随解释变量的分布基本是随机的,图2和图3是云冷杉林和落叶松林蓄积量模型的相对残差随解释变量X35的分布情况,其他8个森林类型的残差图也类似,此处略。

图2 云冷杉林蓄积量的相对残差分布 图3 落叶松林蓄积量的相对残差分布

根据表3中10种森林类型标准蓄积量模型的参数估计值,就可以编制林分材积表,表4为基于激光雷达数据的云冷杉林二元林分材积表。

表4 云冷杉林二元林分材积表

3 讨论与结论

机载激光雷达作为一种现代先进技术手段,为准确估计森林蓄积量,提高森林资源调查工作效率提供了可能[1,15-17]。本文结合我国现阶段对森林资源调查监测的新要求,从森林蓄积量这一核心因子入手,基于1 947个样地的地面实测数据和机载激光雷达数据,利用哑变量建模方法,建立了东北林区10种主要森林类型的蓄积量估计模型,模型的确定系数R2在0.44～0.81之间,与Hollaus等[9]、González-Ferreiro等[12]和Bottalicoa等[13]建立的蓄积量模型的确定系数基本一致,略低于Holmgren等[8]研建的2个蓄积量模型的确定系数(分别为0.82和0.90),显著高于范凤云[10]利用0.39个点/m2的激光雷达数据所建刺槐林和油松林蓄积量模型,其相关系数R分别为0.696和0.453。森林蓄积量模型R2的高低,主要取决于模型中采用的激光雷达数据变量。激光雷达变量包括不同分位数高度、密度和最大回波高度等,但不同的研究所选择的变量不尽相同[4]。本研究提取了反映激光点云高度和强度的98个变量来估计蓄积量,最终选择的模型解释变量为点云高度平均值X35和点云强度平均值X90,这既是正负相关最大的2个激光雷达变量,也是最符合预期的2个变量,所以我们将基于这2个变量的蓄积量模型定义为标准模型。Bottalicoa等[13]提取了49个激光雷达变量作为森林蓄积量的解释变量,结果发现点云高度平均值和三次幂平均值与蓄积量的相关最高,与本研究的结果是基本一致的。

森林蓄积量模型的实用性,主要取决于2项误差指标MPE和MPSE的大小,前者反映对总体蓄积量的平均估计误差,后者反映对林分水平蓄积量的平均估计误差。本研究所建10个森林类型的蓄积量模型,MPE在2.88%～4.42%之间,均未超过5%；MPSE在16.76%～25.52%之间,其中红松林、榆树林、椴树林、水胡黄林小于20%,另外6个森林类型除樟子松林略超25%以外,其他均小于25%。《森林资源规划设计调查技术规程》[24]对小班调查蓄积量的精度等级分A,B,C三级,要求相对误差分别不超过15%,20%,25%。因此,10个森林类型的蓄积量模型有4个可以满足B级的精度要求,其他6个可以满足或基本满足C级的精度要求。采用传统调查方法是非常费工费力的,而基于激光雷达数据构建蓄积量模型,只需通过信息提取就能得到符合精度要求的蓄积量估计值。因此,利用机载激光雷达开展森林资源规划设计调查,在技术上是完全可行的。

根据本研究的相关结果,可以得出以下结论:1)利用机载激光雷达技术估计森林蓄积量,非线性模型总体上要优于线性模型，从森林资源调查的实用性考虑,二元与多元非线性模型之间无实质性差异；2)与森林蓄积量正负相关最大的2个激光雷达变量是点云高度平均值和点云强度平均值,以这2个变量为基础建立的标准蓄积量模型,可以作为森林蓄积量估计的通用模型；3)基于哑变量建模方法,建立具有不同森林类型参数和相同激光雷达变量的蓄积量模型,是林业建模实践中值得推广的一种可行做法；4)本文所建东北地区10种森林类型的蓄积量回归模型,其预估精度达到《森林资源规划设计调查技术规程》的要求,可以在实践中推广应用。