基于遥感信息和WOFOST模型参数同化的冬小麦单产估算方法研究

陈艳玲, 顾晓鹤, 宫阿都, 胡圣武

(1.北京师范大学环境遥感与数字城市北京市重点实验室，北京 100875；2.北京师范大学环境演变与自然灾害教育部重点实验室，北京 100875；3.北京师范大学地理科学学部，北京 100875; 4.国家农业信息化工程技术研究中心，北京100097; 5.河南理工大学测绘与国土信息工程学院，河南焦作 454000)

中国是农业大国，冬小麦是我国主要的粮食作物之一，其播种面积占我国粮食总播种面积的1/5[1]。冬小麦长势监测的实时性和产量预报的准确性是保障国家粮食安全、经济发展的基础。因此，农作物长势监测和估产研究已成为国内外学者的研究热点。迄今为止，卫星遥感和作物生长模拟两项新技术发展最为迅速，潜在应用价值巨大，同时存在一些有待深入解决和研究的问题。

卫星遥感能及时、准确和客观地反映农作物在时间和空间上的分布信息，在大尺度作物长势动态监测和估产方面优势明显。遥感估产研究已从试验阶段逐步进入到实际应用阶段。但是，地学遥感信息表征的只是地表或作物冠层的瞬时物理状况，尚不能真正揭示作物生长发育和产量形成的内在机理过程、个体生长发育状况及其环境气象条件对其的影响[2]。

近年来，作物生长模型模拟技术发展迅速，它综合考虑外部条件、作物遗传特性等因素的影响，通过给定一系列作物参数和环境参数，定量模拟作物的生长、发育和产量形成过程。迄今为止，以荷兰和美国为代表所建立的作物生长模型多达百种。但是，模型中一些参数获取困难，在外界因素的影响下，模拟作物生长过程的一些重要变量地域变异大，难以保证区域尺度上的模拟精度[3]。

综上所述，遥感监测作物宏观状况，表征环境因素对作物的影响，空间连续性强，而作物模型模拟作物生长过程，揭示原因和本质，具有时间连续性的优势，两者各有优点和不足。近年来，将遥感与作物模型两者优势相结合的研究已成为国内外的研究热点[4-9]，通过遥感数据的辅助不仅可以实时准确地获取区域尺度的模型初始输入参数，同时有利于提高作物生长模型的模拟精度。目前，研究方法主要分为驱动法和同化法。驱动法是初期研究的主要方法，使用简单，但是对反演精度要求高。同化法受遥感数据观测时间、空间以及精度的影响较少，应用越来越广泛。同化方法又可分为两种，一种是直接利用遥感数据作为模型的初始化参数，另一种是以植被指数或者反射率等作为同化变量。如Dente[6]通过调整小麦播种时间、萎蔫点和田间持水量三个参数，使利用ENVISAT ASAR和MERIS数据反演LAI和模拟LAI差异最小，实现雷达数据和作物模型的同化，并使同化后的CERES-Wheat模型预测小麦产量的能力显著提高。Ma等[10]利用SCE-UA算法将MODIS-LAI数据产品同化进WOFOST模型中，通过MODIS-LAI值和模型模拟LAI值之间的差异，调整WOFOST模型中的关键参数，取得了较好的产量预测结果。Anderson等[11]分析了2003-2013年巴西主要作物不同生育时期的蒸发胁迫指数(ESI)与产量的关系，发现开花期和灌浆期的ESI与产量相关性最高。

WOFOST模型是由荷兰瓦赫宁根农业大学和世界粮食研究中心共同开发的具有较大影响力的作物生长模型之一，普适性强，是模拟特定气候和土壤条件下一年生作物生长的动态解释性模型，在农作物生长评估、产量预测、土地定量评价等领域得到广泛应用[12]。本研究利用查找表优化算法实现WOFOST模型和遥感反演LAI信息同化，校正模型敏感性参数，以期提高WOFOST模型区域尺度小麦单产模拟精度。

1 材料与方法

1.1 研究区概况

本试验选取河北省藁城市作为研究区，该区地处河北省西南部，介于东经114°38′45″～114°58′47″ 、北纬37°51′00″～38°18′44″之间。全境属于暖温带半湿润大陆性季风气候，四季分明，冬冷夏热，平均气温12.5 ℃。年平均降水量494 mm，7-8月份降水量最多，约占全年的56.2%。年日照时数2 711.4 h，日照率高达61.2%，无霜期190 d。褐土类和潮土类是藁城市面积最大、最主要的土壤类型。藁城市地处华北平原，地势平坦，气候因素分布均匀，是冬小麦的主产区，因此选取该区作为研究区。

图1 研究区概况及样点分布

1.2 数据采集

1.2.1 气象数据采集

本研究所需要的2009-2010年的气象数据由中国气象科学数据共享服务网下载得到，主要包括每日的最高温、最低温、降水量、风速、日照时数、水汽压等气象要素。WOFOST模型的输入参数是总辐射量，因此利用FAO公式(1)将日照时数转化为总辐射量(Rs)。

Rs=(as+bsn/N)Ra

(1)

式中，as和bs是由联合国粮食及农业组织FAO为不同地区提供的与大气质量有关的参数，由于本研究区位于温带地区，因此分别采用0.18和0.55；n为日照时长，可由石家庄气象站获取；N为可日照时数[公式(2)]；Ra为大气上届入射辐射[公式(3)]。

N=24ωs/π

(2)

Ra=37.6dr(ωssinφsinδ+cosφcosδsinωs)

(3)

式中，dr为日地距离系数；ωs为日面中心的时角；δ为地球赤道平面与太阳和地球中心的连线之间的夹角(太阳赤纬)；φ为测点维度，分别由公式(4)计算得到:

dr=1+0.033cos(0.017 2J)

δ=0.420 9sin(0.017 2J-1.39)

(4)

ωs=arcos(-tanφtanδ)

式中，J为日序。

最后将处理好的气象要素按照WOFOST作物生长模型需要的格式建立对应的数据库文件。

1.2.2 农学数据采集

2009-2010年在河北省藁城市进行冬小麦主要生育期的田间试验，选取22个样本点，样本点的分布如图1所示，田间获取叶面积指数(LAI)和产量数据。

设定取样面积为2 500 cm2(50 cm×50 cm)，采用比叶重法测定LAI，计算公式为：

LAI=(w1+w2+w3)/(w1×A×10 000)×S

式中，w1为标叶干重(g)；w2为15株的余叶重；S表示标叶面积(cm2)；A为取样面积(cm2)；w3为剩余叶片干重。

在成熟期，收取每个样本小区1 m2面积冬小麦，晒干后脱粒测产以获取每个样本点的产量数据。

1.2.3 遥感数据采集

2008年6月我国为环境监测与灾害损失监测发射了两颗环境减灾小卫星HJ-1A和HJ-1B。相比于高时间分辨率的MODIS数据和高空间分辨率的Landsat TM数据，CCD数据兼具高时间分辨率(2 d)和高空间分辨率(30 m)特点，可满足农情监测多时相、大范围的需求，有利于提高大尺度地表状况动态监测的准确性。所以，本研究选取环境减灾小卫星CCD数据作为研究所需遥感数据，分别选取时相为2010年5月20日和5月29日的两景晴空过境影像，裁剪至研究区范围，并进行辐射定标、大气校正、几何配准等处理。

1.3 研究方法

1.3.1 采用植被指数法反演叶面积指数(LAI)

在反映作物长势的指标中，LAI是最常用的参数之一，是衡量植被生长和发育状态的重要因子。由于在可见光和近红外波段的植被光谱反射率与LAI相关性高，因而利用LAI与作物植被指数间的关系反演LAI的方法已在区域和大尺度遥感反演中得到广泛应用。本研究基于宽波段反射率，选取了7个在反演冬小麦LAI效果较好的植被指数[13-18]推算LAI，再分别与实测数据进行回归统计，筛选反演指示性、稳定性和预测能力都较好的植被指数。

本研究从2010年河北省藁城市冬小麦开花期(5月20日)和灌浆期(5月29日)植被指数(表1)的22个地面实测点数据中，随机选择15个实测LAI与对应点的植被指数构建统计模型，用剩余的7个数据进行验证。以决定系数(r2)和均方根误差(RMSE)来评价模型的准确性，筛选出与LAI拟合度最好的模型。

表1 冬小麦 LAI 反演各类植被指数公式Table 1 Formulas of the vegetation index for retrieving LAI of winter wheat

1.3.2 WOFOST模型参数敏感性分析

模型的“区域化”和“本地化”是遥感信息和作物模型耦合的首要问题，最关键的环节是模型参数的优化和校正，模型参数的敏感性分析是解决该问题的有效和常用方法，其中扩展傅立叶振幅灵敏度检验法EFAST(extended fourier amplitude sensitivity test)具有较大的应用潜力[19-20]。为解决WOFOST模型在我国的适用性问题，马玉平等选取泰安、郑州和固城三个站点，对该模型参数进行调整，使其更好地模拟华北地区冬小麦的生长过程[21]；郭建茂等对其也进行了适当改进，以适应山东禹城冬小麦生长模拟[12]。本研究基于前人的研究成果利用敏感性分析专业软件Simlab 2.2-EFAST模块对WOFOST模型中冬小麦26个作物参数进行全局敏感性分析。首先，采用均匀分布的方法对26个参数进行插值，得到5 000×26个样本数据。将采样的参数组输入模型，获得对应的5 000个产量数据。最后，利用EFAST方法对各个参数进行敏感性分析，根据敏感性指数的大小筛选出对产量影响较大的敏感因子。

1.3.3 模型同化

查找表优化算法(look-up table，LUT)是一种简单易实现的算法，在模型反演核函数优化等方面得到了广泛应用。因此，本研究选取该算法进行模型同化。叶面积指数LAI的准确性直接影响作物长势好坏和对产量高低的评价，因此选取LAI作为优化比较对象，优化目标函数定义为：

表2 WOFOST模型参数及波动范围Table 2 Range of input variables for WOFOST model

式中，LAIsi和LAIri分别为WOFOST模型模拟和遥感反演LAI值；N为外部同化LAI数据的个数。

本研究将目标函数J作为一个判别条件，当J值最小，即LAIsi和LAIri间差异最小时，确定待优化参数的取值为“最佳”，完成同化。

本研究总的技术流程如图2所示。

2 结果与分析

2.1 遥感反演LAI结果

将小麦开花期和灌浆期的8个植被指数与LAI进行线性拟合，结果(图3)表明，两个时期EVI的拟合效果均最好。两个时期拟合方程的决定系数r2分别为0.913和0.798，RMSE分别为0.410和0.470；利用独立数据进行检验，r2分别为0.858和0.861，RMSE分别为0.531和0.428(图4)，说明利用EVI的反演效果较好。

2.2 WOFOST模型参数敏感性分析结果

WOFOST模型中冬小麦26个作物参数的一阶敏感性指数和总敏感性指数表现基本一致(图4)。在选取的26个作物参数中， TSUM1、SPAN、SLATB1、SLATB2、TSUM1和TSUM2等6个参数比较敏感，一阶敏感性指数和总敏感性指数均超过0.10，其余参数敏感指数均不足0.10。因此，这6个参数是本研究中WOFOST模型作物参数“本地化”时需要优先考虑调整和优化，剩余参数可以认为是对产量形成敏感性小或者不敏感，由实测计算、查阅文献获得或者直接利用WOFOST模型默认值进行固定。

图2 总体流程图

图3 开花期(a)和灌浆期(b)的EVI与LAI的建模结果

2.3 模型同化结果

将筛选出来的6个敏感性因子插值生成5 000×6个样本数据组，分别代入模型生成5 000个相应的LAI模拟值。根据查找表优化算法，从WOFOST模型模拟的LAI值中筛选出与遥感反演的LAI值差异最小的数据组。对LAI的模型模拟值与遥感反演值进行拟合(图6)，开花期和灌浆期的r2分别达到了0.912和0.839，说明模型模拟的LAI与遥感反演的LAI差异较小。再将同化后的模型模拟LAI、遥感反演LAI与实测LAI拟合，同化后的模型两个时期的模拟精度均高于遥感反演精度(图6和图7)，说明加入遥感信息后WOFOST模型模拟的LAI更能反映冬小麦生长的真实状况。

根据上述确定的敏感性参数的最佳取值，得到WOFOST模型的模拟产量数据。将同化前后的模拟产量、真实产量进行对比，模型同化前对产量的模拟精度较低(图8a)，同化后模拟精度明显提高，r2由原来的0.707提高至0.914，RMSE由426.74 kg·hm-2降低至253.67 kg·hm-2，表明模型同化效果较好，也说明利用查找表优化算法进行遥感反演LAI与WOFOST作物生长模型同化对作物单产模拟的可行性。

图4 开花期(a)和灌浆期(b)的反演结果

图5 WOFOST模型参数一阶和全局敏感性参数

图6 开花期(a)和灌浆期(b)模型模拟LAI和遥感反演LAI拟合精度

图7 开花期(a)和灌浆期(b)模型模拟LAI和实测LAI拟合精度

图8 开花期(a)和灌浆期(b)遥感反演LAI和实测LAI拟合精度

图9 同化前(a)和同化后(b)模型模拟产量和实测产量的相关性

3 讨 论

遥感信息与作物模型同化是近年来进行区域作物产量估算的重要方法之一。但多数研究仅局限于田间尺度，在区域或国家尺度的研究较少。本研究以 LAI为变量，借助查找表同化算法将其与WOFOST作物模型进行耦合，通过调整敏感性参数，实现WOFOST产量预测。在所选取的7个植被指数中，EVI是预测LAI效果最好的指数。NDVI是最常用指示植被覆盖度及植被生长状况的植被指数之一。研究表明，大气会降低红光和近红外反射信号的对比度，散射、上行路径辐射中大气的贡献导致红光信号的增强，而近红外信号由于散射和水汽吸收等大气衰减作用而减弱，最终导致NDVI信号的下降和地表植被量被低估。而EVI是一个优化的植被指数，通过加入蓝色波段以增强植被信号，矫正土壤背景和气溶胶散射的影响，对冠层结构变化冠层类型、植被相和冠层结构更加敏感，提高对高生物量区域的敏感性。

在本研究中，所选的冬小麦26个作物参数对产量的全局敏感性指数和一阶敏感性指数分析结果的趋势基本一致。全局敏感性分析综合考虑参数自身及参数间的相互作用对模拟结果的影响，得到的是各参数变化对模拟结果影响的总贡献，相较于局部敏感性分析而言，筛选和提取的敏感参数更有利于模型的改进和优化。在选取的26个作物参数中6个参数比较敏感，敏感性指数均超过0.10，这与马玉平[21]、郭建茂等[12]的研究结果高度一致。其中，DVS=0.5时的比叶面积SLATB2是对冬小麦产量形成影响最大的参数，可解释作物产量变化方差的49.0%；另外SLATB1的敏感性指数也达到了0.237。比叶面积SLATB是叶片单面面积与其干重之比，一般与受光照强弱程度有关，反映作物对不同环境的适应特征。已有研究表明，SLATB与叶片的净光合速率呈正相关，SLATB1和SLATB2的敏感性指数较大，可以变相表征出光合作用较强，有利于LAI的增长和干物质的形成，这就解释了比叶面积对产量形成最为敏感的原因。其次是出苗到开花期的有效积温TSUM1的敏感指数，其值为0.271，而TSUM2的改变对产量形成的影响不显著。这是因为WOFOST模型是利用实际累计温度与TSUM1的比值表征冬小麦的生长发育进程，与LAI的动态变化和干物质的分配等密切相关。另外，叶片衰老指数SPAN的敏感指数也达到了0.191，它是计算叶龄的重要参数，叶龄是表征作物生长发育进程的外部形态诊断指标，是影响作物光合作用的内部因素之一。20 ℃条件下的单叶光能初始利用效率EFFTB3和30 ℃条件下对最大CO2同化率的影响系数TMPF4也比较敏感，敏感性指数分别达到了0.126和0.149。这两个是影响光合速率的重要参数，光合速率的大小决定生物量的形成，温度影响叶片中酶的活性，进而影响光合速率的大小，因此这两个参数对产量的形成较为敏感。

气象条件、土壤状况、田间管理等外部因素以及小麦品种遗传参数的变化会引起相应的产量变化。通过敏感性分析筛选出这些因素中对产量影响较大的参数。光照条件、温度以及光合效率是影响生物量形成的主要因素。上述研究结果表明，6个敏感性因子中，比叶面积SLATB1和SLATB2表征光照强弱，TSUM1反映温度的影响，EFFTB3和TMPF4代表光合速率大小。叶片衰老指数SPAN表征作物的生长发育进程，影响叶片光合能力。除敏感性参数的影响之外，WOFOST模型其他参数的取值大小对估产结果也有影响。本研究只考虑了作物参数的影响，未考虑土壤状况的影响，计划在下一步工作中完善。

本研究将遥感监测信息同化到模型中，通过优化算法调整敏感性参数，使其更接近遥感监测的“客观”状态，WOFOST模型模拟的LAI更能反映冬小麦生长的“真实”状况，使得模拟产量更接近“真实值”，模型模拟河北省藁城市的小麦产量与实测产量的r2高达0.914。LAI作为中间参数，同化后的拟合精度高于遥感监测精度，说明同化目的达到。为同化方法应用于区域尺度作物单产估测积累了理论基础。冬小麦的主要生育时期包括返青期、拔节期、孕穗期、开花期和灌浆期等，本研究只考虑了开花期和灌浆期的影响，值得进一步的研究。

4 结 论

本研究通过对WOFOST模型参数进行适应性调整，建立了适合河北省藁城市冬小麦生长模型，将遥感信息和WOFOST模型相结合，以LAI为状态变量，利用遥感信息校准模型中6个敏感性作物参数，模拟研究区冬小麦产量。研究表明，同化后的区域产量精度相较于未同化有显著提高，说明基于查找表算法的遥感信息与作物模型的同化在区域产量预测方面具有广阔的应用前景。