1 引言
干旱是一种复杂的、 多属性的自然现象， 常表现出频率高、 持续时间长、 影响范围广等特点， 对国民经济尤其是农业生产造成严重影响（王兴等， 2023）。研究表明， 气象灾害造成的经济损失占自然灾害损失的八成以上， 而干旱造成的损失占到气象灾害损失的一半以上（Hao et al， 2015； 张强等， 2022）。干旱事件造成如此巨大损失的主要原因是人们对干旱发展行为认识不足， 缺乏及时准备和有效的缓解措施（Yan et al， 2018）， 因此， 及时、 准确地监测干旱情况， 对于抗旱准备和减少风险至关重要（Sheffield et al， 2014； 宋艳玲等， 2019； 谢五三等， 2021）。
遥感技术具有宏观、 快速、 经济等特点， 能弥补站点监测在空间尺度上的不足， 是实现区域或全球尺度旱情监测的重要手段（Rhee et al， 2010； Zhu et al， 2017）。目前， 国内外学者从不同角度提出了诸多干旱遥感监测模型， 为实现国家综合减灾与农业旱情监测奠定了基础（黄友昕等， 2015； 董婷等， 2015）。干旱常伴随着土壤水分、 植被水分、 冠层温度等多种地表特征的变化（江笑薇等， 2019； 刘二华等， 2020）， 不同地表特征变化会形成不同的遥感信号， 使得不同遥感干旱指数在表征干旱时存在一定的差异（孙灏等， 2012）。因而， 目前干旱监测研究更加趋向于多源信息融合， 以弥补数据自身的局限性和单个干旱指数在区域旱情监测中的不足（米前川等， 2022）。孙灏等（2012）比较分析了13种典型干旱遥感监测指数， 按指数间的相关特征， 将其划分为土壤水分变化、 冠层温度变化、 植被水分变化和作物形态及绿度变化四大类， 建立了遥感干旱指数分类体系。黄友昕等（2015）重点阐述了四类遥感干旱监测指标自身的适用范围， 并在此基础上指出今后干旱遥感指标适用性研究中需要重点解决的6个关键问题。温庆志等（2019）选择四类指标中最能反映淮河流域旱情特征的指数， 通过线性回归模型， 构建了适用于淮河流域的综合遥感干旱监测模型， 并分析了淮河流域干旱的时空分布规律。杜灵通等（2014）在孙灏等提出的干旱遥感监测分类体系基础上， 采用分类回归树构建了多因子综合干旱监测模型。多致灾因子的协同分析， 切实可行地为评估区域干旱的时空分布特征提供了一种方法。
随机森林算法是一种基于分类与回归树的机器学习方法（Breiman et al， 2001）， 具有学习速度快、 运算速度快、 在处理大数据集上高效、 精度高且不易过拟合等优点。与传统线性回归、 分类回归树和神经网络算法相比， 随机森林算法更加稳定， 对噪声和异常值也有很好的容忍性。因此， 本文选择能够反映冠层温度变化、 植被绿度变化、 土壤水分变化和植被水分变化的最优遥感指数， 采用随机森林回归模型， 构建了综合遥感干旱监测模型， 并将该模型应用于西南地区旱情监测业务中， 以评估RFR模型在丘陵山地区域旱情监测的适用性。
2 研究区概况
研究区位于我国西南部， 包含四川省、 云南省、 贵州省和重庆市三省一市， 地理位置介于97°27′E -110°18′E和20°53′N -34°14′N， 属亚热带季风气候和高原山地气候（孙昭萱等， 2022）。区域内地势起伏较大， 呈西高东低， 北高南低的分布趋势， 海拔介于16~7143 m， 平均海拔为1987 m（图1 所示）。地形地貌较为复杂， 由西向东包含了川西高原高山地区、 云贵高原中高山山地丘陵区、 四川盆地及周边山地。受孟加拉湾水汽、 南支槽等影响， 西南地区极易发生极端高温干旱， 是我国旱情高发风险区之一（贾何佳等， 2022）。图1 中“+”和“·”分别为用于构建模型的训练样本和检验模型精度的测试样本站点， 其中训练样本站点个数为173个（占总站点数的70%）， 测试样本站点74个（占总站点数的30%）， 训练样本和测试样本分散分布在研究区内， 但受西南地区地形条件等因素影响， 研究区西北部（川西高原）可用样本站点数明显少于东部地区。文中涉及的地图是基于中华人民共和国自然资源部地图技术审查中心标准地图服务系统下载的审图号为GS（2023）2765号的中国地图制作， 底图无修改。
3 数据来源与方法介绍
3.1 研究数据及处理
3.1.1 遥感数据及处理
为构建遥感干旱监测模型， 本文参照孙灏等（2012）提出的干旱遥感监测指标分类体系， 将15种典型的遥感干旱监测指数划分为4大类（表1 所示）。计算15种遥感干旱监测指数所需的遥感数据主要包括： 2002 -2022年MODIS逐日的地表温度产品（MOD11A1）、 时间分辨率为16 d的MODIS植被指数产品（MOD13A1）， 及2022年MODIS逐日的地表反射率产品（MOD09GA）。其中， MOD11A1是Wan （2014） 采用广义分裂窗算法反演的第六代MODIS LST产品数据， 其空间分辨率为1 km。Wan （2014） 利用北非5个裸土站点数据采用R-based方法对MODIS LST进行验证分析， 结果表明12组验证数据中10组数据LST平均误差在0.6 K以内， 1组数据平均误差小于0.8 K， 另外一组数据LST平均误差小于1.9 K。MOD13A3为MODIS陆地产品组开发的月最大合成植被指数产品， 空间分辨率为1 km， 数据集包含了NDVI、 EVI及各波段反射率等信息。MOD09GA是NASA提供的时间分辨率为1 d， 空间分辨率为500 m， 包含MODIS 7个波段的地表反射率产品。此外， 还包含一个数据质量控制文件， 记录了像元云覆盖、 云阴影、 气溶胶厚度等信息。以上MODIS产品数据均来源于NASA官方网站（http： //ladsweb.nascom.nasa.gov）， 产品存储格式为HDF5， 数据行列号为h26v05、 h26v06、 h27v05和h27v06。本文采用MODIS Reprojection Tool（MRT）工具对MODIS地表反射率、 归一化植被指数和地表温度产品进行批量处理， 统一产品空间分辨率为1 km， 投影坐标为Albers等积圆锥投影， 处理后的文件保存为GeoTIFF格式。
表1 典型干旱遥感监测指数及分类体系Table 1 Classification system of drought remote sensed monitoring index |
| 分类体系 | 典型干旱遥感指数 | 计算公式 | 数据来源 |
|---|---|---|---|
| 冠层温度变化 Canopy temperature change | 温度条件指数TCI | （尹国应等， 2022） | MOD11A1 |
| 温度植被干旱 指数TVDI | （吴天晓等， 2023）
| MOD11A1 MOD13A1 | |
| 植被供水指数VSWI | （宋扬等， 2017） | MOD11A1 MOD09GA | |
| 条件植被温度 指数VTCI | （王蕾等， 2018） | MOD11A1 MOD09GA | |
| 作物形态及绿度变化 Crop form and greenness change | 距平植被指数AVI | （玉院和等， 2018） | MOD13A1 |
| 增强植被指数EVI | （王瑾等， 2022） | MOD09GA | |
| 比值植被指数RVI | （朱高龙等， 2013） | MOD09GA | |
| 归一化植被 指数NDVI | （朱高龙等， 2013） | MOD09GA | |
| 条件植被指数VCI | （张瑾等， 2022） | MOD13A1 | |
| 土壤水分变化 Soil moisture change | 改进的垂直干旱 指数MPDI | （郭浩等， 2017）
| MOD09GA |
| 垂直干旱指数PDI | （陈国茜等， 2014） | MOD09GA | |
| 可见光和短波红外 干旱指数VSDI | （李梦云等， 2016） | MOD09GA | |
| 植被水分变化 Vegetation water content | 全球植被水分 指数GVMI | （王溥等， 2010） | MOD09GA |
| 地表水分指数LSWI | （刘元亮等， 2015） | MOD09GA | |
| 短波红外垂直失水 指数SPSI | （朱琳等， 2011） | MOD09GA |
3.1.2 气象数据及处理
本文采用气象干旱综合指数（MCI）作为参考数据， 用于构建和验证遥感干旱监测模型。MCI指数考虑了60天内有效降水（权重累积降水）、 30天内蒸散（相对湿润度）及90天（季度）降水和150天（半年）降水的综合影响。MCI计算如公式（1） 所示。本文选用了2021 -2022年研究区内247个国家级气象站（图1 ）逐日的最低气温、 最高气温、 平均气温、 风速、 降水数据、 土壤相对湿度等数据， 使用1971 -2022年逐日降水数据计算气候平均值， 数据经过质量控制， 去除数据集中的异常值。其中， 173个站点用于模型训练（占总站点数的70%）， 74个站点用于模型精度评估（占总站点数的30%）， 且保证训练样本和测试样本分散分布在研究区内。
式中： 为季节调节系数， 根据不同季节各地主要农作物生长发育阶段对土壤水分的敏感程度确定， 取值可参考GB/T 20481 -2017标准（2018）； 为近60天标准化权重降水指数， 计算方式见公式（2） 所示：
式中： 为权重累积降水量（单位： mm）， 其计算公式如公式（3） ； N为某一时段的长度； 为贡献参数， 当N为60天时取值为0.85； 为距离当天前第n天的降水量（单位： mm）； 和 分别为近90天和近150天标准化降水指数， 王天等（2022）在其文章中给出了详细的计算方式； 为近30天的相对湿润指数， 其计算公式如公式（4） ； 、 、 和 为权重系数， 其值根据区域差异在GB/T 20481-2017标准（2018）查询得到。杨歆雨等（2022）将MCI划分为MCI>-0.5， 为无旱； ， 为轻旱； ， 为中旱； ， 为重旱； ， 为特旱。
式中： 为某时段相对湿润指数； 为某时段的降水量（单位： mm）； 为某时段的潜在蒸散量， 本文参考GB/T 20481-2017标准（2018）中的Penman-Monteith公式计算得到。
3.2 研究方法
随机森林回归模型（RFR）是由Breiman（2001）提出的一种基于决策树的集成学习算法， 其通过Bootstrap取样法有放回地随机抽取样本和特征， 构建n棵互不关联的决策树， 采用并行方式获得预测结果。每棵决策树都能通过抽取的样本和特征得出一个预测结果， 通过综合所有树的结果取平均值， 得到整个森林的回归预测结果， 这使得RFR拥有更好的预测能力， 对噪声和异常值有很好的容忍度， 能有效地避免过拟合现象（沈润平等， 2016）。
本文利用Python语言构建RFR模型， 来评估干旱与各遥感指标之间的定量关系。如图2 所示， 通过比较MCI指数与四类干旱遥感监测指数的相关性， 选择各类指标中最能反映研究区旱情特征的指数为自变量， MCI为因变量， 利用RFR模型建立遥感干旱监测模型。由于RFR模型预测性能受决策树个数、 决策树最大深度、 最小分离样本数、 最小叶子节点样本、 最大分离特征数， 以及是否进行随机抽样等参数影响（孙兴亮等， 2022）， 因而本文在进行RFR建模时需对以上几个超参数进行调优， 确定各个超参数的最优取值。
4 结果与分析
4.1 遥感干旱指数相关性分析
遥感干旱指数间的相关性表明了其对干旱描述的一致性。图3 为不同遥感干旱指数间的相关系数热力图。从图3 中可以看出， 不同遥感干旱指数在表征干旱特征时存在一定的差异性。部分遥感干旱指数间存在显著的相关性， 但它们与另一部分遥感干旱指数间相关性却较差。例如AVI、 RVI、 NDVI和VCI指数间存在较好的相关性， 其相关系数均优于0.5， 甚至AVI与VCI， NDVI与EVI， 及RVI与NDVI间的相关系数达到了0.90， 0.90和0.99， 但该类指数与VSDI、 SPSI等指数间的相关系数较小。如AVI、 RVI、 NDVI、 VCI与SPSI指数间的相关系数分别为0.03、 0.17、 0.19和0.11。TVDI与TCI和VTCI的相关系数分别为0.29和1.0， 与其他干旱指数间相关系数均小于0.3。LSWI与GVMI、 SPSI间相关系数分别为0.98和0.46， 与其他干旱指数间相关系数均小于0.3（除VSDI外）。以上统计结果较好地说明了遥感干旱指数间存在明显的聚类效应， 类间指数旱情表征一致性较好， 不同类指数间一致性较差（孙灏等， 2012； 温庆志等， 2019）。因而， 在构建综合遥感干旱监测模型时， 结合不同类别遥感干旱指数优势的同时， 还需要考虑遥感干旱指数间的自相关性， 尤其是同类型遥感干旱指数的选择。本文按照不同遥感干旱指数间相关性矩阵和孙灏等（2012）的研究结果， 将15种遥感干旱指数划分为冠层温度变化、 作物形态及绿度变化、 土壤水分变化和植被水分变化四大类， 并从每一类中选择最能反映研究区旱情特征的指数作为RFR综合遥感干旱监测模型的自变量。
4.2 RFR模型建模因子选取分析
箱型图是一种用作显示一组数据分散情况的统计图， 可以直观地探索数据的分布特征。为从四大类干旱遥感监测指数中筛选出每类指标的最优指数， 本文分析了2022年15类干旱遥感监测指数与MCI指数的相关性（图4 ）。冠层温度变化类， TCI、 TVDI、 VSWI和VTCI与MCI指数间相关系数R值的中位数介于0.3~0.4［图4 （a）］， 相比于其他3个指数， TVDI与MCI之间的R值中位数最大， 且R值主要分布在0.3~0.7， 表明TVDI与MCI指数的相关性优于其他三个指数； 从图4 （b）可知， EVI、 RVI与MCI指数的中位数接近且大于其他三个指数， 表明EVI、 RVI与MCI指数的相关性优于其余三个指数， 而对比EVI与RVI结果可知， 尽管EVI中位数要优于RVI， 但EVI上下边缘范围较大， 其相关系数的方差和标准偏差也较大， 综合考虑选择RVI指数作为作物形态与绿度变化类自变量； 土壤水分变化类3个指数与MCI指数间R值均较低， 但与MPDI和VSDI指数相比， PDI与MCI指数间R值分布相对集中， R值间差异相对较小［图4 （c）］， 表明PDI与MCI指数间的相关性要优于MPDI和VSDI； 作物水分变化类结果［图4 （d）］显示， GVMI、 LSWI与MCI指数的中位数接近且大于SPSI， 但与LSWI相比， GVMI箱型图上下边缘距离较短， R值分布相对集中， 说明数据整体方差和标准差相对较小， 在长时间序列上， GVMI与MCI指数的相关性要优于LSWI。综合以上分析， 本文选择TVDI、 RVI、 PDI和GVMI指数分别代表冠层温度变化类、 作物形态与绿度变化类、 土壤水分变化类和作物水分变化类的最优遥感干旱监测指标， 用于输入RFR模型中构建适用于西南地区的遥感干旱监测模型。
4.3 RFR模型构建及精度分析
本文选取TVDI、 RVI、 PDI和GVMI指数作为自变量， MCI指数作为因变量， 分别采用OLS和RFR模型构建遥感干旱监测模型。在RFR模型构建过程中， 选择一组最优的超参数， 不仅能提高模型的学习性能， 还能一定程度上影响模型预测的精度。因此， 本文首先对RFR模型的超参数（决策树个数与最大深度、 最小分离样本数、 最小叶子节点样本、 最大分离特征数， 及是否进行随机抽样）的范围加以确定， 然后在这些范围内确定各个超参数的最终最优取值， 以实现RFR模型超参数的自动调优。
对比图5 （c）和图5 （d）可知， OLS模型构建的综合干旱监测模型与MCI之间的R值为0.76， Bias值为0.02， RMSE值为0.79， 其P值小于0.01， 表明基于OLS模型构建的综合干旱监测指数能较精确地表征区域内旱情状况， 其相关系数通过了显著性检验。RFR模型综合干旱监测指数与MCI指数间的R值为0.90， Bias值为0.02， RMSE值为0.53， P值小于0.01。较高的R值和较低的Bias、 RMSE值表明， RFR模型构建的遥感干旱监测指数在表征区域旱情时比OLS模型更加优秀。此外， 对比图5 （b）和图5 （d）可知， RFR模型训练集拟合精度（R为0.95， RMSE值为0.39）和测试集的预测精度（R为0.94， RMSE为0.41）差异较小， 较好地体现了RFR模型在预测过程中避免过度拟合的能力。
4.4 RFR综合干旱监测模型在西南地区的应用
为检验RFR模型构建的干旱监测模型的应用能力， 本文以西南地区173个气象站观测资料估算的MCI指数为因变量， MODIS数据估算的TVDI、 RVI、 PDI和GVMI指数为自变量， 采用RFR模型定量分析了2022年西南地区旱情时空变化特征（图6 ）， 并将估算结果与MCI插值结果进行对比（图7 ）。对比结果显示， 除7月外， 其余月份MCI指数与遥感估算的旱情在时间和空间分布上较为一致。1 -6月， 西南地区受旱情影响区域较小， 重庆北部、 四川与贵州交界处等局部地区出现轻旱。进入7月， MCI监测结果显示， 重庆东北部、 四川北部地区出现旱情， 干旱等级以中旱为主， 局地有重旱。遥感监测结果显示， 旱情主要以中旱为主， 旱情分布在四川西部和贵州等地。遥感监测结果与MCI指数监测结果差异较大； 8月， 受持续高温天气影响， 重庆大部分地区、 四川东部及四川北部均出现了中旱和重旱， 局地出现特旱。遥感反演结果监测到了四川盆地以东的旱情， 川西地区监测结果与MCI指数插值结果存在较大的差异； 进入9月， 重旱等级以上区域逐步扩大， 四川东部、 重庆中西部和广西北部均受重旱的影响； 10月， 局部地区出现降水， 重庆东北部和四川东部旱情得到缓解， 旱情集中在广西东部地区； 11 -12月， 受长时间无降水的影响， 重庆西南部、 四川南部和广西地区受重旱影响。综合以上结论可知， 受川西地区气象站点个数较少， 站点分布不均， 地形复杂等因素影响（图1 ）， MCI指数与RFR模型估算的遥感干旱监测结果在川西地区存在明显的差异； 在其他区域遥感干旱监测结果与MCI指数时间变化和空间分布较为一致， 能较好地表征区域旱情的时空变化过程， 为防旱抗旱工作提供重要的数据支撑。
图6 2022年基于RFR模型的西南地区遥感干旱监测时空分布图Fig.6 The Spatiotemporal distribution map of drought monitoring in Southwest China based on RFR model in 2022 |
5 结论与讨论
干旱是一种复杂的过程， 传统遥感干旱指数侧重于对单一响应因子进行监测与评估， 缺乏对干旱的综合分析。本文对比分析不同遥感干旱指数与MCI指数的相关性， 筛选出四类遥感干旱监测指数中与MCI指数相关性最好的指数作为自变量， MCI指数作为因变量， 利用RFR模型构建了遥感干旱监测模型。得到以下结论：
（1） 不同遥感干旱指数在表征干旱特征时存在一定的差异性。遥感干旱指数间存在明显的聚类效应， 类间指数旱情表征一致性较好， 不同类间一致性较差。
（2） 遥感干旱监测指数与MCI指数间相关性结果表明， 各类指标中TVDI、 RVI、 PDI和GVMI指数与MCI指数间的相关性较好， 作为自变量构建遥感干旱监测模型。
（3） OLS和RFR模型精度对比显示， RFR模型构建的遥感干旱监测指数在表征区域旱情时比OLS模型更加优秀。OLS模型训练集精度R值为0.78， RMSE值为0.73； RFR模型训练精度R值为0.97， RMSE值为0.33。测试集结果显示， OLS模型R值为0.76， RMSE值为0.79； RFR模型R值为0.90， RMSE值为0.53。
（4） 受川西地区气点个数较少， 站点分布不均， 加之地形复杂等因素影响， MCI指数与遥感干旱监测结果存在明显的差异； 在站点较多区域遥感监测结果与MCI指数时空变化较为一致。
相对于传统的遥感干旱监测模型， 本文创新之处在于选择四类遥感干旱监测指数中最能反映研究区旱情特征的指数作为自变量， MCI指数作为因变量， 利用RFR模型构建遥感干旱监测模型， 用于2022年西南地区旱情时空变化监测。但研究中仍存在一些问题有待解决。（1）本文选择TVDI、 RVI、 PDI和GVMI作为自变量进行建模， 但不同研究区最优遥感干旱指数的选择存在一定差异性， 自变量的选择是灵活可变的， 学者可根据研究区的实际情况进行调整； （2）本文构建的遥感干旱监测模型精度在一定程度上受估算MCI指数的站点个数及站点分布密度影响， 在站点个数较少、 空间分布不均的区域， RFR模型的拟合效果相对较差。为解决这一问题， 可考虑采用模糊数学综合评价模型， 在不使用站点数据作为辅助数据的情况下， 构建遥感干旱监测模型； （3）受西南地区多云雾天气条件影响， 光学遥感数据中存在大量无效值像元， 即使采用均值合成法得到月均值数据， 数据中仍存在部分像元为无效值（图6 ）， 严重限制了光学遥感数据在多云雾地区旱情监测的适用性。后续可采用云下像元地表参数重建技术， 恢复其对云下地表环境的监测评估能力， 从而提高区域旱情监测的时空连续性。