基于Landsat 8遥感影像的阿尔泰山天然林生物量估测

2022-07-15张景路朱雅丽张绘芳地力夏提包尔汉

中南林业科技大学学报 2022年6期

张景路，朱雅丽，张绘芳，刘建，地力夏提·包尔汉

（1.新疆林业科学院现代林业研究所，新疆乌鲁木齐 830000；2.阿尔泰山国有林管理局阿勒泰分局，新疆阿勒泰 836500）

阿尔泰山具有丰富的森林植被资源和生物多样性，是新疆重要的生态屏障之一。阿尔泰山天然林担负着维持区域生物多样性和调节生态平衡的功能，因此需要对天然林加以保护及管理。而森林生物量是区域森林生态系统结构和功能的重要体现，对其进行精确地估测对于研究森林生态系统的固碳能力、生产力以及评估天然林的质量和生态效益具有重要意义[1]。

森林生物量的定量测算具有时空尺度复杂、数据种类和数量庞大、需要跨学科知识储备等特点。我国对于森林生物量的研究始于20世纪70年代末至80年代初，采用的主要方法为基于样地的生物量实地测定，并通过常规统计和尺度转化得到区域植被的总生物量。由于传统的基于样地调查的生物量估测方法具有投入大、破坏性高、难以完整反映区域信息和动态变化等的局限[2-3]，因此随着地球信息技术的日益发展，尤其是RS 的连续动态监测和GIS 的空间数据分析能力的完善，越来越多的国内外学者开始基于地面数据和遥感数据结合对森林生物量进行估测。在国外，Giles 等[4]通过巴西、马来西亚和泰国的陆地卫星TM 数据与实地调查生物量结合建立了基于植被指数、多元回归和前馈神经网络3 种生物量预测模型；Erik 等[5]采用机载激光雷达（LiDRA）和合成孔径雷达（InSAR），结合201 个实地调查样地数据，对挪威北部森林地区的地上生物量进行了估算；Maciej 等[6]采用了单通道InSAR 数据，结合瑞典南部森林的32 个样地的实测数据反演得到林高和林冠密度用以预测生物量。在国内，Wang 等[7]采用SAR-L波段数据和19 个位于日本北海道托马科迈森林分地面真实生物量数据建立了森林生物量的经验回归模型，精度为86%；徐婷等[8]基于Landsat 8 OLI 影像，利用多元逐步回归法建立了江苏省常熟市虞山地区不同类型的森林生物量模型；段祝庚等[9]采用六景Landsat 8-OLI 影像提取的波段反射率及植被指数，结合LiDAR 数据建立了粤西北区域森林指数回归模型；余朝林等[10]采用Landsat TM 影像对浙江省临安、安吉和龙泉3 个毛竹产区分别建立了毛竹林生物量遥感估算模型，包括逐步回归、多元线性和Erf-BP 神经网络等模型；罗庆辉等[11]以Landsat 8-OLI 影像提取的各波段灰度值、不同波段灰度值之间的线性和非线性组合及环境因子为基础，建立了天山雪岭云杉林生物量估测模型并取得了良好的拟合精度。

上述研究主要采用植被生长季单一时点遥感数据进行生物量估测模型反演，未考虑不同时点植被光谱反射变化对模型精度的影响。本研究选取了植被生长季内4 期不同时点的遥感影像数据，通过分析各期遥感影像数据特征变量与地面实测生物量的相关性，选取最优时点影像作为研究数据源，选用岭回归、主成分分析、偏最小二乘法等线性模型将遥感影像数据特征变量与地面实测生物量进行拟合，通过比较分析3 种模型的精度指标选取最优模型来估算研究区生物量，分析其空间分布。本研究为提高山区天然乔木林生物量反演模型精度提供了新思路，同时为阿尔泰山天然乔木林的科学管护提供了理论支持。

1 材料与方法

1.1 研究区概况

研究区位于阿尔泰山国有林管理局布尔津分局和阿勒泰分局所辖区域。地理坐标：86°54′01″～88°38′46″E，47°46′15″～49°10′54″N，海拔900～3 800 m，属温带半干旱大陆性气候，夏季多雨，冬季严寒，年均降水量500～700 mm[12]。阿尔泰山为断块山，有4 级夷平面，母岩多为酸性变质岩和花岗岩，森林和森林草原带的土壤主要是山地灰色森林土，山地棕色针叶林土、山地栗钙土，草甸带的亚高山草甸土及高山草甸土以及荒漠草原和荒漠带的山地棕钙土，灰棕荒漠土，高山顶为高山石漠土等。阿尔泰山山地植被垂直分布明显，1 100 m 以下为山麓草原带；1 100～2 300 m 为森林带，乔木树种主要有西伯利亚落叶松Sibirica Ledeb、西伯利亚冷杉Abies sibirica、云杉Picea asperata、欧洲山杨Populus tremula和疣枝桦Betula platyphylla等，灌木树种主要有绣线菊Spiraea salicifolia、忍冬Lonicera japonica Thunb.、樱桃李Prunus cerasifera Ehrhart、兔儿条Spiraeahypericifolia和圆叶桦Betula rotundifolia[13-16]；2 300 m 以上为山地草甸带与亚高山草甸带。乔木优势树种为西伯利亚落叶松，灌木优势树种为绣线菊、兔儿条和圆叶桦。

1.2 数据来源与处理

1.2.1 样地数据收集与处理

样地数据来自2016年新疆第九次全国森林资源连续清查数据，样方面积（28.3 m×28.3 m）0.08 hm2，共90 个，分布情况如图1所示。数据包括样地号、样地坐标、样地地类、样地地形（海拔、坡向、坡位、坡度）、树种、胸径、树高、郁闭度等。

图1 研究区地理位置及样点分布示意图Fig.1 Geographical location and sample distribution of the study area

本研究中生物量采用相关文献[17-19]中与本研究区气候、地貌相似地区的地上单木生物量方程计算，样地所有单木地上生物量之和即为样地地上生物量，为了研究的方便，将样地地上生物量除以样地面积，换算为每公顷地上生物量即单位生物量。地上单木生物量具体计算公式见表1，式中W地上为各样地地上生物量；D为各样地树种平均胸径。

表1 研究区各优势树种单木生物量模型Table 1 Biomass models of dominant tree species in the study area

得到样地天然乔木林单位生物量信息如表2所示。

表2 样地生物量统计特征Table 2 Statistical characteristics of sample plots biomass

1.2.2 遥感影像数据收集与处理

由于不能获取覆盖研究区的2016年6—9月Landsat 8 遥感影像，且研究区树种主要为慢生树种，一年的生物量变化较小，因此本研究采用2015年6月6日、7月8日、7月24日和9月10日共4 期行列号为144/26 的Landsat 8 遥感影像作为基础数据源，影像分辨率为30 m×30 m，可以完整覆盖研究区域。运用ENVI 5.3 对4 期遥感影像进行波段合成、辐射定标、大气校正、地形校正（地形校正采用ENVI 5.3 扩展模块Topographic correction 模块中的VECA 方法实现）和几何校正等处理，得到研究区所需遥感影像数据，并采用面向对象法进行林地分类，得到研究区天然乔木林空间分布数据，并通过林地变更调查数据以及实地调查对分类结果进行评定[12]。运用ENVI 5.3 软件从处理好的遥感影像中提取Band2、Band3、Band4、Band5、Band6、Band7等6 个单波段数据，归一化植被指数（NDVI）、比值植被指数（RVI）、增强型植被指数（EVI）等3 个植被指数，基于二阶概率统计窗口大小为5×5 像素的均值（Mean,Me）、方差（Variance,Va）、协同性（Homogeneity,Ho）、对比度（Contrast）、相异性（Dissimilarity,Di）、信息熵（Entropy,En）、二阶矩（Second moment,Se）和相关性（Correlation）等8 个纹理信息作为构建生物量反演模型的特征变量。

表3 Landsat 8 数字影像信息Table 3 Landsat 8 digital image information

表4 地形校正模型参数设置Table 4 The parameters setting of terrain correction model

图2 研究区遥感影像示意图Fig.2 Remote sensing image of the study area

为验证地形校正的效果，分别从校正前后的遥感影像中随机提取1 000 个样本点，计算其太阳入射角与各波段值进行相关分析，可见除Band1、Band2 外，其余各波段与太阳入射角的相关性明显减少，地形对反射率的影响被有效减弱。

1.2.3 地形数据收集与处理

本研究采用地形数据来自地理空间数据云平台（http://www.gscloud.cn/）提供的30M 分辨率的DEM 数字高程数据。使用ArcGIS 中Slope 和Aspect 工具生成坡度和坡向，用研究区界限掩膜得到研究区海拔、坡度、坡向数据，作为地形因子变量。

1.3 研究方法

1.3.1 生物量估测模型建立

基于遥感影像的森林生物量反演模型是通过遥感影像所提取的特征变量与样地生物量建立关系的具有确定数学关系的模型[20]。本研究选用岭回归、主成分分析和偏最小二乘回归构建生物量模型（表6），将90 个样本随机分为两部分，其中70 个样本用于建模，20 个样本用于检验。

表5 地形校正前后各波段反射率与太阳入射角相关性Table 5 Correlation between reflectance and solar incidence angle before and after terrain correction

图3 研究区地形数据示意图Fig.3 Topographic data of the study area

表6 估测模型表达式†Table 6 Expression of estimation model

1.3.2 模型精度验证

为了比较不同方法所建立的模型的预测精度差异，本研究选取了4 个评价指标（表7），从不同维度对其进行检验，分别为确定系数（Coefficient of determination，R2）、估计值的标准差（Standard error of the estimate，SEE）、平均系统误差（Mean systematic error，MSE）、模型精度（P）。为直观反映估计值与实测值的偏差，绘制预测值与实测值对比散点图及1∶1 线图。

表7 精度指标公式†Table 7 Precision index formula and description

1.3.3 生物量空间分布

运用ArcGIS 10.2软件的Raster Calculator模块，根据最优生物量模型生成研究区天然乔木林生物量空间分布图，运用ArcGIS 10.2 软件的空间分析工具，将生物量空间分布图分别与研究区海拔、坡向、坡度分布图进行叠加分析，得到研究区天然乔木林生物量在海拔、坡向和坡度上的分布情况。

2 结果与分析

2.1 自变量的选取与检验

2.1.1 特征因子选取

本研究选用经地形校正后的单波段（6 个）、植被指数（3 个）、纹理特征（48 个）、地形因子（3个）共60 个特征变量，以处理完成的遥感影像作为基础数据源，提取实测样地中心点所对应像元及相邻8 像元所对应的特征值，并计算平均值作为最终特征变量。将4 期影像各自对应的60 个特征变量值与样地生物量进行相关性检验，采用Person 相关分析方法筛选出在0.01 水平上显著相关的特征变量（图4）。通过对比可得7月8日遥感影像所提取的与生物量显著相关的特征变量在数量和质量上都优于其他影像，因此本研究选取7月8日遥感影像所提取的与生物量相关性大于0.35的Band2、Band3、Band4、Band7、Band4_Di、

图4 与生物量在0.01 水平上显著相关的特征变量Fig.4 Characteristic variables significantly correlated with biomass at the 0.01 level

Band4_En、Band4_Ho、Band4_Me、Band4_Se、Band4_Va、Band7_En、Band7_Ho、EVI、NDVI、RVI 共15 个特征变量作为建模因子。

2.1.2 自变量检验

由于样本实测数据属于截面数据，可能存在一定程度的异方差问题，因此本研究选取了等级相关系数对样本进行了异方差检验。将因变量和残差绝对值做斯皮尔曼相关分析，得到生物量实测值与残差绝对值相关系数为0.202，显著性0.056，可以认为样本数据不存在异方差问题。

表8 生物量实测值与残差绝对值相关性Table 8 Correlation between the measured biomass and absolute residual value

2.2 生物量估测模型建立

2.2.1 岭回归模型

岭回归是一种用于共线性数据分析的有偏估计回归方法[21]。当方程自变量存在共线性时，单个变量的变化会引起其他变量改变，从而导致计算结果的不稳定，而岭回归能够解决该问题。岭回归在最小二乘法的基础上引入了一个可以保证回归系数稳定的正常数矩阵KI，使得每个变量对计算结果的影响应尽可能小，以得到一个固定的结果[22]。本研究选取与生物量实测值相关系数较高的15 个变量，采用DPS 7.05 软件的岭回归程序模块进行岭回归运算，根据岭迹图判断常数K 取0.2时，模型决定系数（R2）最大值为0.891 4，并且呈极显著性（P＜0.01）各回归系数的估计值基本稳定，得到岭回归模型如下：

表9 岭回归模型系数及显著性Table 9 Coefficient and significance of ridge regression model

2.2.2 主成分分析模型

主成分分析是一种通过降维技术把具有多重共线性的多个变量简化为几个主成分的统计方法[23]。在多变量研究中，由于变量个数过多且彼此之间存在一定的相关性，使得数据信息有一定的重叠，在高维空间中的表现较为杂乱。主成分分析采用降维的方法，构建几个相互独立的新变量来代表原来众多的变量，且尽可能保留原有的信息[24]。本研究使用MATLAB 6.5 对与生物量实测值相关系数较高的15 个变量进行主成分回归分析，计算因子载荷矩阵和得分矩阵，分析可得：前4 个主成分累计贡献率达到93.35%，已包含特征因子中绝大部分信息，其中第一主成分综合反映了除增强型植被指数以外的信息，第二主成分主要反映了波段3、增强型植被指数、归一化植被指数和差值植被指数的信息，第三主成分反映了波段2、波段3、波段4、波段7、波段4 提取的协同性、波段7 提取的协同性的信息，第四主成分反映了波段4 提取的方差信息、波段7 提取的协同性和信息熵的信息。

对所提取的主成分进行回归分析，得到标准化回归方程，再利用载荷矩阵将标准化回归方程转化为线性回归方程得到如下模型，模型相关系数和显著性见表11：

表10 主成分分析荷载矩阵、特征值及百分率Table 10 Load matrix,eigenvalue and percentage of principal component analysis

表11 主成分分析模型系数及显著性Table 11 Coefficient and significance of the principal component analysis model

2.2.3 偏最小二乘法

偏最小二乘法是一种将主成分分析、典型相关分析和多元线性回归分析相结合的数据统计方法[25]。其通过对自变量进行空间转换，提取对因变量解释性强的综合变量，辨别系统中的信息与噪声，克服变量间多重相关性在建模中的影响，使得估测值更接近于真值，相对于主成分分析法对因变量信息有更好的解释能力[26]。本研究选取与生物量实测值相关系数较高的15 个变量，基于交叉有效性原则确定所提取的成分数，参考Press统计量变化趋势图可以看出，提取5 个主成分时，模型效果最好，决定系数(R2)为0.482 7，对这5个主成分进行回归分析并进行线性转换，确定构建偏最小二乘回归模型如下：

表12 偏最小二乘法模型系数及显著性Table 12 Coefficient and significance of the principal component analysis model

2.3 模型精度验证

根据1.3.2 中的精度验证方法对3 种模型精度进行验证。各模型确定系数差异较大，其中岭回归模型最高，为0.918，其次为主成分分析模型（0.824），偏最小二乘法最小（0.626）；估计值的标准差以岭回归模型最佳（16.70 t/hm2），其次为主成分分析模型（20.19 t/hm2），偏最小二乘法模型标准误最大（44.77 t/hm2）；平均系统误差以岭回归模型最佳，为5.39%，其次为主成分分析模型最差，为8.37%，偏最小二乘法模型最差，为-13.79%；岭回归模型的精度最高，达到85.55%，其次为主成分分析模型，偏最小二乘法模型精度最低。将各模型测算出的预测值与实测值进行对比（图5），岭回归模型的预测值与实测值拟合R2最高，为0.927 8。综上，利用岭回归法建立的生物量反演模型各项指标最佳，最适合用于研究区天然乔木林地上生物量的估算。

表13 各模型精度指标对比Table 13 Comparison of accuracy indexes of each model

图5 预测值与实测值对比Fig.5 Comparison of the predicted and measured values

2.4 生物量空间分布格局

研究区单位生物量主要集中在30～200 t/hm2，单位生物量总体上呈现南部低北部高的趋势（图6）。研究区天然乔木林地上单位生物量随着海拔的升高呈现先增大后减小的趋势，在1 801～2 400 m 单位生物量最大（图7），说明这一海拔区域土壤养分、水分和温度更适宜乔木个体的生长；800～1 200 m 单位生物量最小，主要是这一海拔区域降水量比1 201～2 700 m 海拔区域有所减少，土壤水分不能充分满足天然乔木林的生长。在坡度上，单位生物量随着坡度的增加同样呈现先增大后减小的趋势，在缓坡、斜坡最大，险坡最小。主要是坡度影响了土壤养分的分布，坡度越大，土壤养分越容易随降水和地表径流流失，缓坡、斜坡的土壤养分能满足乔木个体生长，险坡土壤养分相对贫瘠不利于乔木生长。在坡向上，整体呈现北坡、西北坡、东北坡（阴坡半阴坡）单位生物量大于南坡、西南坡、东南坡（阳坡、半阳坡），北坡、西北坡单位生物量最大，其次为东北坡，南坡单位生物量最小。在阿尔泰山坡向主要影响降水，大西洋和北冰洋的冷湿气流在北坡、西北坡和东北坡抬升形成丰富降水，土壤水分能满足乔木个体生长；南坡为背风坡，降水量少而蒸发量大，土壤水分不能充分满足乔木个体生长。

图6 天然乔木林地上生物量等级分布示意Fig.6 Level distribution of aboveground biomass of natural forests

图7 天然乔木林地上单位生物量地形分布Fig.7 Topographic distribution of aboveground unit biomass of natural arbor forests

3 结论与讨论

3.1 讨论

3.1.1 生物量反演遥感影像的选择

从不同时点影像所提取的信息看，各期影像所提取的特征变量与实测生物量的相关性差异较大，其中以7月8日影像所提取的特征变量与实测生物量之间相关性最高，7月24日和9月12日次之，6月6日最低。这一现象是由阿勒泰山区的气候和优势树种所决定的。阿尔泰山区的优势树种为落叶松，是一种落叶乔木，据丁晓娟等[27]和孟军贵[28]的研究，落叶松生长与降水和温度具有显著相关关系，由于阿勒泰山属大陆性寒温带气候，7月具有最好的水热条件，是新叶生长和林冠发育的峰值，8—9月逐渐下降。而6月则因为新叶尚在发育，林相枯黄，光谱反映程度欠佳。这一结果与李根柱等[29]对辽东山区长白落叶松叶面积指数月度变化的研究结果相一致。

3.1.2 生物量反演特征变量的选择

从所提取的不同类别的特征变量看，以植被指数与生物量的相关性最大，经地形校正后的单波段数据次之，纹理信息整体较低。其中RVI 和NDVI 相关系数均超过0.85，但是NDVI 与生物量的相关性较RVI 低，这与部分学者[30]的研究不同，这可能是由于研究区植被覆盖度较高，对于高植被区RVI 较NDVI 具有更好的灵敏度，更适合用于生物量的估测。纹理信息数据庞大，与生物量在0.01 水平上显著相关的数量最多，但相关性数值较低，纹理信息数量多、相关性较低这一特点，在李明诗等[31]和吴平等[32]对于区域生物量的研究中也得以证实。纹理信息虽相关性较低，但本研究区天然乔木林针叶林和针阔混交林构成，相对于阔叶林尤其是杂阔类森林，树冠形状相对规整,空间可重复性较高，纹理特征能较为准确地描述森林的空间形态特征，因此在生物量模型的构建研究中选入纹理特征仍然是十分必要的。影响天然乔木林生物量的生境因子包含光照、温度、湿度、CO2浓度、土壤含水量、土壤养分等，如何通过遥感手段获取这些因子进行生物量反演将是后期研究的工作重点。

3.1.3 生物量反演模型的选择

从不同方法所构建的模型的精度看，通过岭回归所构建的生物量模型最适宜研究区。但在面对不同研究区的气候、地形条件和研究对象的差异时，最优的建模方法不尽相同，如吴明山[33]在回归模型在林业的应用研究中得出岭回归对于思茅松多元模型中共线性的消除最为有效合理，而罗庆辉[11]等在天山雪岭云杉生物量不同估测模型对比中，得出偏最小二乘法拥有最佳的预测精度，苟睿坤[22]等对陕西黄龙山林区油松人工林的反演模型对比研究中得到主成分分析得到的预测结果与实际值最为接近。因此在实际研究中具体哪一种方法可以取得最佳精度不能一概而论，需具体问题具体分析。

3.1.4 生物量空间分布格局

研究区天然乔木林地上单位生物量总体呈现出南低北高的格局，即布尔津分局整体高于阿勒泰分局。这一结果与郑栓丽等[34]关于阿尔泰山森林生态系统碳密度与碳储量研究中布尔津高于阿勒泰相一致（图6）。通过分析研究区森林资源二类调查数据可知，布尔津分局成熟林和过熟林所占比例高于阿勒泰分局，因此，布尔津分局单位生物量高于阿勒泰分局。从气候方面分析，水分和温度是影响山区天然乔木林生长的主要因子[12]，来自于大西洋、北冰洋的冷湿气流首先从西北到达布尔津分局，然后沿山体流向东南部的阿勒泰分局，使得布尔津分局区域内的降水量高于阿勒泰分局，进而使得布尔津分局的天然乔木林较阿勒泰分局生长更快，单位生物量更高。影响天然乔木林生物量空间分布的因素不仅包含地形因子，还有气候因子、土壤因子和人为活动等因素，下一步将综合各因素分析天然乔木林生物量的分布。