基于DSSAT作物模型的中美大豆主产区单产模拟与验证

徐春萌，田芷源，陈 威，刘佳佳，白 洁

（1. 北京佳格天地科技有限公司，北京 100190；2. 中国科学院南京土壤研究所，南京 210008；3. 中国农业科学院农业信息研究所，北京 100081； 4. 农业农村部农业信息服务技术重点实验室，北京 100081）

0 引 言

鉴于中国与国外的粮食贸易规模不断加大，对世界主要农作物产量估测的需求日益突出，对估产方法高精度与时效性的要求越来越高[1]，开展客观、及时、大面积的农作物产量预测研究，可以减少对人工统计数据的依赖性，对政府决策具有重要意义。农作物单位产量（单产）的预测研究方法，包括基于统计回归模型和相似年份聚类法的纯气象研究手段，建立植被指数与作物单产关系的纯遥感估产方法，以及逐步发展形成的作物模型同化遥感数据的集合预测方法[2-4]，尽管研究已提出多个作物模型[5-7]。但由于作物模型是基于单个采样点进行开发的，在大面积区域应用时仍面临模型区域化处理的问题。卫星遥感技术可进行大范围地表监测，并获取其中的植被指数[8-10]。遥感-作物模型同化是将遥感数据与基于单点尺度的作物模型进行互补，在考虑大豆生育期生长规律的基础上，借助遥感植被指数反演叶面积，融合生育期的气象时空分布数据，实现大面积范围内的大豆单产估测。

农业技术转让决策支持系统（the Decision Support System for Agro-technology Transfer，DSSAT）是由美国乔治亚大学组织开发的[11]，并于20世纪末引入中国，在中国不同地区进行了适用性验证[12-14]。DSSAT模型可以对众多作物种类进行生长模拟，包括玉米、大豆、棉花等重要粮食经济作物[15]。该作物模型采用绿色叶绿素植被指数（Green Chlorophyll Vegetation Index，GCVI）作为模拟植被生长情况的指标。相比于常用的归一化植被指数（Normalized Difference Vegetation Index，NDVI），GCVI具有当叶面积很高时不会发生指标饱和的优势[16]。

在国内，科研工作者在中国以及国外地区（比如美国Corn Belt）开展了使用遥感信息及作物模型获取作物生物量及产量的探索研究[17-19]，积累了一定的技术方法。随着作物估产需求的提升，近年来越来越多的研究工作关注于同化遥感信息获取详细产量数据（如市县级、地块级别）的应用方法[20-21]。目前针对主粮作物（玉米和小麦）的试验较多，而对于其他重要粮食作物，比如大豆估产模型的研究相对缺乏。考虑到近年中国大豆进口需求的迅速增长，而美国是世界上重要的大豆生产国和中国大豆的进口来源[22]，本研究选取大豆作为作物估产研究对象，将中国吉林省与美国爱荷华州作为研究区域，在统一的估产技术路线基础上，对中美两国代表产区的

大豆单产进行研究，比较方法在不同地区的适用性。本研究通过使用大豆产量统计数据和DSSAT作物模型，并融合长时间序列的气象要素与MODIS中分辨率遥感影像，开展中美代表产区的大豆单产估测。通过在大面积范围内建立大豆单产估测模型，以期为全球重点区域的作物产量分析提供技术参考。

1 材料与方法

1.1 研究区概况

爱荷华州位于美国中北部黄金带，在北纬40°～44°，西经90°～97°之间，其耕地资源丰富，是美国大豆、玉米、小麦等经济及粮食作物的主产区之一。全州土地面积14.47万km2，大豆种植面积约4万km2。绝大部分地区地形平坦，从西北向南地势略微降低，水资源丰富[23]。爱荷华州属于温带大陆性气候，全年温差较大，降雨多集中在夏季且空间分布不均衡，大豆生育期间降雨总量为190～580 mm，日平均太阳辐射量为17～24 MJ/m2。

吉林省位于中国东北部黑土区，在北纬41°～47°，东经122°～131°之间，是中国重要的粮食生产地区之一。全省土地面积18.74万km2，大豆种植面积约2 000 km2。由东南向西北地势降低，东部为山地，中西部为平原地貌，河流众多。吉林属于温带大陆性季风气候，四季分明，雨热同期，大豆生育期间降雨总量为150～900 mm，日平均太阳辐射量为14～23 MJ/m2。

1.2 数据来源

气象参数的模拟资料采用NASA世界能源预测项目（Prediction of Worldwide Energy Resources，简称POWER）所反演的降雨量、日最高温度以及太阳辐射数据，原始气象观测资料来源于POWER数据网的区域数据资料日值集（https://power.larc.nasa.gov/data-access-viewer/），空间分辨率为0.5°×0.5°。本研究使用ArcGIS10.3软件以格点数据为基础进行计算，分别获得爱荷华州和吉林省2008－2017年间每年6－8月的累积降雨量（mm），每年6－8月的日平均太阳辐射量（MJ/m2），以及每年8月的日最高温度平均值（℃），并采用双线性空间插值法将空间尺度细化至500 m×500 m。

爱荷华州的大豆分类数据来自美国农业部（United States Department of Agriculture，USDA）每年发布的农田数据层（Cropland Data Layers，CDL），其空间分辨率为30 m×30 m；使用基于2017年Sentinel-2遥感影像，辅助结合高分1号数据，提取中国吉林省的大豆种植分布数据，其空间分辨率为10 m×10 m，分类结果通过对全省各市县3 802个实地采样点进行验证，精度大于87%，以上数据由北京佳格天地科技有限公司与吉林省农委合作完成。研究中假设吉林省2008－2017年大豆种植区域不变，对其他年份进行统一的大豆分布图掩膜处理。以上大豆分类数据均重采样至500 m×500 m分辨率。

植被指数使用遥感影像MODIS MOD09A1数据集计算GCVI值，其空间分辨率为500 m×500 m。GCVI的计算式如下：

式中NIR和G指近红外光波段（光谱范围841～876 nm）、绿光波段（光谱范围545～565 nm）反射率。

土壤参数来自于国际土壤资料与信息中心（https://www.isric.org/）土壤数据集。根据DSSAT模型中大豆模型SOYGRO需要的基础土壤数据进行提取，包括土层深度、各层质地（黏粒、粉粒、砂粒含量）、pH值、有机质含量、全氮、容重以及持水性能等指标。

田间管理信息选取研究区域典型大豆品种作为研究对象，通过收集美国农业部发布的爱荷华州大豆种植数据（https://quickstats.nass.usda.gov/）与吉林省6个气象站点（榆树、双阳、敦化、辽源、桦甸、延吉）记录的生育期数据，确定大豆各个生育期的日期。氮肥施用水平从50～300 kg/hm2不等。由于2个研究区域普遍不具备灌溉条件，因此在进行大范围估产的情况下将灌溉统一默认为无。

2008－2017年美国爱荷华州大豆单产数据来源于美国农业部调查的爱荷华州大豆产量历史数据（https://quickstats.nass.usda.gov/），以县为统计单位；中国吉林省大豆产量数据来源于吉林省统计局公布的统计年鉴（http://tjj.jl.gov.cn/tjsj/tjnj/），包含各市县的大豆播种面积和总产量，由此计算大豆单产数值。

1.3 单产模拟方法

1.3.1 DSSAT模型调参

本研究选取DSSAT模型中的大豆模型SOYGRO作为模拟工具[24]。按照DSSAT模型的标准与格式，使用历史年份数据（2008－2017年）分别建立气象文件、土壤与田间管理文件、生育期及产量文件，并通过该模型中的最大似然估计（Generalized Likelihood Uncertainty Estimation，GLUE）模块对大豆遗传参数进行率定，分别对爱荷华州和吉林省构建了区域尺度的大豆模型参数，结果如表1所示。

利用历史年份（2008－2017年）数据，通过将大豆生长发育过程中的关键决定参数进行任意组合，包括种植日期、气象数据、氮肥水平播种密度、品种参数等，进而构建多场景的大豆生长情况，进行逐日模拟并绘制叶面积指数（Leaf Area Index, LAI）曲线（模拟起始日期为每年第100天，结束日期为第365天）。通过分析不同情境下LAI随时间的变化曲线，设置了2个28 d大豆生长窗口期：生育早期（每年第201天至第228天）及生育晚期（第229天至第256天）。利用大豆LAI与GCVI的换算关系，如式（2）所示[25]，可将LAI随时间变化曲线转化为GCVI随时间变化的曲线。

结合历史年份（2008—2017年）的气象要素及遥感GCVI指数，获得研究区域内的大豆单产模拟向量系数表。向量系数表由14列（生育日期、气象要素、植被指数等相关参数），784行（28 d × 28 d）数据构成，可通过耦合大豆生长窗口期内的气象要素及遥感植被信息计算得到大豆单位面积产量的估测值。为评价向量系数表的准确性，将2008—2017年气象站点的历史调查数据作为单产真值与单产模拟值（通过向量系数表计算所得）进行拟合，结果显示两者之间的线性回归决定系数R2大于0.80。以上操作是将大豆模型SOYGRO链接到Python编译器，通过Python调用 SOYGRO的功能库完成。

表1 爱荷华州和吉林省的作物模型参数Table 1 Parameters of crop model in Iowa State and Jilin Province

1.3.2 大豆单产模拟

将美国爱荷华州和中国吉林省作为研究区域，分别获取2008—2017年基于MODIS影像反演的逐年GCVI植被指数，以及大豆生育期的气象要素，假设土壤性质和田间管理模式不变，利用研究区域所对应的大豆向量系数表，采用式（3）回归模型对大豆单产进行逐年模拟[26]。

式中Y为大豆单产，kg/hm2；W为生育期内的天气特征（6－8月累积降雨量（mm），6－8月日平均太阳辐射量（MJ/m2），以及8月日最高温度（℃）向量，RMd为生育期内的遥感反演植被指数GCVI向量，系数β为向量系数表中大豆生育早期（第201天至第228天）与生育晚期（第229天至第256天）某一日期组合下的参数。

利用上述回归模型，以谷歌地球引擎（Google Earth Engine）作为遥感数据计算云平台，输入历史年份（2008－2017年）的气象要素及遥感GCVI指数，分别进行2008－2017年美国爱荷华州和中国吉林省每年的大豆单产估测，输出单产模拟值的空间分辨率为500 m×500 m。

1.3.3 大豆单产精度验证

大豆单产精度验证所使用的数据集为2008－2017年爱荷华州各县（n=99）及吉林省各县（n=32～42）模拟大豆单产的平均值，验证集为中美两国政府发布的统计数据。对模拟大豆单产结果的误差分析，采用了3项计算指标进行评价：平均误差（Mean Percentage Error,MPE），均方根误差（Root Mean Square Error, RMSE），平均偏差（Mean Bias Error, BSE）。计算式如（4）～（6）所示。

式中n为统计县个数；mi为作物模型所得到的大豆单产，kg/hm2；di为调查数据代表的大豆单产，kg/hm2。

同时，关于大豆产量的分县数据，在各年份作物模型的产量值与调查数据之间分别建立了线性回归模型，使用相关系数（r）验证两者之间的相关性。

2 结果与分析

2.1 美国爱荷华州大豆单产估测

美国爱荷华州2008－2017典型年份的大豆单产估测值空间分布如图1所示。从图中可以看出，爱荷华州大豆产量的空间变异性在不同年份存在差异：单产较高的大豆（＞4 000 kg/hm2）在2008年分布于东部和中部地区，在2011年主要位于东北地区，在2016年则分布于东部及西南地区；单产较低的大豆（＜2 000 kg/hm2）主要分布于南部地区，但在2013年北部地区的大豆单产也较低。

美国爱荷华州2008－2017年大豆单产的年际变化如图2所示，模拟结果表明，大豆单产在2012年最低，仅为（2 139±193）kg/hm2，而在2010年大豆单产最高，为（4 766±970）kg/hm2，其余年份存在一定的产量波动。调查统计数据表明，爱荷华州大豆单产的最低值为2012年的（3 001±430）kg/hm2，而大豆单产最高值为2016年的（3 991±230）kg/hm2。从10 a间的单产波动来看，模拟值与调查值均呈现出大豆单产在2009年上升，于2011及2012年间下降，随后在2013至2016年上升，并在2017年再次下降的特征。美国爱荷华州大豆单产模拟值与调查数据的年际变化趋势较为一致，但作物模型估测的美国爱荷华州大豆单产的年际变化与统计调查数据相比，增减产的变化幅度较大，表现为单产预测值在高产年份比调查值更高，但在低产年份比调查值更低（图2）。

美国爱荷华州2008－2017年大豆单产估产的验证结果如表2所示。

表2 美国爱荷华州2008－2017年大豆单产模拟值与调查值误差比较Table 2 Errors of estimated soybean yields compared to investigated values from 2008 to 2017 in Iowa State, USA

综合来看，美国爱荷华州大豆模拟的单产区间略大于调查值，2008－2017大豆单产的模拟值为2 139～4 766 kg/hm2，而调查单产为3 002～3 991 kg/hm2（图2），大豆单产估产的平均误差为8.7%～40.7%，估测值的均方根误差为369.5～1 618.2 kg/hm2，平均偏差为-863.7～1373.9 kg/hm2（表2）。其中，2010年模型估测的误差较大，大豆单产的估测值大于单产调查值1 373.8 kg/hm2；而2011、2014年估测精确度较高，单产的平均偏差略高于100 kg/hm2。从10 a平均来看，2008－2017年美国爱荷华州大豆单产估产的验证结果为：平均误差为16.8%，均方根误差为762.8 kg/hm2，平均偏差为107.2 kg/hm2，且估测值比调查值偏大（表2）。

对比2008－2017年美国爱荷华州大豆单产模拟数据与调查数据，在县级尺度上对估产结果的平均值进行了精度验证，典型年份线性回归模型的拟合结果如图3所示。遥感信息结合DSSAT作物模型对美国爱荷华州2008－2017年大豆单产估测值与统计调查的产量值之间的相关系数（r），从0.57到0.78不等。其中，相关程度较好的为2013年和2008年，最差的年份为2016年。整体来看，各县的大豆单产模拟值与调查值相关程度较好。当爱荷华州大豆单产较高时，回归模型的相关系数较小，当单产较低时，回归模型的相关系数较大，这可能与模型所适用的大豆产量区间有关。

2.2 中国吉林省大豆单产估测

中国吉林省2008－2017典型年份大豆单产估测值的空间分布如图4所示。从图中可以看出，吉林省大豆产量的空间变异性在不同年份略有差异。单产较高的大豆（＞3 000 kg/hm2）在2012年分布在北部、中部和东部地区，在其它年份则主要分布于东部地区；而单产较低的大豆（＜2 000 kg/hm2）除2012年以外主要分布于西部地区，部分东部地区单产也较低（图4）。

中国吉林省2008－2017年大豆单产的年际变化如图5所示。其中，吉林省大豆单产的估测数据显示，2009年的大豆单产最低，仅为（1 653±544）kg/hm2，2017年大豆单产最高，为（2 766±933）kg/hm2。从10 a变化趋势来看，中国吉林省大豆单产模拟值在2009年显著下降，然后在2010年单产回升，直到2013年以后单产再次下降，并在2017年出现单产回升（图5）。吉林省的调查统计数据表明，大豆单产最低值为2009年的（1 997±667）kg/hm2，最高值为2013年的（2 797±930）kg/hm2。其中，大豆单产调查值在2009年下降，然后在2010年回升至稳定，直到2013年单产略微上升，并于2016年开始下降（图5）。综合来看，作物模型估测中国吉林省大豆单产的最低年份与调查值最低年份相吻合（2009年），且二者反映的高产年份相同（2008、2013、2015年）。除2017年模型认为产量上升以外，估测与调查所得大豆单产的年际变化波动趋势几乎一致。

中国吉林省2008－2017年运用遥感与作物模型计算所得的大豆单产估测值与调查值验证结果如表3所示。综合来看，2008－2017年中国吉林省大豆单产的模拟值为1 653～2 766 kg/hm2，其产量区间非常接近于调查产量1 997～2 797 kg/hm2（图5），各年间大豆单产估产的平均误差为24.5%～64.4%，估测单产的均方根误差为637.8～1 869.7 kg/hm2，单产平均偏差为-655.2～534.1 kg/hm2（表3）。其中，2015、2017年的大豆单产模型估测误差较大，而2008年的模型估测精确度较高，大豆单产的平均偏差仅为36.0 kg/hm2。总结2008－2017年的10 a平均值，大豆估产在中国吉林省的验证结果为：单产的平均误差为36.3%，均方根误差为1 088.4 kg/hm2，平均偏差为-237.9 kg/hm2，且估测值比调查值偏小（表3）。

表3 中国吉林省2008－2017年大豆产量模拟值与调查值误差比较Table 3 Errors of estimated soybean yields compared to investigated values from 2008 to 2017 in Jilin Province, China

以县为统计单元，对比了中国吉林省2008－2017年间大豆单产的模拟值与调查数据，并对大豆估产结果在各县的平均值进行了精度验证，其中典型年份线性回归模型的拟合结果如图6所示。遥感信息融合DSSAT作物模型对中国吉林省2008－2017年大豆单产估测值与统计调查值之间的相关系数（r）在0.5～0.6，其中r在2012年最高，且在2009、2016、2017年呈现出较好的相关性。整体来看，中国吉林省大豆单产模拟值与调查值的拟合程度较稳定，但其回归模型的相关系数相比于美国爱荷华州偏低。值得注意的是，吉林省模型的回归方程均位于1∶1线右下方，表现为大豆高产时模拟单产比调查单产小，低产时模拟单产比调查单产高（图6）。

3 讨 论

本文应用DSSAT作物模型同化遥感数据对中美两国大豆主产区进行了单产估测模拟研究。结果表明，大豆单产模拟值与调查统计值之间存在一定的相关性，且增减产的年际变化趋势相似；二者在单产的绝对数值上存在一定偏差。本研究的优势在于，通过将遥感信息融入DSSAT作物模型，首次将大豆估产模型由离散的试验点推广至连续的区域尺度，并建立起了中美两国模拟大豆单产的技术路线。

本研究使用作物模型估测单产的方法存在一定的误差，这是由于模拟不同年份作物产量时，不仅需要考虑每年气象因素以及植被长势发生变化，而且产量还受其他环境变量的影响。已有研究表明，在使用DSSAT模型模拟产量时如果仅考虑气象条件的影响，产量模拟值会存在误差[27-28]。研究分析，应用作物模型估算的作物产量为潜在产量，与实际产量之间存在一定的作物产量差[29]。由于地块间管理措施不同，以及天气和土壤状况存在差异，作物产量差具有时空变化特征。本文中大豆单产的估测值与调查值之间存在偏差，且差异在各年份、各地区之间并不相同，可能与地块间的水分管理措施、土壤及天气状况均存在相关性。因此在大面积范围内应用作物产量模型时，建议加强对复杂且相互影响的田间环境变量的考虑，这将是提高农作物单产估测准确度的有效途径。

基于遥感信息与DSSAT作物模型，同时进行中美两国重点地区长时间序列的大豆单产估测与验证工作，具备较强的可操作性与一定的准确性。在比较两国主产区大豆估产结果时，2008－2017年美国爱荷华州10 a平均的单产估测误差要小于吉林省，县域尺度上美国爱荷华州单产模拟值与真实值的相关系数（r）优于吉林省，这与两国的大豆产量区间、调查数据来源不同等因素有关。研究认为，使用该技术路线预测中美两国大豆单产并进行产量数据对比时，应注意大豆高产时单产模拟值偏小，大豆低产时单产模拟值偏大的问题。对其他作物的研究结果表明，参数本地化后的DSSAT模型在站点尺度上已达到模拟值与实测值较好的一致性[27]，而本文使用的大豆估产技术路线可以用作借鉴方法帮助其他农作物实现从站点到区域的产量模拟。

本研究对大豆的模拟误差，与预测其他主粮作物的误差结果相似，其中模拟玉米及小麦单产的平均误差为11%～16%，均方根误差为1 588～1 642 kg/hm2，平均偏差为245～388 kg/hm2[13,19-20]。前人利用纯遥感方法进行大豆估产的结果表明，NDVI不是唯一与大豆单产有关的变量，两者之间的相关性较低[30]。本研究的DSSAT作物模型中加入了气象因素对于产量影响的考虑，结合不同生育时期的遥感指数，进而提高了大豆估产精度与预测单产年际变化的稳定性。

应用遥感信息及作物模型的技术方法，在用于预测当季作物的潜在产量仍然具有挑战性。如何根据天气预报与植被生长指标的实时变化，在作物收获之前获取产量预测值，并以此调整田间管理措施，以及为可能出现的气象灾害制定应对策略[31-32]，是需要进一步研究解决的问题。基于本文使用的区域估产技术路线，建议首先根据多年的历史产量数据获取稳定的本地作物模型参数，并不断积累天气及其他环境变量如何影响作物单产变化的数据分析案例，进而实现为种植者与决策者做出及时科学预报的目标。其次在世界范围内，高分辨率对地观测系统的建立正在逐步加强，尤其是商业卫星的加入，将使得运用高精度遥感影像成为开展作物产量估测工作的发展趋势。除本文使用的MODIS影像以外，随着Sentinel-2，Planet等高分辨率遥感影像的日益积累[33-34]，将提升对全球农作物的监测能力，为获得更高精度的作物分布以及产量数据提供细节丰富的影像信息。这将有助于本文所构建的DSSAT模型结合遥感信息以及历史产量数据开展大豆单产估测的技术路线，得到进一步地推广应用。

4 结 论

本文在获取植被指数、气象要素、大豆种植空间分布以及种植情况调查数据的基础上，通过DSSAT模型，在中国吉林省与美国爱荷华州分别对2008－2017年大豆单产进行估测。研究以500 m×500 m栅格单元为估产目标，在一定程度上细化了作物单产的空间变异性，提高了估产结果的空间分辨率。结果表明，基于长时序的气象要素与MODIS遥感影像，结合种植环境因子所构建的估产技术路线，能够反映大豆单产的空间差异性，以及单产在年际间的波动趋势。估产结果实现了对中美两国代表性农产区在大面积范围内的大豆单产模拟，可以为开展全球农作物的单产预估提供参考。

1）2008－2017年间与大豆单产的调查值相比，应用遥感信息融合作物模型计算大豆单产估测值在爱荷华州大豆单产估测时的平均误差为16.8%，均方根误差为762.8 kg/hm2，平均偏差为107.2 kg/hm2；吉林省的平均误差为36.3%，均方根误差为1 088.4 kg/hm2，平均偏差为-237.9 kg/hm2。通过县域尺度的精度验证可知，本文利用作物模型同化遥感数据对大豆单位面积产量进行估测的技术路线具有一定的准确性，且在爱荷华州的模型估测精度要优于吉林省。

2）大豆单产在2008－2017年间出现了一定的年际波动：在中国吉林省，模型估测大豆单产出现低产及高产的年份与调查数据相吻合，其增减产的波动趋势与调查数据的波动相似；在美国爱荷华州，作物模型估测显示的大豆增减产的趋势与统计调查数据基本一致。