基于非参数方法的黑河地区森林地上生物量估算及时空格局

马 骁 王 华 何 瑞 高厚兴 徐 梦

(东北林业大学林学院,黑龙江 哈尔滨 150040)

本研究以黑龙江省黑河地区两期(2005年和2010年)Landsat TM5遥感影像和森林资源连续清查固定样地数据为研究对象，分别提取纹理、地形和气象因子等特征变量，依据%IncMSE指标筛选出前六个重要性变量，以RF模型估算森林AGB并进行精度评价与时空格局分析，为更大尺度估算森林AGB提供森林可持续经营决策依据。

1 研究区域与数据

1.1 研究区概况

黑河地区位于黑龙江省北部(47°42′～51°03′N，124°45′～129°18′E)，面积约为6.9万平方公里，包括爱辉区、北安、五大连池、嫩江、孙吴县和逊克县。有林地面积291.2万公顷，辖区处于中高纬度地带，为寒温带大陆性季风气候，年均气温-1.3～0.4℃，年平均降水量491～540mm，日照时数2562～2677h，无霜期90～120d[6]。

1.2 遥感数据及预处理

Landsat TM5数据下载于地理空间数据云网站(http://www.gscloud.cn/)。覆盖研究区域的两期共7景影像均选在植被生长季(即5-9月)，与样地调查时间基本一致。由于受云量影响，部分选用相近的年份。首先用ENVI5.3对数据进行辐射定标和大气校正。由于影像获取时间不同，存在着较大的辐射差异[1]，因此，本研究在图像镶嵌前增设这种相对辐射校正处理，之后对多景数据沿道路、河流、山脉等处理接边线完成镶嵌。

1.3 固定样地数据

野外调查数据为黑龙江省黑河地区2005年和2010年两期国家森林资源连续清查样地数据(图1)。其中，2005年为995块样地，2010年为1341块样地(样地面积为0.06hm2)。利用简单随机抽样方法，抽取75%的样地作为训练样本，剩余25%作为检验样本。

图1 研究区位置及样地分布图

1.4 气象及地形数据

本研究考虑将气象和地形因子为估算森林AGB的驱动因子参与到估算模型中。气象数据来源于中国气象科学数据共享服务网(http://data.cma.cn/)，为2005年和2010年生长季(5-9月)日值数据集，包括日照(0.1h)、20～20h降水量(0.1mm)及平均气温(0.1℃)，分别取均值后进行反距离权重(Inverse Distance Weight, IDW)插值，分别获得研究区气象因子空间分布。地形数据为ASTER GDEM，空间分辨率30m，获取于地理空间数据云网站(http: //www. gscloud. cn/)，之后利用ArcGIS计算坡度、坡向。

2 研究方法

2.1 样地AGB计算

根据董利虎[7-8]东北地区主要树种的生物量异速方程分别对树干、树枝和树叶部分的生物量进行估算，相加后得到单木AGB，公式为：

w=α·Db

式中：w为生物量；D为胸径；α、b为模型参数，各种树种模型参数详见文献[7-8]，针对东北地区林分类型的模型计算不同林分类型的AGB，利用固定样地每木检尺数据按树种来计算每株树的AGB，将样地内单木生物量相加并除以样地面积，得到单位面积AGB。

2.2 特征变量的提取

特征变量包括植被指数、纹理因子、主成分分析(Principal Component Analysis, PCA)、气象因子和地形因子共五类。对预处理后的遥感影像提取13个植被指数(详见表1)。本研究用原始波段进行主成分分析以及计算基于灰度共生矩阵(Gray-level co-occurrence matrix, GLCM)的8个纹理因子和5个窗口大小(3×3、5×5、7×7、9×9和11×11)下的纹理因子，方向均为0°、步长为1，公式详见表2。之后利用ArcGIS 10.5的“多值提取到点”工具分别提取2005年和2010年固定样地点的13个植被指数、6个原始波段的不同窗口下8个纹理因子、6个PCA图层、3个生长季气象因子(平均日照、平均降水量、平均气温)及3个地形因子(DEM、坡度、坡向)。

表1 植被指数计算公式

表2 纹理因子计算公式

2.3 随机森林算法

随机森林(Random Forest, RF)算法是许多分类与回归树的集合，是一种非常有效的预测工具。本研究使用R统计软件的Random Forest包实现。由于多余的特征变量会降低模型效率，采用importance函数进行特征变量筛选，其中表现特征变量影响力的指标有均方差增量百分数%IncMSE和节点纯度影响力(IncNodePurity)，本研究选用%IncMSE评价RF模型中重要性特征变量，其表示为变量对精度平均减少地贡献程度。特征变量值越大，则对模型精度贡献越大、重要性越高。

2.4 多元逐步回归模型

为了同非参数的RF算法进行比较，本研究也采用多元逐步回归模型(Multiple Stepwise Regression, MSR)进行森林AGB的 估算。为了筛选出适用于MSR模型的变量，降低变量间的多重共线性。选择在Pearson相关系数显著性检验基础上，引入方差膨胀因子(Variance Inflation Factor, VIF)，剔除VIF值大于10的变量，再进行逐步回归，通过设置标准化残差阈值来剔除异常样本点，同时保证剔除样本点数量小于训练样本数量的10%。将检验样本带入最终模型中，求出预测值并计算模型的精度。

2.5 精度评价

3 结果与分析

3.1 RF算法中的参数与特征变量筛选

为了突出RF算法中MSE随Tree的数量变化趋势，取501～3500之间的值观察其变化趋势。随着Tree数量的增加，MSE逐渐降低并趋于平稳。当ntree分别设置为3000和6000时，模型MSE趋近于平稳。由此可见，仅依靠提高Tree的数量难以提高RF算法精度，考虑计算的效率，本研究取3000作为ntree的最终值。之后以%IncMSE进行重要性排序，分别取前六个变量作为最终建模所用变量，如图2所示。

图2 经重要性排序后保留参与RF建模的特征变量及重要性排名 (a)2005年;(b)2010年

3.2 MSR模型中的特征变量筛选

取RF回归精度最佳特征变量组合进行MSR，通过相关系数显著性检验和VIF指数的筛选出参与MSR的特征变量及残差分析，得到适用于MSR方程的特征变量，如表3所示。

表3 逐步回归模型所用变量

3.3 模型估算结果及精度评价

3.3.1 RF中是否引入纹理因子比较

图3 两期检验精度

3.3.2 RF与MSR的比较分析

表4 不同模型的精度检验结果

3.4 黑河地区两期森林AGB分布

通过对模型的精度检验，将两期最优的RF模型(3×3窗口纹理+非纹理因子)用于森林AGB估算。利用R统计软件实现全域估算，将输出结果在ArcGIS 10.5中进行成图，获得两期森林AGB分布图，如图4所示。对结果进行统计，取值范围均在0～140(t/hm2)之间。

3.5 黑河地区两期森林AGB时空变化分析

将2010年与2005年的森林AGB分布图相减，从而获得黑河地区森林AGB时空变化分布图，如图5所示。利用变化矩阵，进一步观察不同等级的AGB变化情况，结合直方图将整体划分为低：0～40(t/hm2)、中：40～80(t/hm2)和高：大于80(t/hm2)三个等级(图6)。由于两期参与建模变量中均有高程变量，考虑两期变化空间分布情况，结合行政区分析变化与高程之间的关系。利用ArcGIS 10.5中“以表格显示分区统计”功能，获取不同行政区内不同高程区间对应的森林AGB变化均值。

总体上分析黑河地区的森林AGB呈现增长趋势，2005年平均森林AGB为40.64(t/hm2)，到2010年增加至44.75(t/hm2)。从图6可以看出，这种增长趋势并非均匀分布在整个黑河地区，其东部和北部部分林地为主要增长区域，而西北、中部和南部靠近平原区域的森林AGB呈现降低趋势。结合变化矩阵的结果分析，虽然有41.8%的高生物量分布区域流向了中低分布区域，但是2005年高生物量和中生物量区域分别占总面积的为9.5%和54.5%，因此，主要增长集中在高生物量区域。通过分析行政区域分析森林AGB变化均值与高程之间的关系，中海拔的森林AGB所受扰动较小，除嫩江县外均呈现较好的增长趋势，在平原分布较多的五大连池市、北安市及嫩江县，低海拔林区所受扰动较大，降低较为明显，这也与实际情况相符。

图4 黑河地区两期森林AGB分布图 (a)2005年;(b)2010年

图5 黑河地区森林AGB时空变化分布图

图6 两期森林AGB变化矩阵

4 结论

基于非参数方法(RF算法)实现对黑河地区两期Landsat TM5遥感数据进行森林AGB估算及变化分析，在研究中使用多种植被指数、PCA、地形、气象及不同窗口下的纹理因子，并对同一期内影像镶嵌前进行相对辐射校正处理。通过对比分析不同组合和MSR精度，引入一定计算窗口纹理因子能提高模型精度，RF算法建模表现更优。总体上黑河地区森林AGB呈现增长趋势，依据变化矩阵高森林AGB区域增长为较为明显，中海拔区域林地为增长较明显且所受扰动相对较小。本研究中仍然存在估算饱和现象，变化分析中尚未探究引起变化的扰动因素，这些还需进一步探索。