基于随机森林的旱灾损失模型构建与应用
——以河南省雨养小麦为例*

侯陈瑶，朱秀芳*，孙丹一，徐 昆，刘 莹

(1. 环境演变与自然灾害教育部重点实验室，北京师范大学 地理科学学部，北京100875；2. 北京师范大学 地理科学学部遥感科学与工程研究院，北京 100875)

干旱是全球范围内各种气象灾害中分布最广、最为频繁、持续时间最长、影响最为严重的气象灾害，每年由干旱造成的经济损失远超过其余气象灾害[1]。近几十年气候变化导致干旱呈愈演愈烈的趋势，并且模式预估21世纪全球干旱风险将进一步加剧[2]。与此同时，灾害管理已从传统的应急管理过渡到风险管理，而相应的灾害监测、预警与评估机制尚未全面建立，急需加强灾害监测和损失预评估方面的研究工作[3]。

近些年来，在大数据时代的背景下，数据挖掘的算法作为定量精确评估灾害风险和快速评估灾害损失的方法而受到高度重视。例如林江豪等收集广东省造成直接损失的台风数据，利用BP神经网络和空间向量模型VSM构建台风灾害经济损失模型[4]。吴小君等依托流域灾害系统理论，结合江西山洪特点和历史山洪灾害调查数据，建立基于随机森林算法的山洪灾害风险评估模型[5]。李远远等以云南省怒江州泸水县为研究区综合考虑多种地质灾害影响因子，采用确定系数法计算各因子的敏感性值，建立基于确定系数法和支持向量机(CF-SVM)的地质灾害易发性评估模型[6]。数据挖掘算法中的随机森林方法由Breiman于2001年提出[7]，其作为一种较新的机器学习算法具有预测精度高且不易产生过拟合等优点[8]。相比神经网络和线性回归，其表现更稳定，对噪声和异常值有很好的容忍性[9]。该模型已经应用于洪涝风险分析[10]、滑坡危险性评价[11]、泥石流易发性评价[12]、遥感干旱监测[8]等方面，但在干旱灾害损失模型构建方向还鲜少研究。

河南省是我国粮食大省，迄今为止，冬小麦粮食产量已经连续12年位列全国之首，为保证我国粮食安全做出了突出的贡献[13]。干旱对于河南省冬小麦的产量的影响较大，尤其在雨养区，即粮食作物主要依靠天然降水作为水源的地区。综上考虑，本文以河南省为研究区，以雨养区冬小麦为研究对象，将河南省雨养区的冬小麦产量分解为趋势产量和气象产量，计算出减产率。与此同时，以多时间尺度干旱指数SPEI为输入，利用游程理论和河南省冬小麦生育期资料确定生育期内干旱事件，在此基础上利用随机森林建立农业旱灾损失函数，定量评估旱灾引起的产量损失，以期为冬小麦产量损失估算等提供依据。

1 研究区与数据

1.1 研究区

河南省位于我国中东部、黄河中下游，界于31.23°～36.22°N和110.21°～116.39°E之间[14]。处于南北气候过渡带，处于暖温带和亚热带之间，冷热变化和干湿状况主要受季风影响。省内年平均降水量约为500～900 mm，南部及西部山地较多，自东南向西北逐渐减少。降水时空分配不均衡，年内降水的60%集中在7-9月，年际变化大[15]。全省年平均气温一般在12℃～16℃之间，可满足一般作物的两年三熟或一年两熟的生长发育需要。本文选取河南省县一级单位为研究区域，1990-2015年为研究时间段。河南省共174个县区与市辖区，在历年的河南省统计年鉴中，共有106个县有历年冬小麦粮食产量数据，此106个县即为本文的研究区。

图1 研究区地理位置

1.2 研究数据

所用研究数据包括：①来源于1991年到2016年河南省统计年鉴的1990-2015年每个县区的冬小麦播种面积和总产量，该数据主要用于计算各县小麦的单位面积产量；②来源于2000年的河南省统计年鉴的县区有效灌溉面积数据，该数据用于提取以雨养种植为主的县；③来自中国科学院资源环境科学数据中心的2000年中国土地利用类型空间分布图(http://www.resdc.cn)，该数据用于提取2000年河南省耕地空间分布；④来源于全球标准化降水蒸散指数(Standardized Precipitation Evapotranspiration Index, SPEI)数据库(http://sac.csic.es/spei/database.html)的SPEIbase v2.5数据集。该数据集提供了全球1901-2015年1～48个月时间尺度的0.5°分辨率的SPEI月值数据，基于该数据集本文提取了1990-2015年间河南省境内冬小麦生育期内不同时间尺度(1～8月尺度)的SPEI数据集，分别记作SPE-1，SPEI-2，…SPEI-8该数据集一方面用于小麦旱灾损失模型的构建，另一方面用于河南省农业旱灾特征分析。

2 研究方法与技术路线

本文技术路线如图2所示，基本步骤包括：①将实际产量分解为趋势产量和气象产量，计算产量损失率；②根据有效灌溉面积占耕地面积的比例，设定阈值提取以“雨养”为主的县；③计算河南省各县耕地空间分布内的不同时间尺度的SPEI均值；④基于游程理论以及河南省冬小麦生育期对不同时间尺度的SPEI，进行河南省干旱事件的识别和特征分析；筛选影响冬小麦生长的干旱事件作为参与建模的样本数据，利用随机森林的方法构建冬小麦灾损模型，并进行评价和精度验证。

图2 技术路线

2.1 河南省干旱特征分析

SPEI用降水与潜在蒸散之差来代替标准化降水指数中的降水异常，同时考虑了气温对于干旱的影响，因此，SPEI具有对气温敏感以及多时间尺度的特征[16-17]，适合于月以上尺度气候状况的干旱监测。基于SPEI，运用游程理论进行河南省干旱事件的识别和干旱特征变量的提取[18]。具体做法是：以SPEI的分级标准为参考，选取-0.5作为截断水平进行阈值设定，记录SPEI序列中一个或者多个时间内出现连续小于阈值的负游程，进而计算干旱出现的次数，以及每次干旱事件的历时(即负游程的长度)和强度(干旱历时与截断水平包含的面积)。

2.2 县级不同时间尺度干旱指数的计算

本研究所用的SPEI为0.5°×0.5°栅格数据，首先将其重采样为1 km×1 km，然后基于2000年河南省耕地空间分布图提取出河南省耕地范围上的SPEI干旱数据，为了获得每个县区耕地的干旱指数，分别针对每一个县区，找到其耕地范围内的所有干旱指数的栅格求其均值，得到县级的干旱指数。河南省冬小麦通常于9月下旬播种，次年6月收获，考虑到冬小麦在刚播种时期和收获时期受需水量较小，本文定义研究时间段为10月到第2年的5月，共8个月。考虑实际应用需求，结合SPEI具有多时间尺度的特征，分别采用生育期对应月份的SPEI-1、SPEI-2、SPEI-3数值一共24个因子作为建模的输入变量。本文使用不同时间尺度及月份的干旱指数，以分析干旱对河南省冬小麦产量的影响。

2.3 县级冬小麦减产率的计算

粮食作物产量受到多种因素的相互作用，包括自然因素与非自然因素，目前学界将这些因素根据其影响机制划分为农业技术条件、气象条件和随机噪声三大类[19]。对产量进行模拟时，主要把作物产量分解为趋势产量、气象产量和随机产量三部分[20]，由于随机产量的不确定性，通常将其忽略不计[21]。因此，实际产量可以分解为趋势产量和气象产量之和，使用冬小麦单产来衡量产量大小：

yt=Yt+Yc。

(1)

式中：yt为冬小麦实际单产(kg/hm2),Yt为趋势单产(kg/hm2)，Yc为气候单产(kg/hm2)。

在农业生产的过程中，农业技术的改进对粮食产量起到积极作用，主要表现在冬小麦单产的时间序列呈现缓慢上升的趋势，所对应的产量分量即为趋势产量。河南省冬小麦产量受物种属性、自然条件的约束，其理论增长趋势较为符合有增长上限曲线的趋势模型，其中应用最为广泛的是龚珀兹曲线——关于拐点不对称的S形曲线，所以本文采用龚珀兹曲线方程对河南省冬小麦的趋势产量进行拟合：

Yt=Le-ae-bt，(a>0,b>0)。

(2)

式中：Yt为龚珀兹曲线拟合的趋势产量，t表示年序，初始值为1，代表1990年，a和b为待估计参数。

当龚珀兹曲线拟合的趋势产量与实际产量之间的残差序列存在自相关时，不满足计量经济学关于随机扰动项的经典假设，会导致拟合模型结果失真，所以使用自回归移动平均模型(Autoregressive-moving-average model，ARMA)对残差序列进行模拟，以进一步提取变化规律，最后将ARMA对残差的预测结果与龚珀兹曲线得到的趋势产量相加得到最终的趋势产量，以提高趋势产量拟合精度。由龚珀兹曲线和ARMA模型构成混合时间序列模型：

ut=Yt+uc;

φ(B)uc=θ(B)εt。

(3)

式中：ut为混合时间序列预测的趋势产量，Yt为龚珀兹曲线拟合的趋势产量，uc是龚珀兹拟合残差的趋势校正值。φ(B)=1φ1B-φ2B2-…-φpBp，θ(B)=1-θ1B-θ2B2-…-θqBq，B为滞后算子，εt为白噪声序列。

假设研究时段内实际产量未达到趋势产量的年份为气候减产年，于是定义每年的干旱减产率为

rt=max{0,(ut-yt)/ut}。

(4)

式中：rt为第t年的减产率，ut为第t年的趋势产量，yt为第t年的实际产量。其中在农业气象研究中，定义减产率小于-3%为减产年[22]。

2.4 旱灾损失函数构建

由于灌溉为主的县区受干旱的影响较小，而雨养区的冬小麦产量更易因气候变化而减产，所以本文在建立旱灾损失函数时仅考虑雨养为主的县区。假定耕地分为“雨养”和“灌溉”两种，每个县的有效灌溉面积与雨养区面积占比之和为1。给定雨养区的阈值，某个县的雨养区面积大于该阈值即认为该县为以雨养为主的县区。本文设定阈值为0.6，将县区划分成以“雨养”和“灌溉”为主的县，在雨养为主的县区进行建模。综合考虑河南省冬小麦生育期结合游程理论中识别出冬小麦生长阶段的干旱事件，在0.6的阈值下，以公式(4)计算出的减产率为y，以生育期内各月份的SPEI-1、SPEI-2和SPEI-3一共24个干旱指数为自变量x进行随机森林回归建模，利用网格寻参法选择灾损模型的最优参数。

2.5 模型精度验证

本文通过平均绝对值误差(Mean Absolute Error，)和均方根误差(Root-Mean-Square error， )来进行模型的精度评价，具体见公式(5)和公式(6)。

以公式(4)计算出的小麦减产率为真值，计算通过旱灾损失模型估算的小麦减产率和真值之间的平均绝对值误差，以此来验证旱灾损失函数的精度。

(5)

(6)

式中：ri和Ri分别表示第i个县的小麦减产率的真值和估计值，n表示0.6阈值下满足条件的样本点的个数，具体来说阈值0.6所对应的样本点的个数为209。

3 研究结果

3.1 河南省干旱特征分析

干旱是一个逐渐持续累积的过程，某个时期的干旱与当时的降水量、前期降水以及土壤墒情和地表蒸散发等多种要素共同作用有关。图3为1990-2015年间河南省不同时间尺度SPEI的强度变化，其中SPEI-1主要表现月值变化，SPEI-3主要表现季节变化，SPEI-8可以表征河南省整体冬小麦的全生育期的变化。由图3可以看出河南省耕地范围上不同时间尺度的SPEI总体呈现干湿一致的趋势。其中，1992年、1997年、1999年、2000年、2002年与2011-2014年识别出的干旱事件，与河南省实际发生干旱情况相符，都较为严重。其中2014年各干旱指标表征此次旱灾持续时间长，2000年前后和2011年各干旱指标表征发生干旱的整体强度大。

图4反映了河南省冬小麦全生育期干旱历时、次数和强度特征。识别出的干旱历时和干旱强度年均值的分布具有一致性。豫西、豫西南、豫东北地区的干旱历时和干旱强度年均值明显大于河南省其它地方。其中三门峡市、洛阳市、南阳市、驻马店市和濮阳市年均值较大；开封市、周口市和商丘市年均值较小。就干旱发生次数来说年平均干旱次数在豫西北地区整体数值最高，豫东地区其次。洛阳市、济源市、焦作市和周口市数值最高。豫东南地区整体次数数值最小，信阳市最小。

图3 1990-2015年间河南省不同时间尺度SPEI的强度变化图

3.2 冬小麦减产率计算

由于ARMA建模的条件为时间序列平稳，研究区内的4个县区的单产数据在使用混合时间序列建模时，不满足平稳性条件，对其进行剔除，使用剩余102个县区进行进一步的建模。利用龚珀兹曲线-ARMA混合时间序列模型对河南省研究区内各个县区的冬小麦产量进行趋势拟合，得到每个县区的趋势产量，然后利用趋势产量与实际产量计算出每年的减产率。其中龚珀兹曲线的参数估计较为复杂，由于样本量较少，三和值法不适用。当龚珀兹曲线的拟合对象的上限值 可以由一些已知条件事先确定时，能够采用最小二乘法估计参数。对曲线模型两边同时取对数，得到线性化变换为：

(7)

图4 冬小麦全生育期干旱指数(SPEI8-5)空间特征

表1 郸城1990-2015年冬小麦产量

以河南省郸城为例，计算时先观察郸城冬小麦历史产量时间序列，选择不同的曲线上限值代入模型进行试估计，选择拟合优度R2最高时对应的上限值作为龚珀兹曲线趋势模型上限值。经测试，选择520作为郸城冬小麦单产上限值，即L值。模型参数通过显著性检验，且DW检验表明，模型残差序列存在一阶自相关。残差序列单位根检验平稳，可以建立ARMA模型。根据残差序列自相关分析图，建立ARMA(1,1)。ARMA模型残差序列为白噪声，模型可用于拟合。将龚珀兹曲线模型预测值与ARMA(1,1)混合模型结果叠加得到1990-2015年趋势产量值，进而利用公式(4)求出减产率见表1。使用同样的方法求得其余101个县的历年的减产率。

3.3 模型构建与精度验证

利用混合时间序列模型对河南省102个县区1990-2015年历年冬小麦单产进行拟合，得到趋势产量及对应的减产率。为保证样本数量，在0.6的雨养占比的阈值下，利用游程理论和河南省冬小麦生育期对多尺度SPEI进行影响冬小麦生长的样本筛选，最后利用Python中的Scikit-learn下的随机森林回归模型进行建模,利用网格寻参法寻找最优参数，最终确定弱学习器最大迭代次数n_estimators为55,最大特征数max_features为6。

研究中设定训练集和测试集的样本数据比例为8:2。精度验证结果见表2，训练集的真实减产率与多尺度SPEI之间的相关系数达到0.929，测试集相关系数达到0.703。训练集的平均绝对值误差为0.031，测试集平均绝对值误差为0.069。训练集的均方根误差为0.048，测试集的均方根误差为0.083。

表2 训练集和测试集模型精度分析

图5反映训练样本和测试样本的绝对误差区间不同，样本个数也有所差异。当绝对误差小于0.01时，训练样本和测试样本的个数最大。由训练样本和测试样本的累积频率曲线可以看出80%的训练样本的绝对误差在0.05之前，60%的测试样本的绝对误差集中在0.05之前。

图5 训练集与测试集绝对误差直方图与累积频率图

4 结论和讨论

本文以河南省雨养为主的县区的冬小麦为研究对象，利用混合时间序列模型计算出减产率，而后利用游程理论和河南省生育期对多时间尺度的SPEI进行影响冬小麦生长的干旱样本数据的筛选，并进行河南省干旱特征分析，最后建立减产率与筛选出的SPEI之间的旱灾损失模型，使用平均绝对值误差和均方根误差对灾损模型进行了精度验证。结果表明，河南省耕地范围上不同时间尺度的SPEI总体呈现干湿一致的趋势。河南省冬小麦全生育期识别出的干旱历时和干旱强度年均值的分布具有一致性。在雨养比阈值为0.6时构建出的随机森林灾损模型的验证结果显示平均绝对误差为0.069，均方根误差为0.083，能够满足产量损失预评估的要求，可以为产量损失快速评估提供参考依据。

在建模的样本筛选方面，本文重点关注的是雨养小麦的灾损模型，雨养样本的选择对于模型的构建必然存在影响，而河南省各个县都有效灌溉面积，为此我们设定了阈值来区分雨养为主和灌溉为主的样本。然而阈值设定过低受灌溉小麦的影响，模型拟合度可能会过低，阈值设定过高，有可能会使得参与建模的样本量不足。最终本文设定了0.6为阈值。

在农业生产中，只要前期天气条件为农作物提供充分的降水等条件，土壤水分满足作物生长需要，则认为该地区一般不易发生干旱[23]。本文利用游程理论识别干旱事件时关注农作物的生育期，只要识别出有一次干旱事件发生在生育期内就作为备选样本，在进行模型构建时，选择对真实减产率大于3%的数据和生育期内全部月份的干旱指数SPEI进行建模。但是在作物生长季内可能发生多次不连续的干旱事件，最终的产量损失是多次干旱事件的累加结果，基于历史数据无法准确的验证不同生育期的灾损模型。此外，作物对不同程度的干旱胁迫的响应是一个复杂物理化学过程[24]。本文识别的干旱长度最小为1个月，不能反映历时小于1个月的干旱对冬小麦产量的影响。本文在模型的输入方面，主要考虑的是在季节尺度以下冬小麦生育期不同月份的SPEI的值，但SPEI并不是构成引起冬小麦产量损失的唯一因素。影响冬小麦产量损失的过程和因素是综合且复杂，不仅仅与致灾因子有关，还与受灾地区的地理条件、社会经济以及防灾抗灾能力有关。因此在进一步的工作中，可以尝试细分生育期内不同月份的不同尺度的SPEI值，再结合土壤肥力，人为投入等更加多元的要素构建灾损模型，使研究更进一步深入。