基于随机森林回归的玉米单产估测

王鹏新 齐 璇 李 俐 王 蕾 许连香

(1.中国农业大学信息与电气工程学院， 北京 100083; 2.农业农村部农业灾害遥感重点实验室， 北京 100083;3.中国农业大学土地科学与技术学院， 北京 100083)

0 引言

作物长势的动态监测及产量的准确估测，能够为农业经营者的田间管理和国家粮食政策的制定提供有效支撑[1-2]。近年来，随着遥感技术的迅速发展，大范围、多维空间的作物长势监测和产量估测成为可能。目前，经验回归模型是作物产量估测的常用方法之一[3]。

经验回归模型通常选取与作物产量密切相关的特征参数进行估产。在此类研究中，植被指数(Vegetation index, VI)应用广泛[4]。任建强等[5]以美国玉米为研究对象，以各州为估产区，通过筛选的归一化植被指数(Normalized difference vegetation index, NDVI)与玉米单产间的最佳估产模型对2011年各州玉米单产进行了估算，并推算全国玉米单产，结果表明，全国玉米单产的相对误差仅为2.12%。王恺宁等[6]选取Landsat 8 OLI卫星遥感数据，计算冬小麦灌浆期归一化植被指数、比值植被指数(Ratio vegetation index, RVI)、绿度植被指数(Greenness vegetation index, GVI)和增强植被指数(Enhanced vegetation index, EVI)4种植被指数，并与冬小麦单产建立单植被指数和多植被指数的神经网络和SVM模型，结果表明，多植被指数SVM模型的估产精度高于神经网络模型。LIAQAT等[7]以巴基斯坦整个印度河流域为研究区域，通过多种植被指数，如土壤调整植被指数(Soil adjusted vegetation index, SAVI) 和改良土壤调整植被指数(Modified soil adjusted vegetation index, MSAVI)等，与小麦单产建立逐步回归模型，结果表明SAVI与小麦单产的决定系数R2和Pearson相关系数分别为0.74和0.88。然而，作物单产除与植被指数相关外，还与土壤含水率和生长状态密切相关[8]。因此，可通过综合作物生长过程中的水分胁迫指标和生长状态指标提高作物单产估测精度。其中，条件植被温度指数(Vegetation temperature condition index, VTCI)是基于归一化植被指数和地表温度(Land surface temperature，LST)的散点图呈三角形的基础上提出的[9]，可用于定量化地表征作物水分胁迫信息，并已成功应用于陕西省关中地区干旱监测及冬小麦单产估测预测等[10-12]。叶面积指数(Leaf area index，LAI)可表征植物的生长状态和光合作用能力，是作物长势监测及单产估测的重要指标[13]。此外，不同生育时期发生水分胁迫对作物单产的影响程度不同[14]，可通过赋予不同生育时期特征变量不同的权重，构建综合特征参数进行作物单产估测以提高估测精度。王鹏新等[15]利用重采样粒子滤波算法同化VTCI和LAI，并基于组合熵的方法构建加权VTCI和LAI与冬小麦单产的线性回归模型，结果表明不同管理模式下影响冬小麦单产的主要因子不同。

随机森林(Random forest，RF)回归模型是一种流行的机器学习模型，具有抗过拟合和预测精度高的特点[16-17]。应用随机森林回归估测作物单产(尤其是通过综合指数估测单产)的研究相对较少。因此，本文以河北省中部平原地区为研究区域，选取条件植被温度指数和叶面积指数为特征变量，通过随机森林回归算法获取玉米主要生育时期各个特征变量的权重，进而构建加权特征变量与玉米单产间的回归模型，以期为作物长势监测及单产估测提供新思路。

1 材料与方法

1.1 研究区概况

河北省中部平原处于东经114°32′～117°36′，北纬36°57′～39°50′之间(图1)，包括石家庄市、保定市、廊坊市、衡水市和沧州市的部分或全部地区，包含53个县(区)。该区域属暖温带大陆性季风气候，四季分明，降水集中，是华北平原重要的农业生产区之一。该地区年降水量在350～700 mm之间，且时空分布不均，降水主要集中在夏季，占全年的65%～70%，降水量由南向北逐渐减少。冬小麦-夏玉米轮作是该地区的主要耕作制度，该地区夏玉米出苗到拔节期一般在7月上旬至7月中旬、拔节到抽雄期在7月下旬至8月上旬、抽雄到乳熟期在8月中旬至9月上旬、乳熟到成熟期在9月中旬至9月下旬。通过王鹏新等[18]提出的基于时间序列叶面积指数傅里叶变换的作物种植区域提取方法提取了2010—2018年研究区域玉米种植区。

图1 研究区域位置及玉米种植区(2010年)Fig.1 Location of study area and planting area of maize(2010)

1.2 数据获取与处理

1.2.1时间序列VTCI和LAI生成

选取2010—2018年每年7—9月Aqua-MODIS日地表温度产品MYD11A1及日地表反射率产品MYD09GA，经MRT预处理后获得研究区域日LST和日NDVI产品，应用最大值合成技术生成每年7—9月旬时间尺度的NDVI和LST最大值合成产品，基于多年某一旬的NDVI和LST最大值合成产品，运用最大值合成技术分别生成多年的旬NDVI和LST最大值合成产品，基于每年7—9月旬LST最大值合成产品，运用最小值合成技术生成多年的旬LST最大-最小值合成产品。VTCI取值范围为0～1，其值越接近0，表明越干旱，作物受水分胁迫程度越重，其值越接近1，表明越湿润，作物受水分胁迫程度越轻或不受水分胁迫，VTCI计算公式为

(1)

其中

Lmax(Ni)=a+bNi

(2)

(3)

式中L(Ni)——在研究区域内，某一像素的NDVI值为Ni时的地表温度

Lmax(Ni)、Lmin(Ni)——研究区域当NDVI值为Ni时所有像素地表温度最大值和最小值

a、b、a′、b′——待定系数，由研究区域LST和NDVI的散点图近似得到

选取研究区域2010—2018年每年7—9月MODIS叶面积指数产品MCD15A3H，该产品是基于Terra和Aqua卫星上的MODIS传感器获得的，与MOD15A2和MYD15A2产品相比，MCD15A3H产品既有较高的时间分辨率(4 d)又有较高的空间分辨率(500 m)，有利于作物长势监测及产量估测。利用MRT对产品进行预处理得到研究区域叶面积指数产品，原始叶面积指数产品由于云和大气等因素的影响存在数据骤降的现象，因此通过上包络线S-G(Savitzky-Golay)滤波对原始叶面积指数产品进行平滑处理[18]，经上包络线S-G滤波平滑处理后的叶面积指数更加符合玉米生长情况。为使LAI与VTCI具有相同的时间尺度，将玉米各旬所包含的多时相LAI的最大值作为各旬的LAI值，并对上包络线S-G滤波后的LAI进行归一化处理，最大值为7，最小值为0。

1.2.2VTCI和LAI计算

依据玉米4个主要生育时期的划分，将玉米各生育时期包含的多旬VTCI和LAI的平均值作为该生育时期的VTCI和LAI值，如将7月上旬至7月中旬VTCI的平均值作为出苗到拔节期的VTCI值。再叠加研究区域行政边界图，将各县(区)包含的所有像素的VTCI和LAI的平均值作为该县(区)的VTCI和LAI值。以此类推，计算得到研究区域2010—2018年各县(区)玉米各生育时期的VTCI和LAI值。

1.2.3玉米单产数据的来源及异常数据处理

通过查阅《河北农村统计年鉴》得到研究区域各县(区)2010—2016年玉米播种面积和总产量数据，玉米单产由总产量和播种面积计算得到。

将VTCI和LAI与玉米单产进行回归分析的残差的置信区间在[-4 000,4 000] kg/hm2以外的单产数据视为异常数据，在构建模型时将其剔除。

1.3 研究方法

随机森林回归对噪声数据集容忍度较高，对高维数据集具有良好的预测能力[19-20]。它是由一组没有联系的回归决策树{h(x,θk)，k=1,2,…,K}构成的K棵集成决策树，表示为

(4)

式中x——各县(区)玉米4个生育时期VTCI或LAI值及玉米单产数据

K——决策树的数量

θk——独立同分布随机向量

为了提高模型的预测精度并防止出现过拟合情况，以随机森林回归算法结合袋装法得到训练样本子集，并结合随机子空间法得到节点分裂特征[21]。

(1)袋装法通过有放回地随机抽样，从原始样本数据集中重复抽样得到K个与原始样本数据集相等的训练样本N，每个训练样本构成一棵决策树。每次进行Bootstrap重抽样时，未被抽中的样本的概率为(1-1/N)N，当N趋向于无穷大时，未被抽中样本的概率越接近1/e，约为0.368，即原始样本中有36.8%的数据未被抽中，这些数据被称为袋外数据(Out of bag, OOB)，因其未参与回归树的构建，故可用来估计预测袋外数据误差(OOB误差)及评估自变量对因变量的影响程度。另外，基于OOB预测误差可以检验模型的泛化能力，不需再使用测试集检验模型的精度。通过袋装法得到的K个训练样本都不相同，保证了回归树的差异性。

(2)随机子空间法通过袋装法得到K棵回归树后，每个分裂节点随机抽取所有变量(特征)中的Mtry个变量(特征)作为当前节点分裂的特征子集，根据分类回归树(Classification and regression tree，CART)方法在特征子集中选择最优分裂方式进行分裂。通过随机子空间法得到的回归树具有随机性和独立性。在随机森林回归中，树的数量K和随机选择的节点分裂变量(特征)Mtry决定着模型的预测能力。

图3 OOB误差随回归树数量的变化曲线Fig.3 Changing curves of OOB errors with number of regression trees

基于随机森林回归估测玉米单产的流程如图2所示。

图2 基于随机森林回归估测玉米单产的流程图Fig.2 Flow chart for estimating maize yield based on random forest regression

(1)将研究区域各县(区)2010—2016年玉米4个生育时期的VTCI或LAI值及玉米单产数据作为原始样本(共357组数据)输入模型，通过Bootstrap重抽样得到K个训练样本子集并生成K棵回归树。VTCI和LAI估测玉米单产的OOB误差随树的数量K变化曲线如图3所示，可以看出，当K为500时，OOB误差趋于平稳，故将K设为500。

(3)每棵回归树自上向下分裂生长，直到到达某个叶子节点输出估测值，所有回归树构成随机森林。将所有回归树输出的玉米单产求平均值即可得到最终的玉米单产估测结果。

随机森林回归模型不但能精确地估测玉米单产，而且还可给出各个变量的重要性评分，即玉米4个生育时期VTCI或LAI对玉米单产的影响程度。基于基尼系数和基于OOB误差是常用的变量重要性评分的统计量，本研究中基于OOB估测误差得到各变量的重要性。若xj(j=1,2,3,4)为输入变量，则在第k棵树上的重要性Ik为随机置换变量前后袋外数据估测误差的差值[22]。其计算公式为

(5)

变量xj在整个随机森林中的重要性得分为

(6)

式中NOOB——袋外数据样本数

f(xn)——袋外数据中第n个样本值

fk(xn)、fk(x′n)——随机置换变量前后第k棵树上的袋外数据第n个样本的估测值

I(·)——判别函数，当f(xn)=fk(xn)或f(xn)=fk(x′n)时，取值为1，否则为0

由于随机性的引入，模型每次给出的变量重要性评分略有差异，故将10次运行结果的平均值进行归一化处理，作为各个变量的权重。

1.4 估产模型构建

通过随机森林回归方法确定玉米主要生育时期VTCI和LAI的权重，计算2010—2018年各县(区)加权VTCI和LAI。对2010—2016年(除2012年，用来进行精度验证)加权VTCI和LAI与玉米单产进行回归分析，选取拟合程度最优的回归模型对2012年各县(区)的玉米单产进行估测及精度验证，并基于该模型逐像素估测2010—2018年研究区域的玉米单产。

2 结果与分析

2.1 玉米各生育时期的权重

基于随机森林回归模型运行10次输出的各变量重要性的平均值进行归一化处理，得到玉米各生育时期VTCI和LAI的权重(表1)。可以看出，玉米拔节—抽雄期和抽雄—乳熟期的VTCI权重较大，说明受水分胁迫时对玉米单产的影响程度相对较大，主要是因为这两个时期对水分胁迫较敏感，抽雄期前后发生水分胁迫会导致幼穗发育不良，果穗偏小，雄穗在抽出2～3 d后失去散粉能力，甚至有的雄穗不能抽出，或抽穗时间延迟，导致秃尖增长，造成不同程度的玉米产量下降，水分胁迫较重的会造成雌穗部分不育甚至空秆。苗期—拔节期和乳熟—成熟期的VTCI权重相对较小，说明发生水分胁迫对玉米单产的影响较小，主要是苗期发生一定程度的水分胁迫会使根向下生长，有利于玉米植株后期的生长发育，且后期有充足水分时能够弥补之前减少的生长量，乳熟期之后穗粒已经形成，受水分影响不大[23]。LAI对玉米单产的影响以抽雄—乳熟期和乳熟—成熟期较大，苗期—拔节期和拔节—抽雄期较小，表明生长前期LAI与玉米产量的相关性不大，主要是因为光合作用的产物用来进行以根系和叶片为中心的营养生长，抽雄期时LAI达到最大，玉米进入以果穗为中心的生殖生长阶段，LAI与产量的相关性开始增大，这与姚小英等[24]的研究结果较一致。

表1 玉米各生育时期的权重结果Tab.1 Weight results of each growth stage of maize

2.2 单产估测模型选择

将随机森林回归方法计算得到的2010—2016年(除2012年)加权VTCI和LAI与玉米单产基于县域尺度进行线性回归分析，建立不同变量的单产估测模型(表2)。结果表明，基于随机森林回归的加权VTCI和玉米单产的相关性最低(R2=0.001)，且没有通过显著性检验；加权LAI与玉米单产的相关性次之(R2=0.296)；加权VTCI和LAI与玉米单产的相关性最高(R2=0.303)，模型达极显著水平(P<0.001)，表明VTCI和LAI与玉米单产呈显著的正相关关系。因此，基于双变量估产模型的精度高于单变量模型的精度。基于随机森林回归双变量估产模型估测玉米单产时，玉米单产受LAI影响较大，VTCI影响较小，原因可能是研究区域受人为因素的影响较大，当发生水分胁迫时，通过及时灌溉缓解了当地旱情，致使玉米单产对VTCI不敏感。综上所述，基于随机森林回归的双变量估产模型精度最高，可用于估测研究区域2012年各县(区)的玉米单产。

表2 加权VTCI和LAI与玉米单产间的线性回归分析Tab.2 Linear regression analysis between weighted VTCI and LAI and maize yields

2.3 2012年各县(区)玉米单产估测的精度评价

基于随机森林回归双变量估产模型及2012年加权VTCI和LAI对各县(区)玉米单产进行估测(表3)。玉米估测单产与实际单产的相对误差以清苑区最低，为0.35%，以海兴县最高，为37.10%。其中，31个县(区)玉米估测单产与实际单产的相对误差在10%以下，7个县(区)在10%～15%，15个县(区)在15%以上，53个县(区)的平均相对误差为9.85%，均方根误差为824.77 kg/hm2。个别县(区)如海兴县、盐山县的相对误差较大，原因可能是海兴县、盐山县濒临渤海，土壤盐渍化严重，农业生产条件较差，农田水利设施建设和机械化水平较低，不适宜种植经济作物，种植冬小麦和夏玉米是仅有的选择。近年来当地已采取改造重盐碱地的相关措施使玉米单产有所提高，但是玉米生产仍处于较低水平，玉米单产被高估，从而使估测单产与实际单产的相对误差较大。个别县(区)如正定县、藁城区和新乐市实际玉米单产较高，估测单产偏低，玉米单产被低估，原因可能是这几个县(区)是国家粮食丰产科技工程河北省项目区的核心区，田间管理及时，玉米单产受人为因素影响较大。

表3 2012年各县(区)玉米估测单产Tab.3 Estimated yields of maize in each county (district) in 2012

为了进一步验证随机森林回归双变量估产模型的精度，基于2012年各县(区)玉米实际单产与估测单产进行线性回归分析。结果表明，估测单产与实际单产间呈显著的正相关关系(P<0.001)，R2达到0.540；估测单产与实际单产的均方根误差为631.64 kg/hm2，进一步说明基于随机森林回归双变量估产模型的精度较高，可用于研究区域玉米单产估测。

2.4 玉米估测单产的时空变化规律

图4 基于随机森林回归的玉米单产估测结果Fig.4 Estimate results of maize yields based on random forest regression

基于随机森林回归双变量估产模型逐像素估测2010—2018年研究区域玉米单产(图4)，并逐像素统计玉米估测单产。结果表明，2010、2012、2013年玉米估测单产相差不大，西部地区(包括石家庄市和保定市)玉米估测单产在6 600 kg/hm2左右，东部地区(包括沧州市)在6 100 kg/hm2左右，南部地区(包括衡水市)在6 800 kg/hm2左右，北部地区(包括廊坊市)在6 200 kg/hm2左右；2011年玉米估测单产略高于2010年；2014年玉米估测单产略低于2013年；2015、2016、2017年玉米估测单产略高于2014年，西部地区在6 900 kg/hm2左右，东部地区在6 100 kg/hm2左右，南部地区在7 000 kg/hm2左右，北部地区在6 100 kg/hm2左右，2018年西部地区和南部地区玉米单产在7 000 kg/hm2左右，东部地区和北部地区在6 500 kg/hm2左右。以2017年为例，西部地区玉米估测单产为6 868 kg/hm2，东部地区为6 051 kg/hm2，南部地区为6 833 kg/hm2，北部地区为6 045 kg/hm2。研究年份间2011年玉米单产最高，2014年玉米单产较低，原因可能是2011年降水量充沛，玉米单产高于常年，2014年玉米生育期内发生阶段性干旱且局部地区旱情较重，玉米单产下降。

3 讨论

课题组在陕西关中平原的冬小麦干旱监测及单产估测中采用客观赋权法如熵值法确定VTCI的权重[25]，构建的加权VTCI和冬小麦单产的回归模型精度较高，但熵值法基于指标的差异程度确定指标权重，异常数据对权重影响较大，且可能使权重与实际相背，因此确定冬小麦主要生育时期VTCI的权重与实际水分胁迫对冬小麦单产的影响程度不符。在河北省中部平原地区应用随机森林回归确定玉米主要生育时期VTCI和LAI的权重，结果表明随机森林回归确定的VTCI权重以拔节—抽雄期、抽雄—乳熟期的权重较大，根据实际水分胁迫对玉米单产的影响程度[26]可以看出,基于随机森林回归的权重结果更加符合实际情况。主要因为干旱对玉米单产的影响具有非线性的特征，随机森林回归模型对于非平衡数据比较稳健，不易受到异常值的干扰，能有效处理非线性问题。虽然基于随机森林回归确定的玉米主要生育时期VTCI和LAI的权重较合理，但是未考虑农学先验知识，可通过结合主观赋权法如改进的层次分析法进一步修正随机森林回归得到的权重，使权重更加符合实际情况。另外水分胁迫也会影响玉米的生长状态，即VTCI和LAI之间可能存在多元共线性的问题，而随机森林回归模型对多元共线性不敏感，可以很好地预测多个变量的作用，因此随机森林回归模型的精度较高。

影响玉米单产的因素有很多，除了受到水分胁迫和生长状态的影响外，还受到其他因素如温度、洪涝灾害、田间管理、玉米品种等的影响。杨笛[27]通过模拟气候变化、肥料、种植面积、灌溉和品种5个驱动因子对黄淮海夏玉米区玉米单产的影响，表明肥料和品种在玉米增产中的作用和地位随时间在提高，种植面积的增长及灌溉系数的减少不利于玉米增产。通过查阅《河北农村统计年鉴》可以看出,研究年份间灌溉和肥料的使用较多，这可能与研究区域玉米高产有一定的联系。另外，个别年份发生灾害如2016年研究区域部分县(区)玉米苗期发生雹灾，影响玉米出苗，7月又出现涝灾和病虫害使玉米单产略有下降。这些因素对玉米单产的影响不容忽视，综合考虑与玉米单产相关性较大的因素是今后研究的重点。

4 结论

(1)通过随机森林回归确定玉米主要生育时期VTCI和LAI的权重，构建加权VTCI和LAI与玉米单产的单变量和双变量估产模型。结果表明，基于随机森林回归的双变量估产模型精度最高。

(2)基于随机森林回归双变量估产模型估测2010—2018年研究区域玉米单产，结果表明，玉米估测单产在空间上的分布特征为西部地区最高、北部和南部次之、东部最低，年际间的分布特征为在波动中呈先减少后增加的趋势。