基于神经网络和不同立地质量的森林蓄积量遥感估测

2019-04-08，，，

山东科技大学学报(自然科学版) 2019年2期

，，，

(1.山东科技大学测绘科学与工程学院，山东青岛 266590；2.国家海洋信息中心，天津 300171)

森林被誉为“地球之肺”，在全球生态系统的平衡中发挥着重要作用，是人类以及多种物种赖以生存和发展的基础[1]。森林蓄积量是衡量森林健康与否的重要指标，也是政府掌握国家森林资源状况和制定计划采伐、森林经营管理措施的重要依据。传统的蓄积量测定方法是一、二类森林调查，人力、物力投入巨大，而且调查周期长，难以掌握森林的动态发展状况。20世纪70年代以来，遥感技术如雨后春笋，迅速应用于各个领域，利用遥感技术进行森林蓄积量的定量估测也取得了一系列的重要进展[2-4]。

基于遥感技术对森林蓄积量进行估算，主要是利用遥感影像数据以及少量实地调查数据，采用数学方法构建森林蓄积量的估算模型，从而使得模型可以直接应用在森林蓄积量估算中[5]。国外很早就有遥感数据定量估算森林蓄积量的应用。早在1988年，Nelson等[6]就利用Laser数据估算了森林的蓄积量和生物量。Gemmell等[7]基于专题制图仪(thematic mapper, TM)数据，研究了各波段影像以及其他实测数据对估测森林蓄积量的影响度，为TM影像数据在定量估测蓄积量方面提供了理论依据。Fazakas等[8]使用TM影像数据结合最近邻法估算了森林的蓄积量。国内蓄积量的遥感估算相对国外起步较晚，刘琼阁等[1]提出了基于偏最小二乘法的森林蓄积量遥感估测模型，程武学等[9]对基于小班的各种蓄积量预测的合理性进行了诊断，证实了其合理性。近年来，得益于神经网络超强的非线性处理能力，国内外研究人员将其成功地运用到森林蓄积量的估测中，取得了一定的成果。王臣立等[10]运用3层的后向传输(back propagation, BP)神经网络建立了基于植被指数的热带森林蓄积量估算模型，同等条件下较回归分析模型的预测效果更好。许炜敏等[11]构建了结构为10∶3∶1的杉木林蓄积量BP神经网络模型。吴达胜等[12]建立了涵盖地形、地貌、气候气象、土壤、林分结构特征及光谱特征的森林资源蓄积量预测自变量因子集，并应用改进的BP神经网络模型对森林资源蓄积量进行了预测，所取得的模型充分利用了多源数据，且具备很强的泛化能力。

森林地位级表是森林管理中的常用数表之一。早在上世纪90年代，林昌庚等[13]就提出了适合于中国森林调查数据的地位级表编制方法，周春国等[14-15]对其进行了完善，提出了最佳导线曲线的测定方法。研究表明，区分立地质量等级进行预测建模可以显著提高森林蓄积量的估算精度。王艳婷等[16]运用岭估计法，分不同地类分别建立森林蓄积量估测方程，结果表明，当影响蓄积量估测的自变量间存在复共线性时，分地位级建立模型的精度明显优于统一建模的精度。刘俊等[5,17]建立了基于不同立地质量的杉木和松树的森林蓄积量遥感估测模型，采用了更为准确的地位级划分方法，但在建模时只利用了遥感影像主成分分析的第一主成分进行线性回归建模，并没有充分利用遥感影像和遥感辅助数据，由于森林蓄积量受到多方面因素的影响，单一变量难以准确估算蓄积量，因此估算精度不高。

本研究以2009年TM影像数据和凉水自然保护区森林二类小班数据为数据源，采用地位级法划分林分的立地质量等级，结合回归分析与神经网络分析算法，对红松树种以不区分立地质量和区分不同立地质量等级分别建立森林蓄积量的遥感估测模型，旨在为森林蓄积量的遥感监测提供一种新的思路。

1 研究基础

1.1 研究区概况

凉水自然保护区位于黑龙江省伊春市带岭区，小兴安岭山脉的东南段。地理坐标为东经128°47′8″—128°57′19″，北纬47°6′49″—47°16′10″，东西宽13.0 km，南北长17.0 km，总面积为12 133 km2(图1)。该地区深受海洋性气候的影响，具有明显的温带大陆性季风气候特征，总的特点是冬长夏短，冬季严寒干燥，夏季湿凉多雨。年平均气温-0.3 ℃，年平均最高气温7.5 ℃，年平均最低气温-6.6 ℃；年平均降水量676 mm，年平均相对湿度78%，年平均蒸发量805 mm。土壤类型主要有暗棕壤、沼泽土、草甸土、泥炭土4类。该地区森林资源丰富，其中红松是这里的优势树种。

图1 研究区地理位置与TM影像数据图Fig.1 Geographical location and TM image data of the study area

1.2 数据获取

1.2.1 遥感影像数据

实验使用Landsat TM遥感影像，数据在美国地质勘探局(United States Geological Survey，USGS)官网(https://earthexplorer.usgs.gov)上获取，综合考虑时相和云量等因素，选取成像时间为2009年5月29日覆盖凉水自然保护区的TM遥感影像。利用ENVI5.3对影像进行辐射定标、大气校正、裁剪与镶嵌等预处理。

1.2.2 高程数据

凉水自然保护区位于小兴安岭地区，地形复杂多样，海拔290～550 m，以山地、丘陵地形为主。实验运用全球数字高程模型(global digital elevation model，GDEM)数据，对应的分辨率精度为30 m。

1.2.3 森林资源二类调查数据

获取凉水自然保护区2009年森林资源二类调查小班矢量数据，属性表中包括小班林种、起源、地类、公顷株数、郁闭度、高度、自然度、林分高、优势树高、优势径阶、枯倒树种、下木名称、地权使用权、地被总盖度、地被名称、地貌、部位、坡向、坡度、树种组成、林分类型、树种年龄、树种直径、树高、林层结构、龄级、龄组、小班每公顷蓄积量、小班总蓄积量等。

1.3 研究技术路线

运用2009年凉水自然保护区二类森林资源调查数据和同年度该区TM影像作为数据源，首先，对数据进行预处理后，根据前人研究经验和方法，寻找最合适的立地指数导向曲线来划分森林各小班数据的立地质量等级。其次，筛选出共线性低且与森林蓄积量相关性高的遥感因子变量作为模型的自变量，对自变量和蓄积量之间分别进行基于回归分析和基于神经网络分析的模拟，并进行精度检验，筛选出最合适的森林蓄积量估测模型。研究技术路线如图2所示。

图2 技术路线图Fig.2 Technology roadmap

图3 最佳地位级导向曲线Fig.3 The best sue level guide curve

2 地位级划分与因子选取

2.1 立地质量等级评价

采用地位级法进行立地质量的评价。地位级表以林分平均高与平均年龄的关系为依据, 是对同龄纯林分地区分树种编制、反映林地生产力高低的数表, 是森林经营及森林调查中主要的常用数表之一[10]。按相同年龄时，林分高度的变动程度划分为若干个等级，通常为5～7级，以罗马数字Ⅰ、Ⅱ、Ⅲ…等符号依次表示，将每一等级所对应的各个树龄时的平均高度列成表，即为地位级表[10]。编制地位级表的关键就是导向曲线的选定，前人[10-12]提出了多种计算导向曲线的函数模型与方法，但各种导向曲线模型的选择因数据而异。基于凉水自然保护区森林资源二类小班调查数据，对优势树种为红松且郁闭度0.4以上的小班，依据平均树高和平均年龄建立一条代表中等立地等级条件下的地位级导向曲线[11](图3)。通过不同数学模型比较选出对数与逆函数结合的曲线作为红松地位级导向曲线的最优模型，最终所得方程为：

H=28.866-291.720/A-0.524lnA。

(1)

其中，H代表小班的平均高，A代表小班的平均年龄。

目前的研究中，通常以导线曲线为基础，采用比例法确定地位等级的上下界线，以导向曲线的1.35倍为上界线，0.55倍为下界线[10]。用f(A)表示导向曲线方程，则上下界模型分别为1.35f(A)和0.55f(A)。把杉木分为5个地位级数，各等分高度间隔为0.55f(A)、0.75f(A)、0.9f(A)、1.05f(A)、1.35f(A)，按照以上分类办法，得到地位级表，如表1所示。

结合表1，将地位级为Ⅰ、 Ⅱ的归为立地质量差，地位级为Ⅲ的归为立地质量中等，将地位级为Ⅳ、Ⅴ的归为立地质量优。以此数据作为立地质量划分的结果分别用来拟合不同立地质量的森林蓄积量的估算模型。

2.2 影响因子定义

根据以往的研究，依据单变量或者多变量的线性组合，可以不同程度地突出某些需要的信息以及抑制无关的信息[13]。本次调查数据共有小班431个，删除优势树种不是红松和郁闭度低于0.4的小班数据，并利用标准差分析法剔除明显的异常数据。剔除异常样本后各级统计数据如表2。

表1 红松地位级Tab.1 The site lovel of Korean Pine

表2 不同立地质量等级的建模及验证样本数Tab.2 Modeling and validation samples of different site qualities

查阅相关文献后，在前人研究基础上，从TM影像中选取可能影响蓄积量的遥感信息参数，如TM影像的第一波段(B1)、第二波段(B2)、第三波段(B3)、第四波段(B4)、第五波段(B5)、第七波段(B7)；多波段组合(B1/B3、B2/B4、B3/B5、B3×B4/B7、B4×B5/B7、(B1-B7)/B4、(B3+B4+B5)/(B1+B2+B7)，下文分别用B13、B24、B35、B347、B457、B174、B345127表示)；植被指数如归一化植被指数(normalized difference vegetation index,NDVI)、比值植被指数(ratio vegetation index,RVI)、增强型植被指数(enhanced vegetation index,EVI)、差值植被指数(difference vegetation index,DVI)；主成分分析的第一分量(principal componet analysis,PCA)；经过遥感影像数据反演的土壤湿度(soil moisture，SM)、地表温度(land suface temperature，LST)；经过研究区数字高程数据提取的坡向(Aspect)、坡度(Slope)、海拔高度(ASL)等23个变量作为因子，如表3所示。

表3 变量因子表Tab.3 Variable factor table

自变量的选择是建立模型的首要任务。分别计算表3选取的因子在不同立地质量下与获取的蓄积量之间的相关系数R和他们的方差扩大因子(variance inflation factor,VIF)，所得结果如表4所示。

如表4可见，除去地形和地表温度、湿度之外，森林的蓄积量与以上各个因子的相关系数均较高，达到了0.01的显著性水平。无论是否划分立地质量等级，NDVI与蓄积量的相关性均较高，其次是TM影像的第3波段、遥感图像主成分分析的第一主成分等变量，说明NDVI以及遥感影像的红色波段是监测森林蓄积量的敏感波段。然而，红色波段(B3)、NDVI变量的方差扩大因子(VIF)很高，说明其自变量观察值之间存在严重的共线性。各个变量之间的相关性也比较高，会造成模型的冗余[17]。因此在建模之前必须进行因子的筛选，确定进入模型的自变量。

2.3 因子筛选方法

2.3.1 所有子集回归法

表4 不同立地质量等级的变量因子相关信息Tab.4 Information about the selected variable factors of different site qualities

注：1.**表示其达到了0.01的显著性水平，*表示其达到了0.05的显著性水平。2.VIF值表示自变量观察值之间复共线性程度，VIF越大代表共线性越显著。一般来说：当0

应用回归分析处理实际问题时，逐步回归是一种最常用的选择回归变量的方法。但由于这种方法不计算所有可能性，因此得到的变量子集只是局部“最优”，可能会遗漏全局最优的变量子集。因此，本研究选用所有子集回归法，即若有m个可供选择的变量，那么所有可能的回归方程就有2m个，然后利用选元准则挑出最优的方程。常见的选元准则有Cp统计量、自由度调整复决定系数和赤池信息量等。本研究的目的是建立预测模型，故运用使预测值的均方根误差尽量小的Cp统计量[18-19]。马洛斯(Mallows)[20]从预测的角度提出了用Cp统计量来选择自变量。设n为样本数，SSEp表示p个变量模型的残差平方和，则

(2)

使用SAS软件进行最优子集的选择，取Cp统计量作为准则进行所有子集回归，得到森林蓄积量预测的最优子集，如表5所示。

表5 各立地质量等级的最优子集Tab.5 Optimal subset of each site quality grades

2.3.2 基于平均贡献值(MIV)的因子筛选

平均贡献值(mean impact value，MIV)方法选择神经网络输入变量，用于确定输入神经元对应变量的影响大小，其符号代表相关的方向，绝对值大小代表影响的相对重要性。具体过程为：首先，在网络训练终止后，将训练样本P中每一个自变量特征在其原值基础上分别加、减10%构成新的样本P1和P2，将P1、P2分别作为仿真样本利用已经建成的网络进行仿真，得到仿真结果A1和A2，求出A1和A2的差值，即得出该自变量对于因变量的MIV值；其次，根据MIV绝对值的大小为各自变量排序，得到各自变量对网络输出影响相对重要性的位次表；最后，根据位次采用逐步剔除法实现变量的筛选。由于BP神经网络的权重是随机得到的，每次训练生成的MIV值不同，故选取10次计算的均值作为变量的MIV值。使用MATLAB的神经网络工具，构建一个简单的神经网络，并实现各因子的MIV值计算，结果如表6所示。

表6 各因子对网络的平均贡献值(MIV)Tab.6 Mean impact value (MIV) of each factor on the network

3 蓄积量估测模型的建立

3.1 基于回归分析的森林蓄积量估测模型

首先进行单产量的回归，以小班的计蓄积量为因变量、小班提取的TM影像第一主成分遥感因子为自变量，对红松区区分立地质量等级和不区分立地质量两种办法分别建立蓄积量估算模型。根据散点图分布，选取线性函数和幂函数两种模型。各个模型的拟合结果如表7所示。

综合对比可见，无论是何种建模方法，模型的P值都小于0.01，通过了0.01的显著性水平检验，且均通过了F检验，说明森林的蓄积量和TM影像的第一主成分之间存在较好的回归关系。其中，幂函数模型的决定系数均在0.6以上，最高为0.707，较普通线性模型有着显著的提高。另外，基于不同立地质量的模型较不区分立地质量的模型效果略有提升，但提升并不显著。

可见，同等条件下，相较于单产量回归模型，多元线性模型的R2有明显的提高。特别是基于不同立地质量等级进行建模时，模型的拟合效果显著提升，R2分别从0.707、0.620、0.601提升到0.819、0.720、0.712，这一现象表明基于不同立地质量等级分别对森林蓄积量进行建模中的有效性，且遥感因子与森林蓄积量之间存在复杂的关系。为进一步挖掘数据中存在的非线性关系，引进了处理非线性问题表现良好的BP神经网络技术。

表7 单产量因子的回归拟合结果Tab.7 Fitting result of each variable factor

表8 森林蓄积量的多元回归结果Tab.8 Multiple regression result of forest stock volume

3.2 基于BP神经网络的蓄积量预测模型

BP神经网络是当前应用最广泛的神经网络模型之一，是一种按照误差逆向传播算法训练的多层前馈神经网络，具有任意复杂的模式分类能力和优良的多维函数映射能力，可以很好地处理非线性问题[21]。

利用神经网络建模的过程分为数据预处理、网络训练、调参、网络再训练、模型评价等5个阶段。首先，将数据进行归一化以消除量纲的影响；然后进行网络结构的设计，网络共3层，分别是输入层、隐含层和输出层；使用控制变量法分别对比不同传递函数、不同训练函数、不同神经元数等各种网络参数下网络的表现，经多次训练、调参，得到各立地质量等级的网络结构及其在训练集上的系数。其中，模型的传递函数均选择Purelin函数，立地质量为差的数据因样本较少，选择Trainlm为训练函数，其他类别使用的训练函数为Trainrp。模型以网络拟合的均方误差(mean squared error, MSE)为目标损失函数。各模型的网络结构均已在表9中说明。模型的迭代次数设置为1 000次，每个实验重复10次，取10次训练的MSE的均值及其标准差，如表9所示。

可见，利用神经网络技术对不同立地质量的数据进行建模时，模型的均方误差可以收敛于更小的水平，模型的拟合性能得到显著改善。为全面对比分析模型效果，下面进行模型对比验证。

表9 不同立地质量的神经网络模型及训练参数Tab.9 Neural network model and training parameters of different site qualities

3.3 预测模型对比与验证

采用随机预留的验证数据验证模型的精度，将验证集数据应用到模型，得到森林蓄积量的估测值。将估测值与实际蓄积量之间进行线性拟合，得到的R2、斜率及其总体估计精度作为精度评价指标，完成模型的精度检验。其中，曲线斜率越接近1代表预测值和实测误差越小，总体估计精度(P)的计算公式为：

(3)

检验结果如图4和表10所示。综合对比以上验证结果，可见：

1) 模型预测值与实测值之间的R2均达0.7以上，总体估测精度在78%以上，使用遥感影像可以较好地预测森林蓄积量，其中，多自变量模型的预测效果要远好于单产量模型；

2) BP神经网络模型表现要好于单产量回归和多元线性回归模型，其预测值和实测值的R2均达到0.9及以上，曲线斜率也更接近于1，总体预测精度也大大提升，达到了95%以上；

3) 基于不同立地质量等级进行建模，无论采用何种建模方式，较不区分立地质量时，预测精度、R2都有显著的提升。特别是BP神经网络模型的总体预测精度在97%以上，R2达到了0.94以上，达到了很高的预测水平。由此可见基于不同立地质量建立森林蓄积量估算模型的思路是正确的，在同等条件下可以使模型性能得到显著的提升。