基于多源遥感数据和随机森林的综合旱情指标构建

董 婷 任 东 邵 攀 孟令奎

(1.三峡大学计算机与信息学院， 宜昌 443002； 2.武汉大学遥感信息工程学院， 武汉 430079)

0 引言

干旱是全球最常见的自然灾害之一，持续时间长，波及范围广[1-2]。受全球气候变暖影响，近年来干旱灾害呈上升趋势[3]。旱灾不仅影响农业生产，造成直接经济损失，而且给社会发展和人民生活带来巨大影响[4-5]。研究干旱时空演变特征、提高预警预报水平和能力具有重要的现实意义。遥感技术具有快速、经济、宏观、数据连续等优势，是当前实现大面积农业旱情监测的重要手段[6-7]。

植被长势与土壤水分密切相关，通过植被指数比较不同时期植被长势情况，是遥感监测旱情状态的重要途径，常用植被指数包括归一化植被指数(Normalized difference vegetation index, NDVI)[8]、增强型植被指数(Enhanced vegetation index, EVI)[9]、植被状态指数(Vegetation condition index, VCI)[10]等。地表温度(Land surface temperature, LST)也是旱情监测需考虑的一个重要因素，KOGAN[11]基于地表温度建立温度状态指数(Temperature condition index, TCI)，并将VCI和TCI进行线性加权，提出植被健康指数(Vegetation health index, VHI)。农业干旱是因长时间降水异常短缺造成土壤缺水，导致作物产量减少甚至无收的一种气象灾害[12]，气候降水和土壤湿度是农业干旱的重要关联因子。DU等[13]在构建综合旱情指数时考虑降水因子，利用TRMM降水数据建立降水状态指数(Precipitation condition index, PCI)，效果显著。ZHANG等[14]对AMSR-E土壤湿度数据进行归一化，得到土壤湿度状态指数(Soil moisture condition index, SMCI)；通过对SMCI、TCI和PCI进行线性加权构建微波集成旱情指数(Microwave integrated drought index, MIDI)，实验表明，由MIDI反演得到的旱情空间分布与实测值有较好的一致性。

随着长时序遥感数据的积累，除了简单线性加权方法外，利用机器学习方法构建旱情监测模型已成为实施遥感旱情监测的另一主要途径[15-17]。然而，当前基于机器学习的旱情监测方法多采用单一分类器，由于农业旱情关联因子的多样性和不确定性，单分类器旱情监测方法的监测效果不够理想。集成学习通过将多个分类器联合，能够减少单分类器的泛化误差[18-19]，为提高遥感旱情监测精度提供一种可能途径。随机森林是BREIMAN[20]提出的一种多决策树集成学习方法，具有较高的分类精度和拟合能力[21-22]，已在土地覆盖分类[23-25]、生物量估计[26-27]等环境遥感领域得到广泛应用，但在旱情监测方面的相关研究较少。

本文以实测标准化降水蒸散指数作为旱情参考指标，利用随机森林方法构建一种综合旱情监测模型。首先通过相关性分析对不同遥感旱情单因子指数的旱情监测效果进行对比，选取特定时间的建模因子，然后利用随机森林回归算法对多源旱情相关信息进行融合，建立一种多因子集成的旱情状态指数(Integrated drought condition index, IDCI-RF)，最后利用IDCI-RF对中国北部区域大范围旱情空间分布进行评估，以期为实施旱灾状态监测、制定防灾减灾决策提供参考和依据。

1 研究区与数据

1.1 研究区概况

本文针对大空间范围农业旱情监测展开研究，选择北方地区作为研究区(图1)，处于34°38′～45°46′N、108°7′～130°44′E，覆盖范围包括吉林省、辽宁省、北京市、天津市、河北省、山西省、山东省，以及内蒙古自治区、黑龙江省、陕西省、河南省部分区域。

图1 研究区概况Fig.1 Maps of study area

依据全球MODIS土地覆盖分类产品MCD12C1的分类结果，研究区主要土地覆盖类型为草原、作物区和林地，其中，草原地区主要分布在内蒙古自治区、山西省北部和河北省北部部分地区，山东省以作物覆盖为主，作物区还分布在河北省南部、吉林省和辽宁省部分区域(图1b)。根据柯本气候分类法，研究区的气候环境主要是干带草原气候和冬干冷温气候，不同气候类型区域处的降水与气温存在较大差异(图1c)。受东亚季风气候影响，研究区降水不均，易受旱涝等自然灾害影响。

1.2 气象数据

从气象数据共享服务网(http:∥cdc.cma.gov.cn/)获取研究区内所有气象站点1961—2012年间月降水量以及平均气温，通过对所有站点气象数据整理和分析，确定132个包含所有数据记录的气象站(图1a)。分别计算各省份所有站点的月降水量和平均气温算术平均值，以评估该省份的区域旱情状况。

VICENTE-SERRANO等[28]提出标准化降水蒸散指数(Standardized precipitation evapotranspiration index, SPEI)。SPEI结合了帕默尔干旱指数PDSI对蒸散的响应以及标准化降水指数SPI多时间尺度等优点，其监测效果已在全球区域的应用中得到了很好的验证[29-30]。本文利用站点长时序月降水量和平均气温计算132个站点处3种时间尺度(1个月、3个月和6个月)的SPEI值，采用THORNTHWAITE方法[31]估算潜在蒸发量PET。

1.3 遥感数据及预处理

遥感数据包括MODIS植被指数产品MOD13C2和地表温度产品MOD11C3(2003—2012年)，分别提供归一化植被指数NDVI、增强型植被指数EVI以及地表温度LST。利用质量控制文件对数据中噪声点进行掩膜，并对NDVI和LST数据进行归一化处理得到植被状态指数VCI和温度状态指数TCI。利用空间数据聚合方法对5种指数(NDVI、EVI、VCI、LST、TCI)进行重采样，使空间分辨率变为0.25°。

遥感降水量数据采用2003—2012年的热带降雨测量卫星TRMM 3B43数据，首先将降水速率转换为单月降水总量数据，同时计算出1个月、3个月和6个月的累积降水量，记为PRCP1、PRCP3和PRCP6，利用3种累积降水量建立3种降水状态指数，记为PCI1、PCI3和PCI6。

遥感土壤湿度数据采用欧空局气候变化计划项目提供的全球土壤湿度数据(2003—2012年)，包括3类数据集：主动微波数据集、被动微波数据集以及主、被动微波集成数据集。计算3类土壤湿度数据月均值作为月土壤湿度信息(记为SM-A、SM-P以及SM-M)，在此基础上进行研究区裁剪和归一化处理得到土壤湿度状态指数，记为SMCI-A、SMCI-P以及SMCI-M。

归一化方法如表1所示，归一化后的指数取值范围是[0,1]，其中0和1分别代表最干旱和最湿润状态。

本文还获取土地覆盖类型、地表高程和气候分类数据作为模型输入参数。其中土地覆盖类型采用2003—2012年的MODIS土地覆盖类型产品MCD12C1 IGBP数据；地表高程数据采用西部环境与生态科学数据中心提供的1 km分辨率中国数字高程模型DEM数据集；气候分类信息主要依据柯本气候分类法，使用PEEL等[32]提供的Koppen气候分类图作为研究区气候分类参考数据，其空间分辨率为0.1°。利用众数算法将土地覆盖类型和气候分类数据空间分辨率重采样为0.25°，利用3次卷积法将地表高程数据空间分辨率重采样为0.25°。

表1 干旱指数计算方法Tab.1 Remote sensing drought index formulas

2 研究方法

2.1 技术流程

利用随机森林回归算法集成多源遥感数据进行农业旱情监测研究，技术流程如图2所示。

首先利用Pearson相关性分析方法对比17种不同遥感旱情指数在不同时期的旱情监测效果，确定最优建模因子；然后利用随机森林回归方法对不同旱情因子进行学习并对实测指数SPEI进行拟合得到回归模型；基于回归模型对研究区2003—2012年4—10月期间所有遥感影像进行处理，得到多因子集成旱情指数IDCI-RF；最后对IDCI-RF指数进行验证与评价。验证与评价方法包括：对IDCI-RF指数与实测指数进行散点图分析和相关性验证，利用验证样本集评估IDCI-RF的泛化能力；与Bagging和Cubist回归方法进行对比，检验IDCI-RF指数旱情监测精度；讨论IDCI-RF指数在不同空间尺度的旱情监测适用性，同时对比IDCI-RF旱情监测图与站点实测指数分布情况，探讨IDCI-RF指数在大空间范围内的旱情监测能力。

图2 总体技术流程Fig.2 Flow chart of technical process

2.2 随机森林回归方法

随机森林方法是建立在分类回归树基础上的多决策树集成学习方法，在构建过程中随机选择训练样本子集和特征变量，使得该方法具有较高的抗噪能力，能够有效避免分类回归树方法的过拟合问题[20]。

利用R语言实现随机森林回归，基于该方法构建实测指数反演模型，将样本按1∶1的比例随机划分，得到训练样本集和验证样本集，为了验证模型的可靠性，本文重复进行5次随机分组，依次记为样本组1～5。其中，训练样本用于建模，验证样本用于评估模型精度以及泛化能力等。

2.3 模型评价方法

采用均方误差(Mean square error, MSE)作为确定随机森林模型参数取值的统计依据；采用决定系数(Coefficient of determination,R2)、均方根误差(Root mean square error, RMSE)和平均绝对误差(Mean absolute error, MAE)作为评价和验证模型精度的指标。决定系数R2反映数据拟合程度，取值越大，说明拟合值与实测值的相关性越好；均方误差MSE、均方根误差RMSE和平均绝对误差MAE可判断拟合值与实测值的差异，取值越小，说明拟合值越接近实测值。

3 结果与分析

3.1 旱情指数敏感性分析及建模因子确定

在不同旱情阶段植被、降水、气温、土壤等不同旱情关联因子的表现不同。分别提取2003—2012年间研究区内各站点位置处的17种遥感指数，指数包括NDVI、VCI、EVI、LST、TCI、PRCP1、PRCP3、PRCP6、PCI1、PCI3、PCI6、SM-A、SM-P、SM-M、SMCI-A、SMCI-P和SMCI-M，按照不同月份将各遥感指数分别与3种时间尺度(1个月、3个月和6个月)的标准化降水蒸散指数SPEI进行Pearson相关性分析，依据相关性结果确定建模因子。

图3(图中*表示P>0.01)为相关性分析结果。由图3可知，与实测数据的相关性，除了植被指数(NDVI、EVI、VCI)在植被生长季初期与末期以及土壤湿度指数SM-P在5月外，其他的相关性均通过0.01显著性检验。

图3 遥感指数与不同时间尺度SPEI的相关性分析结果(4—10月)Fig.3 Correlations between remote sensing drought indices and in situ different time scales SPEI from April to October

3种植被指数中，VCI与不同时间尺度SPEI的相关性都强于NDVI和EVI；TCI与实测指数相关性也优于LST。对于降水指数，与归一化前降水指数相比，归一化后的降水指数(PCI3和PCI6)与SPEI-3、SPEI-6的相关性更高。此外，由对应月份累积起来的降水指数与其相对应时间尺度SPEI的相关性强于其与其他时间尺度SPEI相关性，如由3个月累积降水量建立起来的归一化降水指数PCI3与SPEI-3的相关性都高于其与SPEI-1和SPEI-6的相关性，PCI1与SPEI-1的相关性也高于其与另外两种时间尺度SPEI的相关性。土壤湿度类指数也表现出类似的规律，归一化后的指数与3种尺度SPEI的相关性均优于归一化前的指数。

由图3可知，从4月到7月再到10月，植被类指数与不同时间尺度SPEI的相关性经历先增大后减小的过程，以SPEI-6为例，在4月，VCI与SPEI-6相关性未通过0.01显著性检验，而在6月，相关性增强(r为0.31)，到7月植被生长旺盛时期，相关系数达到最大(r为0.4)，8月开始，VCI与SPEI-6的相关性又逐渐减弱。这说明以植被指数建立起来的旱情模型在植被生长较为旺盛的阶段监测结果更具参考价值，与JI等[33]的研究结果一致。温度类指数在不同时间也表现出不同的相关性，其在植被生长季初期与实测指数的相关系数比在植被生长季中后期高，如TCI与3种不同时间尺度SPEI的相关性均在6月达到最高。不同类型的土壤湿度指数与实测指数相关性表现也不同，以SPEI-3为例，在6—10月期间，主动土壤湿度指数与SPEI-3的相关性明显优于被动土壤湿度指数和主被动合成土壤湿度指数，说明主动土壤湿度指数适合在植被覆盖度较高的区域进行旱情监测，YUAN等[34]的研究工作也得出相类似的结论。

不同时间尺度的实测标准化降水蒸散指数SPEI可以分析不同的干旱类型，小时间尺度数据适合研究土壤含水率和河道流量，大时间尺度数据可分析地下水变化情况[35-36]。本文重点研究农业旱情，IDCI-RF综合指数的构建以SPEI-3为依据。基于各类遥感指数与SPEI-3的相关性分析结果确定不同月份的模型输入数据，如表2所示。此外，旱情的发生发展与下垫面等环境因素也相关，IDCI-RF除考虑能表征植被、地表温度、降水和土壤湿度状态的4类遥感指数外，还考虑地表高程、土地覆盖类型和气候分类信息3类辅助参数。

3.2 RF参数优化及回归模型构建

表2 4—10月随机森林模型输入参数Tab.2 Inputting parameters of random forest model from April to October

图4 不同模型参数与模型袋外误差的关系Fig.4 Relationship between changing parameters and out-of-bag error

经验值不一定能获得最优模型，本文在RF回归模型构建中，为获得最优mtry，采取逐一增加变量的方法，即对于所有mtry取值分别建立回归模型，找到对应模型误差最小的mtry值作为其最优参数值。当分枝创建中自变量数目mtry最优参数值确定后，通过调节ntree，即决策树的数量，基于可视化分析观察模型误差变化趋势以确定ntree最优参数值。图4分别对应6月5组不同训练模型得到的不同模型参数与模型袋外数据集均方误差MSE的关系，可以发现，当mtry取值依次为2、3、3、2、3时，对应各自模型的MSE最小；当ntree达到1 000时，模型MSE都趋于稳定，因此，6月由5组训练集训练的回归模型中ntree确定为1 000，mtry分别设置为2、3、3、2、3，同样方法得到另外30组模型参数。

最终得到7个月的实测指数回归模型，其中每月基于不同训练集依次对应5个回归模型，模型拟合结果称为多因子集成旱情状态指数IDCI-RF。基于以上回归模型，得到验证样本IDCI-RF值，与验证样本SPEI-3实测值进行对比，评估回归模型泛化能力。图5为各月份第3组验证样本集的实测指数SPEI-3与其拟合数据IDCI-RF的散点图，可以看出所有月份的相关系数都不小于0.75，r依次为0.81、0.81、0.76、0.79、0.75、0.79和0.84，所有相关关系都通过0.005显著性水平检验，说明由随机森林算法得到的拟合模型能很好地估测实测指数，可用于实际旱情监测应用中。

图5 4—10月验证样本实测数据SPEI-3与拟合数据IDCI-RF散点图分析Fig.5 Scatter plots and correlation coefficient r values between IDCI-RF and SPEI-3 of validation samples from April to October

3.3 不同回归方法结果对比

为检验随机森林算法对旱情指数的拟合精度，选择Cubist和Bagging两种集成算法与RF进行对比，Cubist和Bagging模型输入数据可以是连续值和离散值，在土地利用分类、环境监测等领域得到广泛应用[37-38]。分别利用Cubist和Bagging模型对不同月份的5组训练集进行学习得到回归模型，模型拟合结果依次记为IDCI-Cubist和IDCI-Bagging。表3为4—10月不同模型估测的实测指数统计结果。

由表3可以看出，由RF构建的35组回归模型中，IDCI-RF与实测SPEI-3的决定系数R2在0.49～0.71之间，其中，4月、9月和10月的R2都在0.6以上；由Cubist和Bagging回归得到的IDCI与实测值的R2范围依次为0.46～0.69和0.48～0.69，低于RF模型，说明RF估测的拟合值与实测值的相关性更好。3种算法各自建立的35组回归模型中，RMSE和MAE最大值分别为0.69和0.54(RF)、0.72和0.55(Cubist)、0.71和0.56(Bagging)，说明RF估测的拟合值与实测值的差异更小。

从不同月份的不同样本组统计得到的模型精度来看，除了10月(样本组1和样本组4)以及8月(样本组1)以外，其他所有的回归模型中，由RF得到的IDCI-RF与实测值决定系数R2在3种算法中都能达到最高，RF模型表现出更强的预测能力。本文也针对不同月份求得5个样本组的统计指标平均值(表3)，4—10月IDCI-RF与实测值的平均R2依次为0.64、0.63、0.57、0.59、0.54、0.62和0.68，均优于Cubist和Bagging拟合指数，IDCI-RF对应的平均误差指标也最小。分析表明与Cubist和Bagging方法相比，由RF算法得到的IDCI-RF与SPEI-3的一致性更强，用IDCI-RF指数估测SPEI-3指数合理。

3.4 IDCI-RF在不同空间尺度上的敏感性分析

根据各月份5组不同样本集得到的统计评价指标结果，将由第3、3、1、3、1、1、3个样本组构建的回归模型依次确定为4—10月的最优模型。分别提取研究区内各省份所有站点的由最优模型得到的IDCI-RF以及实测指数SPEI-3值，利用散点图分析和统计指标验证方法对IDCI-RF指数在不同省级区域的监测精度进行评估。

图6为7个主要省份(河北省、山西省、内蒙古自治区、辽宁省、吉林省、山东省和陕西省)2003—2012年植被生长季IDCI-RF指数与实测指数的散点图以及统计值验证结果，可以看出IDCI-RF指数在不同省份都能很好地拟合实测指数，所有相关关系都通过0.005显著性检验，决定系数R2都在0.7以上，其中在内蒙古自治区(n=1 557)和辽宁省(n=1 039)(图6c、6d)，IDCI-RF与SPEI-3的相关性最强(R2=0.82)。除吉林省和陕西省外，其他省的实测值与IDCI-RF值的误差指标均在0.5(RMSE)和0.35(MAE)以下。

表3 不同模型实测指数回归统计结果Tab.3 Statistics results for different regression models of in situ reference data

图6 IDCI-RF指数与实测指数SPEI-3的散点图Fig.6 Scatter plots between IDCI-RF and SPEI-3

为验证站点尺度IDCI-RF指数效果，选取研究区部分站点，针对站点处IDCI-RF与实测指数的变化趋势展开分析，选取的站点分散在研究区各区域，能较好地代表研究区不同土地覆盖类型。图7为21个站点在2005—2012年7月时的IDCI-RF与实测SPEI-3变化趋势。

从图7可以看出，IDCI-RF指数能较好地拟合实测指数，但也存在少部分站点处的IDCI-RF与SPEI-3变化不一致，如53480站点的实测SPEI-3在2009年7月至2010年7月呈上升趋势，但其在2010年7月的IDCI-RF值比2009年同期稍低(图7d)；54311站点处SPEI-3在2007年7月到2008年7月呈上升趋势，而IDCI-RF指数在相应时间的变化趋势是略下降的(图7p)；54826站点处的实测指数在2007年7月到2009年7月呈先下降后上升趋势，但IDCI-RF指数表现的变化趋势与实际相反(图7u)。尽管存在少数不一致现象，但在绝大部分情况下，IDCI-RF指数的变化规律能与实测指数SPEI-3保持一致，还有部分站点处的IDCI-RF值与实测值接近(图7a、7e、7o、7r、7t)，说明IDCI-RF指数能较精准地拟合实测指数，可用于判断实际旱情状态变化趋势。

图7 不同站点IDCI-RF与SPEI-3的变化趋势(2005—2012年7月)Fig.7 Year-to-year changes in IDCI-RF and SPEI-3 in representative meteorological sites for July from 2005 to 2012

3.5 旱情空间分布格局

基于研究区2003—2012年多源遥感数据集和3类辅助参数，利用不同月份最优回归模型得到研究区所有位置处IDCI-RF值，分析区域内旱情空间分布特征。以干旱年2009年为例，图8为该年度4—10月的IDCI-RF监测图，为评估IDCI-RF旱情监测精度，不同站点处的实测指数SPEI-3值也依据点的大小进行区分，站点位置处点越大，说明实测指数值越小，旱情越严重，反之亦然。

图8 2009年4—10月IDCI-RF空间分布图以及实测指数SPEI-3分布情况Fig.8 Seasonal changes in drought detected by IDCI-RF and SPEI-3 from April to October in 2009

由图8可知，内蒙古自治区在2009年8—10月旱情有加重趋势，尤其是研究区内内蒙古自治区北部区域，这与实测SPEI-3空间分布一致，该区域站点位置处点变大，说明旱情越来越严重(图8e～8g)。IDCI-RF旱情图显示2009年6月陕西省北部出现严重旱情，从站点大小看出，站点处SPEI-3值较小(图8c)，也说明旱情严重。2009年4月辽宁省所有站点处SPEI-3指数都大于-0.5，IDCI-RF监测图反映的2009年4月辽宁省旱情与实测数据一致(图8a)；从7月开始，IDCI-RF监测图反映出辽宁省出现旱情并表现出加重趋势，SPEI-3也反映出同样的旱情形势(图8d～8f)。

从以上分析可知，由IDCI-RF监测图反映的研究区旱情状态与实测指数SPEI-3分布特征有较好的一致性，2009年4—10月的IDCI-RF旱情空间分布图能较精确地指出研究区干旱受灾区域，也能在一定程度上描绘出不同区域旱情程度及变化趋势，说明IDCI-RF指数可应用于大范围区域农业旱情状态监测业务中。

4 讨论

农业干旱是一种异常复杂的自然灾害，以往的遥感旱情监测指数多侧重于考虑植被、温度、土壤等单一因素，这类指数仅是旱情某方面、某个层次特征的表征，影响监测效果。随着历史数据的积累，机器学习技术被应用于旱情监测中，目前研究多数采用单一分类器对旱情特征进行挖掘。

选择随机森林方法，对多源遥感旱情信息进行集成构建综合旱情指数，取得较优的监测效果。原因为：综合考虑旱情关联因子的多样性，选取特征集包括能表征植被状态、地表温度、气候降水、土壤湿度的多种旱情关联因子，以及3种生态环境因素(地表高程、土地覆盖类型、气候类型)；选择的建模方法随机森林依据分类回归树构建，是一种集成学习方法，与单分类器相比，能够实现信息优势互补，获得更高的拟合精度。下一步尝试将其他集成学习方法应用到遥感旱情监测中。

采用的遥感数据空间分辨率较低(0.25°×0.25°)，混合像元对旱情监测效果有一定影响，后期工作中将采用更高空间分辨率的遥感影像作为基础数据源来开展旱灾监测研究。

5 结论

(1) 利用RF、Cubist和Bagging 3种方法分别建立回归模型，利用RF构建的IDCI-RF与实测SPEI-3的平均决定系数R2为0.54～0.68，优于另外两种方法，RF模型表现出更强的预测能力。

(2) 研究区主要省份的IDCI-RF指数与实测

SPEI-3的决定系数R2都在0.7以上，大部分站点的IDCI-RF变化规律与实测SPEI-3保持一致，IDCI-RF值与实测值接近；由IDCI-RF监测图反映的旱情状态与SPEI-3空间分布特征吻合度较高，IDCI-RF监测图能较准确地指出干旱受灾区域，表明IDCI-RF指数在实际大范围旱情监测中具有应用潜力。