0 引言
【研究意义】冬小麦是中国主要粮食作物之一,及时、准确地估算冬小麦产量,不仅能加强冬小麦生产管理,进一步发挥其生产潜力,同时能为国家调整粮食储备,政府有关部门制定宏观经济计划和粮食政策提供有利的数据支撑[1]。【前人研究进展】目前农作物估产方法总体上可以分为统计预测、农学预测、气象预测和遥感预测四大类[2-6] 。统计预测方法通过建立影响作物单产的统计因子与农作物单产之间的关系模型来预测粮食产量[6]。统计预测方法整体考虑产量的各种影响因子,不考虑作物产量形成的复杂机理,模型一般简单;但是也因其没有明确反映作物生长发育的过程,导致区域外推性不高[7-11]。农学预测方法是通过模拟作物生长的整个过程,建立作物生长状况与作物产量构成要素之间的关系,进而实现产量预测;不同作物产量构成要素不同,如冬小麦的产量构成要素包括有效穗数、每穗粒数和千粒重。该模型机理比较复杂[12],对输入数据要求高。气象预测方法是建立气象数据与作物产量之间的统计关系,并在考虑其他因素的基础上建立统计模型,最终利用统计模型实现作物单产估算的一类模型的总称[13-14]。该模型基本满足农作物产量估算精度需求,但也存在着空间外推和空间适用范围界定难的问题[6]。遥感预测模型是根据生物学原理,在分析收集农作物光谱特征的基础上,通过光谱来获取作物的生长信息,并在作物光谱与产量之间建立关联[6]。随着遥感技术的发展,尤其高分辨率影像的出现,给农作物的估产以及农作物的生长状态监测提供了技术支持,遥感预测方法已成为大面积估产的主流研究方法[15-16]。该类方法一般通过遥感植被指数对农作物产量进行预测[15-24]。综合考虑不同估产方法的优缺点,采用单一的估产方法弊端较大,在精度、可行性等方面存在较大的不确定性[2-16]。近年来在农作物估产领域,利用多种方法构建模型提高农作物估产精度成为热点[25-36]。王培娟等[37]将农学、遥感预测方法相结合用于华北平原冬小麦估产研究,余福水等[38]将气象、农学和遥感预测方法相结合对黄淮海平原冬小麦产量估计,MIRSCHEL等[39]将统计、农学、气象预测方法和专家知识数据库结合,对德国东部地区的多种作物产量估算,HUANG等[40]将农学、遥感预测方法结合,利用改进后的WOFOST模型对美国西部地区的小麦估产探究,均取得较好的估产精度,但模型机理复杂、且数据不易获取。【本研究切入点】该领域仍存在一些不足。比如在大尺度估产方法研究中,针对市、县级尺度的研究较少,而高分辨率像元级尺度的大面积估产的相关研究更少。已有的各种估产方法,优缺点明显,高精度的估产方法普遍比较复杂,应用不方便。【拟解决的关键问题】本文在分析国内外农作物估产方法的基础上,将传统估产方法和遥感估产方法相结合,提出了一种由趋势单产、遥感修正单产和随机误差项三部分组成的混合估产模型,以期对冬小麦的高精度估产相关研究提供借鉴,对大面积冬小麦估产提供指导。1 材料与方法
1.1 研究区概况
北京市地处115.42—117.5°E,39.43—41.05°N,位于华北平原的西北边缘。历年物候数据和统计数据显示,北京市冬小麦一般在10 月初种植,来年6月中旬收割。根据常年统计数据分析,北京市种植冬小麦的区县有12个,其中房山区,大兴区,通州区和顺义区是北京市四大冬小麦主产区[18]。1.2 研究数据
本研究所用到的数据如表1所示,包括2015年冬小麦关键期的三景环境卫星的遥感影像,用来获取多时相NDVI;2015年冬小麦种植面积栅格图,用来获取研究区冬小麦种植范围; 2015年的42个冬小麦实测地块单产数据和近30年北京市各县区冬小麦单产数据（该数据来源于统计年鉴）其中1985—2014年数据用于构建模型,2015年数据用于验证模型精度。Table 1
表1
表1估产应用数据
Table 1The data of yield estimation applications
| 数据名称 Data name | 数据时间 Data time | 分辨率 Resolution | 用途 Purpose |
|---|---|---|---|
| 环境卫星影像 HJ imagery | 2015-04-22 2015-05-07 2015-05-19 | 30 m | 获得多时相NDVI Get multi-phase NDVI |
| 冬小麦种植面积栅格图 The raster map of winter wheat area | 2015 | 10 m | 获得冬小麦种植范围 Get the range of winter wheat area |
| 实测地块单产数据 The yield data of sampling field | 2015 | — | 计算抽样村单产及参与建模 Calculate the yield of sampling yield and build model |
| 市、区县级单产数据 Yield data at city and country levels | 1985-2014 | — | 构建模型 Build model |
| 市、区县级单产数据 Yield data at city and country levels | 2015 | — | 验证模型精度 Verify model accuracy |
| 北京市及区县行政边界 Administrative boundary of Beijing (include city and country) | — | — | 研究区位置 Location of study area |
新窗口打开
冬小麦种植面积栅格图利用高分一号遥感影像2 m全色波段和16 m多光谱波段影像融合,结合实地考察和目视解译进行监督分类,利用最近邻法重采样成10 m分辨率。
实测地块单产数据是北京市统计局提供,通过选取样本地块、冬小麦收割、脱皮、晒干、称重几个主要步骤,最后得出样本地块的实测单产数据。
1.3 混合估产模型
混合估产模型由趋势单产模型、遥感修正单产模型、随机误差项3部分构成,见公式（1）。其具体构建流程如图1所示,主要包括3个步骤：（1）基于1985—2014年的冬小麦亩产数据,利用多项式回归的方法,构建冬小麦的传统趋势单产模型;（2）采用ARIMA模型,对所建的趋势单产模型的时间序列误差进行修正;（3）基于实测地块单产数据和对应地块的NDVI数据（由环境卫星影像获得）,通过多元回归构建遥感单产修正模型。
显示原图|下载原图ZIP|生成PPT图1估产模型流程图
-->Fig. 1Flow chart of yield estimation model
-->
Y= Yt +Ys+U0 （1）
式中,Y为本文混合模型预测单产,Yt为趋势单产模型,Ys为遥感修正单产模型,U0为单产误差项。Yt表达式如下所示：
Yt = Yc + YARIMA + U1 （2）
式中,Yc表示传统趋势单产模型,YARIMA为利用ARIMA模型来修正的趋势单产的系数项,U1为趋势产量与真实产量之间的差值。
1.3.1 传统趋势单产模型 传统趋势单产模型（Yc）的表达式如下
Ycx =aX3+bX2+cX+d+U2x （3）
式中,Ycx表示第X年的趋势单产,单位：t·hm-2,X表示年份,利用近30年（1985—2014年）的统计数据构建模型;a、b、c、d为参数;U2x 表示第X年传统趋势模型产生的随机误差,X是1985到2014年的单产时间序列,U2x 构成了传统趋势模型的时序随机误差项U2 。
1.3.2 传统趋势单产时序误差修正模型 ARIMA（p,d,q）模型也称为差分自回归移动平均模型,p 为自回归项,q 为移动平均项数,d 为时间序列成为平稳时所做的差分次数。
本文针对传统趋势单产时序误差（U2）利用ARIMA模型进行模拟分析,构成最终的趋势产量模型,具体步骤如下：
（1）平稳化处理。采用差分的方法,对时间序列进行平稳化处理,以得到一个平稳随机序列,进而确定差分次数d;（2）采用AIC准则进行定阶[41]。确定所建立的ARIMA模型的阶数p 和q。这一过程以自相关系数和偏自相关系数分析为基础。（3）用最小二乘法对ARIMA模型分析,计算模型参数值;（4）模型的检验。对于得到的拟合模型需要进行随机性检验,即检验模型的残差序列是否为白噪声序列。若残差序列不是白噪声序列,模型需重新构建。在得到最佳的模型后,根据时间序列的历史数据,运用ARIMA模型来进行预测。
本研究利用ARIMA模型对趋势单产误差时间序列分析,预测2015年的传统趋势单产的误差。
1.3.3 趋势单产遥感修正模型 趋势单产误差（U1）利用遥感数据对趋势单产模型进行修正,遥感修正单产模型（Ys）表达式如下：
Ysm = Dmn×Sn+Em+U3 （4）
式中,Ysm 表示实测地块单产数据（真值）与趋势产量之间的差值,单位：t·hm-2,Sn为遥感因子,本文选用作物关键期的NDVI（实测地块内对应像元NDVI之和）,Dmn、Em为参数矩阵,U3 则为遥感修正模型中的随机误差项,n表示实测选取的样本地块的数目,m表示选取冬小麦关键期的影像的数目。研究表明[18,42],作物单产模拟的关键时期是产量形成的关键期,而作物的水分临界期是农作物产量形成的关键时期,大多数农作物的水分临界期都位于花芽分化的旺盛生长期,其中小麦位于孕穗到抽穗期。本研究分别利用冬小麦关键期3到5月一景遥感影像（2015年4月22日）、两景遥感影像（2015年4月22日和2015年5月7日）和三景遥感影像（2015年4月22日、2015年5月7日和2015年5月19日）的NDVI构建3个回归模型,分析3个模型对趋势单产模型的修正精度。由于实测地块只有42个,样本较少,本研究选择对北京市整体趋势单产做遥感修正。
1.3.4 检验方法 本文采用统计学中的平均偏差AD、相对误差ME、决定系数R2 3个指标评价冬小麦估产结果。
式中,n表示选择样本个数,Yi表示第i个的地区的实测单产,Fyi表示第i个地区的预测单产。
式中,Fy为预测单产,Yz为实测单产。
式中,Fyi、Yi分别表示冬小麦预测单产、实测单产。
2 结果
2.1 趋势产量模型构建
2.1.1 传统趋势单产模型 利用1985—2014年北京市以及各区县冬小麦单产数据,对北京市以及各区县的冬小麦亩产做了统计分析,分别对北京市及各个区县冬小麦亩产趋势利用多项式回归模型模拟,结果如下表2所示,模型含义与式（3）相同,Yc为传统趋势模型预测单产,单位：t·hm-2,X代表年份,当X为2015时,可以预测出2015年的单产,与2015年统计数据对比。Table 2
表2
表2北京市冬小麦单产传统趋势模型
Table 2The traditional trend model of winter wheat yield in Beijing
| 地区 Region | 传统趋势模型 Traditional trend model (Yc) | R2 |
|---|---|---|
| 北京市 Beijing | Yc = 0.0487X 3 - 292.2X 2 + 584934X - 4E+08 | 0.8838** |
| 海淀区 Haidian | Yc = 0.0971X 3 - 582.95X 2 + 1E+06X - 8E+08 | 0.688** |
| 房山区 Fangshan | Yc = 0.0566X 3 - 339.95X 2 + 680278X - 5E+08 | 0.7492** |
| 通州区 Tongzhou | Yc= 0.0675X 3 - 405.39X 2 + 811361X - 5E+08 | 0.8762** |
| 怀柔区 Huairou | Yc = 0.0553X 3 - 332.32X 2 + 665846X - 4E+08 | 0.891** |
| 顺义区 Shunyi | Yc = 0.0584X 3 - 350.57X 2 + 701417X - 5E+08 | 0.7571** |
| 大兴区 Daxing | Yc = 0.048X 3 - 288.5X 2 + 577696X - 4E+08 | 0.8779** |
| 昌平区 Changping | Yc = 0.0678X 3 - 407.07X 2 + 815061X - 5E+08 | 0.8469** |
| 密云县 Miyun | Yc = 0.0328X 3 - 196.84X 2 + 394352X - 3E+08 | 0.6443** |
| 平谷区 Pinggu | Yc = 0.0484X 3 - 290.55X 2 + 581791X - 4E+08 | 0.7558** |
| 丰台区 Fengtai | Yc = -0.7881X 2 + 3149X - 3E+06 | 0.7144** |
| 朝阳区 Chaoyang | Yc= -0.8462X 2 + 3382.8X - 3E+06 | 0.8343** |
新窗口打开
2.1.2 ARIMA修正模型 通过针对传统趋势单产模型冬小麦亩产预测值与真实值的差值,构建时间序列,利用统计软件R实现北京市以及各个区县的ARIMA模型,满足AIC准则下,模型中的主要参数如表3所示。
Table 3
表3
表3ARIMA修正模型的参数
Table 3The parameters of ARIMA correction model
| 地区 Region | p | d | q |
|---|---|---|---|
| 北京全市 Beijing | 2 | 0 | 0 |
| 房山区 Fangshan | 4 | 2 | 1 |
| 通州区 Tongzhou | 0 | 0 | 0 |
| 顺义区 Shunyi | 2 | 2 | 1 |
| 昌平区 Changping | 6 | 0 | 0 |
| 大兴区 Daxing | 0 | 0 | 0 |
| 平谷区 Pinggu | 1 | 2 | 1 |
| 怀柔区 Huairou | 2 | 2 | 1 |
| 密云县 Miyun | 8 | 1 | 5 |
新窗口打开
2.2 遥感修正单产模型
趋势单产模型预测冬小麦亩产值与野外实测冬小麦单产值之间的差值作为因变量,以冬小麦关键期遥感影像的NDVI值作为自变量,建立的回归模型。利用一景遥感影像构建的模型中,即m=1,n=42,系数如下所示：
D1,42=[102.3559]
E1=-71.3108
F检验结果为2.064,经查F检验表F0.05(1,40)= 4.085,F<F0.05(1,40),说明4月22日的NDVI与冬小麦产量在0.95水平下存在显著一元线性统计关系。
利用两景遥感影像构建的模型中,即m=2,n=42,系数如下所示：
D2,42=[-161.281 293.2849]
E2=-89.4066
F检验结果为2.759,经查F检验表F0.05(2,39)= 3.238,F<F0.05(2,39),说明4月22日和5月7日的NDVI与冬小麦产量在0.95水平下存在显著二元线性统计关系。
利用三景遥感影像构建的模型中,即m=3,n=42,系数如下所示：
D3,42=[-160.8122 333.0112 -40.0059]
E3=-88.7954
F检验结果为1.812,经查F检验表F0.05(3,38)= 2.852,F<F0.05(3,38),说明冬小麦关键期三景遥感影像的NDVI与产量在0.95水平下存在显著多元线性统计关系。
2.3 应用模型的估产结果
2.3.1 传统趋势单产模型 北京市各区县传统趋势模型由其相关系数R2（表2）可以看出,构建的多项式回归模型与真实的冬小麦单产数据相关性较高,冬小麦的单产变化符合一定的统计规律,且通过了显著性检验。2.3.2 遥感修正模型 利用一景、两景和三景影像的NDVI对趋势单产模型修正,将修正之后的混合估产模型估产结果与2015年北京市统计年鉴中冬小麦单产数据进行对比。结果表明,利用三景影像修正结果（精度平均提升3.55%）优于一景影像修正（精度平均提升2.70%）和两景影像修正（精度平均提升1.47%）。在市级尺度和县级尺度上,3种遥感修正模型修正后的混合估产模型估产精度均超过90%。由于利用三景影像修正结果最稳定,对于估产精度提升幅度最大,故本研究选择三景遥感影像修正模型来构建最终的混合估产模型。
2.3.3 三种模型估产结果对比 以统计数据作为真实值,将传统趋势产量模型、ARIMA模型修正后的趋势产量模型和混合估产模型进行对比,结果见表4。真实值为2015年北京市统计年鉴中冬小麦单产数据。
Table 4
表4
表4冬小麦估产模型结果
Table 4The estimated values of winter wheat yield by different estimation models
| 地区Region | 1* | 2* | 3* | 4* | 5* | 6* | 7* | 8* | 9* | 10* |
|---|---|---|---|---|---|---|---|---|---|---|
| 北京市 Beijing | 5.3568 | 5.6378 | 0.9475 | -0.1027 | 5.5351 | 0.9667 | -0.2432 | 5.2919 | 0.0649 | 0.9879 |
| 海淀区 Haidian | 5.0973 | 5.2162 | 0.9767 | — | 5.2162 | 0.9767 | — | 5.2162 | 0.1189 | 0.9767 |
| 房山区 Fangshan | 4.8323 | 5.8813 | 0.7829 | -0.2681 | 5.6131 | 0.8384 | -0.3137 | 5.2994 | 0.4672 | 0.9033 |
| 通州区 Tongzhou | 6.1304 | 6.3359 | 0.9665 | 0.0060 | 6.3419 | 0.9655 | -0.1961 | 6.1458 | 0.0155 | 0.9974 |
| 密云县 Miyun | 4.2048 | 4.5472 | 0.9186 | 0.0669 | 4.6141 | 0.9027 | -0.3544 | 4.2597 | 0.0549 | 0.9869 |
| 平谷区 Pinggu | 4.8644 | 5.6155 | 0.8456 | -0.1833 | 5.4322 | 0.8833 | -0.4487 | 4.9834 | 0.1191 | 0.9755 |
| 大兴区 Daxing | 5.1090 | 5.5772 | 0.9084 | -0.0323 | 5.5448 | 0.9147 | -0.4121 | 5.1327 | 0.0237 | 0.9954 |
| 昌平区 Changping | 4.0850 | 4.2028 | 0.9711 | -0.1061 | 4.0967 | 0.9971 | 0.0496 | 4.1463 | 0.0614 | 0.9850 |
| 顺义区 Shunyi | 5.5596 | 6.0121 | 0.9186 | -0.1500 | 5.8621 | 0.9456 | -0.0535 | 5.8086 | 0.2490 | 0.9552 |
| 怀柔区 Huairou | 4.2200 | 4.4198 | 0.9527 | -0.3828 | 4.0369 | 0.9566 | 0.3758 | 4.4127 | 0.1927 | 0.9543 |
新窗口打开
在市级尺度上,传统趋势模型、ARIMA修正后的趋势模型、混合估产模型精度均达到93%以上。而在县级尺度上,传统趋势模型在房山区（精度78.29%）、平谷区（精度84.56%）精度偏低,引入ARIMA模型对传统趋势单产模型进行修正,修正后的6个县估产精度提高,两个县估产精度降低,其中怀柔县精度提升最大,为9.07个百分点,但仍存在房山区（精度83.84%）、平谷区（精度88.33%）精度未达到90%,整体上平均精度提高1.59个百分点,也说明了本研究引入ARIMA模型修正的必要性。传统趋势单产模型和ARIMA修正后的趋势单产模型仅考虑统计学中的意义,能反映出产量变化的整体趋势,却很难体现出地域性的差异,并且由突发自然灾害或其他原因引发的产量变化很难用传统统计模型预测,故引入遥感修正单产模型,汲取遥感影像具有时效性和空间分辨率高的优点,构建混合估产模型,与真实值（2015年统计数据）比较,最终精度均达到90%以上,混合估产模型与遥感模型修正前的趋势模型对比,平均估产精度提高3.55%。
总之,混合估产模型在市级尺度上北京市冬小麦估产精度达到98.7%,除房山的估产精度为90.3%之外,其他区县的精度均达到95%以上,满足县级尺度估产的要求,整体精度比较高。其中海淀区由于统计年鉴数据不全,没有完整的单产数据,无法进行ARIMA模型和遥感修正单产模型修正处理。
混合估产模型与传统趋势单产模型和趋势单产模型相对误差进行对比,如图2所示,整体误差大小排序为传统趋势单产模型相对误差?趋势单产模型相对误差?混合估产模型相对误差,反映出混合估产模型精度最高;并且在利用ARIMA模型修正后的趋势产量模型,相比传统趋势单产模型在整体估产精度上有一定程度提高;而利用遥感修正模型修正后的混合估产模型,整体估产精度进一步提升,房山、平谷地区的估产精度明显提高。
3 讨论
统计预测模型能够稳定反映冬小麦单产变化,但推广性差,且由于统计数据多以行政单元为基础,基于统计预测模型估算出的作物产量往往尺度较低,只能反映行政单元尺度（如市、县）作物产量的整体水平;遥感预测模型能够及时快速地获得大面积作物信息,并可以预测出像元尺度的产量,但其预测精度不稳定。
显示原图|下载原图ZIP|生成PPT图2相对误差对比
-->Fig. 2A contrastive study of relative error
-->
本文提出的混合估产模型综合考虑了统计预测模型和遥感模型特点,利用长时间序列冬小麦单产统计数据,充分考虑冬小麦在时间尺度上的变化趋势,估算趋势单产,得到预测年份较低尺度作物产量的估计值,进而利用估测年份实测地块样方作物产量与遥感影像反映作物长势的植被指数NDVI进行回归建模,植被指数NDVI能反映不同空间地域冬小麦长势的差异,即冬小麦在空间尺度上的变化,进一步修正趋势单产,得到最终的产量估计值。该模型是时间尺度和空间尺度相结合的估产模型。
与已有的主流估产方法相比,该方法扬长避短,相比于统计预测模型,该方法能体现出地域性的差异,可以反映出突发自燃灾害或其他原因引发的产量变化。相比于遥感预测模型,我们的预测方法由于融合了统计学的方法,能更好的反映出产量变化的整体趋势,使得估产精度更稳定。相比农学估产模型来说,简单易操作且估产精度较高。此外,本文方法可以得到不同尺度（北京市市级尺度、各县区级尺度、地块级尺度和像元级尺度）产量估算结果,例如图3为利用本文提出的混合估产模型对顺义地区2015年冬小麦产量的估产结果。这里要特别说明的是由于没有像元级的实测单产,本文没有进行像元级别的单产估测结果精度评价。
此外,本文实测地块只有42个,实测地块数量太少在一定程度上影响着模型的精度,可通过增加每个区县实测地块数量,针对每个区县做遥感修正模型,来提高模型精度。
显示原图|下载原图ZIP|生成PPT图3混合估产模型估产结果
-->Fig. 3The results by using the hybrid yield estimation model
-->
4 结论
4.1 本文所提出的混合估产模型精度要高于传统趋势单产模型,混合估产模型在市级尺度和区县级尺度均取得较高精度。具体来说,混合估产模型对北京市的冬小麦单产预测精度达到98.7%,各区县估产精度均超过90%,除房山（90.3%）外,各县单产预测相对精度均超过95%;而传统趋势单产模型对北京市的冬小麦单产预测精度达到94.75%,但在区县尺度上,传统估产模型预测精度较低,对房山区的估产精度不足80%。
4.2 在趋势修正单产模型的建立过程中,引入ARIMA模型可以提高传统趋势单产模型的精度。修正后的趋势单产模型冬小麦单产预测精度平均提高了1.59%。
4.3 本文建立的遥感修正模型,利用三景遥感影像修正结果最优,此方法使冬小麦估产精度整体提升3.55%,尤其是房山、平谷等区县,精度明显提升。
本文所提出的混合估产模型在市级尺度和县级尺度上预测冬小麦单产均取得较高精度,充分考虑冬小麦时间尺度和空间尺度上的变化,对农作物估产有一定的指导意义。The authors have declared that no competing interests exist.
