基于Shapley值组合预测的玉米单产估测

王鹏新 乔 琛 李 俐 周西嘉 许连香 胡亚京

(1.中国农业大学信息与电气工程学院， 北京 100083； 2.中国农业大学土地科学与技术学院， 北京 100083)

0 引言

粮食的生产与安全对提升我国农业的经营管理水平、完善农作物的种植结构和确保粮食安全等具有重要意义[1-3]，因此及时准确估测农作物的产量有利于维持国家稳定和促进经济发展。

干旱作为影响农作物生长发育及产量的重要因素已受到广泛的关注[4-5]。王鹏新等[6]在归一化植被指数(Normalized difference vegetation index，NDVI)和地表温度(Land surface temperature，LST)的散点图呈三角形区域分布的基础上，提出了条件植被温度指数(Vegetation temperature condition index，VTCI)的干旱监测方法，可用于反映农作物生长过程中的水分亏缺信息。孙威等[7]对VTCI冷热边界的确定方法进行了完善，并验证了利用 VTCI 进行干旱监测的可行性。农作物产量不仅受到水分胁迫的影响，还与作物的生长状态密切相关[8]，叶面积指数(Leaf area index，LAI)通过对作物的光合速率和干物质的积累量的反映能较好地表征农作物的生长状态[9]。

遥感技术对比传统统计调查方法，凭借其覆盖范围广、重访周期短等独特优势被广泛应用于作物长势监测及产量估测，同时对大规模的农业生产调查、评价、监测和管理具有独特的作用[10]。随着大数据技术的发展，遥感大数据将推动农业遥感估产的发展[11]。目前，遥感估产方法中常采用的统计模型、机理模型和半机理模型等均能够较好地对作物进行估产[12]。但由于在实际应用中机理和半机理模型存在需要输入较多的参数问题，因此机理和半机理模型存在一定的局限性。而统计模型其估产的精确度依赖于选取遥感影像的时相，对作物的生长和产量形成的机理解释性不强[12]，因此在实际应用中同样具有一定局限性。在作物生长过程中，经常受到各种因素影响，同时，这些因素在作物不同生育时期产生不同的影响，即使采用相同的估算方式也会得到不同的结果与估测精度。虽然其估产精度不同，但是这些模型可以提供不同的有用信息，如将一些方法丢弃，就会失去有用信息，从而影响估产的精度[13]。因此，可以将单个估测模型提供的有用信息进行组合，以提高估产的精度。组合预测(Combination forecasting，CF)将不同预测模型进行有效组合，可视为对无限逼近真实数据生成过程的有效补充[14]。

BATES等[15]首次提出了组合预测理论，这一预测理论在提高精度的同时更充分利用了预测样本所表达的信息，受到了国内外广大学者重视。组合预测根据各单一模型的信息贡献度，进而基于计算得到的单一模型权重构建组合预测模型，从而实现减少预测误差、提高预测精度的目标[16]。特别是在时间序列中分析真实数据生成过程中，通常具有区制转换(Regime shift)或参数漂移等特性。组合预测方法的引入，可减少由参数错误或模型错误带来的预测误差[17]，甚至在单一预测结果存在有偏性的情况下，通过组合能产生具有无偏性的预测结果[18]。因此，组合模型具有普遍性的优点，最终的预测结果更接近实际数据。根据组合预测方法综合手段的不同，可分为权重综合法和区域综合法两种。与区域综合法相比，权重综合法得到了广泛的应用。目前，组合预测现在已经在多个领域内得到了应用[19-20]，但是在农业领域还少有报道。

本文以河北中部平原为研究区域，选取与玉米长势和产量密切相关的VTCI和LAI为特征变量，采用极限梯度提升树(Extreme gradient boosting，XGBoost)和随机森林(Random forest，RF)两种机器学习算法模型分别估测研究区域的玉米单产，并借鉴经济学中合作博弈论Shapley值利益分配理论，通过组合预测模型中的权重合成思想，以单一预测模型的均方误差为基础确定单一模型权重得到组合预测模型，以期为玉米长势监测和产量估测提供新方法。

1 材料与方法

1.1 研究区概况

河北中部平原属黄淮海平原组成部分，是我国重要的玉米种植和生产基地[21]，位于114°32′～117°36′E与36°57′～39°50′N之间，包括石家庄等5个市的53个县(区)，其土地面积约为5.30×104km2。河北中部平原地处典型温带大陆性季风气候区，四季分明、雨热同期。河北省中部平原年降水量350～800 mm。降水时空分布不均，南方的降水量比北方高，夏季降水量高于冬季，丰水年降水量与枯水年降水量相差较大。该地区的耕种管理制度为一年两熟，结合该地区玉米实际生长状况，本文将河北中部平原玉米生长划分为4个生育时期：7月上旬至7月中旬的出苗-拔节期、7月下旬至8月上旬的拔节-抽穗期、8月中旬至9月上旬的抽穗-乳熟期和9月中旬至9月下旬的乳熟-成熟期。根据王鹏新等[22]提出的作物分类方法，进而获得研究区域2010—2018年玉米种植区分布图，其中河北中部平原2017年玉米种植区分布图如图1所示。

1.2 数据获取与处理

1.2.1VTCI时间序列的生成

选取2010—2018年每年7—9月空间分辨率和时间分辨率分别为1 000 m、1 d的MODIS日地表温度产品(MYD11A1)及日地表反射率产品(MYD09GA)。利用MRT对日地表温度和日地表反射率产品进行预处理之后，得到研究区域的日LST和日NDVI产品，运用最大值合成法生成旬LST与NDVI的最大值合成产品；基于生成的多年某一旬的NDVI和LST最大值合成产品再次使用最大值合成技术生成多年的旬NDVI和LST的最大值合成产品；基于生成的多年某一旬的LST最大值合成产品利用最小值合成技术逐像素提取最小值，计算得到多年旬LST的最大—最小值合成产品，并以此通过计算生成河北中部平原旬VTCI时间序列数据。VTCI的计算公式为[6,23]

(1)

(2)

LSTmax(NDVIi)=a+bNDVIi

(3)

LSTmin(NDVIi)=a′+b′NDVIi

(4)

式中NDVI——归一化植被指数

ρNIR、ρred——近红外、红光波段的反射率

LST——地表温度

LSTmax(NDVIi)、LSTmin(NDVIi)——当NDVIi为某一特定值时所有像素地表温度最大值和最小值，即热边界和冷边界

LST(NDVIi)——研究区域内某一像素的NDVI值为NDVIi时的地表温度

a、b、a′、b′——LST和NDVI散点图近似得到的待定系数

1.2.2LAI时间序列的生成

选取2010—2018年7—9月的MODIS叶面积指数产品(MCD15A3H)，其空间分辨率和时间分辨率分别为500 m、4 d。利用上包络线Savitzky-Golay滤波对经MRT处理后得到的原始LAI产品进行平滑处理以消除云层覆盖、大气溶胶等因素引起的数据骤降的现象[23]，滤波处理后的叶面积指数变化趋于平稳且更加符合玉米的生长物侯特征，解决了原始数据存在的质量问题。由于LAI与VTCI时间分辨率和取值范围不同，因此对滤波后的叶面积指数进行归一化处理，并通过观察多时相MODIS的MCD15A3H原始数据相元统计直方图将取值范围设置为0～7，进而计算得到研究区域2010—2018年玉米主要生育时期的LAI时间序列数据。

1.2.3玉米生育时期VTCI和LAI的计算及异常点数据处理

基于研究区玉米生育时期的划分结果，将生育时期内所包含的多旬VTCI和LAI平均值作为研究区域玉米该生育时期逐像素的VTCI和LAI值；通过叠加河北中部平原各县(区)行政边界图，将各县(区)所包含的玉米生育时期逐像素VTCI和LAI值的平均值作为生育时期县(区)尺度的VTCI和LAI值。同时在构建回归模型时，剔除加权VTCI与LAI残差置信区间在[-4 000,4 000] kg/km2以外的异常点数据。

1.3 研究方法

1.3.1极限梯度提升算法

极限梯度提升算法(XGBoost)是一种基于梯度提升的决策树(Gradient boosting decision tree)集成算法[24]，其基分类器主要包含分类和回归树(Classification and regression tree，CART)。本文玉米单产估测是回归问题，所以其基模型选择为回归树。含有K棵决策树XGBoost的树集成模型定义为

(5)

xi——样本所对应的特征变量VTCI与LAI

fk——第k个决策树的预测函数

树集成优化模型可以定义为

(6)

(7)

l(yi,i)=(yi-i)2

式中l(yi,i)——损失函数，即均方误差

Ω(f)——正则化项

γ——复杂度参数

T——树中叶子节点个数

λ——固定系数

ω——叶子节点量化权重向量

求解式(6)的优化问题，即求解CART树的结构。通过保留训练好的前t-1轮树模型不变，在第t轮时添加一个新的预测函数，迭代计算得到最终预测结果[25]

(8)

ft(xi)——第t轮加入的新的预测函数

XGBoost模型可以通过输出玉米4个生育时期VTCI或LAI的特征变量重要性来评估不同特征变量对玉米产量的影响程度。在XGBoost中常用基于增益、覆盖度、频率的特征重要性指标进行特征重要性评价，其中基于增益(gain)表示各特征变量对XGBoost模型中每棵树采取每个特征变量的贡献而计算的模型相对贡献度为

(9)

式中G——VTCI与LAI基于增益获取的特征重要性得分值

GL、GR——左、右子树中的梯度

HL、HR——左、右子树中的二阶梯度

c——给定的临界值

1.3.2随机森林算法

随机森林(RF)是一种基于袋装法(bagging)理论实现的集成学习算法，其基模型是决策树模型。随机森林具有极好的准确率，同时可以评估各特征在回归问题上的重要性。随机森林将若干没有联系的回归决策树{y(x,θk),k=1,2，…,K}构成K棵集成决策树[26]，可以定义为

(10)

式中y(x)——估测研究区域玉米的单产

x——特征变量的输入，即VTCI与LAI

文献[26]基于RF算法确定了河北中部平原玉米各个生育时期的权重，进而构建了加权VTCI和LAI与玉米单产之间的双参数估产模型，结果显示，双参数估产模型精度相对较高，达到极显著的水平(P<0.001)，且基于随机森林算法确定的玉米各生育时期权重均较为合理。因此，本文引用其基于随机森林算法确定的河北中部平原玉米各生育时期的VTCI和LAI权重进行估产模型的构建。

1.3.3组合预测模型

XGBoost算法与RF算法均为统计理论的机器学习方法，但二者各有优势。XGBoost算法通过引入正则化项在很大程度上避免了过拟合问题的出现，同时采用优化算法降低了问题计算复杂度，在处理大规模数据集方面具有明显的优势，已被广泛应用于不同领域的研究[27-29]。相对于XGBoost，RF更适合处理高维数据，对噪声数据的容忍性较高。通过结合XGBoost算法和RF算法的优点，参考经济学中合作博弈论的Shapley值理论，通过计算确定组合预测模型总误差在单一预测模型之间的分布，并在此基础上确定各单一预测模型的权重，具体步骤如下：

(1)求取Shapley值

(11)

(12)

式中i——预测模型序号，即XGBoost与RF预测模型序号

E′i——第i个预测模型的Shapley值，即XGBoost与RF预测模型计算得到的均方误差

s——包含模型XGBoost与RF的所有集合

|s|——预测模型个数

E(s)——集合s的收益

E(s-{i})——集合s中去除成员i后的收益

(2)计算权重

组合预测模型中若某单一预测模型所获得的误差分配值越大，表示预测精度越低，在组合预测模型中的权重就越小。基于此原则，预测模型i的权重λi定义为

(13)

最终的组合预测模型可描述为

(14)

式中Y——组合预测模型估测的玉米单产

Yi——模型XGBoost与模型RF估测的玉米单产

以通过计算得到的2010—2017年玉米4个生育时期的VTCI和LAI数据作为特征变量，相对应的玉米单产数据作为目标变量(每年53组数据，共424组数据)，以2010—2017年(除2012年)数据作为训练集合，2012年数据作为测试集合。由于特征变量与目标变量之间差值较大，特征变量值的范围在0～1之间，而目标变量值在3 000～9 000 kg/km2之间。由于特征变量和目标变量处于不同区间的值域差异可能对结果造成不同的影响，取对数后不会改变数据的性质和相对关系，对训练数据集和测试数据集中的目标变量进行取对数处理，并通过对估产模型输出的估测数据进行取对数还原，得到最终的估测单产数据。采用均方根误差(RMSE)、平均相对误差(MRE)和决定系数(R2)等指标对模型的精度进行评价。

2 结果与分析

2.1 估产模型构建

构建XGBoost估产模型与RF估产模型，借鉴Shapley值理论，计算得到XGBoost与RF估产模型的Shapley值。根据Shapley值确定XGBoost估产模型与RF估产模型权重，进而完成对组合预测模型的构建。

2.1.1XGBoost与RF的估产模型构建

基于XGBoost与RF算法输出的特征重要性值，通过归一化计算得到研究区域玉米各生育时期VTCI与LAI的权重(表1)。可以看出，玉米生育后期的LAI权重高于玉米生育前期LAI的权重，表明玉米生育前期的LAI对玉米产量影响程度较弱，玉米生育后期的LAI对玉米产量影响程度较强。原因可能是玉米LAI变化规律呈前期增长缓慢、中期增长快速、后期下降缓慢的偏峰曲线，其中玉米在出苗-拔节期和拔节-抽雄期主要以分化茎叶的营养生长为主，此阶段玉米生长迅速，叶片迅速增多增大；以抽雄-乳熟期为界，玉米进入以生殖生长为主的生育后期，光合产物的分配模式主要以果穗为中心，是玉米产量形成的重要时期。因此玉米生育后期的LAI对玉米产量的影响程度强于玉米生育前期。拔节-乳熟期的VTCI权重高于出苗-拔节期和乳熟-成熟期的权重，表明拔节-乳熟期的VTCI对玉米产量的影响程度较强，出苗-拔节期和乳熟-成熟期的VTCI对玉米产量的影响程度较弱。原因可能是在拔节-抽雄期、抽雄-乳熟期玉米进入营养生长的高峰期，此阶段玉米生长迅速，对土壤中的水分吸收也最为急迫，若发生水分胁迫将减缓玉米根茎叶的生长发育，降低光合作用对玉米干物质积累速率，减少蛋白质和有机质的合成，造成玉米粒重明显降低，最终影响玉米产量，因此拔节-抽雄期、抽雄-乳熟期的VTCI对玉米产量的影响程度较强。而在出苗-拔节期，玉米由于植株矮小，对水分的需求量较少，在乳熟-成熟期，玉米处于生育后期，生长变缓，对一定的水分胁迫表现出一定的忍受力，因此出苗-拔节期、乳熟-成熟期VTCI对玉米产量的影响程度较弱。综上所述，基于XGBoost算法与RF算法确定的研究区域玉米各生育时期的权重均较为合理。

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

分别将XGBoost估产模型与RF估产模型得到的玉米各生育时期的特征权重进行加权VTCI和LAI计算，进而构建基于加权VTCI和LAI与玉米单产之间的回归估测模型，将2010—2017年(除2012年)数据代入XGBoost估产模型与RF估产模型，进行可视化分析(图2)。可以看出，XGBoost估产模型R2为0.31，均方根误差为940.91 kg/km2，RF估产模型R2为0.30，均方根误差为947.50 kg/km2，XGBoost估产模型与RF估产模型均通过显著性检验(P<0.001)。结果表明，XGBoost估产模型与RF估产模型玉米单产估测精度均较高，可用于研究区域各县(区)玉米的单产估测。

2.1.2组合预测模型的构建

为进一步提高玉米单产估测精度，进行组合预测模型的构建。需要确定XGBoost与RF估产模型的权重，基于2010—2017年(除2012年)由单一估产模型输出的各县(区)的估测单产数据(共371组)，分别计算XGBoost估产模型、RF估产模型和组合预测模型的均方误差，将各模型的均方误差代入式(11)，将得到的XGBoost与RF估产模型的Shapley值代入式(13)，则可以得到XGBoost与RF估产模型的权重分别为0.48与0.52，则组合预测模型为

Y=0.48YXGBoost+0.52YRF

(15)

通过计算，组合预测模型的R2为0.32，模型通过显著性检验且达显著水平(P<0.001)。同时，组合预测模型均方根误差为939.81 kg/km2，与单一XGBoost估产模型与RF估产模型相比，组合预测模型决定系数与均方根误差均得到提升。结果表明，组合模型的玉米单产估测精度优于单一估产模型，可用于河北中部平原的玉米单产估测。

2.2 玉米估测单产的精度评价

将2012年数据分别代入XGBoost估产模型、RF估产模型与组合预测模型，将玉米的的估测单产与实际单产数据进行数据可视化分析(图3)。对比分析，可以看出XGBoost估产模型、RF估产模型和组合预测模型的估测单产与实际单产均呈显著的正相关关系。同时，R2均不小于0.50，可以反映实际单产的波动均有50%以上能被估测单产的波动所描述，即玉米的实际单产与估测单产之间的误差较小。其中，组合预测模型R2达到0.52，均方根误差为831.14 kg/km2，平均相对误差为9.86%，对比单一XGBoost估产模型与RF估产模型，组合预测模型的决定系数、均方根误差与平均相对误差均得到提升。进一步表明，基于Shapley值的组合预测模型的估产精度优于单一估产模型。

2.3 区域玉米单产的估测与时空变化分析

将组合估产模型应用于河北中部平原2010—2018年逐像素玉米单产估测(图4)。结果表明，研究区域玉米的估测单产随年际变化呈现先减少后增加的波动变化趋势。其中，2010—2014年玉米估测单产整体呈现逐年下降的趋势，并在2014年玉米平均单产达到最低，平均单产在6 000 kg/km2左右，原因可能是2014年河北中部平原降水较少且发生阶段性干旱，导致玉米单产减少严重；2010、2011、2013年平均估测单产相差不大，均在6 350 kg/km2左右。对比2010—2014年玉米估测单产年际变化趋势，2015—2018年的玉米单产整体呈逐年上升的趋势，其中2018年玉米平均估测单产最高，约为6 500 kg/km2，2016年和2017年玉米平均估测单产次之，约为6 400 kg/km2，2015年玉米平均估测单产达到最低，约为6 200 kg/km2。

研究区域玉米的单产估测空间上呈现西部地区玉米估测单产最高，南部和北部地区玉米单产次之，东部地区玉米单产最低的分布特征。河北中部平原的北部地区中，2012年玉米的估测单产最高，2014年玉米估测单产最低，分别约为6 000 kg/km2和4 500 kg/km2，2015、2016、2017年玉米估测单产相差不大，均约为5 000 kg/km2，2011年和2018年玉米的估测单产略低于2012年；东部地区中，2016年玉米的估测单产最高，2014年玉米的估测单产最低，分别约5 000 kg/km2和4 500 kg/km2，2016—2018年东部地区的玉米估测单产相差不大，均约为4 700 kg/km2。南部地区玉米的估测单产分别在2012年和2014年达到最高和最低，分别约6 500、4 500 kg/km2，其余年份玉米估测单产相差不大。西部地区中，2012、2015年玉米的估测单产较高且相差不大，均约为6 700 kg/km2，2010、2014年玉米估测单产较低，均约为4 500 kg/km2。

3 讨论

经过对不同估产模型特点的研究分析发现，在玉米单产估测过程中，不同估产模型对特征因素的信息提取方式不完全一样。同时以往单一估产模型均注重特征因素对产量的影响，而缺少考虑估产模型对特征因素的信息提取以及对估测产量的影响。本文在考虑VTCI和LAI作物长势监测指标的基础上，也充分考虑了估产模型对特征因素的信息提取以及对估测产量的影响，利用RF对噪声数据的容忍性较高以及XGBoost较低的计算复杂度等优点，借鉴经济学中Shapley值理论计算得到XGBoost和RF估产模型的权重，进而构建组合估产模型，实现了单一估产模型的综合利用。在今后的研究中，可以在Shapley值理论的基础上尝试加入其他估产模型，以期进一步全面地提取特征因素的信息，避免单一估产模型在估产过程中造成有用信息的浪费，从而实现玉米单产估测精度的进一步提升。

本文虽然选取了与玉米长势和单产密切相关的VTCI和LAI作为特征变量，但玉米的生长发育除了受到水分胁迫和生长状态的影响之外，还受到其他自然因素和人为因素的影响，例如温度、土壤肥力、田间管理水平等因素。因此在基于VTCI与LAI作物生长指标的基础上，未来研究应进一步综合考虑与玉米单产相关性较大的其他因素。此外，基于Shapley值构建的组合预测模型对河北中部平原各县(区)玉米单产估测精度虽较高，但缺少对农学先验知识的考虑，今后研究中可以通过主观赋权法进一步修正XGBoost估产模型和RF估产模型权重，使模型权重更优。

4 结论

(1)通过借鉴组合预测模型思想，利用经济学合作博弈论Shapley值利益分配理论确定了XGBoost估产模型与RF估产模型的权重，进而构建组合预测模型。结果表明，基于Shapley值的组合预测模型估测精度较高，且优于单一估产模型的精度。

(2)将组合预测模型应用于研究区域2010—2018年逐像素的玉米单产估测。结果表明，从研究区域玉米单产估测的年际变化看，河北中部平原2010—2018年各县(区)的玉米估测单产波动变化，但总体呈先减少后增加的趋势；从玉米单产估测的空间分布看，整体呈西部地区最高，南部和北部地区次之，东部地区最低的特征。