基于SWAP-IES 的旱区春小麦长势和产量模拟

金建新 ，丁一民 ，孙振源 ，朱 磊 ※

（1. 宁夏大学土木与水利工程学院，银川 750021；2. 宁夏回族自治区黄河水联网数字治水重点实验室，银川 750021；3. 宁夏农林科学院农业资源与环境研究所，银川 750002）

0 引 言

宁夏引黄灌区是中国西北干旱地区重要的小麦生产基地，在全球气候变化和局部冲突的背景下，小麦对保障粮食安全的重要性日益突出，近年来播种面积大幅增加。及时、准确地对小麦长势和产量进行监测，是动态调整田间管理措施，确保小麦稳产增产的重要基础[1]。

基于过程机理的作物生长模型能够根据区域气象条件、土壤性质和管理措施等，定量描述单点尺度作物生长发育和产量形成过程，是优化管理、辅助决策、生长评估的重要工具。目前已从单一模型发展到综合模型阶段，国内外主流的模型主要有农业生产系统模拟器（agricultural production systems simulator，APSIM）、作物-环境资源综合系统（crop environment resource synthesis，CERES）、土壤-水-大气-植物模型（soil-water-atmosphereplant，SWAP）、侵蚀-土地生产力影响评估模型（erosionproductivity impact calculator，EPIC）、世界粮食作物研究模型（world food studies，WOFOST）等[2]。其中，SWAP模型在农田蒸散发、土壤水分运动的模拟中具有明显优势，被广泛用于旱区作物非充分灌溉的模拟研究[3-5]。但品种、土壤的空间异质性导致作物生长模型难以在不同区域上进行准确模拟。同时，模型的结构误差和简化假定（例如不考虑虫害、病害等）也会给模拟结果带来一定的不确定性[6]。借助数据同化方法，将不同观测手段获取的作物状态信息融合到不断演进的作物模型中，校正模型模拟方向，实现作物长势的连续动态监测和准确预估，是提高模型模拟精度的重要手段[7]。

目前，主要的数据同化方法有强迫法、参数优化法和顺序同化法等[8-10]。集合卡尔曼滤波算法（ensemble kalman filter，EnKF）作为一种序列数据同化方法，是作物生长模型数据同化研究中最常用的数据同化算法[11-12]。利用EnKF 算法在WOFOST 模型中融入遥感获得的叶面积指数信息，可将冬小麦估产的R2从无同化的0.06 提高到0.41[13]，在华北平原尺度上冬小麦估产的R2从无同化的0.28 提高到0.65[14]。利用APSIM-EnKF 同化模拟框架也可有效提高小麦估产精度，相对于同化前，小麦估产的相对误差减小了8.7%[15]。与EnKF 算法相比，迭代集合平滑器算法（iterative ensemble smoother，IES）可以同时同化所有可用的观测信息，在强非线性问题中表现得更好[10,16-17]。目前IES 算法主要用在地下水同化模拟、水文地质参数反演等领域，并取得了较好的模拟效果[18-19]。利用IES 算法结合SWAP 模型构建春小麦同化模拟系统的研究还有待开展。

除同化算法外，同化变量观测时期的选择也是影响模拟精度的重要因素。以往作物模型同化研究大多基于Sentinel[20]、Landsat[21]、MODIS[7]等卫星遥感数据获取的作物和土壤状态信息，上述数据具有周期固定、易受不利天气影响等缺点，为获取较好的模拟结果，通常利用作物生育期内所有可用观测数据。近年来，随着低空遥感、机器视觉等信息技术的进步，可以利用无人机、地面相机等手段更加便捷灵活地获取多种农情观测数据，已被广泛应用于田块尺度作物长势监测研究中[22-23]。尽管上述观测手段较为灵活，但频繁观测仍然需要较高的观测成本[24]。因此，亟需开展不同生育期观测信息有效性的评估研究，以节约劳动力和观测成本。此外，水分胁迫是干旱区作物生长面临的主要制约因素之一[25]，评估同化模拟系统在不同水分处理下的模拟精度，对于干旱区作物生长监测和产量估算具有重要意义，但相关研究还较为缺乏。

因此，本文的主要研究目的如下：1）利用IES 算法和SWAP 模型构建SWAP-IES 春小麦同化模拟系统；2）在不同水分处理条件下评估同化叶面积指数、土壤水分及其组合在旱区春小麦生长模拟和估产中的作用，以期为利用无人机、地面相机等观测设备开展不同尺度下的旱区春小麦不同水分管理条件下生长监测、产量预估和精准管理提供技术参考。

1 材料与方法

1.1 田间试验及数据收集

试验于2019—2020 年在宁夏永宁县望洪镇进行，地理位置为106°11′E，38°10′W，海拔1 125 m，属典型的大陆性干旱半干旱气候，四季分明，降雨稀少且多集中在夏季和秋季，冬季寒冷干燥，春旱频发，多年平均降雨量198.19 mm，多年平均蒸发量1 186.73 mm，土壤质地为砂壤土，耕层平均容重为1.45 g/cm3，平均田间持水量为24.38%。

供试春小麦品种为宁春55 号，为灌水量单因素试验，以农户常规灌水量为对照（T5），以450 m3/hm2为递减梯度设置5 个灌水处理，各处理灌水次数均为7 次，具体灌水量见表1。小区长10 m、宽6 m、面积60 m2，每个处理3 次重复，共计15 个小区。小麦播种方式为匀播，采用滴灌供水，滴灌带间距50 cm。施肥总量为氮（N）195 kg/hm2、磷（P2O5）120 kg/hm2、钾（K2O）60 kg/hm2，其中30%氮肥和70%磷钾肥基施，剩余养分在苗期、拔节期和灌浆期随水追施，其他田间操作均保持一致。

表1 各处理灌水定额Table 1 Irrigation quota of each treatmentm3·hm-2

试验收集的数据包括以下部分：

1）土壤含水率：在每个小区中间位置的滴灌带旁布设一根TDR 土壤水分测定管，在小麦生育期内每隔10 d测定一次土壤含水率，土壤剖面100 cm 深度按照20 cm分层，共计5 层，灌水前后和降雨后加测。

2）株高和生物量：利用卷尺测量株高，每个小区随机选择10 株小麦，从地表量至生长点；量完后全部挖出带回实验室洗净，烘箱调至105 ℃杀青30 min 后，在80 ℃温度下烘至恒质量，利用分析天平称质量，株高和生物量每个生育期测定一次。

3）叶面积指数：利用LAI-2 200 植物冠层分析仪（LICOR，美国）测定，在小麦每个生育期，选择晴朗天气的16:00—18:00 背光测量，每个小区测定3 次取算术平均值。

4）产量：在小麦成熟期，每个小区随机选择一个1 m2的样方，剪掉样方内的所有麦穗带回实验室，搓掉称质量，各重复取平均值，换算成每公顷的产量。

5）气象数据：利用MC-QX 小型气象站每30 min自动记录一次试验区附近气象资料，包括太阳辐射、温湿度、2 m 高度风速、降雨、相对湿度等。

1.2 SWAP 模型及春小麦基础参数

SWAP 模型是荷兰瓦赫宁根大学开发的一种水分驱动的模型，自1978 年发布以来，经过不断的补充和完善，已广泛应用于优化灌溉制度、估产、农田管理、农业区域评价等[26]，输入参数主要包括气象参数、土壤参数、作物参数等。参考王学梅等[27-28]在干旱区小麦上SWAP模型参数敏感性分析的结果，选择敏感性较强的5 个作物参数作为同化系统中待校正的参数。其他敏感性相对较低的参数通过田间试验、模型默认或参考他人研究成果的方式确定，具体参数值见表2。

表2 SWAP 模型中春小麦主要基础作物参数值及来源Table 2 Parameter values and sources of the main basic crops of spring wheat in the SWAP model

2019—2020 年的田间试验观测数据一方面用于验证同化模拟系统在不同水分条件下模拟的准确性，另一方面用于确定SWAP 模型中敏感性相对较低的基础参数，作为SWAP-IES 同化模拟系统中旱区春小麦的基本信息。利用田间试验数据率定的基础作物参数和方法如下：

1）开花到成熟时的积温（TSUMAM）：利用气温数据和观测的物候期结合有效积温公式计算得到；2）出苗时的生物量（TDWI）：在出苗时取样烘干称质量测定；3）出苗时的叶面积指数（LAIEM）：在出苗时测定样穴内全部植株的叶片面积换算得到； 4）比叶面积（SLA）：通过测定全生育期内叶片面积和干质量计算得到； 5）最大CO2同化速率：通过测定特征生育期光响应曲线计算得到； 6）叶面积指数日最大相对增速（RGRLAI）：利用全生育期叶面积指数和积温计算得到。

对于土壤参数，将试验区0～100 cm 土壤分为3 层，土壤水力参数根据实测的土壤机械组成、干容重及水分特征曲线数据，利用RETC（retention curve program for unsaturated soils）软件拟合得到，具体参数值见表3。

表3 试验区不同深度土壤水力参数Table 3 Hydraulic parameters of soil at different depths in the experimental area

1.3 SWAP-IES 同化模拟系统

以LAI 和SW 为同化的观测变量，利用IES 数据同化算法，不断调整SWAP 模型中的关键敏感参数，逐渐减小模型的模拟误差，直至达到预设的收敛标准（图1）。

图1 土壤-水-大气-植物模型-迭代集合平滑器同化模拟系统流程图Fig.1 Flow chart of soil-water-atmosphere- plant(SWAP)-iterative ensemble smoother(IES) data assimilation system

IES 算法是一种能同时利用所有观测数据的连续型数据同化方法，其同化思想来源于后验概率密度函数采样的最优化方法，它利用样本集合估计参数-状态的敏感性信息，使用高斯牛顿法迭代更新参数[29]。假设模型输入与输出之间关系可以表示为

式中dobs是实测值向量，F是正向模型，m是未知参数向量，ε为符合均值为0 协方差为CD=E[ε×εT] 的观测误差向量。然后IES 可以分以下几个步骤实施：

第一步：从参数的先验分布中产生Ne个样本构成初始样本集合M0：

第二步：在第l（l=1,2,...）次迭代中，给定所有时刻的历史观测数据，参数样本Ml可以根据式（3）完成更新：

式中βl为调整迭代步长的参数，CM=△M0(△M0)T/(Ne-1)代表参数的先验协方差，其在整个迭代过程中保持不变。△M0表示矩阵M0与其均值的偏差，Gl为基于集合平均的敏感性矩阵。

第三步：重复步骤二，直到达到IES 预设的收敛标准。

1.4 数据同化模拟方案

实践中春小麦品种特性往往具有较大的不确定性，是制约作物生长模型准确模拟作物生长过程的关键因素。因此，本文用最敏感的5 个模型参数来表征未知的品种特性信息，利用田间试验数据评估SWAP-IES 系统是否能够通过吸收观测数据校正模型中未知的关键敏感参数，实现作物长势和产量的准确模拟。根据SWAP 模型敏感性分析结果[27-28]，选择TSUMEA、CVO、SPAN、EFF、CVL 5 个关键敏感的作物参数作为待校正参数，根据文献 [30-31] 及实测数据得到待校正参数的先验均值与上下限，如表4 所示。

表4 待校正参数的先验均值及取值范围Table 4 The priori mean and value range of the parameters to be corrected

本文共设置了12 个同化模拟情景（表5）。情景1不同化观测数据，仅利用待校正参数的先验值进行模拟，情景2、情景3 和情景4 使用了所有6 次观测数据，用来评估吸收不同观测变量及其组合对同化结果的影响。情景5～情景12 按生育期分别添加观测数据，用来评估观测时期对同化结果的影响。模拟中假定土壤在水平方向上是匀质的，作物参数中除待校正的参数外其余参数在不同处理下均保持不变。土壤垂直方向水分运动范围假设为1 m，将垂直方向土壤剖面离散为50 个网格，每个网格2 cm，每个网格压力水头的初始值统一设置为-100 cm，上边界条件为大气边界，下边界条件为自由排水边界。参数TSUMEA、CVO、SPAN、EFF、CVL 在同化过程中以0.04、0.08、0.10、0.24、0.08 为方差生成先验值集合。用于数据同化的LAI 和SW 等观测变量为大田试验的实测数据，利用高斯噪声将观测变量分别以0.10 和0.01 为方差进行扰动，生成观测值样本集合。情景1～情景4 在模型中同化观测数据后优化得到的不确定性参数取值见表6。

表5 同化模拟情景设置Table 5 Assimilation simulation scenario settings

表6 不同处理各同化情景下待校正参数的取值Table 6 The values of parameters to be corrected under different assimilation scenarios of different treatments

1.5 模拟精度评价指标

采用决定系数（R2）、均方根误差（SRMSE）、标准均方根误差（SNRMSE）和相对误差（SRE）等指标对模拟结果进行评价，其中R2越接近1，SRMSE和SRE越小表示模拟精度越高，SNRMSE≤10%表示精度极高，10%＜SNRMSE≤20%表示高精度，20%＜SNRMSE≤30%表示中等精度，SNRMSE＞30 %表示精度较低[32]。

2 结果与分析

2.1 同化不同观测变量对土壤水分模拟的影响

图2 为SWAP-IES 系统同化不同观测变量对2019和2020 年T1～T5 处理土壤水分模拟的影响。不同化任何观测数据时，模型对土壤水分模拟的精度最低（R2=0.48，SRMSE=0.058 cm3/cm3，SNRMSE=43.35%）。同化观测变量后，模型对土壤水分的模拟均达到中等以上精度，其中同化LAI 时模拟的R2=0.73，SRMSE=0.043 cm3/cm3，SNRMSE=29.26%，为中等精度。同化LAI+SW 时，对土壤水分的模拟精度较同化LAI 有所提高，但仍然为中等精度（SNRMSE＞20%）。仅同化SW 时SWAP-IES 系统对土壤水分模拟的精度最高，R2=0.87，SRMSE=0.029 cm3/cm3，SNRMSE=19.08%，表现为高精度。这说明同化土壤水分可显著提高SWAP-IES 系统对土壤水分的模拟精度，为后期准确模拟作物水分胁迫状况、腾发量等奠定基础。

图2 同化不同观测变量对土壤水分模拟的影响（2019—2020）Fig.2 Effects of assimilation of different observation variables on soil moisture simulation（2019—2020）

2.2 各处理同化不同观测变量对生长指标模拟的影响

以2019 年T1 为例来说明SWAP-IES 系统同化不同观测变量时各待校正参数样本对春小麦LAI、生物量和株高的模拟（图3）。不同化任何观测数据时，各样本对LAI、生物量和株高模拟的平均值和实测值的SRMSE分别为0.61 m2/m2、1 550.19 kg/hm2和7.89 cm。同化SW时，各样本对3 指标模拟的平均值和实测值的SRMSE分别为0.52 m2/m2、984.33 kg/hm2和6.08 cm（图3a）。同化LAI 时对3 指标模拟的SRMSE分别降低为0.30 m2/m2、733.94 kg/hm2和2.68 cm（图3b）。同化LAI+SW 时对3 指标模拟的SRMSE分别为0.32 m2/m2、437.63 kg/hm2和2.28 cm（图3c）。SNRMSE和R2也类似，同化SW 时对各指标模拟的SNRMSE较无同化降低12.8%～54.47%，R2提高17.24%～33.33%。同化LAI 时，各指标模拟的SNRMSE较无同化下降了50.07%～79.89%，R2提高了42.18%～66.67%。同化LAI+SW 时，各指标模拟的SNRMSE较无同化降低45.9%～82.88%，R2提高48.28%～68.4%。说明在T1 处理中仅同化LAI 时模型对LAI 的模拟精度最高，而同化LAI+SW 对生物量和株高的模拟效果最好。

图3 SWAP-IES 系统同化不同观测变量对小麦生长指标的模拟（2019 年T1 处理）Fig.3 Simulation of wheat growth indicators by assimilating different observation variables in the SWAP-IES system (T1 in 2019)

可见，在SWAP-IES 系统中仅利用参数样本的先验值无法准确模拟春小麦生长过程，同化不同观测变量，均能改善春小麦生长指标的模拟效果，降低各样本模拟的不确定性。

表7 为SWAP-IES 系统同化不同观测变量对2019—2020 年各处理生长指标模拟的R2、SRMSE和SNRMSE值。不同化任何观测数据时，模型对LAI、株高和生物量模拟的R2介于0.31～0.67，SRMSE分别为0.63～0.96 m2/m2、5.90～7.99 cm 和1 289.00～1 808.00 kg/hm2，SNRMSE均大于20%，各处理均处于中等精度以下。同化SW 时模型对LAI、株高和生物量模拟的R2介于0.49～0.83，SNRMSE介于10.47%～27.96%，SRMSE分别介于0.38～0.69 m2/m2、4.30～6.28 cm 和789.00～1 134.00 kg/hm2，仍有部分指标处于中等精度以下。同化LAI 时模型对LAI 模拟的精度显著提高，R2介于0.76～0.96 之间，SRMSE介于0.22～0.48 m2/m2，SNRMSE介于7.66%～19.38%，各处理均达到高精度水平。同化LAI+SW 能显著改善生物量的模拟结果，R2可达到0.73～0.92，特别是T1～T3 处理，R2较无同化最大可提高112.5%。SRMSE介于384.00～705.00 kg/hm2，SNRMSE介于6.52%～14.02%，模拟结果均达到高精度以上。株高无明显规律，但2 a 中各处理同化LAI 和LAI+SW对株高的模拟均达到高精度或极高精度，R2为0.71～0.96，SRMSE为2.30～5.16 cm，SNRMSE为3.25%～10.84%。可见选择不同的观测变量，SWAP-IES 系统对春小麦生长指标模拟结果有所差异，需根据研究目的选择适宜的观测变量进行同化模拟。

表7 SWAP-IES 系统同化不同观测变量对春小麦生长指标模拟的R2、SRMSE 和SNRMSETable 7 R2，SNRMSE and SNRMSE values of simulating spring wheat growth index by using SWAP-IES system to assimilate different observation variables

2.3 同化不同观测变量对估产的影响

图4 为SWAP-IES 系统同化不同观测变量对2019 和2020 年T1～T5 处理春小麦估产的结果。无同化时模型对春小麦产量预测的精度较低，R2为0.45，SRMSE为1 259.34 kg/hm2，SNRMSE为28.38%，这是因为参数先验值及模型本身的不确定性等导致的。同化SW 改善了模型对春小麦产量的预测效果，R2提高到0.66，SRMSE减小了355.57 kg/hm2，SNRMSE减小了6.73%，但是改善效果有限，模拟结果仍为中精度。同化LAI 时，对产量模拟的精度进一步提高，R2为0.79，SRMSE为698.49 kg/hm2，SNRMSE为16.74 %，达到高精度水平。同化LAI+SW 时，模型对春小麦产量模拟的精度最高，R2为0.87，SRMSE为561.54 kg/hm2，SNRMSE为13.45%。无同化情景下估产的SRE介于10.89%～40.34%，各处理均比同化情景大。T1、T2 和T3 处理下2 a 均为同化LAI+SW 时估产的SRE最小，T1 两年分别为6.13%和8.38%，T2 分别为3.87%和4.69%，T3 分别为6.11%和7.4%。T4 和T5 处理同化LAI 时模型估产的SRE均在10%以内（4.05%～9.17%），小于其他情形。可见在SWAP-IES 同化模拟系统中融入LAI 观测数据后，其估产精度显著提高，继续添加SW，估产能力会进一步提高，特别在中度或重度水分胁迫处理下（T1～T3）。

图4 同化不同观测变量的估产效果（2019—2020 年）Fig.4 Yield estimation effect of assimilating different observation variable（2019—2020）

2.4 观测变量的观测时间对估产的影响

图5 为同化数据观测生育期对估产精度的影响。单个生育期观测数据中同化开花期的LAI+SW 对提高SWAPIES 系统估产精度的作用最大（R2从无同化的0.45 提高到0.74，SRMSE和SNRMSE分别从无同化的1 259.34 kg/hm2、28.38%降低到676.19 kg/hm2、17.11%），其次为同化拔节期和孕穗期的观测数据。同化多个生育期的观测数据可显著提高模型的估产精度，同化拔节期和开花期的LAI+SW观测数据时估产的R2=0.79，SRMSE=610.80 kg/hm2，SNRMSE=15.12%，同化拔节期、孕穗期和开花期3 个生育期的LAI+SW 观测数据时估产的R2=0.83，SRMSE=598.99 kg/hm2，SNRMSE=14.46%。说明SWAP-IES 系统在同化模拟过程中，应重点选择开花期、拔节期等关键生育期的观测数据进行同化模拟，以实现提高模拟精度、减少田间观测强度的目的。

图5 同化数据观测生育期对SWAP-IES 系统估产精度的影响Fig.5 Impacts of assimilating data of growth period on the yield estimation accuracy of the SWAP-IES system

3 讨 论

本文将SWAP 模型和IES 数据同化算法相结合，构建了SWAP-IES 同化模拟系统，以春小麦大田试验为基础，评价了观测变量的选择对准确模拟春小麦长势和产量的价值。作物生长及农田土壤环境各指标间相互影响，在作物生长模型中同化任意一个指标，均可提高与该指标相关的其他指标的模拟精度，特别是LAI 和生物量等关联性较强的指标[33]。本文中同化LAI 显著提高了SWAPIES 系统对春小麦LAI、生物量和产量的模拟精度，R2达到0.76～0.96，在无水分胁迫或轻度水分胁迫下改善效果更为显著，这和DONG 等[34]在SAFY 模型中同化LAI对加拿大东北部冬小麦模拟的结论一致。本文中T4 和T5 处理同化LAI 时估产的相对误差均在10%以内，说明小麦受轻度或不受水分胁迫时，LAI 就能反映大部分作物生长和农田环境信息，仅同化LAI 就可实现农田尺度的准确估产，这和LI 等[35-36]的研究结果一致。借助遥感信息将作物生长模型扩展到区域尺度，同化LAI 后在县域尺度估产的误差在14.8%以内[37-38]，说明在正常生长条件下，同化LAI 后在农田和区域尺度均能准确模拟小麦生长过程和产量形成。

在重度水分胁迫条件下，仅同化LAI 无法实现准确估产，本文T1 和T2 处理在仅同化LAI 时估产的相对误差为5.94%～15.32%，显著低于同化LAI+SW 的情景，这可能是因为小麦在水分胁迫下，作物会产生耐旱的保护机制，而模型无法准确描述这一生理变化，造成对耗水量模拟误差，导致模拟的产量偏小。同化LAI+SW 后，修正了作物生长发育和土壤水平衡过程，对土壤水分的模拟精度显著提高，进而提高了估产效果，说明所构建的SWAP-IES 系统中LAI 为春小麦准确估产的必须观测变量，在水分胁迫较重时还需要增加SW 作为观测变量，这与PAN 等[39-40]研究结果一致。除了土壤水分外，增加腾发量[41]等其他观测数据，同样能提高水分胁迫逆境等条件下模型对作物产量的模拟精度。作物不同生育期的生长和土壤水分信息对估产的价值不同[42]，在同化模拟系统中融入更多有效信息后，模拟精度逐步提高，但是过多的观测数据会增加观测和计算成本，因此本文评价了观测变量的观测时期对估产的影响，寻求合理的观测生育期，以提高同化模拟的效率。受作物品种、土壤、气象等空间异质性制约，目前作物生长模型同化系统在区域估产方面往往精度有限，如何提高遥感反演观测变量的准确度，科学选择观测变量和同化时序，以及模型从验证向预估转变等是未来研究的重点。

4 结 论

本文融合土壤-水-大气-植物模型（soil-water-atmosphere-plant，SWAP）和迭代集合平滑器算法（iterative ensemble smoother，IES）构建了SWAP-IES 作物同化模拟系统，利用宁夏引黄灌区春小麦灌水量试验多时序叶面积指数（leaf area index，LAI）和表层土壤水分（soil water content，SW）的观测数据，评价了同化不同观测变量在春小麦长势和产量模拟中的价值。主要结论如下：

1）同化LAI+SW 和SW 时，对土壤含水率模拟的R2从无同化的0.48 分别提高到0.78 和0.87，表明同化变量中增加SW 能显著改善模型对土壤水分的模拟效果，为不同水分处理下准确模拟春小麦生长过程和产量形成奠定基础。

2）同化LAI+SW 时估产效果最好，R2从无同化的0.45 提高到0.87。无同化时模型估产的误差最大，相对误差介于10.89%～40.34%，T1～T3 处理同化LAI+SW 时估产效果最好，相对误差为3.87%～8.38%，T4～T5 处理同化LAI 时估产效果最好，相对误差为4.05%～9.17%。

3）单个生育期中开花期的观测数据对提高SWAPIES 系统估产精度的作用最大，同化开花期LAI+SW 时估产的R2可从无同化的0.45 提高到0.74，其次为拔节期。同时同化拔节期和开花期2 个生育期LAI+SW 观测数据时估产的R2可达到0.79。