基于全局敏感性分析的WOFOST模型参数优化

谢松涯，张宝忠

(1.中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室，北京 100038；2.国家节水灌溉北京工程技术研究中心，北京 100048)

0 引 言

玉米作为中国的3大粮食作物之一，在全国分布广泛。根据2016年发布的《中国统计年鉴年》，我国玉米播种面积占全国农作物播种面积的22.06%，占粮食作物播种面积的32.52%，总产量为21 955.2 万t，占粮食总产量的35.63%[1]。

近年来作物生长模型与田间试验的结合，在水肥管理、作物种植体系优化方面得到了极大发展[2,3]。世界范围内广泛运用的作物生长模型主要有，澳大利亚的APSIM模型，美国的DSSAT、EPIC模型和荷兰的WOFOST、SWAP模型等。美国的作物模型实用性和可操作性较强，模型的各个模块由不同单位研发，从不同方面模拟环境因子对作物生长影响，具有模拟思路多的优点，同时也存在模块基本假设不同而导致模拟精度不同的现象；澳大利亚的作物模型以土壤特性变化为模拟中心，可通过选用不同模块对比不同方法的优劣；荷兰的作物模型具有机理性、通用性强的优点，提供了作物产量、土壤水分与养分的多种模拟思路[4-6]，据此本文选择荷兰作物生长模型WOFOST为研究对象。

作物模型在运用时需要对模型参数进行率定和适用性验证，由于作物模型具有参数众多的特点，难以对所有参数进行调整率定，据此利用敏感性分析的方法对模型的参数进行分析，进而筛选出需要率定的参数[7,8]。敏感性分析方法主要包括局部和全局敏感性分析，局部敏感性分析研究单个参数变化对模型结果的影响，全局敏感性则分析多个参数变化和参数之间相互作用对模型结果的影响[9,10]，全局敏感性分析方法主要包括Morris法[11]、多元回归法[12]、基于方差分解理论的Sobol方法[13]、傅里叶幅度灵敏度检验法(FAST)[14]和扩展的傅里叶幅度灵敏度检验法(EFAST)[7]，这些方法在水文模型和作物生长模型的敏感性分析中均有广泛的运用[8,15,16]。EFAST方法是Seltelli等人结合Sobal法的优点对FAST法改进的全局敏感性分析法，具有需样本数少、计算稳定和计算效率高的特点[8]，

本文将运用WOFOST模型对北京地区夏玉米的生长过程和产量进行模拟，运用扩展的傅里叶灵敏度检验法(EFAST)对WOFOST的作物参数进行全局敏感性分析，分析作物参数对于夏玉米的叶面积指数、干物质量和产量3个输出变量的影响，利用调参软件PEST对敏感性参数优化，以进一步优化WOFOST模型在北京地区夏玉米整个生长过程的动态模拟。

1 模型敏感性分析与优化方法

1.1 WOFOST模型

WOFOST模型是荷兰瓦赫宁农业大学和世界粮食研究中心共同开发研制的作物生长模拟模型。近几十年来，WOFOST模型已经在学术和农业界得到了广泛的运用[17,18]。WOFOST模型以气象数据为驱动，通过调整土壤、管理和作物参数来控制和调整作物的生长过程。模拟内容主要包括呼吸作用、同化作用、蒸腾作用、干物质的分配等。模型气象数据主要包括每日最高气温、最低气温、蒸汽压、降雨量、平均风速、辐射量等。土壤数据主要包括土壤饱和含水率、土壤导水率、田间持水量和凋萎系数等。作物管理数据来源于不同年份的实际田间管理数据，主要包括播种时间、灌水量等。作物参数包括不同生育阶段所需要的积温、不同生育阶段二氧化碳的同化效率、不同器官的同化物转化效率等[19,20]。

1.2 敏感性分析方法

本研究运用的敏感性分析方法为扩展傅里叶幅度检验法(EFAST)，EFAST法是Saltelli等人结合Sobol法和傅立叶幅度敏感性检验法(FAST)的优点所提出的全局敏感性分析方法[7]。EFAST法通过分解模型方差，进而求出各参数和参数之间的相互作用对总方差的贡献量，获得各参数的一阶和多阶敏感性指数。模型总方差可分解为：

(1)

式中：V为模型的总方差；Vi为参数xi通过参数xi作用所贡献的方差；Vij为参数xi通过参数xi、xj所贡献的方差；V12…k为参数xi通过参数x12…k所贡献的方差。

定义参数及参数相互作用的方差与总方差的比值为敏感性指数，反映参数xi对模型输出总方差的直接贡献率，即参数xi的一阶敏感性指数Si可表示如下：

(2)

同理，参数xi的二阶、三阶敏感性指数可表示为：

(3)

(4)

参数xi的总敏感性指数即为各阶敏感性指数之和，表示如下：

ST,i=Si+Sij+Sijm+S1…i…k

(5)

总敏感性指数反映了参数直接贡献率和通过参数之间的相互作用间接对模型输出对总方差的影响。由于WOFOST模型中包含多个参数之间的相互作用，本研究选用全局敏感性分析方法来分析WOFOST模型中参数对输出结果的影响，EFAST方法通过对模型输出方差的分解，可定量地获得参数的一阶和总敏感指数。这就使得EFAST方法可以同时检验多个参数的变化对WOFOST模型结果的影响，并且可分析每一个参数变化对模型结果的直接和间接影响。

敏感性分析借助专业软件Simlab完成。选择49个作物参数作为作物文件的输入值，在夏玉米参数的默认值的基础上取±15%的范围，参数在此范围内均匀分布。采用蒙特卡罗(Monte Carlo)方法，对参数进行7 000次随机取样，同时以7 000次参数随机取样作为WOFOST模型的作物参数文件，对夏玉米产量、干物质量、叶面积指数等结果进行批量输出。最后运用EFAST方法对输出结果进行全局敏感性分析。

1.3 PEST模型优化参数及验证方案

PEST为参数优化软件，利用外部模型的输出文件和软件内部输入文件进行数据交换连接，从而达到参数优化的目的。PEST软件的核心是求解目标函数的最小值，软件运用了Gauss-Marquardt-Levevberg(GML)算法求解，该算法是基于牛顿法和梯度下降法的一种非线性优化方法，能够在多维的参数空间内优化模型输入参数，迭代逐步逼近目标函数最小值，具有快速收敛，运行次数少的特点[21,22]。目标函数方程为外部模型多个输出变量的计算值与实际观测值的带权重最小二乘差异函数Ψ，其公式为：

(6)

式中：F为一系列作物参数；M(ti)为在时间i的观测值；S(F，ti)为在模型时间i的模拟值；wi为观测值的权重系数；n为观测值的个数。

本研究包括三类观测对比值，分别为叶面积指数、干物质量和作物产量，根据目标函数代表模型输出值的种类和模型模拟的要求，设计不同的权重系数wi，使叶面积指数、干物质量和作物产量模拟误差的比例控制在一定的范围内。PEST工具极大缩短了调参时间，同时降低了人为主观因素的影响，实现了自动化调参功能。

模型的验证采用观测值和模拟值的相对误差(RE，relative error)、一致性指数(index of agreement，d)[]和标准均方根误差(normalized root mean square error，nRMSE)[]三个指标来评价，它们均可以反映模型值与实测值之间的相对差异程度，同时也是无量纲统计量，可以进行不同变量之间的比较。当RE越靠近0、d越接近1说明模拟精度越高, 当nRMSE<10%时，此时模型效果为优；当10%≤nRMSE<20%时，模拟效果为良；当20%≤nRMSE<30%时，模拟效果为中等；当nRMSE>30%时，模拟效果为差。相对误差RE、一致性指数d和标准均方根误差nRMSE的计算式分别为：

(7)

(8)

(9)

式中：Si为第i个模拟值；Mi为第i个观测值；Sm、Mm为分别为模拟值和观测值的均值；n为样本数。

2 试验区概况与指标获取

2.1 试验区概况

本研究的试验数据来源于中国水利水电科学研究院大兴试验基地，该试验基地位于北京大兴南30 km处，地处 39°37.25′(N)，116°25.51′(E)。试验区位于华北平原地区，属于半干旱大陆性季风气候，多年平均气温是12 ℃，平均降雨量为556 mm，年均日照数为2 620 h。降雨量小于蒸发量，最小蒸发量为980 mm，最大蒸发量为1 100 mm，该地区最主要的作物种植模式是冬小麦-夏玉米连作模式，在正常年份冬小麦需补充灌溉，以保证作物对水分的需求，夏玉米在平水年以上，生长期内通常不需要补充灌溉。

2.2 监测指标与获取方法

本文所选择的田间试验对象为夏玉米(纪元16)。夏玉米为2016年6月16日播种，2016年9月27日收获，播种密度约为4.54 万株/hm2。试验的气象数据来自于试验站的气象站，试验区的土壤主要为沙壤土。土壤主要性质见表1。

研究涉及到的观测数据主要为夏玉米的产量、干物质量和叶面积指数，夏玉米的叶面积指数每7～10 d左右测一次。玉米的干物质量每10 d左右测定一次。

3 结果与讨论

3.1 全局敏感性分析

3.1.1 产量的参数敏感性分析

对于产量，一阶敏感性指数Si>0.1的前6位的参数依次为二氧化碳同化率在12 ℃下的矫正因子(TMNFTB12)、最大光合速率在30 ℃下的校正因子(TMPFTB30)、开花到成熟之间的积温(TSUM2)、贮存器官生长同化物转化效率(CVO)、生育期为2(DVS=2)时的可见光的消散系数(KDIFTB2)、温度为0 ℃时的单叶片初始量子效率(EFFTB0)，其余参数一阶敏感性指数均小于0.1；全局敏感性指数(ST,i>0.10)前六位参数的分布与一阶敏感性指数的分布一致，其值分别为0.434、0.345、0.294、0.258、0.200、0.144。萌芽温度的下限(TBASEM)、温度为40 ℃时单片叶初始量子效率(EFFTB40)、生育期为2时的老化矫正因子(RFSETB2)的值均超过0.1，其值分别为0.138、0.119、0.114。其余参数全局敏感性指数均小于0.1。具体敏感性参数分布见图1。

3.1.2 干物质量随时间变化的参数敏感性

干物质量(TAGP)在整个生育期不断的变化，分析参数的敏感性随时间变化将有助于分析各个参数在不同生育期的作用，选取不同生育期一阶敏感性指数Si>0.05、全局敏感性分ST,i>0.1的参数，在此基础上综合取全局敏感性指数较大的前10个作物参数作为分析对象，分别为：CO2同化率在12 ℃下的矫正因子(TMNFTB12)、在35 ℃时叶面积的生命周期(SPAN)、生育期为0时的比叶面积(SLATB0)、生育期为1.1时的根干物质的分配系数(FRTB11)、30 ℃下积温的增长量(DTSMTB30)、茎同化物转换效率(CVS)、根的同化物转换效率(CVR)、储存器官的同化物转换效率(CVO)、叶片的同化物转换效率(CVL)和生育期为1时的单叶片CO2同化效率(AMAXTB100)，干物质量随时间变化的一阶以及全局敏感性参数变化如图2所示。

图1 产量的全局和一阶敏感性分析结果Fig.1 First order and total sensitivity analysis result of summer maize yield

图2 干物质量随时间变化的全局和一阶敏感性分析结果Fig.2 First order and total sensitivity analysis results of summer maize TAGP

可以看出，选取参数的一阶敏感性指数和全局敏感系指数的变化趋势基本趋于一致，但全局敏感性分析，不仅包含某一参数对模型结果的贡献率，同时也包含不同参数之间的交互作用。CO2同化率在12℃下的矫正因子(TMNFTB12)在夏玉米的整个生长周期，都表现出很大的敏感性，在出苗后的73 d内一阶和全局敏感性指数都保持稳定的增长，但在73 d后至整个生育期的结束，TMNFTB12的敏感性指数出现大幅度的下降。在夏玉米生长前期，SLATB0、CVR、CVS的一阶、全局敏感系指数都有较高的数值，但在整个生育期都有不同程度的下降，全局敏感系指数在后期有所上升，但上升幅度较小，均未超过0.1。在夏玉米生长前期，CVO、SPAN的一阶、全局敏感系指数数值较小，但在生育的中后期，CVO的一阶以及全局敏感系指数逐渐上升至较高的数值，SPAN的全局敏感系指数也有较大幅度的上升。

3.1.3 叶面积指数随时间变化的参数敏感性

选取不同生育期一阶敏感性指数Si>0.05、全局敏感性分析ST,i>0.1的参数，在此基础上取全局敏感性指数较大的前10个作物参数作为分析对象，分别为：根的同化物转换效率(CVR)、生育期为0时的比叶面积(SLATB0)、CO2同化率在12 ℃下的矫正因子(TMNFTB12)、生育期为1时的单叶片CO2同化效率(AMAXTB100)、最大光合速率在0 ℃时的矫正因子(TMPFTB0)、生育期为1时的比叶面积(SLATB1)、茎同化物转换效率(CVS)、水分限制时叶片的相对死亡率(PERDL)、出苗至开花有效积温(TSUM1)和在35 ℃时叶面积的生命周期(SPAN)，叶面积指数随时间变化的一阶以及全局敏感性参数变化如图3所示。

图3 叶面积指数随时间变化的全局和一阶敏感性分析结果Fig.3 First order and total sensitivity analysis results of summer maize LAI

可以看出CVR、SLATB0、TMNFTB12的一阶以及全局敏感性指数在前期有较高的数值，但随着时间的推进，一阶以及全局敏感性指数降低至0.1以下；SPAN及SLATB1参数在生育前期，一阶以及全局敏感性指数均较小，但在出苗72 d后，两个参数的敏感性指数均出现较大幅度的增长。

3.2 参数优化验证结果

3.2.1 优化参数的选择

综合选取CO2同化率在12 ℃下的矫正因子(TMNFTB12)、在35 ℃时叶面积的生命周期(SPAN)、生育期为0时的比叶面积(SLATB0)、生育期为1时的比叶面积(SLATB1)、茎同化物转换效率(CVS)、根的同化物转换效率(CVR)、储存器官的同化物转换效率(CVO)、最大光合速率在30 ℃下的校正因子(TMPFTB30)、生育期为2(DVS=2)时的可见光的消散系数(KDIFTB2)、温度为0 ℃时的单叶片初始量子效率(EFFTB0)等10个对夏玉米产量、干物质量、叶面积指数敏感的参数，运用PEST软件，对2016年的夏玉米数据进行调参，并利用2015-2016年的数据进行模拟验证。

3.2.2 参数优化验证

夏玉米干物质量、叶面积指数的优化及验证结果如图4，由图4可知，2016年模拟干物质量(TAGP)一致性指数d=0.995，接近于1，标准均方根误差nRMSE=12%，在10%～20%之间，模拟结果良好。利用2015年干物质量(TAGP)的实测数据验证调整参数，得到一致性指数d=0.993、标准均方根误差nRMSE=12%，拟合结果良好。2016年模拟叶面积指数(LAI)一致性指数d=0.948，标准均方根误差nRMSE=19%，2015年叶面积指数(LAI)拟合结果一致性指数d=0.919、标准均方根误差nRMSE=18%。2016年模拟产量(WSO)的相对误差RE=0.15%，2015年验证产量的相对误差RE=2.9%。

从模型验证角度看，对2015、2016年的模拟结果与观测值对比，相对误差RE小于5%，一致性指数d均大于0.9，标准均方根误差nRMSE均小于20%。综合上述得出，WOFOST模型调参后对于北京大兴地区的夏玉米的适应性良好。

图4 干物质量、叶面积指数的优化验证结果Fig.4 Amount of dry matter、leaf area index optimization results

4 结 语

本研究运用中国水利水电科学研究院大兴试验基地2016年夏玉米数据，利用扩展的傅里叶幅度检验法(EFAST)对WOFOST模型的49个作物参数进行敏感性分析，选出对产量、干物质量、叶面积指数敏感性较高的10个参数进行调参，并利用2015年的试验数据进行模拟验证。得到的主要结论如下。

(1)对产量敏感的参数主要有：二氧化碳同化率在12 ℃下的矫正因子(TMNFTB12)、最大光合速率在30 ℃下的校正因子(TMPFTB30)、开花到成熟之间的积温(TSUM2)、贮存器官生长同化物转化效率(CVO)、生育期为2(DVS=2)时的可见光的消散系数(KDIFTB2)。

(2)对干物质量敏感的参数主要有：CO2同化率在12 ℃下的矫正因子(TMNFTB12)、在35 ℃时叶面积的生命周期(SPAN)、生育期为0时的比叶面积(SLATB0)、生育期为1.1时的根干物质的分配系数(FRTB11)、30 ℃下积温的增长量(DTSMTB30)。

(3)对叶面积指数敏感的参数主要有：根的同化物转换效率(CVR)、生育期为0时的比叶面积(SLATB0)、CO2同化率在12 ℃下的矫正因子(TMNFTB12)、生育期为1时的单叶片CO2同化效率(AMAXTB100)、最大光合速率在0 ℃时的矫正因子(TMPFTB0)。

(4)2016年调参后的夏玉米模拟值与试验观测值对比，产量的相对误差RE=0.15%，干物质量的一致性指数d=0.995，标准均方根误差nRMSE=12%，叶面积指数的一致性指数d=0.948，标准均方根误差nRMSE=19%。利用2015年数据对参数进行验证，夏玉米的模拟值与试验观测值对比，产量的相对误差RE=2.9%，干物质量的一致性指数d=0.993、标准均方根误差nRMSE=12%，叶面积指数的一致性指数d=0.919、标准均方根误差nRMSE=18%。表明优化后的参数对于北京地区夏玉米模拟效果良好。

敏感性分析方法及调参软件在作物生长模型中的应用减少了人为主观因素的影响，极大的缩短调参时间。

□