基于Landsat 的高山松地上生物量动态变化估测模型构建

廖 易 张加龙 鲍 瑞 许冬凡

（西南林业大学林学院，云南 昆明 650233）

森林生物量是森林生态系统最基本的数量特征[1]，它是森林生态系统长期代谢过程的积累，在一定程度上体现了森林的固碳能力。区域森林生物量的精确估测对该地区的森林资源经营管理有重要指示作用，对区域森林生态系统的研究和碳汇评价有重要意义[2]。森林地上生物量（AGB）是森林生物量的重要组成部分。当前，遥感技术已被广泛运用于森林AGB 的估测，成为森林AGB估测的主要方法[3]。

林分的变化会引起遥感图像光谱特征的变化，而光谱特征的变化对林分生物量值的改变有一定的指示作用。与光学影像相比，虽然激光雷达在捕捉森林结构方面有较大优势，但其采集成本高，计算要求高，难以获取大面积的数据[4]。此外，长时间序列的激光雷达数据有限，难以全面获取森林生长变化[5]。目前大面积、长时间的森林生物量监测系统主要采用Landsat 时间序列等多光谱卫星和实地调查数据[6]。Landsat 存档数据和预处理产品的开放，促进了其在森林监测活动中的研究[7]。当前，森林生物量估测方法应用较多的有线性回归[8]、随机森林[9]、支持向量机[10]等。随着基于单期影像估测森林生物量方法的成熟，出现了大量的估测生物量动态的研究，其中大多采用（定期）变化量建模[11]。

变化量建模相对于静态数据建模有一定的优势。陆驰[12]通过研究证明定期变化量模型提升了预测精度。张加龙等[13-14]对比静态和动态（即定期变化量）的森林生物量建模，也得出动态模型比静态模型效果更好。国外学者Gómez 等[15]将变化过程分为状态变量（静态）和过程变量（动态）构建模型，结果表明过程变量比状态变量对森林AGB 具有更强的预测能力。Puliti 等[16]使用双时相（AGB 变化量）和单时数据估测AGB 变化量，结果表明双时相方法可以得到更精确的AGB 变化值。Boisvenue 等[17]以树木的年龄和截距拟合变化量模型，并结合遥感信息的变化估测生物量，结果表明该模型能够更好的估计生物量及其变化。这些研究虽然都构建了动态变化模型，但其中大多只采用定期变化量建模，Gómez等虽然采用变化率进行研究，但并未比较定期变化量和变化率的效果，且这些研究的建模效果仍然有一定的提升空间。

故本研究拟基于1987—2017 年的森林资源连续清查样地和Landsat TM/OLI 数据，采用多元线性回归和随机森林建模方法，在年平均变化量和定期变化量建模研究的基础上加入变化率进行香格里拉市高山松AGB 模型构建，以期运用不同类型的遥感因子变化量，来提高高山松森林AGB 动态变化遥感估测的精度，为森林AGB 的动态估测研究提供技术支持。

1 研究区概况

研究区香格里拉市位于云南省西北部，地处北纬26°52′～28°52′，东经99°20′～100°19′，总面积为11 613 km2。该区域年平均气温为6.3 ℃，年平均降水量为651.1 mm，日照充足、太阳辐射强，属寒温带季风气候。境内水系发达，水资源丰富，全境河流均属金沙江水系。区域内地形起伏较大，最高点海拔5 545 m，最低点海拔1 503 m，平均海拔3 459 m。香格里拉地处“三江并流”的中心地带，有丰富的自然和生物资源，且风光独特，是西南地区重要自然资源储存地。境内有包括褐土、红壤和高山草甸土等在内的11 个森林土壤类型，森林覆盖率较高，达到75%[18]，高山松（Pinus densata）、云南松（Pinus yunnanensis）与云冷杉（Picea asperata）为优势树种[19]。

2 材料与方法

2.1 数据来源与处理

2.1.1 样地数据

研究所采用的样地数据为1987—2017 年共7 期的国家森林资源连续清查固定样地数据，各年份样地数量分别为19、22、23、16、16、17、23，共计136 块，其中包含部分复测样地。各年份样地胸径、树高统计结果见表1。样地大小为28.28 m × 28.28 m，面积为0.08 hm2，测树因子包含胸径、树高和株数等。

表1 样地测树因子统计情况Table 1 Statistical table of tree measuring factors in sample plots

2.1.2 遥感数据

研究获取的Landsat 影像数据共7 期21 景，从美国地质勘探局（http://glovis.usgs.gov/）网站获取。所有影像的获取年份与该地区连清高山松固定样地调查数据相对应，影像经过辐射定标和FLAASH 大气校正预处理。因研究区域内地形起伏较大，使用坡度匹配法进行地形校正。参照已有的SPOT-5 影像（北京1954 坐标系）对选取的Landsat 影像进行几何校正，误差控制在1 个像元内。在ENVI 5.1 软件中对各部分影像进行拼接，然后使用香格里拉市行政边界数据进行裁剪，分别得到各年份的香格里拉市影像。

2.2 遥感因子提取

提取5 种类型573 种遥感因子用于筛选自变量，详见表2。因子提取在ENVI 5.1 中进行，提取因子时以样地中心点为基础，落入Landsat 影像的对应位置并提取因子值。

表2 遥感因子Table 2 Remote sensing factors

2.3 高山松样地地上生物量计算

根据单木生物量模型[20]来计算高山松单木生物量（W），该模型为项目组研究成果，是利用实测的116 株高山松，选用胸径（D）、树高（H）作为模型的自变量建立，模型见公式（1）。使用连清数据中每块高山松样地的平均胸径、平均树高代入该单木模型计算平均单木生物量，然后乘以样地高山松株数得到样地生物量的近似值。

2.4 生物量和遥感因子变化量的计算

本次研究中，提出将定期变化量和变化率分别按照5 a 和10 a 进行研究。年平均和定期变化量计算时，剔除未进行重复观测和变化量为负值的数据，变化率数据采用7 个年份的共同样地计算。计算变化量后采用皮尔逊相关性法进行遥感因子筛选，每种变化量均选择10 个相关性最强且显著的遥感因子用于建模。计算得到的生物量变化量统计值见表3。

表3 地上生物量变化量统计表Table 3 Statistical table of variations of AGB

2.4.1 定期变化值

定期变化量是指从1987 年、1992 年、1997 年、2002 年、2007 年、2012 年、2017 年 的 样 地 数 据中选取2 个年份，首先筛选出这2 个年份的重复观测样地，用较大年份的AGB 和遥感因子值减去较小年份的对应值即得到变化量值。本研究计算了5 a 和10 a 的定期变化量，计算公式如下：

式中： Δb为定期变化量，b为高山松AGB 或对应的遥感因子值，n和m为用来计算变化量的2 个不同年份。

2.4.2 年平均变化值

年平均变化量包括5、10、15、20、25、30 a变化的年平均值，将定期变化量值比上对应的变化年份即可得到年平均变化量。公式如下：

式 中： Δg为年平均变化量，Δb为定期变化量，n和m为用来计算变化量的2 个不同年份。使用该公式计算AGB 及对应的遥感因子变化值。

2.4.3 变化率值

当 Δx→0时，Δy/Δx有极限，则称函数y=f(x)在点x0处可导，这个极限叫做f(x)在x0处的导数，这个导数即为该点处的瞬时变化率[21]，简称变化率，通常记为f′(x0)，表达式为：

将高山松AGB 取重复观测样地数据按年份进行排列，以年份和生物量分别作为横纵坐标进行拟合并生成曲线，对曲线每个年份处求导得到此处导数作为变化率，即表示生物量在该年份处5 a 或10 a 的瞬时变化速率。本研究中得到变化率建模数据的过程中需进行2 次相关性筛选，第1 次筛选出20 个与原始样地AGB 相关性最高的遥感因子，计算其与对应AGB 的变化率后，进行第2 次相关性筛选，得到10 个相关性最高的遥感因子变化率用于建模。

2.5 建模方法

本研究随机抽取80%的数据采用多元线性回归和随机森林方法建模，剩余20%数据用于检验。本研究中同一种变化类型中多元线性回归（MLR）与随机森林（RF）模型输入的建模因子相同。

2.5.1 多元线性回归

当多个自变量与因变量之间是线性关系时，可以进行多元线性回归分析。多元线性回归模型的基本公式为：

式中：β0为 常数项，β1～βn表示回归系数，ϵ为随机误差，xn表 示自变量，在本研究中代表遥感因子，y表示因变量，在本研究中代表AGB。使用方差膨胀因子（VIF）[22]来克服变量之间的共线性问题。

2.5.2 随机森林（RF）

随机森林[23]是一个集成学习模型，其中含有多个决策树，常被应用于多领域的数据分类和非参数回归[24]。本研究中随机森林模型算法的实现基于Python 语言“Sklearn”包中提供的“Random Forest Regressor”算法。建模过程中需对模型参数进行调参，以得到较好的模型结果。主要进行调整的参数为：模型的最大迭代次数n_estimators，决策树的最大深度max_depth，叶节点的最小样本数min_samples_leaf，内部节点再划分所需最小样本数min_samples_split。调参时根据样本量因素，min_samples_leaf 和min_samples_split设为默认值。

2.6 模型精度评价方法

本研究采用的精度评价指标为决定系数（R2）、均方根误差（RMSE）及相对均方根误差（rRMSE），计算公式如下：

式中：yi为真实值，为模型回归值，为真实值均值，n为样本个数。

3 结果与分析

3.1 建模因子分析

本研究通过皮尔逊相关性分析筛选用于建模的遥感因子，由于因子较多，每种变化量选出10 个因子用于建模。5 种变化量最终选出的遥感因子见表4。由表4 可知，采用RF 方法进行5 种类型的变化量建模，选入的50 个（每种变化量选入10 个）建模因子中48 个为纹理因子，其中SK 有16 个，出现频率高，HO、EN、DI 各有6个，ME、SM 各有4 个、VA 有3 个、CO 有2 个、CC 只有1 个，另有植被指数因子ND53 有2 个，表明纹理因子在生物量建模中具有较强的敏感性。

表4 5 种变化量的建模因子Table 4 Modeling factors of 5 variables

3.1.1 多元线性回归模型

采用筛选出的遥感因子变化量进行MLR 建模，得到各类型的MLR 模型结果见表5。3 类5 种变化量的MLR 模型均通过方差检验。由表5可知，最终选入模型的因子均为纹理因子，这表明纹理因子对MLR 的变化率建模较为敏感。

表5 5 种变化量的MLR 模型Table 5 MLR model with 5 type variations

3.1.2 随机森林模型

调用Python 语言“Sklearn”包中提供的“Random Forest Regressor”算法建模，得到RF的最佳模型参数见表6。随机森林会将输入的因子根据重要性进行排序，并输出各因子的建模贡献度值。本研究中5 种变化量类型的遥感因子贡献度（%）见图1。图中5 种变化量建模中贡献度最高的因子分别为R17B5ME、R11B7SK、ND53、R19B7HO 及R11B5SK，其中R19B7HO 在所有因子中贡献值最高，达50.26%。

表6 RF 模型最佳参数Table 6 Optimal parameters of RF model

图1 5 种变化量类型的RF 建模因子贡献度Fig. 1 Contribution degree of RF modeling factors for 5 variation types

3.2 模型精度评价

3 类变化量数据分别采用MLR 与RF 方法构建AGB 变化估测模型，得出结果的对比见表7，表中的模型编号及其变化量类型与表5、表6 相对应。

表7 3 类5 种变化量的MLR 与RF 建模精度评价Table 7 Accuracy evaluation of MLR and RF modeling based on 3 types and 5 kinds

由表7 对比所有模型的拟合效果，发现模型3、2、1、5、4、8、7、6、10、9 的R2依次增大，其中MLR 模型4 的R2最高（0.465），RF 模型9 的R2最高（0.956）。采用RF 建模，模型9 的拟合RMSE 最小，为0.664；采用MLR 建模，模型3 的RMSE 最小，为2.318。MLR 模型拟合的rRMSE 值最小为48.742%（模型4）；RF 模型拟合的rRMSE 值最小为35.883%（模型9）。以上结果表明模型9 的拟合效果最好，即RF 模型的5 a 变化率模型拟合效果最好。MLR 和RF 模型中，各变化量的建模拟合效果从劣到较优依次为年平均变化量、10 a 定期变化量、5 a 定期变化量、10 a 和5 a 变化率。表7 中，模型3 和模型8 的RMSE 在同种模型（MLR 或RF）中较小，这是变化的时间尺度造成的，模型3 和8 均为年平均变化量建模，表示地上生物量一年的平均变化量，而其他4 种变化量表示的是地上生物量至少5 年的变化，因此其变化数值相对较小，最终导致模型的RMSE 也较小。

模型预测检验结果与拟合效果对应，即模型3、2、1、5、4、8、7、6、10、9 的检验效果依次变优。RF 模型8、7、6、10、9 的RMSE 分别为2.082、20.547、17.244、4.422、2.285，MLR 模型3、2、1、5、4 的RMSE 分 别为2.342、21.491、21.277、5.301、4.865。预测效果的rRMSE 最小为42.718%（模型9）。表明5 a 变化率RF 建模的预测效果最好。

经对比，同一变化量类型的RF 建模效果优于MLR 建模效果，其中同种变化量类型MLR 与RF 模型的R2的最大差值为0.695，最小差值为0.491，且5 种变化量中5 a 变化率的RF 和MLR建模效果均为最优。综上，认为其符合一定规律，即3 类5 种变化量中5 a 变化率建模效果最好，年平均变化量最差，具体优劣排序为年平均变化量＜10 a 定期变化量＜5 a 定期变化量＜10 a 变化率＜5 a 变化率。

3.3 地上生物量反演

根据模型精度评价结果，采用最优的5 a 变化率RF 模型反演，先反演出各年份的变化率，再通过变化率公式反推得出AGB，结合已有的高山松分布面积范围数据[14]得出香格里拉市各年份变化率和AGB 空间分布图（图2、图3），高山松AGB 统计见表8。由表8 可知，1987—1997 年香格里拉市高山松AGB 总量逐渐增长，最高为929.02 万t，1997—2017 年AGB 总量逐渐减少，最低为872.25 万t。

图2 各年份变化率分布图Fig. 2 Distribution of change rate in each year

图3 各年份高山松AGB 分布图Fig. 3 AGB distribution of P. densata in each year

表8 各年份高山松AGB 反演结果Table 8 Result of the estimating AGB of P. densata in each year

4 结论与讨论

本研究以1987—2017 年每5 年1 期的香格里拉市高山松固定样地为基础，提取对应样地的遥感因子，经过处理得到3 类5 种变化量数据，分别为定期变化量（5、10 a 定期变化量）、年平均变化量、变化率（5、10 a 变化率），分别以这5 种变化量进行多元线性回归（MLR）和随机森林（RF）建模并比较各模型的精度。结果表明：

1）在同种变化量类型中，RF 模型的效果均优于对应的MLR 模型，其中RF 模型和MLR 模型的R2最大差值达0.705，RMSE 最多降低了11.236，可见非参数模型在遥感森林AGB 的动态变化估测中相较于参数模型具有一定的优势。

2）所有类型的变化量中，无论采用MLR 方法还是RF 方法建模，变化率（5、10 a 变化率）的模型效果均优于年平均变化量和定期变化量，其中5 a 变化率的RF 建模效果最好，其MLR 和RF 模 型 的R2分 别 为0.454 和0.956，RMSE 分 别为2.543 和0.664。采用变化率建模能够提高模型的精度。

3）纹理因子在本研究的5 类遥感因子中为出现频次最多的因子。本研究共筛选出相关性高的50 个建模因子，其中48 个为纹理因子，因此其具有最强的敏感性。在各种9 种纹理因子中，偏度（SK）出现次数最多，最具敏感性，相关性（CC）出现次数最少，最不具敏感性。

本研究在前人的基础上进一步提出采用年平均变化量与瞬时变化率建模，对比了3 种变化量（变化率、年平均变化量、定期变化量）模型的建模效果，得出结果表明瞬时变化率模型相较于年平均变化量和定期变化量模型在估测精度上有一定的提升，可以为地上生物量的精确估测提供一定的参考。

本研究样地数据时间跨度大，调查技术手段的改变也可能造成调查结果的不一致性，最终可能会对模型的评价结果产生一定的偏差。变化量模型精度较高，有利于区域AGB 的动态变化研究，在今后的研究中可以进一步挖掘更多类型的变化量用以提升模型精度。研究区香格里拉市地处高原地区且地形起伏明显，海拔、坡度、坡向等地形因素可能对AGB 有较大影响，但本研究暂未对该方面进行探究，今后的研究可以从此方面着手以探究地形对AGB 变化量模型的影响。