基于随机森林偏差校正的农业干旱遥感监测模型研究

刘 冀 张 特 魏 榕 张 茜 刘艳丽 董晓华

(1.三峡大学水利与环境学院， 宜昌 443002； 2.三峡库区生态环境教育部工程研究中心， 宜昌 443002；3.南京水利科学研究院水文水资源与水利工程科学国家重点实验室， 南京 210098；4.水利部应对气候变化研究中心， 南京 210029)

0 引言

干旱是一种普遍的自然灾害，对农业、生态环境、社会经济产生重大影响[1]。在全球气候变暖的大背景下，未来全球干旱频率将逐渐上升[2]。淮河流域是我国重要的农业基地，亦是我国的干旱多发区[3]。因此，针对淮河流域开展农业干旱监测模型研究对抗旱工作的开展具有重要的现实意义。

帕尔默指数(Palmer drought severity index，PDSI)、标准化降水指数(Standardized precipitation index，SPI)、标准化降水蒸散指数(Standardized precipitation evapotranspiration index，SPEI)被广泛应用于干旱监测，并取得了较好的监测效果[4-6]，尤其是SPEI指数与土壤水分之间有较强的相关性[7-8]。干旱是一种区域性现象，这些指数通过站点观测数据计算得出，仅能准确反映站点附近的旱情。传统方法通过空间插值实现区域干旱监测，其结果受地理因素与插值方法等影响，存在较大的不确定性[9]。遥感技术能够提供降水量、地表温度以及植被长势等空间分布信息，为实现区域农业干旱监测提供了有力支撑。一些研究尝试应用多源遥感数据构建综合干旱监测指数，并证实了多源遥感数据在干旱监测中的应用潜力,但这些研究采用权重方法结合多源遥感干旱因子，难以描述不同干旱因子间的非线性关系[10-12]。近年来，一些学者采用支持向量机、人工神经网络、分类回归树、深度人工神经网络等机器学习方法，融合多源遥感信息构建模型，预测旱情空间分布，为实现区域干旱监测提供了新思路[13-14]。

随机森林(Random forest，RF)是基于集成思想的机器学习算法，较单一分类器的误差更加稳定，具有较强的抗过拟合能力[15]。董婷等[16]采用RF、Cubist和Bagging 3种方法在黄淮地区构建了综合干旱指数，结果表明，RF算法的拟合能力更强。PARK等[9]研究表明，相比增强回归树、Cubist方法，RF算法对旱情预测更加准确。RF算法基于多棵决策树的平均结果进行预测，其结果相对准确可信，但同时也会导致一定的偏差，尤其是对极值的预测能力较弱[17]。对于干旱灾害而言，极端情况往往带来更大的损失，更应予以充分重视。相关研究表明，对RF算法进行偏差校正能有效提升预测结果的准确性[18]。

本文在综合比较多种RF偏差校正方法适用性的基础上，合理选取最优偏差校正方法，构建基于多源遥感数据的淮河流域农业干旱监测模型，并分析模型在农业干旱监测中的适用性，以期提高空间旱情监测的准确性，为抗旱工作的开展提供科学支撑。

1 研究区概况

淮河流域位于我国东部(111°55′～121°25′E， 30°55′～36°36′N)，地跨湖北、河南、安徽、江苏及山东五省，是我国重要的农业基地。流域面积约270 000 km2，其中约有2/3的面积为平原和1/3的面积为丘陵山地。淮河流域地跨我国南北气候带，淮河以北属暖温带区，淮河以南为北亚热带区，流域内气候温和，年平均降雨量约1 040 mm，年平均气温11～16℃，年内降雨集中，汛期降水占全年降水量的60%～80%，其独特的气候条件和地理位置，导致流域内的旱涝灾害频发。流域位置、高程及土地利用概况如图1所示。

图1 淮河流域地理位置、高程及土地利用情况Fig.1 Geographic location, elevation and land use type of Huaihe River Basin

2 数据与方法

2.1 研究数据

研究数据主要包括1970—2014年的气象数据以及2001—2014年的遥感数据两部分。气象数据主要包括淮河流域内31个气象站的逐日气象数据，含14个农业气象站的逐旬土壤相对湿度数据，统一处理至月尺度。遥感数据包括TRMM(Tropical rainfall measure mission)3B43降雨数据、MODIS(Moderate-resolution imaging spectroradiometer)的MOD11A3地表温度数据和MOD13A3归一化植被指数(Normalized vegetation index，NDVI)、STRM(Shuttle radar topography mission)高程数据。此外还包括土壤有效持水量(Available water capacity of soil，AWC)栅格数据，统一处理至1 km空间分辨率。各项数据来源信息如表1所示。

表1 研究数据来源及说明Tab.1 Sources and instructions of data

2.2 研究方法

2.2.1SPEI指数

SPEI指数采用降水与潜在蒸发的亏缺程度来量化干旱程度[19]，计算式为

(1)

Pi——降雨量,mm

PETi——潜在蒸发量,mm

k——时间尺度

i——分析时段起始时间

j——分析时段终止时间

SPEI指数可有效评估不同时间尺度的干旱情况，其中3个月尺度的SPEI指数与土壤湿度有较好的相关性[8]。通常采用Thornthwaite公式计算PET，但该公式受温度影响较大，因此本文采用Penman-Monteith公式计算的参考蒸发量ET0代替PET，并计算1970—2014年各站点SPEI3指数。当SPEI<-2时，定义为极端干旱(极旱)；当-2.0≤SPEI<-1.5时，定义为严重干旱(重旱)；当-1.5≤SPEI<-1.0时，定义为中等干旱(中旱)；当SPEI≥-1.0时，定义为无旱。

2.2.2遥感指数

对各像元1、3个月的平均降水量、地表温度、NDVI进行归一化处理，即可计算出不同时间尺度(1、3个月)的降水状态指数(Precipitation condition index，PCI)、地表温度状态指数(Temperature condition index，TCI)以及植被状态指数(Vegetation condition index，VCI)，计算方法如表2所示。

表2 遥感指数计算方法Tab.2 Computation of different remote sensing drought monitoring indexes

2.2.3随机森林算法

随机森林是一种基于决策树的集成算法，在回归及分类问题中应用广泛。该方法通过Bootstrap自助法在原数据集中有放回地随机抽取n个样本集，每次抽取2/3样本容量的数据作为袋内数据，建立n棵决策树构建随机森林，并以n棵决策树回归结果的平均值进行预测。该方法的优点在于训练速度相对较快且不需要进交叉验证，同时抽样随机性和特征选择随机性使得随机森林不易陷入过拟合[20]。随机森林方法中的两个重要参数分别为决策树的数量ntree及随机选取特征变量的数量mtry，算法详见文献[17]。本文采用R语言中的Random forest程序包(http:∥cran.at.r-project.org/src/contrib/Archive/randomForest/randomForest_4.6-7.tar.gz)进行计算。

2.2.4偏差校正方法

为提高RF模型精度，采用4种偏差校正方法对RF的回归结果进行校正。

(1)简单线性回归法

简单线性回归法(Simple linear regression, SLR)主要通过建立袋内数据预测值与实际值间的线性关系式，并应用该线性关系式对袋外样本数据的预测值进行修正，以达到偏差校正效果。计算公式为

yobs=a+bpre

(2)

式中yobs——实际值pre——预测值

a、b——系数

(2)偏差估算法

ZHANG等[21]提出的偏差估算法(Bias corrected, BC)以预测值和实际值之间的残差为因变量，以原样本数据及预测值为自变量，通过RF构建估计残差的模型，在原始预测结果上进行偏差校正。计算式为

r=yobs-pre

(3)

=RFres(X,r)

(4)

cor=pre+

(5)

式中r——实际残差——预测残差

RFres——残差计算函数

X——样本数据

(3)残差旋转法

SONG[18]提出基于残差旋转(Residual rotation, RR)的偏差校正方法，能有效降低原始RF的均方误差。该方法通过建立预测值与残差间的线性关系，将拟合直线旋转至与X轴重合位置，通过旋转矩阵估计残差对原始预测值进行校正。

(4)最优角度残差旋转法

最优角度残差旋转法(Best-angle residual rotation，BRR)基于RR法，对旋转角进行寻优，使整体均方误差最小。例如，将拟合曲线旋转α可与X轴重合，在(α-β，α+β)区间寻找一个最优角度，采用该角度下的旋转矩阵对预测值进行校正。

2.2.5模型精度评估方法

以站点SPEI3为因变量，以对应站点的PCI1、TCI1、VCI1、PCI3、TCI3、VCI3(1、3分别代表1、3个月时间尺度)、高程、AWC共8个变量为自变量，建立2001—2014年作物生长季(4—10月)的RF模型，采用上述4种偏差校正方法对RF预测结果进行校正，并选取最优的校正方法构建随机森林偏差校正干旱监测模型(Bias-corrected random forest drought condition，BRFDC)。采用均方根误差(Root mean square error，RMSE)最小原则确定RF模型参数mtry，并选取决定系数(Determination coefficient，R2)、RMSE以及干旱等级预测准确度评估BRFDC模型精度，综合比较后确定最优校正方法，最后通过站点土壤相对湿度数据及流域典型干旱事件记录对BRFDC模型的监测能力进行评估。此外，考虑到淮河流域特殊的气候分区特征，采用C均值模糊聚类法将淮河流域划分为两个气候区分别研究，计算流程如图2所示。

图2 计算流程图Fig.2 Flowchart of procedures used

3 结果与分析

3.1 流域分区及RF参数优化

以各气象站年降雨量与参考蒸发量的差值(P-ET0)为依据，将淮河流域各气象站划归为北(区域1)、南(区域2)两个分区，站点聚类结果如图1所示，各分区气象站特征如表3所示。由表3可知，区域1、区域2的多年平均降雨量分别为755、1 018 mm，其气候差异主要体现在年降雨量，而年平均气温T和ET0的差异较小。

表3 不同区域气候特征Tab.3 Climatic characteristics of different regions

本文构建的RF模型包含8个变量，从1～8逐一取值以确定最优参数mtry，而ntree参数统一设置为1 000。将不同参数mtry的RF模型独立运行100次，各参数取值下的平均RMSE如图3所示。由图3可知，当RMSE最小时，区域1中4—10月各模型的mtry取值分别为4、7、4、3、5、4、5，区域2的mtry取值分别为5、5、7、6、5、7、5。

图3 不同区域各参数取值对模型RMSE的影响Fig.3 RMSE of different parameter values on model of region 1 and region 2

3.2 偏差校正方法精度分析

采用SLR、BC、RR及BRR共4种方法对RF模型监测结果进行偏差校正，不同方法校正后监测结果与站点SPEI3的散点图如图4所示。由图4可知，区域1中RF监测值与站点SPEI3的R2和RMSE分别为0.719和0.534，经偏差校正后的R2为0.719～0.897，RMSE为0.335～0.533；区域2中RF监测值与站点SPEI3的R2和RMSE分别为0.751和0.502，经偏差校正后的R2为0.751～0.874，RMSE为0.362～0.501。由此可见，偏差校正后较校正前精度均有一定程度上的提高，其中SLR法的校正效果并不明显，而BC、RR以及BRR法对模型精度有显著提升。整体而言，BRR方法的校正效果最优，校正后各区域的R2分别为0.897和0.874，RMSE分别为0.335和0.362。

图4 不同区域各校正方法的精度评估Fig.4 Scatter diagrams of drought index estimated by different bias-correcting methods and SPEI for three-month time scale

不同校正方法对干旱等级监测的准确率如表4所示。由表4可知，RF对区域1中旱、重旱和极旱的监测准确率分别为45.1%、40.2%和0，经校正后的中旱、重旱和极旱监测准确率分别为44.4%～69.3%、39.2%～74.2%和0～50%；RF对区域2中旱、重旱和极旱的监测准确率分别为36.6%、18.6%和0，经校正后的中旱、重旱和极旱监测准确率分别为39.2%～64.1%、18.6%～43.3%、0～33.3%。整体而言，SLR法的校正效果较弱，而BC、RR及BRR法则有效提升了干旱等级监测的准确率，其中BC方法的校正效果最优，尤其是对于极端干旱监测具有显著的优势。综合考虑R2、RMSE以及干旱等级监测准确率的评估结果，最终采用BC方法构建BRFDC模型。

表4 各校正方法在不同区域对各干旱等级监测准确率Tab.4 Accuracy of different drought conditions monitored by all bias-correcting methods drought in two regions %

3.3 干旱监测能力评估

为评估BRFDC模型在农业干旱监测中的适用性，采用农业气象站相对土壤湿度及流域典型干旱事件对BRFDC模型的监测能力进行验证。

采用BRFDC模型计算指数、站点SPEI3分别与农业气象站10、20 cm的土壤相对湿度进行相关性分析，相关系数如表5所示。由表5可知，多数站点的SPEI3、BRFDC模型计算指数与10、20 cm的土壤湿度显著相关(P<0.01)，二者都能有效表征土壤墒情的变化。相比SPEI3，共有10个站点10 cm土壤相对湿度与BRFDC模型计算指数相关性更强，此外，9个站点20 cm土壤相对湿度与BRFDC模型计算指数相关性更显著。由于土壤水分是农业干旱的决定性因素，因此BRFDC模型计算的指数较SPEI3更适宜监测农业干旱。土壤水分变化受诸多因素影响，SPEI3仅采用气象要素进行计算，在表征土壤水变化过程上略有不足，而BRFDC模型综合考虑降水、地表温度、植被长势及环境因素，其计算指数更贴近于真实土壤水分的变化过程。在构建农业干旱监测模型时，需综合考虑各种潜在影响因素。

表5 BRFDC模型模拟干旱指数与不同深度土壤相对湿度相关性分析Tab.5 Correlation analysis between drought index simulated by BRFDC and soil relative moisture in different depths

根据文献[22]记载，2001年淮河流域发生了多年罕见的春、夏、秋连旱。河南省3—5月无有效降水，8月信阳地区创50多年降水最低纪录；安徽省5—9月降水量较常年偏少50%～60%；江苏省淮北地区5—6月降水量偏少，7—8月降水与往年基本持平，旱情得到缓解，9—11月降水较同期偏少81%，旱情进一步发展；山东省淮河流域汛期降水偏少，且主要集中在7—8月。

采用反距离平方(Inverse distance weighting，IDW)法对站点SPEI3进行空间插值(简称IDW模型)，并应用RF模型和BRFDC模型模拟了2001年5—10月淮河流域旱情空间分布，结果如图5所示。由图5可知，3种方法对旱情发展过程的监测结果与真实旱情基本一致，并且在空间上具有较好的一致性。IDW模型通过简单的数学方式即可实现空间旱情监测，能大致反映干旱的时空特征，但受站点个数及插值方法限制，无站点区域的旱情监测结果存在不确定性。例如9月流域西北部出现中度旱情，由于该区域仅有开封站一个站点，而未能监测到该区域的旱情；10月流域南部旱情减轻，而该区域仅霍山站一个站点，对旱情等级有所误判，而RF和BRFDC模型的监测结果则有效反映了这些区域的旱情。RF模型以多源遥感数据为输入，考虑了降水、地表温度、植被状态、高程、AWC等因素的空间异质性，其空间监测结果更加可靠。通常极端干旱将带来更大的损失，需要重视，但RF对6、9月流域南部以及10月流域中部极端干旱的监测效果并不理想，而BRFDC模型则准确识别了流域内的极端旱情。

图5 淮河流域2001年5—10月IDW、RF、BRFDC模型干旱监测结果Fig.5 Drought maps of Huaihe River Basin from May to October in 2001

整体上看，BRFDC模型与土壤湿度具有较好的相关性，并且能有效反映淮河流域2001年5—10月的干旱事件的时空演变过程，对于极端旱情的监测结果也相对准确，表明该模型可用于淮河流域农业干旱监测。

3.4 讨论

农业干旱涉及降水、植被、温度以及环境因素影响，其过程非常复杂。结合多源遥感数据进行监测，能有效解释农业干旱的复杂性。相比原始RF而言，本文构建的BRFDC模型具有更准确的干旱监测能力，对极端干旱监测的准确性从0提升到33.3%～50%，模型在淮河流域农业干旱监测中取得了较好的效果。BRFDC模型依据免费的遥感数据即可实现区域农业干旱监测，能有效减少开展干旱监测业务的成本，为缺资料地区提供可靠的旱情空间分布信息。通常降水对干旱程度起决定性作用，而本文仅采用了TRMM卫星的降水产品，其误差势必会对模型精度产生一定影响，提高降水产品的精度能显著改善监测结果的准确性。由于不同的卫星降水产品在不同区域、时间的精度有所差异[23]，融合多种卫星降水产品是提高干旱监测准确性的重要途径。此外，本文在构建模型时，考虑了降水、地表温度、植被长势、高程及AWC等影响因子，但农业干旱的复杂程度远不止如此。在后续的研究中，仍有一系列其他遥感数据需要考虑，以期提高模型的监测性能，例如MODIS蒸发数据[24]、SMOS(Soil moisture and ocean salinity)土壤水分数据[25]以及土地利用类型等数据。

4 结论

(1)BRR法校正后，监测结果的R2和RMSE分别为0.897、0.874和0.335、0.362，优于其他方法；BC法校正后，预测结果的R2和RMSE分别为0.870、0.854和0.365、0.387，同时显著提高了模型对极端干旱监测的准确率，达到33.3%～50.0%。

(2)BRFDC模型计算指数与站点实测10、20 cm的土壤相对湿度具有显著相关性，且BRFDC模型比SPEI3更适于农业干旱监测。

(3)相比IDW空间插值法及原始RF模型，BRFDC模型模拟结果能更加准确地反映淮河流域2001年典型干旱事件的时空演变过程。