水氮数据联合驱动的冬小麦数据同化研究

李金敏，韩景晔，余丹阳，史良胜

（武汉大学水资源与水电工程科学国家重点实验室，武汉430072）

0 引 言

中国是世界上最大的小麦生产国和消费国，小麦也是中国第3大农作物，仅次于水稻和玉米[1,2]。因此对冬小麦进行准确的生长模拟和产量预测是十分重要的。作物模型从系统科学的观点出发，运用数学物理方法和计算机技术，对作物生长过程中的光合、呼吸、蒸腾等重要生理生态过程进行定量描述和预测[3]，是作物生长发育及其对环境影响评估的重要工具[4]。

在过去的半个世纪中，研究者们开发出很多作物模型来模拟冬小麦的生长过程，比如AquaCrop[5,6]，DSSAT[7,8]，WOFOST[9]，SWAP[10]和STICS[11]等。然而，作物模型的精度通常会受到模型结构、参数和输入等因素引起的不确定性的影响，导致模拟结果和实际情况有偏差[11]。

数据同化是一种集成观测和模型这2种基本科学研究手段的重要方法，通过融合多源观测和带有不确定性的模型，从而改善对作物生长过程的模拟[13,14]。目前，数据同化方法在作物生长模拟和产量估计研究中受到广泛关注[15-19]。Wang 等结合遥感数据反演的LAI和叶片氮累积量（LNA），利用粒子群优化算法对作物生长模型RiceGrow 的播种日期、播种密度和施氮量等参数进行初始化参数校准，提高了水稻产量的预测精度[20]。Huang 等采用集合均方根滤波（EnSRF）方法，将遥感反演的LAI 和LNA作为观测加入到WheatGrow 模型数据同化系统中，成功地估计了冬小麦LAI和LNA，提高了产量的估计精度[21]。然而，这些作物模型数据同化系统大多关注作物叶面积指数（LAI），对土壤和作物水氮状态的数据同化研究相对较少。实际上，水和氮是作物生长的2大关键因素，对于土壤和作物水氮状态的准确估计是十分重要的[22-26]。因此，本研究中尝试利用在数据同化系统中同时更新土壤含水量（SM）和作物叶片氮累积量（LNA），探究水氮数据联合驱动更新对作物生长模拟的影响。

本研究中的冬小麦生长模型为SWAP-WOFOST 模型，该模型采用WOFOST 的作物生长模块，兼顾土壤-作物水氮平衡，已经成功地应用于作物的生长模拟和产量估计研究中[17，18，27]。基 于SWAP-WOFOST 模 型 构 建 虚 拟 观 测 模 拟（OSSE）系统，采用集合卡尔曼滤波（EnKF）方法进行顺序性数据同化，围绕SWAP-WOFOST-EnKF 系统开展系统研究：土壤含水量和叶片氮累积量的联合更新对模拟结果的影响、观测数据的价值分析和不同数据同化更新策略对模拟结果的影响。

1 材料与方法

1.1 SWAP模型

SWAP（Soil-Water-Atmosphere-Plant）模型是荷兰瓦赫宁恩大学开发的一个垂向一维农业水文模型，它可以用来模拟非饱和区域内土壤水、溶质和热量的运移与作物生长之间的关系。SWAP 采用Richards 方程模拟土壤水运动，通过对流、弥散和扩散方程来描述溶质运移转化等过程。SWAP模型中的作物生长部分采用WOFOST 模块，通过量化光合作用、呼吸消耗等过程，并考虑水肥胁迫等环境影响，形成的干物质在不同的器官间进行分配，而这一切都受到作物生育期（DVS）的影响，DVS在种子萌发时为0，在开花时取1，最终成熟时DVS=2。通过耦合土壤氮模块，计算出土壤中氨态氮、硝态氮的动态平衡。土壤中氮素通过作物根系吸收进入作物体内，并根据作物氮需求进行氮素分配，通过计算氮素胁迫因子NNI来影响作物光合、叶生长、氮素分配等过程。

1.2 数据同化方法

集合卡尔曼滤波（EnKF）是一种顺序性数据同化算法，可以通过不断获取观测数据对模型状态和参数进行持续更新。由于其具有较低的计算成本，并且具有一定程度非线性问题适应能力，在农业水文领域有着广泛的应用[16,28-30]。

进行同化的状态向量通常由模型参数、状态和观测共同组成：

式中：yt为t时刻的状态向量；为模型参数，是模型状态变量，如土壤含水量和作物叶面积指数等；是观测值，上标T表示矩阵转置。

模型观测可以表示为：

因此有：

式中：H为观测算子，将状态向量矩阵转换为观测值相对应的模型预报值。

采用EnKF对模型状态向量进行更新：

式中：j为样本序号为数据同化更新后的状态向量为模型预测的状态向量；d t,j为扰动后的观测值；Kt是卡尔曼增益，记为如下公式。

式中：为模型预测状态向量的协方差矩阵；为状态向量与观测向量的互协方差；为观测向量的自协方差；为观测误差。

1.3 虚拟观测模拟系统

本研究中，采用虚拟观测模拟（Observing Simulation System Experiment， OSSE）系统完成相关研究内容。OSSE 系统共包括N+1次样本模拟，其中的1个样本模型结果作为参照值，余下的N个样本由带有偏差的先验参数抽样后模拟得来，采用参照值来评价模拟效果。在本研究中，样本量N=50，与de Wit等和Curnel等的研究保持一致[31,32]。

1.3.1 SWAP-WOFOST 模型输入数据

气象数据来源于武汉大学灌溉排水试验场的气象站，包括气温（最高和最低温度）、太阳辐射、降雨、空气水汽压和风速，数据频率为1 d，模拟生长时间从2017年的11月18日至2018年的5月10日，共174 d。模型上边界为大气边界，下边界采用隔水边界。施肥情况为在作物播种当天施加氮素100 kg/hm2。

1.3.2 不确定性参数

在进行数据同化模拟时，需要通过扰动不确定性参数来生成样本，因此敏感参数的选择是非常重要的一步。本研究中采用Morris筛选法[33]，主要分析模型输出的产量（Yield）、LAI和叶片氮累积量（LNA）对50个参数的敏感性。Morris方法首先计算参数的基效应：

式中：EEi为第i个模型参数的基效应数值；x1，x2，…，xm为m个模型参数的取值；Δ为参数增量；f为模型状态输出。

通过计算基效应EEi绝对值的均值来得到敏感性指标μ*，μ*越大，表示参数越敏感。本研究采用的样本数量为1 000，模型参数为50 个，共需计算51 000 次运算，计算的敏感性结果见图1。由图1 可知，模型的50 个参数并不都是敏感参数，只有部分参数在模型运算过程中起到重要的影响。其中土壤水参数ALFA、NPAR、OSAT，土壤氮参数TCSF_N 和作物参数EFF、CVL、FS_090 对于模型输出的产量、LAI和LNA都有影响，尤其是TCSF_N，敏感性指标数值很高；有些参数只对产量敏感，比如CVO 和FO_110。结合HU 等的研究，本文中选取了ALFA、NPAR、TSUMEA、EFF、CVL、CVO、TCSF_N等关键参数[16]。其中新加入的TCSF_N 为蒸腾浓度流系数，与作物根系吸收土壤氮的速率有关，可以反映作物吸收氮的能力，从而影响作物生长。假设所有不确定性参数服从对数正态分布，初始均值和方差为：ln（ALFA）～（-3.0，0.10），ln（NPAR）～（0.53，0.05），ln（TSUMEA）～（-6.91，0.03），ln（EFF）～（-0.60，0.08），ln（CVL）～（-0.29，0.03），ln（CVO）～（-0.22，0.03），ln（TCSF_N）～（-0.36，0.10）。

1.3.3 算例设计

本研究基于SWAP-WOFOST 模型构建冬小麦作物生长的数据同化系统，利用多种观测联合估计作物水氮状态。选择合适的状态变量是进行数据同化模拟的关键。传统研究大都将土壤含水量、LAI作为状态变量，本研究在此基础上加入了叶片氮累积量（LNA）这一状态变量，进行水氮状态的联合更新。其中观测时间为作物生长的第10到第160 d，观测频率为10 d一次，全生育期一共16次观测。土壤水观测为0～10 cm的平均土壤含水量（SM），观测误差为±0.04 m3/m3，LAI的观测误差为±0.4 m2/m2，LNA的观测误差为±4 kg/hm2。

算例设计见表1。其中，C1 为利用传统的土壤含水量和LAI观测进行数据同化，状态变量为负压（H）、土壤含水量（SM）、叶干重（LV）、叶龄和叶片氮累积量（LNA）。在C1 的基础上，C2 的观测中加入叶片氮累积量（LNA），试图探索加入LNA观测对数据同化系统的影响。C2 与C3 的对比是为了研究2种同化策略对土壤水分剖面和作物生长模拟的价值，其中SSPE 指的是同时更新作物模型的状态和参数，而USO 只更新模型状态。在C4中，观测误差由固定值改变为观测值的10%，通过和C2 结果对比，研究观测误差对数据同化系统模拟的影响。在C5 中，通过将观测频率从10 d 增加到20 d，从而探究不同观测频率对作物生长模拟的影响。

表1 数据同化算例设计Tab.1 Case design of data assimilation

2 结果与分析

2.1 更新LNA的重要性

图2展示了算例C1和C2模拟的LAI、LNA、作物产量以及5、20 和40 cm 处土壤含水量随时间的变化，其中横坐标为作物播种后天数，“Reference”表示OSSE 系统模拟的真值，“Open-loop”表示不进行数据同化的确定性模型结果。由图2（a）可以发现，通过在算例C1 中加入LAI观测可以实现较好的LAI模拟，但是仍有一定改善空间；算例C2 中加入LNA观测后，作物LAI的模拟更加接近真实值，表明加入作物叶片氮信息有助于LAI的更新，但是总体来看，C1和C2对于LAI的模拟结果相差不大。由图2（b）中可知，算例C1 中由于仅加入LAI的观测，不能很好地模拟出作物叶片氮累积量的信息，尤其是在作物生长的中间阶段（第100～130 d）。相比于C1，C2可以更好地模拟出作物生长后期的LNA变化，从而可以更加精准地模拟作物的产量。然而，无论是C1 还是C2，对于作物早期的状态都没有能够很好地重现，这主要是由于数据同化过程中给定的观测误差较大，尤其在早期作物LAI和LNA都比较小的情况下，给出较大的观测误差会影响状态向量的更新，因此对于作物早期状态的估计都存在一定的改善空间。由图2（c）可以发现，C1 中对于作物产量的估计效果很差，而加入LNA观测的C2算例能够较好地重现产量的积累过程。

图2（d）～图2（f）展示了算例C1 和C2 模拟的5、20 和40 cm 处的土壤含水量随时间的变化，可以发现C1 和C2 都能较好地重现出土壤含水量的剖面分布。整体来看，C1 对于土壤含水量的模拟效果优于C2，可能是由于加入LNA观测后，引入了新的观测误差，导致C2 中对于土壤含水量的模拟精度略低于C1。

对比算例C1 和C2 对于4 种状态变量（SM、LAI、LNA和作物产量）的更新结果可以发现，仅加入LAI这一个作物状态的观测可以改善模型对LAI的模拟，但是对于作物产量的模拟效果不佳，主要是由于没有很好地更新模型的氮素信息；在模型能够准确模拟水分状态的情况下，通过加入LNA的观测，可以有效提升作物对氮素状态的模拟，从而提升作物LAI和产量的估计精度。

2.2 SSPE和USO数据同化策略对比研究

图3 展示了不同数据同化策略下算例C2 和C3 模拟的LAI、LNA、作物产量以及5、20 和40 cm 处土壤含水量随时间的变化。相比于算例C1，可以发现C2 和C3 都能很好地重现作物LAI和LNA随时间的变化，主要是因为C2 和C3 中都加入了LNA的观测，能够很好捕捉作物氮素信息。然而，相比于同时更新模型状态和参数的C2，仅更新状态的C3在作物产量的估计上效果不佳。

图3（d）～图3（f）展示了算例C2 和C3 模拟的5、20 和40 cm 处的土壤含水量随时间的变化，可以发现不更新参数的C3 算例不能很好地模拟5 cm 处和20 cm 处的土壤含水量，对于土壤水分状态的模拟存在偏差。

对比算例C2 和C3 对于4 种状态变量（土壤含水量、LAI、LNA和作物产量）的更新结果可以发现，只更新模型状态的C3 算例尽管能够很好地模拟作物状态，但是由于缺乏对土壤水分状态的正确估计，导致其产量估计依旧出现很大的误差；而同时更新状态和参数的SSPE 方法可以很好重现土壤—作物水分和氮素的变化，因此能够更好地实现作物的产量估计。图4展示了在SSPE方法中土壤水参数NPAR、作物参数EFF和氮素参数TCSF_N 的变化过程，可以发现3 类参数在数据同化模拟过程中都会往参照值靠近，不确定性逐渐减小，即SSPE方法可以同时减小模型参数和过程的不确定性，从而获得更加准确的模拟效果。

对比C1、C2和C3这3种算例可以发现，C1很好地重现了土壤含水量的变化，但是不能准确模拟作物氮素状态；而C3很好地重现了作物叶片氮素累积的变化，但是不能准确模拟土壤的水分状态，因此这2个算例对于最终作物产量的估计都存在一定偏差。C2 由于同时更新了土壤—作物的水氮状态，有效提升了产量估计的精度。因此，准确的土壤—作物水氮状态的模拟对于作物生长的模拟和产量估计至关重要。

2.3 观测误差的影响

图5 展示了在不同观测误差下，算例C2 和C4 模拟的LAI、LNA和作物产量随时间的变化。C2 采用固定的观测误差值，SM为±0.04 m3/m3，LAI为±0.4 m2/m2，LNA为±4 kg/hm2，而C4采用相对的观测误差，即任一时刻的观测误差为对应观测值的10%。由图5可知，相比于C2，C4可以获得更加精准的LAI和LNA模拟效果，尤其在作物生长的早期，由于C2 采用固定的观测误差，而作物生长早期LAI和LNA值都较小，观测误差就相对较大，因此同化效果不佳；C4 中采用相对观测误差，在整个模拟过程中都能较好地重现作物生长的状态。尽管C4能够完成更好的作物LAI和LNA模拟，C2 的模拟效果也不算差，而且C2估计的最终产量比C2更接近参照值。考虑到观测误差通常取决于观测手段和观测仪器，采用固定观测误差可能更符合实际应用情景。

2.4 观测频率的影响

图6 展示了在不同观测频率下，算例C2 和C5 模拟的LAI、LNA和作物产量随时间的变化。其中C2 采用10 d 一次的观测频率，而C5 为20 d 一次。由图6 可知，减小观测频率的C5 在LAI和LNA模拟过程中波动比C2 要大一些，但是在2 种算例中，最终产量的估计却是一致的，都非常接近参照值。这表明，适当地减小观测频率并不会降低模型的预测能力，数据同化系统在不同观测频率下具有一定的稳健性。对于冬小麦而言，间隔20 d 一次的观测几乎可以覆盖所有的生育阶段，因此从权衡模型表现和观测成本的角度来看，每个生育期1次观测即可满足本研究中的数据同化模拟要求。

3 结 论

传统的数据同化系统只有土壤水和LAI的观测，无法有效估计作物氮素状态，通过加入叶片氮累积量的观测，可以很好地重现作物叶片氮素变化过程。相比于只更新状态的USO方法，同时更新模型状态和参数的SSPE 方法在更新土壤含水量方面具有更加重要的价值。水分和氮素状态的错误估计都会引起作物产量模拟存在偏差，因此水氮数据的联合驱动对于作物生长模拟和产量估计是十分必要的。此外，本研究中的数据同化系统在不同观测误差和观测频率下都具有良好的稳健性。