基于同化叶面积指数和条件植被温度指数的冬小麦单产估测

张树誉，孙辉涛，王鹏新，景毅刚，李 俐

(1.陕西省气象局， 陕西 西安 710014； 2.中国农业大学信息与电气工程学院， 北京 100083)

传统的作物估产主要采用统计调查、农业气象预报等方法，成本偏高、效率较低，难以实现区域尺度高精度作物产量的估测[1]。目前，作物生长模型已成功应用于单点尺度的作物生长发育过程模拟和产量估测，但由于难以获取能够准确描述大范围非均匀性地表、近地表环境空间的模型输入参数，使其难以应用于区域尺度[2]；而卫星遥感具有及时、宏观、信息量大等特点，并能够定量地监测和描述作物在区域尺度的生长状况和反映环境因子对作物生长的综合影响[3]。数据同化方法能够有效融合作物生长模型和卫星遥感观测信息，对现代化农业的发展有着十分重要的意义，已经成功应用于区域产量估测[4-5]。

同化变量(参数)的选取对同化结果精度至关重要，基于多变量数据同化策略能够综合多因子对作物产量共同作用的影响，已经成为农业数据同化的研究趋势[6-7]。叶面积指数(LAI)是评估作物籽粒潜在产量的重要指标[8]，土壤水分(SM)由于控制着作物水分的胁迫信息，与作物产量变化密不可分[9]。Ines等[10]采用改进的集合卡尔曼滤波(EnKF)算法对CERES-Maize模型模拟和遥感数据观测的LAI和SM实施同化，结果表明，同时同化LAI和SM双变量的模型估测结果精度更高。Huang等[11]采用四维变分(4DVAR)算法对SWAP模型模拟和遥感数据观测的LAI和ET实施同化，结果表明，与单独同化LAI或ET单变量相比，同时同化LAI和EI双变量明显提高了冬小麦单产的估测精度。

目前，获取同化所需的观测土壤水分的方式主要是微波遥感技术[10，12]，但由于其在高植被覆盖度条件下对土壤浅层水分的反演精度较低[13]，因此会影响同化结果的精度。王鹏新等[14]在归一化植被指数(NDVI)和地表温度(LST)的散点图呈三角形区域分布的基础上，提出了条件植被温度指数(VTCI)的干旱监测方法，并在陕西省关中平原冬小麦的干旱监测、预测及其影响评估等研究中得到了广泛应用[15-18]。以往的研究表明，VTCI与土壤浅层水分具有显著相关性[19]，适合高植被覆盖度下的干旱监测[20]。基于此，本研究以LAI和VTCI为同化系统状态变量，采用粒子滤波(PF)同化算法对CERES-Wheat模型模拟和遥感观测的LAI和VTCI实施同化，运用熵的组合赋权方法获取加权观测和加权同化变量值，结合样点实测单产构建基于观测和同化变量的估产模型，并分别评价其单点和区域尺度的估测结果精度，以期更准确地进行作物长势监测、干旱监测预测及其影响评估。

1 材料与方法

1.1 试验区概况

关中平原位于陕西省中部，素有“八百里秦川”之称，其行政区域包括西安、宝鸡、咸阳、渭南、铜川5个省辖市和杨凌国家农业高新技术产业示范区。区域内地势平坦，土层深厚，蓄水性好，但是区域内水资源相对不足，干旱是造成粮食减产的主要原因[21]。作物种植模式在灌溉区域主要为冬小麦和夏玉米轮作，而在旱作区域多为夏季休闲式的冬小麦连作。冬小麦播种时间一般为10月上、中旬，播种5～6 d后出苗，4月下旬进入抽穗期，乳熟期为5月中、下旬，5月下旬或6月上旬收获[22]。研究选择2008—2014年冬小麦生长季内，关中平原12个典型的小麦种植区作为试验样点(图1)，其中眉县常兴镇、扶风县城北、三原县鲁桥镇和临渭区蔺店镇为灌溉样点，其余为旱作(雨养)样点。

图1研究区域

Fig.1 Study area

1.2 数据来源及预处理

1.2.1 CERES-Wheat模型模拟LAI和VTCI 在CERES系列模型中，CERES-Wheat模型是面向小麦生长和发育过程的决策等级模型。CERES-Wheat模型的输入数据主要包括气象数据、土壤数据、田间管理数据和小麦品种遗传特性参数。除小麦品种遗传特性参数外，均由实地调研和实测获取。冬小麦的遗传特性参数控制着其生长发育进程，直接关系到植株物理形态的发育与作物产量的形成，因此模型在实际使用前先结合当地实际测量的LAI、生物量、单产、收获日期等对这些参数进行本地化标定[22-23]。

同化所需模拟LAI和VTCI时间序列数据的获取。运行标定后的模型得到样点模拟LAI时间序列数据。研究所需的样点模拟VTCI不能由模型直接模拟得到，但以往的研究结果表明，该地区多年旬尺度的VTCI与土壤浅层(0～20 cm)水分具有显著相关性[19]，因此本文将观测VTCI与模拟浅层水分数据进行线性回归分析，建立经验关系，获取样点模拟VTCI时间序列数据。

1.2.2 MODIS遥感数据反演LAI和VTCI 基于Aqua-MODIS日地表反射率产品(MYD09GA)和日地表温度产品(MYD11A1)，利用MODIS重投影工具MRT(MODIS re-projection tool)对原始数据进行裁剪、投影转换、重采样、格式转换等预处理操作，得到研究区日LST和日地表反射率数据。利用日地表反射率数据计算日NDVI，应用最大值合成技术分别生成旬NDVI和旬LST最大值合成产品，并依据VTCI计算方法获取旬VTCI[14,24]。本研究LAI选取2008—2014年冬小麦主要生育期内所有MCD15A3产品，时间分辨率较高，每4 d合成一次，空间分辨率为1 km，更有利于农作物长势和物候的监测。利用MRT将所有数据转换为统一的Lambert等面积投影，对投影后的数据作裁剪处理，得到关中平原的遥感观测LAI。以样点所在像素为中心3像素×3像素内所有像素的LAI和VTCI的均值作为同化所需观测LAI和VTCI时间序列数据。

1.3 研究方法

(1)

(2)

根据敏感性分析结果，本研究将粒子数设置为200，对于LAI和VTCI，基于田间实测和以往经验，CERES-Wheat模型模拟误差标准差分别设定为0.5和0.1，分别设定遥感观测的13%和3%作为观测误差标准差。

1.3.2 熵的组合赋权方法 将冬小麦越冬后的生育期划分为返青期(3月上旬—3月中旬)、拔节期(3月下旬—4月中旬)、抽穗～灌浆期(4月下旬—5月上旬)和乳熟期(5月中旬—5月下旬)4个生育时期[18]。应用研究区域2008—2014年冬小麦4个生育时期(n)的LAI和VTCI数据构建数据矩阵A(amn)N×4(N为研究选取的总年份数)，计算各个生育时期的熵值hn[26]：

(3)

计算第n个生育时期的差异性系数gn，并对其进行归一化处理得到第n个生育时期的权重wn：

(4)

(5)

1.4 时间尺度的匹配与估产模型的构建

由于研究选用的观测(未同化)和同化VTCI均是以旬为尺度，而观测LAI是4 d合成产品，同化LAI以1 d为步长。为了解决两变量时间尺度不一致的问题，在冬小麦的4个生育时期将观测LAI采用S-G滤波方法[27]插值成以1 d为步长，然后将各生育时期S-G滤波后的观测LAI和经过PF算法的同化LAI累加求和分别作为该生育时期的累积观测LAI和累积同化LAI。取各生育时期内所包含的多旬观测和同化VTCI的均值分别作为该生育时期的观测和同化VTCI。运用熵的组合赋权方法获取各生育时期的加权LAI和VTCI，将观测和加权LAI和VTCI分别与样点实测单产进行一元线性回归分析构建冬小麦单产估测模型(研究选用2012年用于单产估测模型的精度验证，故2012年各样点的相关数据均未参与模型的构建)。

2 结果与分析

2.1 LAI和VTCI的同化结果与分析

通过对2008—2014年12个样点的综合对比分析，以2014年代表性较好的乾县石牛乡旱作样点为例，进行LAI和VTCI的同化结果分析(图2)。关中平原的冬小麦一般在10月上旬播种，之后进入越冬阶段，次年的3月上旬进入返青期，返青期后LAI逐渐增大，至抽穗～灌浆期达到最大。经过PF算法的同化LAI变化趋势(图2a)，可以看出，同化LAI最大值出现时间在4月25日左右，在合理范围内。在曲线上升阶段，模拟LAI变化趋势从3月2日至3月14日左右出现强烈的振荡现象，由于综合了遥感观测LAI变化趋势，同化LAI变得更为平滑；在曲线下降阶段，遥感观测LAI变化趋势从5月17日至5月24日左右出现先缓慢下降再缓慢上升又迅速下降的过程，而同化LAI从5月12日至5月18日左右变化趋势与遥感信息基本一致。由于CERES-Wheat模型是机理模型，在针对研究区进行参数和模型准确标定后，其LAI绝对值更接近于实际情况，根据遥感观测特点，如果不受到严重的云污染影响，遥感LAI变化趋势可直接反映作物LAI的实际变化情况。因此，经过PF算法的同化LAI变化趋势与冬小麦实际生长变化情况更为相符。

研究已经证明VTCI是一种实时的干旱监测方法，VTCI的值越小，干旱程度越严重，土壤水分含量较小，作物受水分胁迫的程度就愈严重；反之，则干旱程度越轻，土壤水分含量较大，作物受水分胁迫程度较轻，其与降水量数据有显著的相关性[19]。为验证同化VTCI是否能较观测VTCI和模拟VTCI更好地与降水量数据结合，将观测VTCI、模拟VTCI和同化VTCI与旬累积降水量数据对比分析(图2b)，可以看出，整体上同化VTCI均能较好地表达模型模拟和遥感观测VTCI，具体表现为：乾县石牛乡3月上旬和5月下旬的降水量仅为6.4 mm和5.1 mm，但观测VTCI分别达到0.76和0.80，属于不旱范围[28]，受到模拟VTCI的影响，同化VTCI分别调整为0.34和0.32，调整幅度较大，结果更能反映样点实际水分胁迫程度。另外3月中旬、4月上旬、4月下旬等同化VTCI在原观测基础上均有不同程度的调整，结果能更好地与旬累积降水量数据结合，更加符合关中平原实际干旱情况。

图2 2014年乾县石牛乡叶面积指数(LAI)和条件植被温度指数(VTCI)变化趋势

Fig.2 Comparison of seasonal changes of leaf area index (LAI) and vegetation temperature condition index (VTCI) in Shiniu of Qian County in 2014

2.2 估产模型的构建与分析

将2008—2014年(除2012年)每年主要生育时期的累积观测LAI和VTCI以及累积同化LAI和VTCI分别运用熵的组合赋权方法计算获取每年加权LAI和VTCI，并将加权变量值分别与样点每年实测单产进行线性回归分析，构建基于观测和同化LAI和VTCI的冬小麦单产估测模型(表1)。结果表明，各估产模型均通过了1%的显著性水平检验，无论是基于观测还是同化变量，基于LAI和VTCI双变量构建的估产模型的估测结果精度均高于单独基于LAI或VTCI单变量构建的结果，且基于同化变量构建的估产模型的估测结果精度明显优于基于观测变量构建的结果。其中，基于观测LAI单变量构建的估产模型的估测结果精度最低(R2=0.279)，可能的主要原因是MODIS遥感观测LAI受云和混合像素的影响会导致整体偏低，与冬小麦实际情况偏差较大。通过对比分析，可以看出，基于同化LAI和VTCI双变量构建的估产模型的估测结果精度最高(R2=0.547)，确定其为最优的单产估测模型。

应用基于同化LAI和VTCI双变量构建的最优单产估测模型对2012年12个研究样点进行单产估测，并将估测单产与实测单产进行一元线性回归分析。结果表明，回归模型的决定系数较高(R2=0.594)，达显著水平(P<0.01)；估测单产与实测单产间间的均方根误差为527.10 kg·hm-2，表明应用基于PF算法的同化变量构建的单产估测模型的单点尺度估产结果误差较低，精度较高。

2.3 区域冬小麦单产的估测与分析

应用最优单产估测模型对2008—2014冬小麦生长年关中平原冬小麦区域单产进行估测(图3)。结果表明，从年际变化上看，关中平原2008—2014年的单产呈现个别年份波动，但整体平稳增长的趋势，其中2013年由于旱情严重而减产较为严重，这与关中平原实际单产及旱情的分布特点较为一致。

表1 基于观测和同化LAI和VTCI的冬小麦单产估测模型

注：Y表示估测的冬小麦单产，单位为kg·hm-2。 Note:Yrepresents the estimated yield (kg·hm-2).

图3 2008—2014年关中平原冬小麦单产的估测结果

Fig.3 Estimated yields of winter wheat in the Guanzhong Plain in the years from 2008 to 2014

通过对区域冬小麦估测单产的空间分布(图3)进行分析，可以看出，关中平原冬小麦单产整体呈现中部单产高于东部和西部，西部单产高于东部的空间分布规律，这与关中平原实际情况相符。综上所述，基于PF算法的同化LAI和VTCI双变量构建的估产模型的区域单产估算结果在时空分布上均与关中平原实际情况一致，说明基于PF算法的同化变量构建的估产模型的精度高，适用于区域冬小麦单产的估测。

3 结论与讨论

粒子滤波(PF)同化算法通过构建一组具有权重的随机样本，以样本均值代替复杂的后验概率积分运算，随着粒子数目的不断增加，这些粒子的概率密度函数逐渐逼近最优贝叶斯估计的效果，可用于非线性非高斯动态系统[29]。研究应用PF算法对CERES-Wheat模型模拟和遥感数据反演的LAI和VTCI实施同化。结果表明，相比于模型模拟和遥感反演的LAI，同化LAI具有良好的时间和空间连续性，在MODIS-LAI偏低的情况下同化LAI有所提升，且同化LAI变化趋势更加符合关中平原冬小麦实际变化情况。同化VTCI能综合表达模拟和观测VTCI的变化趋势，且与旬累积降水量的相关性更高，能更准确地反映冬小麦水分胁迫的程度。

运用熵的组合赋权方法分别构建基于观测LAI和VTCI以及经过PF算法的同化LAI和VTCI的冬小麦单产估测模型。结果表明，基于同化LAI和VTCI构建的估产模型的单产估测精度明显优于基于观测LAI和VTCI构建估产模型的单产估测精度，且基于LAI和VTCI双变量构建估产模型的估测结果优于基于LAI或VTCI单变量构建的结果。采用最优估产模型分别对单点和区域尺度的冬小麦单产进行估测，结果表明，2012年12个样点的估测单产与实测单产间的相关性达显著水平，两者间的均方根误差为527.10 kg·hm-2；区域冬小麦单产估测结果在时空分布上均与关中平原实际情况相符，说明无论是单点还是区域尺度，基于PF算法同化LAI和VTCI双变量构建的估产模型的估测精度均较高。

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

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

[3] Wiegand C L, Richardson A J, Kanemasu E T. Leaf area index estimates for wheat from Landsat and their implications for evapotranspiration and crop modeling[J]. Agronomy Journal, 1979,71(2):336-342.

[4] 王鹏新,孙辉涛,王 蕾,等.基于4D-VAR和条件植被温度指数的冬小麦单产估测[J].农业机械学报,2016,47(3):263-271.

[5] 靳华安,王锦地,柏延臣,等.基于作物生长模型和遥感数据同化的区域玉米产量估算[J].农业工程学报,2012,28(6):162-173.

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

[7] Nearing G S, Crow W T, Thorp K R, et al. Assimilating remote sensing observations of leaf area index and soil moisture for wheat yield estimates: An observing system simulation experiment[J]. Water Resources Research, 2012,48(5):W5525.

[8] Huang J, Sedano F, Huang Y, et al. Assimilating a synthetic kalman filter leaf area index series into the WOFOST model to improve regional winter wheat yield estimation[J]. Agricultural and Forest Meteorology, 2016,216:188-202.

[9] He B, Li X, Quan X, et al. Estimating the aboveground dry biomass of grass by assimilation of retrieved LAI into a crop growth model[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015,8(2):550-561.

[10] Ines A V M, Das N N, Hansen J W, et al. Assimilation of remotely sensed soil moisture and vegetation with a crop simulation model for maize yield prediction[J]. Remote Sensing of Environment, 2013,138(6):149-164.

[11] Huang J, Ma H, Su W, 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.

[12] Zhang H, Qin S, Ma J, et al. Using residual resampling and sensitivity analysis to improve particle filter data assimilation accuracy[J]. IEEE Geoscience and Remote Sensing Letters, 2013,10(6):1404-1408.

[13] Bindlish R, Jackson T J, Gasiewski A J, et al. Soil moisture mapping and AMSR-E validation using the PSR in SMEX02[J]. Remote Sensing of Environment, 2006,103(2):127-139.

[14] 王鹏新,龚健雅,李小文.条件植被温度指数及其在干旱监测中的应用[J].武汉大学学报:信息科学版,2001,26(5):412-418.

[15] 张树誉,李 慧,王鹏新,等.条件植被温差指数干旱监测方法的研究与应用[J].干旱区资源与环境,2014,28(4):118-122.

[16] 王 蕾,王鹏新,田 苗,等.效率系数和一致性指数及其在干旱预测精度评价中的应用[J].干旱地区农业研究,2016,34(1):229-235.

[17] 林 巧,王鹏新,张树誉,等.不同时间尺度条件植被温度指数干旱监测方法的适用性分析[J].干旱区研究,2016,33(1):186-192.

[18] 李 艳,王鹏新,刘峻明,等.基于条件植被温度指数的冬小麦主要生育时期干旱监测效果评价Ⅲ—干旱对冬小麦产量的影响评估[J].干旱地区农业研究,2014,32(5):218-222.

[19] 王 维,刘翔舸,王鹏新,等.条件植被温度指数的四维变分与集合卡尔曼同化方法[J].农业工程学报,2011,27(12):184-190.

[20] 王鹏新,孙辉涛,解 毅,等.基于LAI和VTCI及粒子滤波同化算法的冬小麦单产估测[J].农业机械学报,2016,47(4):248-256.

[21] 何慧娟,卓 静,李红梅,等.基于MOD16产品的陕西关中地区干旱时空分布特征[J].干旱地区农业研究,2016,34(1):236-241.

[22] 刘骁月,王鹏新,张树誉,等.基于作物模型模拟年际生物量变化的冬小麦干旱监测研究[J].干旱地区农业研究,2013,31(1):212-218.

[23] 贺 鹏,王鹏新,解 毅,等.基于动态模拟的冬小麦水分胁迫敏感性研究[J].干旱地区农业研究,2016,34(1):213-219.

[24] 孙 威,王鹏新,韩丽娟,等.条件植被温度指数干旱监测方法的完善[J].农业工程学报,2006,22(2):22-26.

[25] 姜志伟,陈仲新,任建强,等.粒子滤波同化方法在CERES-Wheat作物模型估产中的应用[J].农业工程学报,2012,28(14):138-146.

[26] 苏 涛,王鹏新,刘翔舸,等.基于熵值组合预测和多时相遥感的春玉米估产[J].农业机械学报,2011,42(1):186-192.

[27] Savitzky A, Golay M J E. Smoothing and differentiation of data by simplified least squares procedures[J]. Analytical Chemistry, 1964,36(8):1627-1639.

[28] 张树誉,孙 威,王鹏新.条件植被温度指数干旱监测指标的等级划分[J].干旱区研究,2010,27(4):600-606.

[29] Weerts A H, El Serafy G Y H. Particle filtering and ensemble kalman filtering for state updating with hydrological conceptual rainfall runoff models[J]. Water Resources Research, 2006,42(9):W9403.