基于机器学习算法的九寨沟自然保护区植被生态水储量定量反演研究

周湘山，杨武年，罗 可，张宇航，文 艳，唐晓鹿

（1.中国电建集团成都勘测设计研究院有限公司，四川 成都 610072；2.成都理工大学，四川 成都 610059）

0 前 言

森林生态系统约占地球陆地表面积的33%，被称为“地球之肺”，是地球上最大的陆地生态系统，是全球生物圈中重要的一环，也是地球上重要的基因库、碳贮库、蓄水库和能源库，对维系整个地球的生态平衡有着不可替代的作用［1］。森林生态系统具有一项重要的生态服务功能就是涵养水源［2］，其涵养着陆地90% 以上的淡水资源［3］。林冠层是森林生态系统降水调节过程的第一个环节，也是森林生态系统中水分重新分配的第一个界面层［4］，林冠层对降水的调节作用主要表现为林冠截留和穿透雨两部分；其次是树干径流，是指降雨沿着树干流至树木根部的过程，树干径流量在一次降雨中所占的比率比较低，一般低于5%；枯落物层作为森林生态系统的第二个水文作用层，主要功能是增加土壤的有机质含量、改善土壤的理化性质，提高土壤孔隙率，防止土壤严重板结、增强土壤入渗性能，减少土壤的无效蒸发［5］，使更多的水分能够储存在土壤中供植物生长；土壤层是森林生态系统的第三个水文作用层，在整个森林生态系统中水源涵养能力最强，不仅能够调节降水资源的分配，也是生态系统水分和元素循环的主要储蓄库，是评价森林生态系统水源涵养能力的重要指标。研究不同森林植被类型储水量及其水文功能，是生态修复、重建过程中的重要问题，也是需要迫切解决的基础理论问题。

目前已有研究成果中，基本采用单一遥感影像，通过各类植被指数、纹理指数等地面参数与野外实测数据建立回归关系反演或估算植被含水量。部分研究增加了LIDAR、INSAR、高光谱等多源遥感融合反演的成果，但利用机器学习算法建立回归模型的应用较少。常用的机器学习算法包括：随机森林、支持向量机、人工神经网络等。机器学习算法较以往的线性参数模型有更好的拟合精度，对植被含水量数据需要满足的统计条件较低，能够有效地解决高维数据（变量个数较多、变量类型多样）在建立模型过程中所带来的变量选择困难。随着机器学习算法的发展及森林植被生态水相关遥感研究的深入，利用机器学习算法对森林植被进行分类、植被生态水反演模型构建，将进一步提高植被生态水储量的估算精度和准确性。

1 研究区概况

研究区位于九寨沟漳扎镇，漳扎镇位于九寨沟县南部，原名九寨镇，地处N33°02′N～33°21′N，108°38′E～104°03′E，位于长江系嘉陵江上游白水江源头。地势整体呈南高北低，地貌以高山为主，海拔在1 800～4 800 m之间；气候温和，年平均气温为12.7℃。春季气温较低且温差较大，平均气温位于9～18℃；夏季较凉爽；秋季气候宜人但昼夜温差较大；冬季较寒冷，日照充足，年平均日照约1 600 h。保护区内土壤大部分以半湿润、半干旱暖的温带基带山地淋溶褐土为代表，林下土壤则多为山地棕壤和山地暗棕壤，腐殖质含量较高。保护区内野生动植物资源非常丰富，已发现的野生动物有600余种，其中属于国家重点保护动物有20多种。此外，保护区还是四川省的第二大林区，森林覆盖率超过70%，有四川红杉、白皮杉等多种国家级珍稀保护植物。

2 数据采集及处理

2.1 实测数据分析

2.1.1 野外数据采集

由于九寨沟自然保护区范围的气候条件复杂，云层覆盖较多，考虑不同季节植被的含水量会有较大变化，而不同年份的相同季节植被含水量较之变化不大的因素，尽量保持下载的Sentinel-2B和Landsat8 OLI影像成像时间与所采集样点实测数据的季节相契合，布点方式考虑了植被类型多样和交通便利两种因素。野外工作由A、B两个小组一起完成，共采集54个样地数据。样地采用圆形样地方式，以一点为圆心，10 m距离为半径，对圆内所有胸径大于5 cm的树木进行单木检尺，主要记录样地编号、树种名称、1.3 m高度处胸径和树高。其中针叶林类型14个，其主要树种有油松、云杉、马尾松、冷杉、云杉等；阔叶林类型16个，其主要树种有栎树、桦树、杨树等；混交林类型24个。

2.1.2 生物量计算

本文采用模型法对样地林木生物量进行估算，在样地内测得的数据有胸径和树高，利用已有的生物量方程计算样地内单棵树木干枝叶各部位的生物量。

根据野外采集的胸径D计算得出单棵林木的枝、干、叶生物量（干重），分别为BM枝、BM叶、BM干。由野外采集的林木的枝、干、叶少量样本，可以得到分别对应器官的样本鲜重FM枝、FM叶、FM干，带回室内后将样品放置烘箱中，将叶类样本在70℃条件下烘干，枝干部分样本在80℃条件下烘干，反复称量至恒重，分别得到对应林木器官的干重DM枝、DM叶、DM干。

据生物量（干重）和含水率的关系，可以计算出林木的枝、干、叶植被含水量，并将森林植物样地的含水量统一换算到标准单位吨/公顷（Mg/hm2）。

2.1.3 研究区植被含水量

按照54个样地的实测数据，根据各类树种生物量模型计算出各类树种的干、枝、叶等部位的生物量，通过含水率和含水量计算出各类样地所有树种干、枝、叶和总的含水量。

2.2 遥感影像分类

本文主要提取研究区植被林种分布及面积，选择计算机“解译为主，目视解译为辅”的信息提取方法。将研究区地物类型分为针叶林、阔叶林、针阔混交林、草地、灌木、水域、裸地和雪地，如图1所示。对分类结果进行混淆矩阵精度验证，其中总分类精度为81.73%，Kappa系数为0.795 6，满足分类要求。

图1 研究区土地利用类型分类

2.3 特征波段选取

本文结合植被光谱特性和Sentinel-2B、Landsat8OLI卫星传感器参数参数，选取了Sentinel-2B影像 中 的 Band1、Band2、Band3、Band4、Band5、Band6、Band7、Band8、Band8b、Band9、Band11、Band12和Landsat8 OLI卫星影像中的Band1、Band2、Band3、Band4、Band5、Band6、Band7、Band9等原始波段地表反射率作为可选的特征因子。

2.4 植被指数提取

本研究主要选取了基于代表绿度的垂直植被指数（Perpendicular Vegetation Index）、三角植被指数（Triangle Vegetation Index）、差值 环 境植被指数（Difference Vegetation Index）、归一化绿蓝差异指数（Normalized Green-Blue Difference Index）、绿蓝比值指数（Green-Blue Ratio Index）、可见光抗大气指数（Atmospherically Resistant Vegetation Index）、绿红比值指数（Green-Red Ratio Index）、有效叶面积指数（Specific Leaf Area Vegetation Index）、归一化水分指数（Normalized Different Moisture Index）、归一化红外指数（Normalized Difference Infrared Index）等10种植被指数作为遥感特征因子。

2.5 纹理指数提取

Haralick（1979）定义了8种常用的纹理特征，即：均值、方差、协同性、对比度、相异性、熵、二阶矩、相关性［6］。

3 模型算法

3.1 XGBoost

XGBoost算法是由陈天奇等人［7］提出来的，与一般的梯度boost方法相比，XGBoost对目标函数进行了二阶泰勒展开，并在训练过程中利用二阶导数加快了模型的收敛速度。其核心思想就是不断地添加树，不断地进行特征分裂来生长一棵树，每次添加一个树，其实是学习一个新函数f（x），去拟合上次预测的残差。最后只需要将每棵树对应的分数加起来就是该样本的预测值。

3.2 MARS

多元自适应回归样条（Multivariate Adaptive Regression Splines，MARS）是一种数据分析方法。该方法以样条函数的张量积作为基函数，分为前向过程、后向剪枝过程与模型选取三个步骤。其优势在于能够处理数据量大、维度高的数据，而且计算快捷、模型精确。

3.3 随机森林

随机森林算法是通过集成学习的思想将多棵决策树集成的一种算法，它的基础组成单元是决策树，而随机森林可以看作是若干棵决策树的集成。基本组成单元采用CART算法［8］，其本质属于机器学习的一大分支——集成学习（Ensemble Learning）方法，其工作原理是生成多个分类器或者模型，各自独立地学习和做出预测。

4 模型精度评价

本研究选取交叉验证（Cross Validation，CV）的方式来检查模型的精度，采用决定系数（R-squared，R2）和均方根误差（Root Mean Square Error，RMSE）进行评价。

4.1 Sentinel-2B

采用三种机器学习方法建立生态水反演模型，并通过十倍交叉验证法进行精度验证。其中XGBoost模型决定系数R2为0.45，均方根误差RMSE为54；MARS模型决定系数R2为0.37，均方根误差RMSE为57；随机森林模型决定系数R2为0.40，均方根误差RMSE为56。

4.2 Landsat8 OLI

采用三种机器学习方法建立生态水反演模型，并通过十倍交叉验证法进行精度验证。其中XGBoost模型决定系数R2为0.38，均方根误差RMSE为57；MARS模型决定系数R2为0.35，均方根误差RMSE为59；随机森林模型决定系数R2为0.32，均方根误差RMSE为60。

4.3 联合反演

采用三种机器学习方法建立生态水反演模型，并通过十倍交叉验证法进行精度验证。其中XGBoost模型决定系数R2为0.57，均方根误差RMSE为48；MARS模型决定系数R2为0.46，均方根误差RMSE为53；随机森林决定系数R2为0.43，均方根误差RMSE为55。

4.4 模型评价

根据分析结果得出，通过Sentinel-2B影像和Landsat8 OLI影像分别提取的54个样本点的各类特征波段、植被指数和纹理指数指标进行重要性分析，选取相对重要性和节点纯度高的指标带入XGBoost模型、MARS和RF三种模型进行相关性计算，发现提取指标联合反演的相关性最高，Sentinel-2B影像提取的指标单独反演的相关性其次，Landsat8 OLI影像提取的指标单独反演的相关性最低。同类影像提取指标中XGBoost反演的相关性最高，MARS反演的相关性其次，RF反演的相关性最低。

5 植被生态水定量反演

利用森林生物量方程计算了54个样地中的不同树种的地上生物量，通过不同树种地上生物量和含水率数据计算得到54个样地植被生态水含量。对植被生态水含量反演方法进行研究，选择特征波段、植被指数和纹理指数进行植被生态水含量反演，利用所建立的XGBoost模型植被生态水反演模型进行回归分析，对研究区2019年12月份的植被生态水含量进行反演，反演结果如图2所示。

图2 研究区生态水储量分布

根据本研究区的土地利用类型分类结果，研究区植被生态水储量主要分布在针叶林、阔叶林、混交林和灌木林等四种森林类型，其中针叶林生态水储量为351.94×104Mg，平均值210.33 Mg/hm2；阔叶林生态水储量132.52×104Mg，平均值127.88 Mg/hm2；针阔混交林生态水储量14.48×104Mg，平均值为132.52 Mg/hm2；灌木林生态水储量为25.02×104Mg，平均值为150.54 Mg/hm2。研究区植被生态水储量为523.97×104Mg，平均值为175.53 Mg/hm2。由于本研究基于的是光学遥感影像的定量反演，不具备LIDAR和INSAR等影像数据的穿透性，无法获取枯落物和土壤层的生态水含量。同时基于生物量公式及含水率进行计算，因此研究区植被生态水储量主要代表各森林类型的枝、干、叶等器官的地上植被生态水含量。

6 结论及展望

6.1 结 论

（1）根据特征选择结果，在同一Sentinel-2B变量三种模型下的重要性分析中，纹理指数指标占比最多，其次为特征波段，植被指数占比最低。在同一Landsat8 OLI变量三种模型下的重要性分析中，纹理指数指标占比最多，其次为植被指数。在同一联合变量三种模型下的重要性分析中，纹理指数指标占比最多，其次为特征波段，植被指数占比最低。可以看出三种模型下纹理指数指标的重要性占比均最多，是九寨沟自然保护区及川西高原生态水定量反演的重要参数指标。

（2）根据模型反演结果，同类影像提取指标中XGBoost模型反演的相关性最高，指标构建中联合反演的相关性最高。在XGBoost模型中，结合Landsat8 OLI和Sentinel-2B影像作为预测变量提供了最好的生态水储量估计。通过不同传感器数据集和建模算法的协同作用，为基于遥感的九寨沟自然保护区及川西高原生态水定量反演提供了一种新的方法。

（3）根据植被生态水含量定量反演结果，按照不同植被类型生态水储量总量占比，从大到小依次为针叶林、阔叶林、灌木林和针阔混交林。其中针叶林分布面积最大；生态水储量平均值最大；生态水储量总量最大，占比高达67.17%。

6.2 展 望

（1）对于植被生态水的时序监测，应该以不同时间（季节/年份）、不同植被类型、固定样地的野外数据作分析研究，既简化了野外工作，又方便建立研究区植被生态水数据库。

（2）分析地形因子对植被生态水分布的影响。经野外采集数据时发现不同海拔、坡度条件下，植被生长状况和植被类型都不完全一致，说明海拔、坡度可能是影响植被生态水含量的因素，具体相关性有待进一步研究。

（3）本研究中植被生态水反演采用了光学影像提取地面反演参量建立模型，未采用高光谱影像、SAR影像及其结合的研究。由于高光学影像对植被表层信息捕获能力更强，SAR影像的后向散射信息还包含了枝干之间更为复杂的回波信号，可以作为植被垂直空间信息的补充。后续可以采用高光谱影像、SAR影像开展相关性反演工作，进一步提升植被生态水反演相关性及精度。