王利民，姚保民，刘 佳，杨玲波，杨福刚

基于SWAP模型同化遥感数据的黑龙江南部春玉米产量监测

王利民1，姚保民1，刘 佳1※，杨玲波2，杨福刚1

（1. 中国农业科学院农业资源与农业区划研究所，北京 100081；2. 浙江大学环境与资源学院，杭州 310012）

农作物种植类型的地理分布差异，气候条件差异、土壤环境不同等因素的影响，需要开展农作物生长模型参数区域化、本地化的研究工作；通过改善区域气象数据空间化方法以提升插值精度的研究，也需要得到应用的重视。针对以上问题，该文以SWAP（soil-water-atmosphere-plant model，土壤-水-大气-作物模型）模型为基础，以中国黑龙江省南部地区作为研究区域，以其主要农作物春玉米为目标作物，确定研究春玉米的作物生长模型参数，并综合考虑纬度及海拔对气温的影响情况，研究将协同克里金（coKriging）方法引入作物生长模型气象数据插值获取中，从而提高模型输入参数中气象数据精度，并以叶面积指数（leaf area index，LAI）及蒸散发（evapotranspire，ET）数据作为同化遥感数据源，通过优化玉米灌溉量和出苗日期，获取了研究区2013年的玉米产量空间分布成果，与统计资料结果对比，玉米总产量监测结果的2达到了0.939 4，均方根误差（root mean squared error，RMSE）达到了148 065 t，平均绝对误差（mean absolute error，MAE）为114 335 t。研究区15个县市区的预测单产和统计单产之间的决定系数达到了0.724 5，RMSE为598.5 kg/hm2，MAE为531.5 kg/hm2。研究结果表明，利用SWAP模型，以协同克里金方法获取气象数据空间插值成果作为输入数据，通过同化LAI和ET遥感数据，可以有效进行黑龙江南部区域的玉米产量遥感监测，为区域作物生长及生产力的遥感监测预测提供参考。

作物；遥感；模型；玉米；叶面积指数；蒸散发；协同克里金；同化；SWAP模型

0 引 言

玉米是中国主要的大宗农作物之一，其种植面积常年位居中国农作物种植面积的首位，2016年中国的玉米种植面积达到3 676.8万hm2[1]。近年来，随着中国种植业结构调整政策的不断推进，尤其是在中国“镰刀弯”地区调减玉米种植面积，增加豆类等作物面积成为优化中国种植业结构的重要举措，监测大尺度区域玉米作物的种植及生长情况，对中国农业政策的制定、调整及评估具有重要的意义[2]。然而，当前区域尺度的作物面积及产量信息获取的传统手段——统计调查，已逐渐难以适应新形势下对于农业信息及时、准确、客观获取的要求，利用遥感手段进行农情信息的快速、精确、低成本获取，逐渐成为大尺度区域作物信息获取的重要来源及研究方向[3-5]。采用同化技术结合遥感数据与作物生长模型是实现大尺度、高频率、高精度的农作物产量监测预测的有效途径。

作物生长模型使用严格的数学模型方法，考虑对作物生长具有关键影响的光、温、水、土、肥及田间管理等因素，描述作物从播种到收获的整个发育期的生长过程，并模拟作物不同发育期内的各项生理、生化参数，从而为科学、合理、定量模拟作物生长提供可靠的手段[6]。当前最常用的作物生长模型包括WOFOST（world food studies）、DSSAT（decision support system for agrotechnology transfer）、SWAP（soil-water-atmosphere-plant model）等[7-9]。作物生长模型具有机理性强、精度较高等优势，然而，其主要是针对地块尺度上的作物生长状态模拟，而在大尺度区域上的应用则受到了土壤条件、作物类型、田间管理等参数不同及难以获取等因素的影响而受到制约。而遥感数据具有覆盖范围大、数据时效性好、成本低等一系列优势，通过数据同化技术，将遥感数据与作物生长模型进行耦合，从而实现作物生长模型在大尺度区域上的应用，逐渐成为解决区域尺度作物生长监测的主要技术方向[10-12]，相关的同化算法主要包括以集合卡尔曼滤波（ensemble Kalman filter，EnKF）为代表的更新同化算法和以复合型混合演化算法（shuffled complex evolution method developed at the University of Arizona，SCE-UA）为代表的参数优化同化算法[6]。Cheng等研究了基于WOFOST模型和时序HJ-1卫星数据同化的地块尺度玉米产量估算方法，获得了较高的精度，反演产量的相对误差为-0.99%，均方根误差为530.81kg/hm2[13]；黄健熙等研究使用了EnKF方法，同化MODIS LAI产品和WOFOST模型，实现了衡水市冬小麦产量的估测，均方根误差达880kg/hm2[12]；姜浩则研究了基于SWAP模型及SCE-UA算法同化MODIS LAI数据的2017年华北地区冬小麦产量预测方法，2达到0.511，均方根误差423 kg/hm2[14]。在作物生长模型与遥感数据同化进行作物监测领域内，当前研究及应用较多的为WOFOST模型[15-21]，SWAP模型等其他模型的研究则相对较少[10,14,22]，尤其是针对玉米作物。

当前作物生长模型同化遥感数据的相关研究主要集中在同化算法、同化参数等的研究中[23-27]，而对作物生长模型影响较大的气象因子研究相对较少。大部分作物生长模型所需要输入的气象数据是由地面固定气象站点获取，并使用插值方法获取[28]，仅少部分是通过其他手段如TRMM（tropic rainfall measurement mission）数据获取[9]。而当前作物生长模型应用中，常用的气象数据插值方法主要是反距离权重、普通克里金插值等方法[12,20-21,28-29]。反距离加权插值法根据距离的远近对数据进行插值，具有计算简便、对数据分布无要求等特点；克里金插值法以属性变化是随空间位置变化存在自相关关系为依据，并假设样本点的变异函数是对待估样本取值的无偏、最优估计为前提进行插值，结果平滑性较好，能够明显降低歧异值对样本整体分布的影响[30]。然而，这两种方法都未能充分考虑到气温等数据除了空间自相关性外，与其他因素，如海拔、纬度等存在密切的关系，从而导致其插值结果可能与实际地面温度存在偏离，导致作物生长模型输入数据误差进而影响精度。综合考虑纬度、植被、海拔等多方面因素对气象要素的影响，进而实现气象数据的高精度获取及输入，将对作物生长模型的优化及作物生长监测具有重要的意义。

针对以上需求，本文以中国玉米的主产区之一，东北平原黑龙江省西南部的哈尔滨附近作为主要研究区，进行基于SWAP模型的玉米产量遥感同化监测研究，通过使用协同克里金（coKriging）方法，综合考虑纬度、海拔等气象要素对气温等气象数据的影响，实现SWAP模型的本地化作物参数标定，并使用MODIS ET及LAI数据作为输入遥感数据，通过同化作物出苗日期和灌溉量，实现了研究区玉米产量的监测，并与统计年鉴资料进行对比，进行精度验证，评价SWAP模型同化遥感数据在东北区域玉米产量监测上的可行性和准确性，为大尺度区域作物生长与生产力遥感监测预测提供支撑。

1 研究区概况

研究区位于黑龙江省南部，地处中国东北平原北部、松嫩平原南部，横跨哈尔滨、绥化和伊春3个地级市，总面积约52 457 km2，约占黑龙江省总面积的11%，纬度范围45.1°～47.6°N，经度范围125.3°～129.4°E。研究区大部分为平原地貌，地势低洼，仅在东部存在部分山区，主要河流为松花江和嫩江。研究区四季分明，冬季漫长寒冷，夏季短暂清凉，属于温带大陆性气候区，冬季1月平均气温约-19 ℃，夏季7月平均气温约23 ℃。年平均降水量在500～600 mm左右，降水主要集中在夏季。研究区的土壤类型主要为黑土和草甸土，土壤养分含量丰富，适合于各种农作物生长。研究区是中国商品粮的主要产地之一，主要作物类型为玉米和水稻，其他如大豆、瓜果等也有少量种植。研究区位于中国“镰刀弯”地区中的东北冷凉区，属于需要调减玉米种植面积的主要地区之一，监测研究区玉米种植及生长情况具有重要的意义。2013年研究区的玉米种植面积约为185万 hm2，产量为1 448万 t，约占黑龙江省玉米种植面积和总产量的30%。研究区春玉米的播种时间一般为5月上旬，5月20日左右进入出苗期，并在9月下旬进入成熟期，总生长期约125 d。范围及所包括的县城如图1a所示。

2 研究数据及预处理

2.1 气象数据及其预处理

气象数据从中国气象科学数据共享服务网获得，包括2013年全年东北的逐日最高温、最低温、日照时数、水汽压、风速和降水量等气象要素数据以及各气象台站的地理坐标位置。同时，使用Angstrom公式将原始气象数据中的日照时数数据转化为辐射量数据[31]，以满足SWAP模型对于输入数据的要求。具体的计算公式如下

式中R为太阳辐射量，MJ/(m2·d)；a和b为经验系数，分别设为0.29和0.42；为日照时数，h，即气象站点记录值；为当地当日的最大日照时长，h，计算公式如公式（2）所示；R为大气上界入射辐射量，MJ/(m2·d)，计算公式如公式（3）所示。

式中d为日地距离系数；G为太阳常数，其值为0.082 0 MJ/(m2·min)；为太阳赤纬，（°）；ω为太阳时角，（°）；是气象站点的纬度，（°）。其中，、ω可以从气象站点的纬度及观测日期（儒略历）计算，公式如下

2.2 土壤数据

土壤数据资料参考使用了中国科学院南京土壤研究所制作的1km空间分辨率的中国土壤类型数据[32-33]，，该数据包含了土壤类型（主要黑土及草甸土）、土壤砂粒含量、土壤粉粒含量、土壤黏粒含量、有机质、全氮、全磷、全钾含量等参数，如图2所示：

a. 砂粒a. Coarse sandb. 粉砂b. Siltc. 黏粒c. Clay

d. 有机质d. Organic mattere. 全氮e. Total nitrogenf. 全磷f. Total phoshporusg.全钾g. Total potassium

2.3 地面实测数据

在玉米作物生长的早期和中期，分别在研究区实地获取春玉米作物及当地的土壤等数据。地面调查工作主要包括地物类型样本点照片采集、土壤和作物生长早期光谱曲线采集、叶面积指数采集、土壤含水率数据采集、土壤样本采集等，获取了作物识别样本点600多个，并采集部分样本点位置土壤样本及测量叶面积指数，分析水分、养分含量（表1），并对当地玉米的物候期、种植结构、作物产量等进行了走访调查。

其中，土壤水分测量包含土壤10 cm深的土壤剖面含水率，即相应深度的土壤水含量。叶面积指数使用LAI2200（美国LI-COR）植物冠层分析仪采集获取。

表1 土壤样本养分、含水率及玉米叶面积指数值统计

2.4 基于MODIS遥感数据的LAI、ET产品数据

下载获取2013年5月1日~2013年10月1日之间的MODIS LAI和ET数据产品。其中，叶面积指数使用MODIS LAI标准产品中的MOD15A2产品，为8 d合成产品，空间分辨率为1 000 m。ET数据使用的是MOD16A2产品，为8 d合成，空间分辨率也是1 000 m。数据统一转换为UTM坐标系，图1b给出了2013年第257天的LAI数据[34]。

在MODIS LAI和ET产品数据中，因为云层、水汽、气溶胶等的影响，造成LAI及ET的数据存在缺失及噪音现象，因此需要对数据产品进行去噪处理，剔除异常数据。本研究使用基于上包络线的Savitzky-Golay（S-G）滤波对LAI和ET数据集进行滤波处理，去除噪声，并对缺失数据进行填补，形成完整的数据集。

同时，根据前人研究，由于MODIS影像混合像元等因素的影响，其LAI产品的峰值大小明显小于地面实测数据。本研究通过利用MODIS LAI数据与地面实测数据进行拟合，对LAI的峰值大小进行调整。MODIS原始LAI峰值数据基本在2～3之间，而玉米的叶面积指数实测值基本在4～7之间，存在显著的低估。由于本文仅获取了两次LAI实测值，无法使用Logistic公式进行LAI观测值的调整。考虑到地面实测的LAI数据包括了生长初期LAI数据和基本成熟时期的LAI指数，而在早期，实测LAI和MODIS观测LAI之间差异并不显著，本研究根据前后两次实测LAI与遥感观测LAI比值的变化趋势，进行LAI调整系数的计算，方法如下：

首先，以第一个时间节点的遥感观测LAI数据作为最小值基准，以第二个时间节点的遥感观测LAI数据作为最大值基准，利用、对第个时间点的遥感观测LAI值x进行归一化，得到归一化后的各时间点LAI观测值l；其次，获取第一时间节点的实测LAI数据，和第二时间节点的实测LAI数据；计算实测值倍数与观测值之间的比值/和；最后根据增大倍数次根号函数拟合的原则，调整LAI增大的倍数，得到最终调整后的观测LAI值y。具体公式计算调整后的LAI值如下

式中为根号的次数，本研究使用=5，使LAI调整值变化两端更加平缓。图3为调整后的LAI值，可以看出，调整后的LAI的大小与实际情况更加吻合，有利于进行模型LAI模拟值和观测值的评估。

2.5 作物空间分布数据

使用Landsat-8时序数据，采用随机森林分类方法，基于地面调查及Google Earth高分辨率影像方法选取样本数据，进行研究区的玉米作物提取，结果如图1所示，总体分类精度为93.2%。同时与2013年研究区各县市的玉米种植面积统计资料进行对比（图4），结果表明，各县市的作物面积监测结果和统计结果基本都在1∶1线附近，2达到了0.951 2。作物空间分布数据主要用于SWAP模型作物区域的裁剪及区域总产量的计算，其计算过程此处不做赘述。

图3 MODIS LAI滤波前后及调整后的值对比

图4 玉米种植面积遥感监测结果与统计资料结果对比

3 研究方法

本文研究的总体技术路线如图5所示，下文将对具体技术内容进行介绍。

3.1 基于coKriging方法的空间气象数据获取

在单点尺度的SWAP模型中，输入的各类气象数据其往往可由当地气象站点直接观测数据得到。然而，在区域尺度连续作物模型应用中，为了获取空间连续分布的气象数据，需要将点分布的气象站点观测数据转化为面状连续的气象数据。本文选用协同克里金（coKriging）方法[30,35-36]，通过引入与气温数据存在高度相关性的海拔数据和纬度数据，将离散分布的气温数据转化为面状数据，进而输入SWAP模型中用于作物生长发育过程模拟。

协同克里金方法的公式如下所示

式中λi和λj是变量的权重系数，分别表示x0处主变量Z1(x1i)和协同变量Z2(x2j)的贡献程度，由观测点的距离、变量变异函数及协变异函数决定。n1和n2分别表示在x0附近用于估值的区域化变量个数。x1i、x2j分别是x0附近的用于估值的第i个主变量值和第j个协同变量值。

在本研究中，考虑到研究区气象站点数量相对较少，为了提高精度，使用了东北及华北一部分的全部气象站点数据参与coKriging空间插值计算。同时，主要对SWAP模型中较为重要的日最高气温和日最低气温使用coKriging插值，考虑的相关变量包括了纬度和DEM数据，DEM数据分辨率为30 m，来源为GDEM（global digital elevation model）。而对于其他气象参数如水汽压、日照时数等气象数据，依然使用普通克里金方法进行插值。

3.2 SWAP模型参数标定

本研究使用SWAP模型作为玉米作物产量监测的主要模型。SWAP模型是荷兰Wageningen University and Research Center开发的用于模拟田块尺度的非饱和带水流及盐分运移、热量传递和作物生长与产量的模型[37-38]。该模型为有作物生长条件下的土壤水分运动模拟模型，适用于多层土壤，并能够考虑地下水动态变化的影响，模型中采用的模拟方法包括详细作物生长模拟、简单作物生长模拟两种方式，详细模拟方式可以模拟作物的生长过程，而简单模拟方式仅可以模拟作物的最终产量。本文使用详细作物模拟方法，其以1 d为步长，模拟作物的生理生态过程，包括作物发育期、CO2同化、呼吸作用、作物蒸腾、干物质分配、土壤水分平衡等。对于研究区玉米参数，其中，比叶面积（specific leaf area，SLA）参数参考刘维等[26]研究结果，TSUM0、TSUM1、TSUM2参数则通过作物物候期结合气温数据计算获取，TDWI、LAIEM等参数通过实地调查获取，RGRLAI、SPAN、CVL、CVO、CVR、CVS等参数则主要参考模型默认参数及李琪等[39]、陈思宁等[11]、张素青等[25]研究结果并通过“试错法”获取[40]，其余如DVSI、DVSEND 、RML、RMO、RMR、RMS参数使用SWAP模型玉米默认参数，各参数如表2和表3所示，表中DVS（development stage）指的是作物生长发育阶段，无量纲。

表2 研究区玉米比叶面积参数设置

表3 研究区玉米SWAP模型参数设置

3.3 遥感数据与模型同化方法

使用SWAP模型，通过同化MODIS LAI及ET数据产品，获取区域尺度高精度的玉米作物单产信息。SWAP 模型输入的参数主要有气象、土壤、作物和田间管理等，各参数对作物的发育和产量形成影响的差异不同，该研究仅选择那些对作物产量较为敏感的参数进行优化。针对玉米，主要选择灌溉量和出苗日期作为模型的主要优化参数。一般而言，春玉米的出苗日期在有限区域内较为一致，根据调查及查阅相关资料[41]，该文将研究区春玉米出苗日期初始值设置为5月20日；而研究区的春玉米灌溉制度较为复杂，一般在湿润年份以雨养为主，在一般年份则会在拔节期或抽穗期进行少量灌溉[42]。为简化模型运算，该文设置SWAP模型春玉米灌溉日期为6月30日，灌溉量初始值设置为70 mm。在SWAP运行过程中，将根据同化方法和代价函数，不断尝试修正出苗日期及灌溉量，使得模型模拟LAI及ET结果与遥感数据趋势尽可能一致，并最终实现区域产量的同化反演。

该文中，同化方法选择使用SCE-UA方法，该方法由美国亚利桑那大学Duan等[43-44]于20世纪90年代提出，该算法是下山单纯形算法的发展，采用多个单纯形并行地搜索解空间的策略，这种策略被证明有助于克服下山单纯形算法可能会收敛于局部最小的缺点[45]。SCE-UA算法可以有效地解决高维度参数的全局性优化问题，能快速获取全局最优解并避免陷入局部最优现象，具有良好的全局优化性能和效率[43]。

同化代价函数的选取方面，由于MODIS的LAI产品数据存在较高的低估，且本文校正LAI的方法仅有两个生育期的数据，难以有效保证LAI值的准确性，因此模型同化的代价函数选用了黄健熙等[46]提出的基于一阶差分的代价函数。该方法利用遥感LAI监测数值与模型模拟LAI之间的单调性作为目标函数，通过最小化两者之间的单调性差异，达到优化模型输入参数。这一方法可以避免直接对观测LAI和模拟LAI值的大小的比较，从而避免因为遥感LAI低估可能造成的模型误差。其计算公式如式（11）所示，Dif为第个观测值或模拟值处一阶差分的正负性，则代表LAI或ET，threshold为阈值，取决于遥感数据的精度。获取这一数据后，将所有观测时间的观测值正负性与模拟值正负性进行对比，两者越接近则表明参数优化效果越好[46]。

使用SCE-UA方法不断调整SWAP模型的输入参数（本文为出苗日期及灌溉量），当代价函数达到最小或迭代超过一定次数限制（1 000次）时，所获得的参数即为优化参数。使用优化参数输入模型，最后输出结果即为玉米作物的产量信息。值得注意的是，SWAP模型能在作物成熟前，以1 d间隔输出作物实时产量数据，可用来评价不同时间段作物的长势情况。

3.4 精度验证

精度验证使用黑龙江统计局编著的《黑龙江县（市）农村经济社会统计概要2014》中黑龙江各县市农作物种植面积和产量数据作为基准资料，计算包括春玉米面积识别精度及春玉米产量反演平均绝对误差（mean absolute error，MAE）、均方根误差（root mean square error，RMSE）及决定系数2。

4 结果与分析

4.1 气象数据空间化结果

使用coKriging方法，在气象站点日最高气温、日最低气温基础数据上，考虑DEM及纬度因子，插值获取气温数据的空间分布栅格数据。同时，也使用普通克里金方法，对气温数据进行了插值，并使用舍一法交叉验证对插值的精度进行对比。图6和图7是6月10日最高气温及最低气温普通克里金插值及协同克里金插值的结果，可以看到，两者直接存在较为显著的差异，且协同克里金呈现相对更加复杂细致的气温分布趋势。气温整体呈现从东南至西北逐渐降低的趋势。根据舍一法交叉验证结果表明，最低气温普通克里金插值结果的平均标准误差（average standard error）为1.51 ℃，而对应的协同克里金法的平均标准误差仅为0.31 ℃；最高气温普通克里金插值结果的平均标准误差为1.14 ℃，而对应的协同克里金法的平均标准误差仅为0.30 ℃。结果表明，相比普通克里金方法，考虑DEM和纬度因素的协同克里金方法可以获得更加精确的气温空间化插值数据产品，这对于依赖于气温数据作为驱动作物生长发育期计算模拟的作物生长模型具有重要的意义，有利于提高作物生长模型输入数据的精度，进而提高模型作物长势和反演精度。

图6 2013年6月10日最高气温插值结果

图7 2013年6月10日最低气温插值结果

4.2 玉米长势及产量模拟结果

根据SWAP模型同化LAI和ET数据，获取研究区逐日玉米产量预测结果，分别选取位于生长早期的2013年8月1日和旺盛时期的2013年8月19日的玉米产量预测结果，并使用Natural break方法，将其分为5级，用于分析不同时期研究区春玉米的长势模型模拟结果（图8）。

从图中可以看出，在2013年8月1日，研究区的西南部，主要包括双城以及肇东市南部地区的玉米长势相对最好，而在研究区的中部地区，长势较好，研究区的北部则长势一般。8月1日的长势情况整体呈现从南到北逐渐降低的过程，这主要是由于研究区春玉米物候期差异决定的，受积温差异影响，南部玉米生长相对较早，而北部相对较晚。而在2013年8月19日，随着玉米长势的逐渐变化，研究区的中部玉米长势逐渐变好，长势分级出现西南-东北逐渐降低的趋势。

图8 研究区春玉米长势模型模拟结果

根据SWAP模型同化遥感数据获取研究区玉米的单产分布图，结果如图9所示，可以看出，研究区玉米单产呈现出较强的区域差异。统计玉米平均单产，得到研究区玉米单产均值为7 831 kg/hm2。从图上可以发现，研究区中部地区单产最高，达到了8 800 kg/hm2以上，而在南部的双城北部，哈尔滨市辖区南部，以及中部的呼兰、兰西、北林交界区域玉米单产相对较低，在6 118 ~7 756 kg/hm2之间。分析其主要原因，可能是由于该区域8月份中旬出现的强降雨天气导致了最后阶段玉米产量受到一定的影响。

收集研究区15个县市区2013年玉米产量的统计信息，并根据SWAP模型拟合获取的区域作物产量空间分布图，统计各行政区划内的玉米产量（包括单产和总产），对比统计产量和预测产量，评价模型精度。结果如图10所示，从图中可以看出，春玉米的预测产量与统计产量较为接近，2达到0.939 4，均方根误差RMSE为148 065 t，平均绝对误差MAE为114 335 t。在单产上，研究区15个县市区的模型预测单产与统计结果之间的决定系数达到了0.724 5，RMSE为598.5 kg/hm2，MAE为531.5 kg/hm2。精度评价表明基于SWAP模型同化遥感数据进行黑龙江南部区域尺度玉米作物产量估算的准确性和适应性。

图9 2013年研究区玉米单产SWAP模型同化预测结果空间分布图

图10 2013年研究区各县市玉米产量预测与统计结果对比

5 讨 论

SWAP作物生长模型在单点尺度的作物生长模拟中有着重要的应用，以其高效率、高精度、机理明确等优点，为田间管理、水肥控制等技术研究提供了一种便捷、低成本的技术手段；而遥感数据具有覆盖范围大、获取成本低等优势；通过将面状的遥感数据与传统的单点尺度作物生长模型进行耦合，结合大尺度、区域化的模型参数优化获取，可以实现作物生长模型在区域尺度的应用。该文以SWAP模型区域化应用为出发点，从遥感数据获取、遥感数据预处理、农作物类型空间分布图制作、气象数据空间化、模型参数选择，以及数据同化、结果精度验证等方面进行了完整的研究与分析，并获得了较高的监测精度。

由于农作物类型地理空间分布及其影响因子的复杂性，更大尺度上，更为通用的遥感数据与作物生长模型耦合的数据、方法等方面还有待更为深入的研究。在同化技术所需的遥感数据获取方面，该文选用了MODIS LAI及ET遥感产品数据，考虑到MODIS产品的分辨率相对较低，在精细农业中的应用有一定的局限；在今后的研究中，可进一步结合使用中高分辨率的卫星遥感数据如HJ-1、Landsat等进行作物生长模型同化研究[47-48]；并且MODIS LAI产品在相比实际作物LAI值往往存在一定的低估[49]，如果不选用合适的预处理技术或同化技术，可能会对同化结果造成影响，当前有一些研究逐渐开始尝试结合冠层反射率[50]、冠层温度[51]以及土壤调节植被指数SAVI（soil adjusted vegetion index）等指数[52]来改善类似的缺点。

在气象数据插值方面，当前作物生长模型的气象数据处理过程中，往往直接使用较为简单、快捷的反距离权重插值或者普通克里金插值，仅仅考虑了气象数据的空间自相关性。在大尺度区域应用中，气象数据往往会与其他因子存在较强的相关性[36,53]，如气温与纬度、海拔的关系密切[54-55]，不考虑这一因素，可能导致插值结果与实际的偏离。针对气温这个表示作物生长发育阶段的重要参数，该文在研究中使用了协同克里金插值代替普通克里金插值方法，最低气温插值结果的平均标准误差由1.51 ℃降低到0.31 ℃，最低气温插值结果的平均标准误差由1.14 ℃降低到0.30 ℃，表明该方法可以有效提高气温数据插值的精度，从而为作物生长模型提供精确的输入数据。在后续的研究中，可以尝试结合NDVI、纬度、TRMM（tropical rainfall measuring mission）降雨数据等，使用协同克里金、地理加权回归克里金（Geographically weighted regression Kriging）等方法[49]，进一步提高作物生长模型所需的降水等数据的精度，以达到进一步提高区域作物生长模型对作物产量监测与预测的精度。

6 结 论

本文以中国玉米的主产区之一，黑龙江省南部地区为研究区，对SWAP模型春玉米模拟进行了参数的本地化研究，并使用SCE-UA同化算法，基于MODIS的LAI和ET数据集，使用一阶差分代价函数，进行了作物生长模型与遥感数据的同化。结果表明，研究区内15个县、市、区的玉米产量与统计产量直接基本吻合，总产量的2为0.939 4，均方根误差RMSE为148 065 t，平均绝对误差MAE为114 335 t，单产的决定系数达到了0.724 5，RMSE为598.5 kg/hm2，MAE为531.5 kg/hm2。表明基于SWAP模型同化MODIS数据进行黑龙江南部区域春玉米产量监测的准确性和适用性。

[1]中国统计出版社. 中国农村统计年鉴[M]. 北京：中国统计出版社，2017.

[2]余晓洋，田帅，吕一，等. 基于Shift-Share模型的“镰刀弯”地区粮食增产结构变动分析：兼论玉米种植结构调整[J].中国农业资源与区划，2017，38(5)：14－20. Yu Xiaoyang, Tian Shuai, Lv Yi, et al. Structure variation of grain yield in “sickle bend” region based on shift-share model[J]. Chinese Journal of Agricultural Resources and Regional Planning, 2017, 38(5): 14－20. (in Chinese with English abstract)

[3]刘佳，王利民，姚保民，等. 基于多时相OLI数据的宁夏大尺度水稻面积遥感估算[J]. 农业工程学报，2017，33(15)：200－209. Liu Jia, Wang Limin, Yao Baomin, et al. Ningxia rice area remote sensing estimation on large scale based on multi-temporal OLI data[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(15): 200－209. (in Chinese with English abstract)

[4]Massey Richard, Sankey Temuulen T, Congalton Russell G, et al. MODIS phenology-derived, multi-year distribution of conterminous US crop types[J]. Remote Sensing of Environment, 2017, 198: 490－503.

[5]Zhong Liheng, Gong Peng, Biging Gregory S. Efficient corn and soybean mapping with temporal extendability: A multi-year experiment using Landsat imagery[J]. Remote Sensing of Environment, 2014, 140: 1－13.

[6]黄健熙，黄海，马鸿元，等. 遥感与作物生长模型数据同化应用综述[J]. 农业工程学报，2018，34(21)：144－156. Huang Jianxi, Huang Ha, Ma Hongyuan, et al. Review on data assimilation of remote sensing and crop growth models[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(21): 144－156. (in Chinese with English abstract)

[7]Diepen C A, Wolf J, Keulen H, et al. WOFOST: A simulation model of crop production[J]. Soil Use & Management, 2010, 5(1): 16－24.

[8]Shafiei Mojtaba, Ghahraman Bijan, Saghafian Bahram, et al. Uncertainty assessment of the agro-hydrological SWAP model application at field scale: A case study in a dry region[J]. Agricultural Water Management, 2014, 146: 324－334.

[9]Ovando Gustavo, Sayago Silvina, Bocco Mónica. Evaluating accuracy of DSSAT model for soybean yield estimation using satellite weather data[J]. Isprs Journal of Photogrammetry & Remote Sensing, 2018, 138: 208－217.

[10]Huang Jianxi, Ma Hongyuan, Liu Junming, et al. Regional winter wheat yield estimation by assimilating MODIS ET and LAI products into SWAP model [C]// Agro-Geoinformatics 2013 Second International Conference on, IEEE, 2013. (Agro-Geoinformatics).

[11]陈思宁，赵艳霞，申双和. 基于集合卡尔曼滤波的PyWOFOST模型在东北玉米估产中的适用性验证[J]. 中国农业气象，2012，33(2)：245－253. Chen Sining, Zhao Yanxia, Sheng Shuanghe. Applicability of PyWOFOST model based on ensemble Kalman filter in Simulating maize yield in Northeast China[J]. Chinese Journal of Agrometeorology, 2012, 33(2): 245－253. (in Chinese with English abstract)

[12]黄健熙，武思杰，刘兴权，等. 基于遥感信息与作物模型集合卡尔曼滤波同化的区域冬小麦产量预测[J]. 农业工程学报，2012，28(4)：142－148. Huang Jianxi, Wu Sijie, Liu Xinquan, et al. Regional winter wheat yield forecasting based on assimilation of remote sensing data and crop growth model with Ensemble Kalman method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(4): 142－148. (in Chinese with English abstract)

[13]Cheng Zhiqiang, Meng Jihua, Wang Yiming. Improving spring maize yield estimation at field scale by assimilating time-series HJ-1 CCD data into the WOFOST model using a new method with fast algorithms[J]. Remote Sensing, 2016, 8(4): 303.

[14]姜浩. 基于作物模型同化遥感物候信息的冬小麦估产方法研究[D]. 成都：电子科技大学，2011. Jiang Hao. Study on Winter Wheat Yield Estimation Method Based on Remote Sensing Phenological Information Assimilated by Crop Model[D]. Chengdu: University of Electronic Science and Technology of China, 2011. (in Chinese with English abstract)

[15]李秀芬，马树庆，宫丽娟，等. 基于WOFOST的东北地区玉米生育期气象条件适宜度评价[J]. 中国农业气象，2013，34(1)：43－49. Li Xiufen, Ma Shuqing, Gong Lijuan, et al. Evaluation of meteorological suitability degree during maize growth period based on WOFOST in Northeast China[J]. Chinese Journal of Agrometeorology, 2013, 34(1): 43－49. (in Chinese with English abstract)

[16]栾庆祖，叶彩华，莫志鸿，等. 基于WOFOST模型的玉米干旱损失评估：以北京为例[J]. 中国农业气象，2014，35(3)：311－316. Luan Qingzu, Ye Caihua, Mo Zhihong, et al. Maize yield loss assessment for drought based on WOFOST model: A case study in Beijing[J]. Chinese Journal of Agrometeorology, 2014, 35(3): 311－316. (in Chinese with English abstract)

[17]Cheng Zhiqiang, Meng Jihua, Qiao Yanyou, et al. Preliminary study of soil available nutrient simulation using a modified WOFOST model and time-series remote sensing observations[J]. Remote Sensing, 2018, 10(1): 64.

[18]Ceglar Andrej, Kajfež-Bogataj Lučka. Simulation of maize yield in current and changed climatic conditions: Addressing modelling uncertainties and the importance of bias correction in climate model simulations[J]. European Journal of Agronomy, 2012, 37(1): 83－95.

[19]Huang Jianxi, Tian Liyan, Liang Shunlin, et al. Improving winter wheat yield estimation by assimilation of the leaf area index from Landsat TM and MODIS data into the WOFOST model[J]. Agricultural & Forest Meteorology, 2015, 204: 106－121.

[20]Zhao Yanxia, Chen Sining, Shen Shuanghe. Assimilating remote sensing information with crop model using ensemble Kalman filter for improving LAI monitoring and yield estimation[J]. Ecological Modelling, 2013, 270(2): 30－42.

[21]王丽媛. 遥感数据与作物模型同化的冬小麦估产研究[D]. 杭州：浙江大学，2018. Wang Liyuan. Estimating Winter Wheat Yield by Assimilation of Remote Sensing Data into Crop Growth Model[D]. Hangzhou: Zhejiang University, 2018. (in Chinese with English abstract)

[22]Huang Jianxi, Ma Hongyuan, Su Wei, et al. Jointly assimilating MODIS LAI and ET products into the SWAP model for winter wheat yield estimation[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(8): 4060－4071.

[23]Dong Taifeng, Liu Jiangui, Qian Budong, et al. Estimating winter wheat biomass by assimilating leaf area index derived from fusion of Landsat-8 and MODIS data[J]. International Journal of Applied Earth Observation and Geoinformation, 2016, 49: 63－74.

[24]陈劲松，黄健熙，林珲，等. 基于遥感信息和作物生长模型同化的水稻估产方法研究[J]. 中国科学：信息科学，2010，40(S1)：173－183. Chen Jingsong, Huang Jianxi, Lin Hun, et al. Rice yield estimation by assimilation remote sensing into crop growth model[J]. Science China Information Sciences, 2010, 40(S1): 173－183. (in Chinese with English abstract)

[25]张素青，张建涛，李继蕊，等. WOFOST模型在河南省夏玉米主产区的校准与验证[J]. 河南农业科学，2014，43(8)：152－156. Zhang Suqin, Zhang Jiantao, Li Jirui, et al. Calibration and validation of WOFOST in main maize-producting region in Henan[J]. Journal of Henan Agricultural Sciences, 2014, 43(8): 152－156. (in Chinese with English abstract)

[26]刘维，侯英雨，吴门新，等. WOFOST模型在东北春玉米产区的验证与适应性评价[J]. 气象与环境科学，2017，40(3)：9－15. Liu Wei, Hou Yingyu, Wu Menxing, et al. Validation and adaptability evaluation of WOFOST model in spring maize area of Northeast[J]. Meteorological and Environmental Sciences, 2017, 40(3): 9－15. (in Chinese with English abstract)

[27]Wang Liyuan, Huang Jingfeng, Gao Ping, et al. Estimating winter wheat yield by assimilation of MODIS LAI into WOFOST model with Ensemble Kalman Filter[C]// 2017 6th International Conference on Agro-Geoinformatics, IEEE, 2017.

[28]黄健熙，马鸿元，田丽燕，等. 基于时间序列LAI和ET同化的冬小麦遥感估产方法比较[J]. 农业工程学报，2015，31(4)：197－203. Huang Jianxi, Ma Hongyuan, Tian Liyan, et al. Comparison of remote sensing yield estimation methods for winter wheat based on assimilating time-sequence LAI and ET[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(4): 197－203. (in Chinese with English abstract)

[29]石淑芹. 基于多源数据的吉林省玉米生产力区划研究[D]. 北京：中国农业科学院，2009. Shi Shuqin. Maize Productivity Regionalization in Jilin Province Based on Multiple Source Data[D]. Beijing: Chinese Academy of Agricultural Sciences, 2009. (in Chinese with English abstract)

[30]卢月明. 面向大气污染指数分析的时空协同克里金插值方法[D]. 北京：中国测绘科学研究院，2018. Lu Yueming. Spatio-temporal Cokriging Interpolation for Air Pollution Index Analysis[D]. Beijing: Chinese Academy of Surveying and Mapping, 2018. (in Chinese with English abstract)

[31]Allen Richard G, Pereira Luis S, Raes Dirk, et al. Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56[Z]. FAO, Rome, 1998, 300(9): D5109.

[32]史学正，于东升，高鹏，等. 中国土壤信息系统(SISChina)及其应用基础研究[J]. 土壤，2007，39(3):329－333. Shi Xuezheng, Yu Dongsheng, Gao Peng, et al. Soil information system of China (SISChina) and its application[J]. Soils, 2007, 39(3): 329－333. (in Chinese with English abstract)

[33]国家科技基础条件平台—国家地球系统数据共享服务平台-土壤科学数据中心：http://soil.geodata.cn.

[34]NASA LANDS DAAC[EB/OL]. [2019-05-26]. https://ladsweb. modaps.eosdis.nasa.gov/search/

[35]Guo Qingyan, Zhang Min, Yang Xiandong, et al. Spatial interpolation of monthly mean temperatures based on cokriging method[C]// International Symposium on Computational Science and Computing, Springer. 2018.

[36]Park Sunyurp. Integration of satellite‐measured LST data into cokriging for temperature estimation on tropical and temperate islands[J]. International Journal of Climatology, 2011, 31(11): 1653－1664.

[37]Wageningen University and Research--Soil Water Atmosphere Plant[EB/OL]. [2019-05-26] https://www.swap.alterra.nl/.

[38]Kroes J G, van Dam J C, Bartholomeus R P, et al. SWAP version 4; Theory description and user manual[ZOL]. Wageningen, Wageningen Environmental Research, 2017, Report 2780. [2019-05-26]Available at: http://library.wur.nl/ WebQuery/wurpubs/fulltext/416321

[39]李琪，胡秋丽，朱大威，等. 基于WOFOST模型的吉林省春玉米干旱复水模拟研究[J]. 农业现代化研究，2019，40(1)：155－162. Li Qi, Hu Qiuli, Zhu Dawei, et al. Simulation of drought and rewatering of spring maize based on WOFOST model in Jilin Province[J]. Research of Agricultural Modernization, 2019, 40(1): 155－162. (in Chinese with English abstract)

[40]蔡福，米娜，纪瑞鹏，等. 基于锦州春玉米田间试验的WOFOST模型参数的确定及性能评价[J]. 生态学杂志，2019，38(4)：1238－1248. Cai Fu, Mi Na, Ji Ruipeng, et al. Determination of crop parameters for WOFOST model and its performance evaluation based on field experiment of spring maize in Jinzhou, Liaoning[J]. Chinese Journal of Ecology, 2019, 38(4): 1238－1248. (in Chinese with English abstract)

[41]李正国，杨鹏，唐华俊，等. 近20年来东北三省春玉米物候期变化趋势及其对温度的时空响应[J]. 生态学报，2013，33(18)：5818－5827. Li Zhenguo, Yang Peng, Tang Huajun, et al. Trends of spring maize phenophases and spatio-temporal responses to temperature in three provinces of Northeast China during the past 20 years[J]. Acta Ecologica Sinica, 2013, 33(18): 5818－5827. (in Chinese with English abstract)

[42]肖俊夫，刘小飞，刘祖贵，等. 我国春玉米产区灌溉问题分析与研究[C]//北京，玉米产业技术大会，2008.

[43]Gupta V K. Optimal use of the SCE-UA global optimization method for calibrating watershed models[J]. Journal of Hydrology, 1994, 158(3/4): 265－284.

[44]Duan Qingyun. The Shuffled Complex Evolution (SCE-UA) Method[D]. Tucson: University of Arizona Tucson, 1992.

[45]栾承梅. 流域水文模型参数优化问题研究[D]. 南京：河海大学，2006. Luan Chengmei. Study on Automatic Optimization of Hydrologic Models' Parameters[D]. Nanjing: Hohai University, 2006. (in Chinese with English abstract)

[46]黄健熙，马鸿元，田丽燕，等. 基于时间序列LAI和ET同化的冬小麦遥感估产方法比较[J]. 农业工程学报，2015，31(4)：197－203. Huang Jianxi, Ma Hongyuan, Tian Liyan, et al. Comparison of remote sensing yield estimation methods for winter wheat based on assimilating time-sequence LAI and ET[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(4): 197-203. (in Chinese with English abstract)

[47]陈艳玲，顾晓鹤，宫阿都，等. 基于遥感信息和WOFOST模型参数同化的冬小麦单产估算方法研究[J]. 麦类作物学报，2018，38(9)：121－130. Chen Yanling, Gu Xiaohe, Gong Adou, et al. Estimation of winter wheat assimilation based on remote sensing information and WOFOST crop model[J]. Journal of Triticeae Crops, 2018, 38(9): 121－130. (in Chinese with English abstract)

[48]Huang Jianxi, Tian Liyan, Liang Shunlin, et al. Improving winter wheat yield estimation by assimilation of the leaf area index from Landsat TM and MODIS data into the WOFOST model[J]. Agricultural & Forest Meteorology, 2015, 204: 106－121.

[49]徐胜男. 基于遥感数据修正WOFOST模型的江苏省冬小麦生长模拟研究[D]. 南京：南京信息工程大学，2016. Xu Shennan. Simulation of Winter Wheat in Jiangsu Province with the Modified WOFOST Model Based on Remote Sensing Data[D]. Nanjing: Nanjing University of Information Science & Technology, 2016. (in Chinese with English abstract)

[50]Bouman B A M. Linking physical remote sensing models with crop growth simulation models, applied for sugar beet[J]. International Journal of Remote Sensing, 1992, 13(14): 2565－2581.

[51]王纯枝，宇振荣，辛景峰，等. 基于遥感和作物生长模型的作物产量差估测[J]. 农业工程学报，2005，21(7)：84－89. Wang Chunzhi, Yu Zhenrong, Xin Jingfeng, et al. Yield gap estimation by combining remote sensing and crop growth model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2005, 21(7): 84－89. (in Chinese with English abstract)

[52]吴伶，刘湘南，周博天，等. 多源遥感与作物模型同化模拟作物生长参数时空域连续变化[J]. 遥感学报，2012，16(6)：1173－1191. Wu Ling, Liu Xiangnan, Zhou Botian, et al. Spatial-time continuous changes simulation of crop growth parameters with multi-source remote sensing data and crop growth model[J]. Journal of Remote Sensing, 2012, 16(6): 1173－1191. (in Chinese with English abstract)

[53]Chen Yuanyuan, Huang Jingfeng, Sheng Shaoxue, et al. A new downscaling-integration framework for high-resolution monthly precipitation estimates: Combining rain gauge observations, satellite-derived precipitation data and geographical ancillary data[J]. Remote Sensing of Environment, 2018, 214: 154－172.

[54]程朋，任鸿瑞. 山西高原气温空间插值分析[J]. 河北师范大学学报：自然科学版，2016，40(1)：73－79. Cheng Peng, Ren Hongrui. Spatial interpolation analysis of annual average air temperature in Shanxi Plateau[J]. Journal of Hebei Normal University: Natural Science Edition, 2016, 40(1): 73－79. (in Chinese with English abstract)

[55]Ishida T, Kawashima S. Use of cokriging to estimate surface air temperature from elevation[J]. Theoretical and Applied Climatology, 1993, 47(3): 147－157.

Maize yield monitoring in Southern Heilongjiang based on SWAP model assimilative remote sensing data

Wang Limin1, Yao Baomin1, Liu Jia1※, Yang Lingbo2, Yang Fugang1

(1.,,, 100081; 2.,310012)

Crop yield monitoring and forecast has great impact on food security, ecological environment, and farmers’ incomes. Crop growth and yield monitoring and forecast by using crop growth model has the advantages of clear mechanism, high precision, and high monitoring frequency, but its monitoring scale is usually limited to land block level. Along with the fast development of remote sensing satellite technology, using remote sensing data combined with crop growth model to accurately monitor regional crop growth and yield in a large regional scale has gradually become an important means of regional level and even national level crop growth monitoring. However, due to differences in crop types, climatic conditions, soil conditions and monitoring areas, the regionalization and localization of crop growth model is the major bottleneck of crop growth monitoring by using crop growth model combined with remote sensing data, and it is urgent to conduct targeted studies on the identification of assimilation parameters of crop growth model, pre-processing of meteorological data, and the setting of crop parameters. Based on soil-water-atmosphere-plant model (SWAP), and by taking the major commodity grain production base of China, Northeast China Region as a study region, in this the paper, we conducted a study by taking the major crop of spring maize of the region as its target crop. Firstly, we used Landsat to obtain maize () planting area in the study area, and used it as the basic data for estimating the total maize yield in the study area. The overall accuracy of maize area classification was 93.2%, with² of 0.951 2. By considering the influence of latitude and altitude on temperature, in the study, we used the coKriging method in crop growth model meteorological data interpolation acquisition, so as to improve the precision of input parameters of the model. The result showed that the average standard error of minimum temperature of coKriging method was 0.31 ℃, while that of the Kriging method was 1.51 ℃. The average standard error of maximum temperature of coKriging method was 0.30 ℃, while that of the Kriging method was 1.14 ℃. In the study, leaf area index (LAI) and evapotranspiration (ET) were used as assimilative remote sensing data sources, and we proposed a novel method to adjust the LAI product of MODIS to make it closer to actual value. By optimizing maize irrigation and crop emergence date, we obtained spatial distribution result of maize yield of the study area of 2013. The monitoring result was compared with the statistical data. The² reached 0.939 4, with RMSE of 148 065 t, and MAE of 114 335 t. Moreover, the correlation coefficient of predicted yield and statistical yield reached 0.724 5, with RMSE of 598.5 kg/hm², and MAE of 531.5 kg/ hm²。The study result showed that, using SWAP model, taking meteorological data spatial interpolation results obtained by using coKriging method as input data and assimilation of LAI and ET remote sensing, can effectively conduct corn yield remote sensing monitoring of the study region, which provided reference for the remote sensing monitoring and forecast of crop growth and productivity of the region.

crops; remote sensing; models; corn; leaf area index (LAI); evapotranspiration; coKriging; assimilation; SWAP model

王利民，姚保民，刘 佳，杨玲波，杨福刚. 基于SWAP模型同化遥感数据的黑龙江南部春玉米产量监测[J]. 农业工程学报，2019，35(22)：285－295. doi：10.11975/j.issn.1002-6819.2019.22.034 http://www.tcsae.org

Wang Limin, Yao Baomin, Liu Jia, Yang Lingbo, Yang Fugang. Maize yield monitoring in Southern Heilongjiang based on SWAP model assimilative remote sensing data[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(22): 285－295. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2019.22.034 http://www.tcsae.org

2019-04-19

2019-06-26

国家重点研发计划“粮食作物生长监测诊断与精确栽培技术”课题“作物生长与生产力卫星遥感监测预测”（2016YFD0300603）

王利民，博士,研究员，主要从事农业遥感监测业务运行研究。Email：wanglimin01@caas.cn

刘 佳，研究员，主要从事农业遥感监测业务运行研究。Email：liujia06@caas.cn

10.11975/j.issn.1002-6819.2019.22.034

S252+.9

A

1002-6819(2019)-22-0285-11