EN中文

地球信息科学学报 >

2015 , Vol. 17 >Issue 7: 871 - 882

DOI: https://doi.org/10.3724/SP.J.1047.2015.00871

同化叶面积指数和蒸散发双变量的冬小麦产量估测方法

  • 包姗宁 , 1, 2 ,
  • 曹春香 1 ,
  • 黄健熙 , 3**, * ,
  • 田丽燕 3 ,
  • 马鸿元 3 ,
  • 苏伟 3 ,
  • 倪希亮 1
展开
  • 1. 中国科学院遥感与数字地球研究所 遥感科学国家重点实验室,北京 100101
  • 2. 中国科学院大学,北京 100049
  • 3. 中国农业大学信息与电气工程学院,北京 100083
*通讯作者：黄健熙(1976-),男,副教授,博士生导师,主要从事农业定量遥感等方面的研究。E-mail:

作者简介：包姗宁(1991-),女,硕士生,辽宁阜新人,研究方向为环境健康遥感诊断。E-mail:

收稿日期: 2014-09-19

  要求修回日期: 2015-01-19

  网络出版日期: 2015-07-08

基金资助

国家自然科学基金项目(41371326)

水利部公益性行业科研专项经费项目(1261430110032)

收起

Research on Winter Wheat Yield Estimation Based on Assimilation of Leaf Area Index and Evapotranspiration Data

  • BAO Shanning , 1, 2 ,
  • CAO Chunxiang 1 ,
  • HUANG Jianxi , 3, * ,
  • MA Hongyuan 3 ,
  • TIAN Liyan 3 ,
  • SU Wei 3 ,
  • NI Xiliang 1
Expand
  • 1. State Key Laboratory of Remote Sensing Science, The Institute of Remote Sensing and Digital Earth of Chinese Academy of Sciences, Beijing 100101, China
  • 2. University of Chinese Academy of Sciences, Beijing 100049, China
  • 3. College of Information and Electrical Engineering, China Agricultural University, Beijing 100083, China
*Corresponding author: HUANGjianxi, E-mail:

Received date: 2014-09-19

  Request revised date: 2015-01-19

  Online published: 2015-07-08

Copyright

《地球信息科学学报》编辑部 所有
Fold

摘要

同化遥感信息到作物生长过程模拟模型,是估测区域作物产量的重要方法之一。同化变量的选取对同化结果精度至关重要。本文在标定WOFOST作物模型参数的基础上,优化了WOFOST模型的默认灌溉参数。利用ET和LAI作为同化变量,分别构建了时间序列趋势信息的代价函数和四维变分代价函数;采用SCE-UA算法最小化代价函数, 重新初始化WOFOST模型初始参数——作物初始干物质重、作物35 ℃生命期和灌溉量。最后利用MODIS LAI产品（MCD15A3）、MODIS ET产品（MOD16A2）,同化到作物模型估测产量,并对比分析了水分胁迫模式下同化单变量（ET或LAI）和同化双变量（ET和LAI）的估产精度。结果表明：同化双变量ET和LAI的策略,优于同化单变量LAI或ET,双变量策略的冬小麦产量估测精度为R2=0.432,RMSE=721 kg/hm2;单独同化高精度LAI对提高估产精度具有重要作用,其冬小麦产量估测精度为R2=0.408,RMSE=925 kg/hm2;单独同化ET的趋势信息改善了WOFOST模型模拟水分平衡的参数,但是,产量估测精度（R2=0.013,RMSE=1134 kg/hm2）与模型模拟估测产量精度（R2=0.006,RMSE=1210 kg/hm2）相比改善效果有限。本研究为其他区域的遥感数据与作物模型的双变量数据同化的作物产量估测研究提供了参考价值。

关键词： 蒸散发; 叶面积指数; 同化; WOFOST作物模型; 估产

本文引用格式

包姗宁 , 曹春香 , 黄健熙 , 田丽燕 , 马鸿元 , 苏伟 , 倪希亮 . 同化叶面积指数和蒸散发双变量的冬小麦产量估测方法[J]. 地球信息科学学报, 2015 , 17(7) : 871 -882 . DOI: 10.3724/SP.J.1047.2015.00871

Abstract

Assimilating remote sensing information into crop growth model is an important approachto estimate regional crop yield. The assimilation algorithm and corresponding assimilation variables are the keys of the assimilation system, which greatly impact the accuracy of assimilation results. In thepaper, the default irrigation parameters of WOFOST were optimized firstly with the help of calibrating WOFOST crop model parameters. Then, ET data was chosen as the assimilation variable to build the cost function of time series trends using MODIS ET products (MOD16A2) and WOFOST simulation. And LAI data was assimilatedwith the cost function of four dimensional variational data assimilation method using MODIS LAI (MCD15A3) products and WOFOST simulation. Furthermore, parameters including the crop initial dry matter(TDWI), the lifetime of crop in 35℃(SPAN) and the irrigation(RIRR) were optimizedcontinuously by SCE—UA algorithm, which would stop running the program when the cost function isoptimal. The estimatedcrop yield results were obtained usingfour methods comparatively under water limited mode, including the method thm2t does not assimilate and methods thm2t assimilateET, assimilateLAI and assimilate both ET and LAI.We address the assimilation of double variables to be the methodthm2t ET and LAI arebothm2ssimilated. Finally, the accuracies of yield estimation by assimilating double variables and a single variable under water limited mode were compared and analyzed. The results indicated thm2t the method of assimilating double variables was better than assimilating a single variable, which got the highest accuracy (R2=0.432, RMSE=721 kg/hm2). The method of assimilating high precision LAIsignificantly improved the accuracy of yield estimation (R2=0.408, RMSE=925 kg/hm2). The method of assimilating ET demonstratedbetter performance when the WOFOST model simulates the water balance during crop growing period, but had a limited impacton improving the accuracy of yield estimation (R2=0.013, RMSE=1134 kg/hm2) compared with model simulation (R2=0.006, RMSE=1210 kg/hm2). This research provided a reference for studies in other areas on predicting crop production at regional scale thm2t based on assimilating double variables.

Key words： evapotranspiration; leaf area index; assimilation; WOFOST crop model; yield estimation

1 引言

极端气候频繁发生和水资源缺乏,是影响我国粮食安全和农业可持续发展的主要问题之一。华北平原作为我国的冬小麦主产区,粮食总量占全国小麦总产量的56%左右,但长期以来遭受以干旱为主的农业气象灾害,使产量受到严重影响,因此,开展水分胁迫模式下冬小麦产量估测具有重要意义。
估测作物产量的方法之一是直接应用作物生长过程模拟模型（简称作物模型）,这种模型运用数理公式揭示作物生长发育的动力学机制和物理过程,擅长在给定农业资料的条件下模拟田间尺度作物生长发育过程及其不同水分条件下的预测产量。然而,在区域尺度方面,由于空间范围广、地表环境非均匀,一些宏观资料和模型参数难以获取,致使作物模型的应用受到限制。估测产量的另一种方法是运用遥感植被指数与产量构建的经验模型预测作物产量[1],这种方法利用遥感信息空间连续和时间动态变化的优势,监测作物宏观状态和反映环境因子对作物的综合影响[2-5]。但是,卫星数据受运行周期和天气影响只能提供离散的观测资料,无法揭示作物生长和产量形成的机理,而这正是作物模型的优势。因此,运用数据同化技术结合作物模型和遥感数据,成为估测区域尺度作物产量的一种重要方法（本文简称同化方法）,该方法不仅发挥作物模型和遥感数据各自的优势,并且已成为定量遥感和精确农业研究重点和发展趋势[6-8]。
运用同化方法估测作物产量已经具备一定的研究基础。1999年deWit[9]利用遗传算法同化NOAA卫星数据与作物模型,证明在混合像元影响下同化NOAA卫星数据能提高作物模型模拟产量的结果;2005年Mo等[10]利用NOAA卫星数据反演时间序列叶面积指数（LAI）与SVAT作物模型同化估测华北平原冬小麦产量,估测产量与实测产量的均方根误差为1124 kg/hm2;2006年闫岩等[11]利用MODIS LAI产品和实测LAI同化CERES-Wheat模型改进预测产量结果,与实测产量比较,相对误差分别为-8.55%和1.82%;2008年Dente等[12]通过同化多源遥感反演LAI与CERES-Wheat作物模型模拟LAI,提高小麦区域产量估测精度结果显示,LAI反演精度越高,估测产量精度越高;2010年陈劲松等[13]通过同化环境卫星反演的LAI与WOFOST作物模型反演的LAI,提高了作物模型模拟区域产量的精度,均方根误差为35 g/m2;2011年Tripathy等[14]通过遥感观测NDVI确定播种日期,再通过同化遥感观测LAI与WOFOST作物模型模拟LAI估测冬小麦产量,相对误差小于6%。综上可知,以往研究在同化遥感数据和作物模型方面已经取得丰富成果,但其估产精度仍然偏低。为进一步提高作物同化估产的精度,主要从2个方面进行改进：（1）提高作物参数（如LAI）和模型参数（如出苗日期、灌溉日期等）精度及同化方法估测产量的精度,例如,刘翔舸等[15]利用同化方法提高遥感数据反演LAI精度;罗治勇[17]和何锦[18]等人分别利用SWAP模型模拟农田水分动态变化;Ayse等[19]通过同化蒸散发（ET）模拟灌溉;郭建茂等[20]利用遥感数据反演蒸散发替代作物模型模拟蒸散发,在此基础上同化LAI提高估测产量精度;黄健熙等[21]比较单独同化MODIS ET和单独MODIS LAI 2种方法的估产精度;（2）选取多种变量代替单一变量同化改进估测产量精度,该方法的依据是同化多个变量,可协同提高作物模型估测产量精度,其结果优劣取决于同化变量的选取,本文针对采用同化变量,蒸散发（Evapotranspiration,ET）（也称作陆地表层水汽通量）代表土壤、水面蒸发和植被蒸腾的总和,表征农作物水分胁迫的状况,是地表水量平衡和热量平衡的重要参量[17-18]。同化遥感观测ET到作物模型可优化作物模型,对土壤水分平衡的模拟[19],提高水分胁迫模式下的冬小麦生长过程模拟的精度[20]。
本研究旨在验证同化双变量优于同化单变量估产的理论;对比同化ET和LAI双变量和同化ET或LAI单变量对产量的优化程度;对比同化ET和同化LAI对模型模拟蒸散发和产量的优化程度。主要方法是：首先,标定WOFOST作物模型,设定水分胁迫模式,并优化模型灌溉参数,同时采用SCE-UA数据同化方法集成到作物模型,构建作物估产同化系统;然后,引入MODIS数据产品并运行系统,通过改变代价函数和同化变量形成4种同化方案（WOFOST模型模拟、单独同化MODISET与WOFOST模型模拟ET、单独同化MODIS LAI与WOFOST模型模拟LAI、同化MODISET和MODIS LAI与WOFOST模型模拟ET和LAI双变量）;最后,分析4种同化方案估测产量的结果。

2 研究数据与方法

（1）研究区为保定和衡水的冬小麦主产区,位于河北省中南部（图1）,地形以平原为主,耕地占总面积的60%以上,大部分地区适宜冬小麦生长。气候属温带季风性气候,年日照时数2400~3100 h,年均降水量300~800 mm,1月平均气温在3 ℃以下,7月平均气温18~27 ℃,四季分明。冬小麦作为该地区的主要粮食作物之一,播种期一般在10月上旬,返青期在2月下旬至3月上旬,3月下旬至4月中旬为拔节期,4月下旬至5月上旬为抽穗扬花期,5月中下旬为灌浆期,6月中旬至下旬为收获期。
Fig. 1 Location of the study area

图1 研究区域位置

（2）WOFOST作物生长过程模拟模型（简称WOFOST模型）由荷兰瓦赫宁农业大学和世界粮食研究中心共同研发,用于分析产量的年际变化,产量和土壤条件的关系,不同品种的差异、种植制度、气候变化对产量的影响,区域生产力的限制因素等[22];由于它对作物生长过程的模拟（图2）是通用的,因此,可通过改变作物参数来分析不同的作物[23];此外,WOFOST模型能根据需要分别模拟潜在、水分胁迫、养分胁迫3种生产水平下的作物生长过程。综合模型以上特点,通过搜集冬小麦的作物参数标定WOFOST模型,同时考虑研究区域（保定衡水地区）的冬小麦产量主要受干旱影响,开展了水分胁迫模式的作物生长过程模拟。
Fig. 2 The schematic diagram of WOFOST cropmodel

图2 WOFOST作物模型示意图

（3）气象数据由研究区周围6个国家气象局农业气象观测站（点）提供（http://cdc.cma.gov.cn/）,包含2008-2009年冬小麦生长发育期内模型所需的6个气象要素（最高气温/最低气温、总辐射、水汽压、风速、降水）;所有气象要素数据以1 d为步长;采用反距离权重插值法（IDW）,形成研究区内1 km分辨率大小的栅格格式气象数据。
作物数据根据获取难度的不同主要来源于3种途径：田间观测、经验值和WOFOST模型的缺省值（表1）。
Tab. 1 The table of WOFOST cropmodel parameters

表1 WOFOST作物模型参数表

参数 中文描述 取值 来源
TSUM1 出苗到开花积温 777.0 ℃·d-1 实测计算
TSUM2 开花到成熟积温 672.0 ℃·d-1
DLO 发育最适光长 14 h 文献
DLC 临界光长 8 h
CVL 叶同化物转化效率 0.740 kg·kg-1 实测计算
CVO 贮存器官同化物转化效率 0.830 kg·kg-1
CVR 根同化物转化效率 0.694 kg·kg-1
CVS 茎同化物转化效率 0.740 kg·kg-1
SLATB 比叶面积 DVS=0
0.00224 hm2·kg-1
DVS=2.0
0.00350 hm2·kg-1
SMW 凋萎土壤含水量 0.072 cm3·cm-3 实测
SMFCF 田间持水量 0.325 cm3·cm-3
SM0 饱和含水量 0.506 cm3·cm-3
TBASE 出苗最低温度下限 0℃ 文献
LAIEM 出苗时叶面积指数 0.13 hm2·hm-2 WOFOST
RGRLAI 叶面积指数最大增长速率 0.00817 hm2·hm-2·d-1
KDIFTB 漫射光的消散系数 0.6
EFFTB 单叶片同化CO2的光能利用率 0.47（kg·hm-2·h-1）/（J·m-2·s-1）
AMAXTB 最大CO2同化速率 45 kg·hm-2·h-1
Q10 温度升高10 ℃呼吸作用变化速率 2
RML 叶相对维持呼吸作用 0.03 kg·kg-1·d-1
RMO 贮存器官的维持呼吸作用 0.01 kg·kg-1·d-1
RMR 根的维持呼吸作用 0.015 kg·kg-1·d-1
RMS 茎的维持呼吸作用 0.015 kg·kg-1·d-1
蒸散发数据来自MODISET产品（即MOD16 A2）,MOD16 A2所采用的ET算法自发布产品前后时间段内已经得到广泛的验证,主要包括8 d合成产品（8 d）和月合成产品（Monthly）的ET、LE、ETP全球数据[24]。本研究选用2009年第1~152 d的MOD16 A2产品共19景。该产品经过简单的影 像镶嵌和投影变换后即可直接由影像DN值转化得到ET值（ET=DN×0.1）,本研究提取研究区的8 d合成ET值表示1 km空间分辨率上8 dET估测值的 总和。
叶面积指数采用2009年第1~152 d的MCD15 A3产品共38景。MCD15 A3是MODIS标准产品中的陆地4级产品,与MOD15 A2（Terra-MODIS）和MYD15 A2（Aqua-MODIS）产品不同,MCD15 A3数据产品同时使用Terra和Aqua卫星上的MODIS传感器数据进行反演,时间分辨率更高,每4 d合成一次,空间分辨率1 km,更有利于农作物长势和物候的监测。但是,由于云和大气的影响,MCD15 A3产品带有较大的噪声,由于LAI属于相对变化稳定的状态参数,本研究假设时间序列LAI曲线的局部最小值代表云雾引起的噪声,局部最大值代表真实值,采用Savitzky-Golay滤波（简称S-G滤波）算法对时间序列LAI曲线进行“去噪”处理。利用29个实测样点数据验证S-G滤波MODIS LAI时间序列曲线,发现该曲线的时序变化趋势与实际冬小麦物候相符（R2=0.62,RMSE=4 d）。但是,MODIS数据空间分辨率较低,像元内包含了除冬小麦以外的多种地物类型,造成LAI数值普遍偏低（< 2 m2·m-2）。鉴此,本文结合29个实测样点的冬小麦LAI数据 与对应的MODIS LAI数据构建回归模型,并利用该回归模型纠正整个研究区域的时间序列MODIS LAI产品,从而得到校正后的LAI（记为adjust-LAI）与实测LAI的绝对误差由3.64 m2·m-2减小至1.28 m2·m-2。图3所示为LAI的S-G滤波结果及修正结果（DOY指儒略日,一种连续纪日法）。
Fig. 3 The comparison chart of LAI before and after adjustment

图3 LAI调整前后对比图

3 冬小麦产量估测同化系统构建

3.1 系统平台建成过程

SCE-UA算法集合了控制随机搜索法和竞争进化概念,以及复合体的混合概念[25-27],适用于遥感信息与作物生长模型的同化过程[6]。本研究采用Intel Fortran Complier软件作为系统平台集成WOFOST作物模型和SCE-UA同化算法代码构建冬小麦产量估测同化系统。该系统技术流程如图4所示：首先,处理遥感数据产品,标定WOFOST作物模型并改进作物模型灌溉参数;然后,分别利用模型模拟同化变量和遥感估测同化变量构建代价函数,同化算法通过调整模型初始参数,使代价函数不断迭代更新达到最小值（代价函数达到最小的过程即同化的过程）;最后,输出优化后的产量、蒸散发和叶面积指数。图4中LAI代表叶面积指数,ET代表蒸散发,S-G代表Savitzky-Golay滤波器,RIRR代表灌溉量,TDWI代表作物初始干物质重,SPAN代表35 ℃时生命期（即35 ℃时作物生存天数）。下面将分别详述作物模型标定、灌溉参数优化、同化代价函数设计和模型初始参数选择4个系统构建的重要过程,并运行系统进行单产模拟。
Fig. 4 The flow chart of assimilating ET and LAI jointly

图4 同化ET和LAI双变量技术流程图

3.2 作物模型标定

WOFOST作物模型需要大量的模型参数,模型参数的变化范围和敏感性分析实验结果来自文献[28],依据这些参数的敏感程度和变化范围,本研究的模型标定方案如下：（1）敏感程度较高且变化范围较小的参数和不敏感的作物参数：通过查阅文献或者直接采用WOFOST模型的默认值来确定;（2）敏感程度高且变化范围大的作物参数,利用实测数据计算或查阅文献获得,或依据先验知识等确定参数可能的取值范围,通过“试错法”确定。具体参数见表1[28]（不包含本研究所选用的模型初始参数）。
标定作物模型并运行后得到ET和LAI 2个参数的模拟数据,处理成时间序列曲线,分别与遥感估测ET和LAI曲线（MODISET和MODIS LAI）、实测ET和LAI曲线对比（图5）,其中,ET曲线为8 d合成的数据,代表8 d ET的和,LAI曲线为4 d合成的数据,代表4 d的LAI,为清晰呈现ET趋势变化,利用样条插值法（SPLINE）对ET折线进行平滑处理。从图5中发现存在2个问题：（1）WOFOST模型模拟的ET在关键时间节点上的趋势变化不符合 实测ET的变化;（2）由于尺度效应,MODISET估测值普遍偏低,与实测数据相差较大（MODIS LAI也存在类似情况,但经过回归模型调整后的adjust-LAI数值接近于实测LAI）,与实测点上的ET趋势略有不同（图5（b）89 d ET）,但总体趋势一致, 因此,本文旨在比较同化双变量和单变量对估 产结果精度的影响（文中假设MODIS ET的趋势 可信）。
Fig. 5 The comparison chart of time series curves

图5 引入灌溉前后时序曲线对比图

3.3 灌溉参数优化

模型模拟ET趋势变化偏离实测ET这一现象的主要原因是“WOFOST模型默认灌溉量为零”,因此,本文将标定模型灌溉参数,然而由于实际灌溉量难以精确统计。为估测灌溉参数（包括灌溉日期、灌溉量2项参数）,本文综合研究区域地理、农时、物候等经验信息和气象、蒸散发等实测数据判断灌溉日期,并根据谢静等[29]提出的保定地区最优灌溉量标定模型灌溉量,最终采取灌溉标定方案（表2）。结果显示（图5,引入灌溉WOFOST-ET）,WOFOST模型经标定灌溉参数后提高了模型模拟ET的精度,尤其是冬小麦生长后期（DOY 113-145）的ET趋势变化得到明显改善,与实测情况基本相符。对于模型模拟LAI方面,标定灌溉参数增加了冬小麦生长后期（DOY 113-145）LAI的数值,但趋势变化没有得到改善。
Tab. 2 The defining scheme of the irrigation parameters in WOFOST

表2 标定WOFOST模型灌溉参数的方案

生育阶段 拔节-抽穗 灌浆-成熟
灌溉量（cm） 8.006 9.839
灌溉日期 4.6-4.13 5.23-5.30
儒略日DOY 97-104 137-144
遥感数据日期 第97天 第137天
WOFOST模型参数 RIRR1 RIRR2

3.4 同化代价函数设计

代价函数是农业数据同化系统的核心,用于评价观测资料与模型模拟的拟合程度,寻找使代价函数最小的控制变量,即寻找使拟合程度最好的控制变量。本研究根据MODISET和MODIS LAI的数据特点,分别采用不同的代价函数同化ET和LAI。
单点尺度方面,MODISET数值普遍低于实测值,但趋势变化信息相对比较准确。如果用最小二乘法直接同化MODISET和模型模拟ET则会引起较大误差,导致产量结果不切实际地偏低。为避免该问题,本研究利用MODISET和模型模拟ET的时序变化趋势差异最小来构建代价函数,具体思路是利用一阶差商代表ET时序曲线的变化趋势,以遥感观测ET一阶差商与模型模拟ET一阶差商的一致性表述趋势的相似性。具体的构建步骤为：
（1）计算每个时间节点上ET的一阶差商。如式（1）所示, D i 表示第i时间节点上的一阶差商值; ε 代表ET数据的最小变化量,等于MODISET最小变化量（0.1 mm）。其中,数值-1,0,1分别代表相应时间节点上的ET变化趋势为下降,不变和升高。
（2）判断遥感观测ET和模型模拟ET时序变化趋势的异同,如果相异,代价函数加一,否则不变。如式（2）所示,FET代表同化ET的代价函数值; D obs i 和 D sim i 分别表示第 i 时间节点上遥感观测ET和模型模拟ET的一阶差商值,用于计算代价函数的ET步长为8 d,冬小麦全生育期内共有M个MODISET数据。FET的取值范围是0-M。
D i = - 1 , if ( E T i + 1 - E T i - 1 ) < - 2 ε 0 , if abs ( E T i + 1 - E T i - 1 ) ≤ 2 ε 1 , if ( E T i + 1 - E T i - 1 ) > 2 ε （1）
F ET = F ET + 1 , if ( D obs i ≠ D sim i ) F ET , if ( D obs i = D sim i ) i = 1,2 , … M （2）
F LAI = 1 2 ∑ k = 1 2 X k - X k 0 T B - 1 X k - X k 0 + 1 2 ∑ j = 1 N Y sim j - Y obs j T Q - 1 ( Y sim j - Y obs j ) （3）
MODIS LAI经过实测数据的回归模型校正后精度较高,适合直接用于同化模型模拟LAI,本研究采用四维变分方法构建代价函数。四维变分代价函数主要包括2部分：（1）描述模型初始参数的优化值与初始值的距离 X K - X K 0 ;（2）描述同化变量的遥感估测值和模拟值之间的距离 ( Y sim j - Y obs j ) 。代价函数的具体数学表述形式如式(3)所示, F LAI 代表同化LAI的代价函数值; X K - X K 0 T 代表矩阵 X K - X K 0 的转置;B矩阵是指代模型初始参数的误差协方差矩阵; B - 1 X K - X K 0 是B矩阵的逆; Q ( Y sim j - Y obs j ) 是同化变量LAI遥感估测值的误差协方差矩阵,由实测值和MODIS LAI共同计算得到;N表示冬小麦生育期内的MODIS LAI（步长为4 d）个数。四维变分同化方法的精度主要取决于模型误差协方差矩阵B和观测量误差协方差矩阵Q,确定B矩阵需要考虑模型参数的不确定性和研究区域的经验信息,本研究将B定义为模型初始参数的误差矩阵,优化的初始值见表3,误差矩阵B依据公式（3）经验运行的统计信息标定为 7.8 0 0 0.7 ;Q矩阵表征观测误差,本研究选取研究区29个实测样本点,分别计算7个重要物候阶段中调整后的LAI（adjust-LAI）与实测LAI值的绝对误差的均值标定Q矩阵。计算结果为：返青阶段LAI的平均绝对误差等于0.43 m2·m-2;起身期等于0.67 m2·m-2;拔节期等于0.86 m2·m-2;孕穗期等于0.61 m2·m-2;抽穗期等于0.77 m2·m-2;开花期等于0.93 m2·m-2;成熟期等于0.47 m2·m-2。通过对计算机计算的FLAI值的统计,其取值主要分布于3000-5000之间。
同化双变量LAI和ET的代价函数既要包含以趋势法计算的ET代价,也要包括以四维变分法计算的LAI代价,但是,二者数量级相差较大,不能直接将二者代数和作为总代价,否则,将会造成数量级小的单变量代价比重过低。因此,需要对ET和LAI的代价值进行式（4）所示处理,其中,average ( F LAI ) 和average ( F ET ) 分别代表ET和LAI代价函数值的平均值,通过对所有代价值的分析和经验得到,再以式（5）的形式计算总代价。
F LAI = F LAI average ( F LAI ) F ET = F ET averag e ( F ET ) （4）
F = F ET + F LAI （5）

3.5 模型初始参数选择

蒸散发ET是反映农田水分平衡的重要参数,对灌溉十分敏感,适合作为判别水分胁迫是否发生的指标,因此,本研究同化ET选择灌溉参数作为相应的待优化初始参数。同时根据WOFOST模型敏感性分析结果[28]可知叶面积指数LAI对出苗日期（IDEM）、作物初始干物质重（TDWI）和35 ℃时的生命期（SPAN）敏感度最高,但IDEM与TDWI相关性很强,不适合共同作为待优化初始参数,本研究仅选择TDWI和SPAN,作为同化LAI的待优化初始参数。同化双变量ET和LAI选择RIRR,SPAN和TDWI共同作为待优化初始参数,其中,RIRR分解成RIRR1和RIRR2 2个参数,指代不同日期的灌溉量。同化前TDWI和SPAN初始值参考作物模型默认值,RIRR1和RIRR2初始值参考优化灌溉量（表2）,待优化初始参数初始值如表3所示。
Tab. 3 The initial values of optimization parameters

表3 模型初始参数值

参数名称 TDWI SPAN RIRR1 RIRR2
初始值 210 27 8.006 9.839
单位 kg/hm2 d cm cm

3.6 单产模拟比较

本研究中共设计4组同化方案（表4）对比估产精度：第1组是对照组,不进行同化处理,只记录WOFOST模型直接模拟的结果;第2组同化蒸散发ET,优化2个灌溉日期的灌溉量（RIRR1和RIRR2）,基于趋势法构建代价函数;第3组同化的变量为叶面积指数LAI,优化作物初始干物质重（TDWI）和35 ℃的作物生命期（SPAN）,代价函数构建方法是四维变分法;第4组同化双变量ET和LAI,优化4个模型参数（RIRR1, RIRR2, TDWI, SPAN）,总代价函数是ET代价和LAI代价处理后的代价和。
Tab. 4 Data assimilation schemes

表4 数据同化方案

编号 同化变量 模型初始参数 代价函数
第1组 无 无 无
第2组 ET RIRR1&RIRR2 趋势法
第3组 LAI TDWI&SPAN 四维变分法
第4组 ET&LAI TDWI&SPAN&RIRR1&RIRR2 趋势&四维变分法
系统构建完成后,运行并输出4种同化方案下单点和区域尺度的作物产量、蒸散发和叶面积指数,对比分析估产结果,进一步探索提高估产精度的同化方法。

4 系统估产结果与分析

4.1 模型初始参数

模型初始参数的优化结果如表5所示,其中,第1组是模型模拟结果,4个模型参数未经优化,与初始值相等;第2组只优化灌溉参数RIRR1和RIRR2,2个参数经优化后数值升高,其余参数不变;第3组只优化TDWI和SPAN,优化后的TDWI和SPAN数值变大;第4组优化4个参数,相比于初始值TDWI升高幅度最大,2个灌溉量升高幅度较小,SPAN优化后值变小。
Tab. 5 The final values of optimization parameters

表5 模型初始参数终值表

编号 TDWI（kg/hm2） SPAN（d） RIRR1（mm） RIRR2（mm）
第1组 210 27 8.006 9.839
第2组 210 27 8.150 11.120
第3组 300 33 8.006 9.839
第4组 299.85 22 10.130 9.950

4.2 同化变量分析

图6（a）表示同化前后ET时相曲线的变化,从图中可看出：第1组至第4组ET曲线比MODISET更接近实测ET曲线;第1组ET曲线在DOY73-89,以及DOY97-1052个阶段波形凸起,与实测ET差距较大;第2组相比其他组ET曲线模拟效果得到提高,在DOY97-105阶段波形更平缓,但在DOY73-105阶段的局部趋势变化和数值大小仍然偏离实测ET。提高同化方案主要分为3个方面：（1）选取更高精度的ET遥感产品,或更可靠的ET观测产品;（2）改进同化ET的代价函数,使用其他同化趋势的代价函数提高同化趋势的结果;（3）改进作物模型本身模拟ET的精度,但此方法十分复杂,使用前2种方法更佳。
图6（b）表示LAI实测点同化前后LAI曲线的变化,从图中可看出：第1组曲线显示模型模拟LAI值域范围与实际LAI相差最大;第2组与第1组相比,仅在DOY125-145阶段提高了LAI数值;第3组LAI数值提高最为显著,但在DOY 137之后的曲线趋势变化与实测不符;第4组的LAI同化结果相比于第3组数值略有下降,但在DOY 137之后的曲线趋势更加符合实测趋势。总体上,第3-4组LAI曲线优于其他两组,说明同化LAI使模型模拟LAI的精度得以提高;同时,由于校正后的遥感数据精度很高,因此,同化结果精度相对ET而言精度更高,间接证明ET遥感产品精度是制约模型模拟结果的重要因素。
Fig. 6 The comparison chart of ET and LAI curves before and after data assimilation

图6 数据同化前后对比图

4.3 区域产量对比分析

利用ARCMAP软件的ZONAL工具,将产量结果从像元尺度统计到区域尺度（县级行政单位）,并根据保定市和衡水市2009年冬小麦官方统计产量分析各同化方案模拟的区域产量精度,结果如表6所示,由于本研究在水分胁迫模式下进行同化,预测产量普遍低于实际产量。图7是各组预测产量结果与统计产量的散点图,各组散点的收敛程度和产量数值关系：第1组<第2组<第3组<第4组。
Fig. 7 Scatterplotsofwinterwheatyieldwithm2ssimilation

图7 冬小麦产量散点图

Tab. 6 The statistical table of winter wheat yield (kg/hm2)

表6 冬小麦产量统计表（kg/hm2）

— 均值（μ） 最大值（max） 最小值（min） 决定系数（R2） 均方根误差（RMSE）
统计产量 6022 6428 5592 — —
第1组 4844 5343 4392 0.006 1210
第2组 4918 5392 4498 0.013 1134
第3组 5131 5669 4478 0.408 925
第4组 5355 5863 4769 0.432 721
利用像元级同化产量结果和县级统计产量绘制空间分布图（图8）,显示了44组同化产量和统计产量的空间分布。由图8可得出,官方统计冬小麦产量在整个研究区域的西部和南部偏高,但在研究区域的东北部偏低。按照这一空间分布规律分析4组产量结果发现第4组同化产量分布图与统计产量的空间分布规律吻合程度最高,而且产量空间变异性较大,产量数值最接近统计产量,因此,估产精度最高。而第1组和第2组的西部产量结果偏低,产量的空间差异不明显。第3组的产量结果整体值域偏低,但产量空间差异较为明显,精度略低于第4组的预报产量。
Fig. 8 Spatial distribution of winter wheat yield

图8 冬小麦产量分布图

综合以上：由于MODIS ET精度有限,同化ET方法对模型模拟ET优化程度较小,导致产量结果精度提高不大;由于MODIS LAI经调整后精度较高,同化LAI方法对模型模拟LAI和产量提高较大;即使同化结果受限于MODIS ET精度,但是,同化ET和LAI双变量的方法仍然明显优于仅同化LAI方法。这是由于叶面积和植被蒸腾之间存在制约关系（如图2所示,植被最大蒸腾速率受叶面积直接控制,相对蒸腾速率直接影响光合产物积累,间接影响叶片面积增长,因此,植被蒸腾和叶面积之间不是完全独立的,2变量会随着彼此改变而改变）,即使ET优化程度很小,但在LAI的协同作用下产量结果仍然得到“双重”优化,使精度优于同化单变量的方法。

5 结论

本研究对比分析了单独同化ET、单独同化LAI、同化ET和LAI双变量3种方式的冬小麦产量预报精度结果表明,在水分胁迫模式下,同化ET和LAI双变量能取得比同化ET或LAI单变量更高的精度,说明双变量数据同化有更大的优势和应用潜力,均值为5355 kg/hm2,决定系数R2=0.432,均方根误差RMSE=721 kg/hm2;单独同化LAI也取得较高的估测精度,均值为5131 kg/hm2,R2=0.408,RMSE=925 kg/hm2;单独同化MODISET能显著改善WOFOST模型中的土壤水分敏感的参数,但是对冬小麦产量的提升作用不明显,同化ET的产量精度最低,均值为4918 kg/hm2,R2=0.013,RMSE=1134 kg/hm2,略优于模型模拟产量结果（均值为4844 kg/hm2,R2=0.006,RMSE=1210 kg/hm2）。
同化ET和LAI双变量的冬小麦估产精度仍然偏低的原因：（1）由于MOD16ET产品是主要面向全球尺度生产的,并且ET算法中采用了分辨率较粗的MERRA GMAO气象再分析数据（原始空间分辨率为 0.5 ° × 2 / 3 ° ）,导致MOD16ET产品的区域尺度精度不高,影响同化系统精度;（2）遥感反演参数与作物模型模拟参数之间的时空尺度不匹配和遥感反演参数精度限制了同化模型的精度;（3）代价函数的构建影响了同化结果精度。因此,在今后的研究中,农业同化系统应集中在提高遥感反演农田参数的精度、遥感与作物模型之间的尺度转换模型,以及改进代价函数方面。

The authors have declared that no competing interests exist.

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
Tennakoon S B, Murty V V N, Eiumnoh A. Estimation of cropped area and grain yield of rice using remote sensing data[J]. International Journal of Remote Sensing, 1992,13(3):427-439.

[2]
Duveiller G, Baret F, Defourny P.Monitoring crop growth inter-annual variability from MODIS time series: Performance comparison between crop specific green area index and current global leaf area index products[C]. 2011 6th International Workshop on the Analysis of Multi-temporal Remote Sensing Images (Multi-Temp), IEEE, 2011:21-24.

[3]
Duveiller G, Baret F, Defourny P.Crop specific green area index retrieval from MODIS data at regional scale by controlling pixel-target adequacy[J]. Remote Sensing of Environment, 2011,115(10):2686-2701.

[4]
张艳楠. 遥感估产技术研究现状与展望[J].现代农业科技,2010(7):52-56.

[5]
陶伟国,徐斌,杨秀春.草原产草量遥感估算方法发展趋势及影响因素[J].草业学报,2007,16(2):1-8.

[6]
黄健熙,武思杰,刘兴权,等.基于遥感信息与作物模型集合卡尔曼滤波同化的区域冬小麦产量预测[J].农业工程学报,2012,28(4):142-148.

[7]
马玉平,王石立,张黎,等.基于遥感信息的华北冬小麦区域生长模型及模拟研究[J].气象学报,2005,63(2):204-215.

[8]
黄健熙,李昕璐,刘帝佑,等.顺序同化三种时空分辨率LAI的冬小麦估产对比研究[J].农业机械学报,2015,46(1):241-248.

[9]
De Wit A J W. Application of a genetic algorithm for crop model steering using NOAA-AVHRR data[C]. Proceedings of SPIE 3868, Remote Sensing for Earth Science, Ocean, and Sea Ice Applications, 1999:167-181.

[10]
Mo X, Liu S, Lin Z, et al.Prediction of crop yield, water consumption and water use efficiency with a SVAT——crop growth model using remotely sensed data on the North China Plain[J]. Ecological Modelling, 2005,183(2):301-322.

[11]
闫岩,柳钦火,刘强,等.基于遥感数据与作物生长模型同化的冬小麦长势监测与估产方法研究[J].遥感学报,2006,10(5):805-811.

[12]
Dente L, Satalino G, Mattia F, et al.Assimilation of leaf area index derived from ASAR and MERIS data into CERES——Wheat model to map wheat yield[J]. Remote sensing of Environment, 2008,112(4):1395-1407.

[13]
陈劲松,黄健熙,林珲,等.基于遥感信息和作物模型同化的水稻估产方法研究[J].中国科学:信息科学,2010,40(增刊):173-183.

[14]
Tripathy R, Chaudhari K N, Mukherjee J, et al.Forecasting wheat yield in Punjab state of India by combining crop simulation model WOFOST and remotely sensed inputs[J]. Remote Sensing Letters, 2013,4(1):19-28.

[15]
刘翔舸,刘春红,王鹏新,等.基于卡尔曼滤波的小麦叶面积指数同化方法[J].农业工程学报,2010,26(1):176-181.

[16]
邢雅娟,刘东升,王鹏新.遥感信息与作物生长模型的耦合应用研究进展[J].地球科学进展,2009,24(4):445-449.

[17]
罗治勇. 农田水分动态研究[D].西安:长安大学,2008.

[18]
何锦. 基于SWAP模型的农田动态水分研究[D].西安:长安大学,2006.

[19]
Irmak A, Kamble B.Evapotranspiration data assimilation with genetic algorithms and SWAP model for on-demand irrigation[J]. Irrigation science, 2009,28(1):101-112.

[20]
郭建茂. 基于遥感与作物生长模型的冬小麦生长模拟研究[D].南京:南京信息工程大学,2007.

[21]
黄健熙,马鸿元,田丽燕,等.基于时间序列LAI和ET同化的冬小麦遥感估产方法比较[J].农业工程学报,2015,31(4):197-203.

[22]
谢文霞,王光火,张奇春.WOFOST模型的发展及应用[J].土壤通报,2006,37(1):154-158.

[23]
邬定荣,欧阳竹,赵小敏.作物生长模型WOFOST在华北平原的适用性研究[J].植物生态学报,2003,27(5):594-602.

[24]
Mu Q, Heinsch F A, Zhao M, et al.Development of a global evapotranspiration algorithm based on MODIS and global meteorology data[J]. Remote Sensing of Environment, 2007,111:519-536.

[25]
Duan Q Y, Gupta V K, Sorooshian S.Shuffled complex evolution approach for effective and efficient global minimization[J]. Journal of Optimization Theory and Application, 1993,76(3):501-521.

[26]
Duan Q Y, Sorooshian S, Gupta V K.Effective and efficient global optimization for conceptual rainfall—runoff models[J]. Water Resource Research, 1992,28(4):1015-1031.

[27]
Duan Q Y, Sorooshian S, Gupta V K.Optimal use of the SCE_UA global optimization method for calibrating watershed models[J]. Journal of Hydrology,1994,158:265-284.

[28]
马冠南. 基于MODIS数据与WOFOST作物模型同化的区域冬小麦产量估测方法研究[D].北京:中国农业大学,2012.

[29]
谢静. 保定地区冬小麦水分生产函数及节水灌溉制度研究[D].石家庄:河北农业大学,2011.

Options
文章导航

/

〈 〉