基于全局敏感性分析的WOFOST模型參數(shù)優(yōu)化

謝松涯，張寶忠

(1.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實驗室，北京 100038；2.國家節(jié)水灌溉北京工程技術(shù)研究中心，北京 100048)

0 引 言

玉米作為中國的3大糧食作物之一，在全國分布廣泛。根據(jù)2016年發(fā)布的《中國統(tǒng)計年鑒年》，我國玉米播種面積占全國農(nóng)作物播種面積的22.06%，占糧食作物播種面積的32.52%，總產(chǎn)量為21 955.2 萬t，占糧食總產(chǎn)量的35.63%[1]。

近年來作物生長模型與田間試驗的結(jié)合，在水肥管理、作物種植體系優(yōu)化方面得到了極大發(fā)展[2,3]。世界范圍內(nèi)廣泛運(yùn)用的作物生長模型主要有，澳大利亞的APSIM模型，美國的DSSAT、EPIC模型和荷蘭的WOFOST、SWAP模型等。美國的作物模型實用性和可操作性較強(qiáng)，模型的各個模塊由不同單位研發(fā)，從不同方面模擬環(huán)境因子對作物生長影響，具有模擬思路多的優(yōu)點(diǎn)，同時也存在模塊基本假設(shè)不同而導(dǎo)致模擬精度不同的現(xiàn)象；澳大利亞的作物模型以土壤特性變化為模擬中心，可通過選用不同模塊對比不同方法的優(yōu)劣；荷蘭的作物模型具有機(jī)理性、通用性強(qiáng)的優(yōu)點(diǎn)，提供了作物產(chǎn)量、土壤水分與養(yǎng)分的多種模擬思路[4-6]，據(jù)此本文選擇荷蘭作物生長模型WOFOST為研究對象。

作物模型在運(yùn)用時需要對模型參數(shù)進(jìn)行率定和適用性驗證，由于作物模型具有參數(shù)眾多的特點(diǎn)，難以對所有參數(shù)進(jìn)行調(diào)整率定，據(jù)此利用敏感性分析的方法對模型的參數(shù)進(jìn)行分析，進(jìn)而篩選出需要率定的參數(shù)[7,8]。敏感性分析方法主要包括局部和全局敏感性分析，局部敏感性分析研究單個參數(shù)變化對模型結(jié)果的影響，全局敏感性則分析多個參數(shù)變化和參數(shù)之間相互作用對模型結(jié)果的影響[9,10]，全局敏感性分析方法主要包括Morris法[11]、多元回歸法[12]、基于方差分解理論的Sobol方法[13]、傅里葉幅度靈敏度檢驗法(FAST)[14]和擴(kuò)展的傅里葉幅度靈敏度檢驗法(EFAST)[7]，這些方法在水文模型和作物生長模型的敏感性分析中均有廣泛的運(yùn)用[8,15,16]。EFAST方法是Seltelli等人結(jié)合Sobal法的優(yōu)點(diǎn)對FAST法改進(jìn)的全局敏感性分析法，具有需樣本數(shù)少、計算穩(wěn)定和計算效率高的特點(diǎn)[8]，

本文將運(yùn)用WOFOST模型對北京地區(qū)夏玉米的生長過程和產(chǎn)量進(jìn)行模擬，運(yùn)用擴(kuò)展的傅里葉靈敏度檢驗法(EFAST)對WOFOST的作物參數(shù)進(jìn)行全局敏感性分析，分析作物參數(shù)對于夏玉米的葉面積指數(shù)、干物質(zhì)量和產(chǎn)量3個輸出變量的影響，利用調(diào)參軟件PEST對敏感性參數(shù)優(yōu)化，以進(jìn)一步優(yōu)化WOFOST模型在北京地區(qū)夏玉米整個生長過程的動態(tài)模擬。

1 模型敏感性分析與優(yōu)化方法

1.1 WOFOST模型

WOFOST模型是荷蘭瓦赫寧農(nóng)業(yè)大學(xué)和世界糧食研究中心共同開發(fā)研制的作物生長模擬模型。近幾十年來，WOFOST模型已經(jīng)在學(xué)術(shù)和農(nóng)業(yè)界得到了廣泛的運(yùn)用[17,18]。WOFOST模型以氣象數(shù)據(jù)為驅(qū)動，通過調(diào)整土壤、管理和作物參數(shù)來控制和調(diào)整作物的生長過程。模擬內(nèi)容主要包括呼吸作用、同化作用、蒸騰作用、干物質(zhì)的分配等。模型氣象數(shù)據(jù)主要包括每日最高氣溫、最低氣溫、蒸汽壓、降雨量、平均風(fēng)速、輻射量等。土壤數(shù)據(jù)主要包括土壤飽和含水率、土壤導(dǎo)水率、田間持水量和凋萎系數(shù)等。作物管理數(shù)據(jù)來源于不同年份的實際田間管理數(shù)據(jù)，主要包括播種時間、灌水量等。作物參數(shù)包括不同生育階段所需要的積溫、不同生育階段二氧化碳的同化效率、不同器官的同化物轉(zhuǎn)化效率等[19,20]。

1.2 敏感性分析方法

本研究運(yùn)用的敏感性分析方法為擴(kuò)展傅里葉幅度檢驗法(EFAST)，EFAST法是Saltelli等人結(jié)合Sobol法和傅立葉幅度敏感性檢驗法(FAST)的優(yōu)點(diǎn)所提出的全局敏感性分析方法[7]。EFAST法通過分解模型方差，進(jìn)而求出各參數(shù)和參數(shù)之間的相互作用對總方差的貢獻(xiàn)量，獲得各參數(shù)的一階和多階敏感性指數(shù)。模型總方差可分解為：

(1)

式中：V為模型的總方差；Vi為參數(shù)xi通過參數(shù)xi作用所貢獻(xiàn)的方差；Vij為參數(shù)xi通過參數(shù)xi、xj所貢獻(xiàn)的方差；V12…k為參數(shù)xi通過參數(shù)x12…k所貢獻(xiàn)的方差。

定義參數(shù)及參數(shù)相互作用的方差與總方差的比值為敏感性指數(shù)，反映參數(shù)xi對模型輸出總方差的直接貢獻(xiàn)率，即參數(shù)xi的一階敏感性指數(shù)Si可表示如下：

(2)

同理，參數(shù)xi的二階、三階敏感性指數(shù)可表示為：

(3)

(4)

參數(shù)xi的總敏感性指數(shù)即為各階敏感性指數(shù)之和，表示如下：

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

(5)

總敏感性指數(shù)反映了參數(shù)直接貢獻(xiàn)率和通過參數(shù)之間的相互作用間接對模型輸出對總方差的影響。由于WOFOST模型中包含多個參數(shù)之間的相互作用，本研究選用全局敏感性分析方法來分析WOFOST模型中參數(shù)對輸出結(jié)果的影響，EFAST方法通過對模型輸出方差的分解，可定量地獲得參數(shù)的一階和總敏感指數(shù)。這就使得EFAST方法可以同時檢驗多個參數(shù)的變化對WOFOST模型結(jié)果的影響，并且可分析每一個參數(shù)變化對模型結(jié)果的直接和間接影響。

敏感性分析借助專業(yè)軟件Simlab完成。選擇49個作物參數(shù)作為作物文件的輸入值，在夏玉米參數(shù)的默認(rèn)值的基礎(chǔ)上取±15%的范圍，參數(shù)在此范圍內(nèi)均勻分布。采用蒙特卡羅(Monte Carlo)方法，對參數(shù)進(jìn)行7 000次隨機(jī)取樣，同時以7 000次參數(shù)隨機(jī)取樣作為WOFOST模型的作物參數(shù)文件，對夏玉米產(chǎn)量、干物質(zhì)量、葉面積指數(shù)等結(jié)果進(jìn)行批量輸出。最后運(yùn)用EFAST方法對輸出結(jié)果進(jìn)行全局敏感性分析。

1.3 PEST模型優(yōu)化參數(shù)及驗證方案

PEST為參數(shù)優(yōu)化軟件，利用外部模型的輸出文件和軟件內(nèi)部輸入文件進(jìn)行數(shù)據(jù)交換連接，從而達(dá)到參數(shù)優(yōu)化的目的。PEST軟件的核心是求解目標(biāo)函數(shù)的最小值，軟件運(yùn)用了Gauss-Marquardt-Levevberg(GML)算法求解，該算法是基于牛頓法和梯度下降法的一種非線性優(yōu)化方法，能夠在多維的參數(shù)空間內(nèi)優(yōu)化模型輸入?yún)?shù)，迭代逐步逼近目標(biāo)函數(shù)最小值，具有快速收斂，運(yùn)行次數(shù)少的特點(diǎn)[21,22]。目標(biāo)函數(shù)方程為外部模型多個輸出變量的計算值與實際觀測值的帶權(quán)重最小二乘差異函數(shù)Ψ，其公式為：

(6)

式中：F為一系列作物參數(shù)；M(ti)為在時間i的觀測值；S(F，ti)為在模型時間i的模擬值；wi為觀測值的權(quán)重系數(shù)；n為觀測值的個數(shù)。

本研究包括三類觀測對比值，分別為葉面積指數(shù)、干物質(zhì)量和作物產(chǎn)量，根據(jù)目標(biāo)函數(shù)代表模型輸出值的種類和模型模擬的要求，設(shè)計不同的權(quán)重系數(shù)wi，使葉面積指數(shù)、干物質(zhì)量和作物產(chǎn)量模擬誤差的比例控制在一定的范圍內(nèi)。PEST工具極大縮短了調(diào)參時間，同時降低了人為主觀因素的影響，實現(xiàn)了自動化調(diào)參功能。

模型的驗證采用觀測值和模擬值的相對誤差(RE，relative error)、一致性指數(shù)(index of agreement，d)[]和標(biāo)準(zhǔn)均方根誤差(normalized root mean square error，nRMSE)[]三個指標(biāo)來評價，它們均可以反映模型值與實測值之間的相對差異程度，同時也是無量綱統(tǒng)計量，可以進(jìn)行不同變量之間的比較。當(dāng)RE越靠近0、d越接近1說明模擬精度越高, 當(dāng)nRMSE<10%時，此時模型效果為優(yōu)；當(dāng)10%≤nRMSE<20%時，模擬效果為良；當(dāng)20%≤nRMSE<30%時，模擬效果為中等；當(dāng)nRMSE>30%時，模擬效果為差。相對誤差RE、一致性指數(shù)d和標(biāo)準(zhǔn)均方根誤差nRMSE的計算式分別為：

(7)

(8)

(9)

式中：Si為第i個模擬值；Mi為第i個觀測值；Sm、Mm為分別為模擬值和觀測值的均值；n為樣本數(shù)。

2 試驗區(qū)概況與指標(biāo)獲取

2.1 試驗區(qū)概況

本研究的試驗數(shù)據(jù)來源于中國水利水電科學(xué)研究院大興試驗基地，該試驗基地位于北京大興南30 km處，地處 39°37.25′(N)，116°25.51′(E)。試驗區(qū)位于華北平原地區(qū)，屬于半干旱大陸性季風(fēng)氣候，多年平均氣溫是12 ℃，平均降雨量為556 mm，年均日照數(shù)為2 620 h。降雨量小于蒸發(fā)量，最小蒸發(fā)量為980 mm，最大蒸發(fā)量為1 100 mm，該地區(qū)最主要的作物種植模式是冬小麥-夏玉米連作模式，在正常年份冬小麥需補(bǔ)充灌溉，以保證作物對水分的需求，夏玉米在平水年以上，生長期內(nèi)通常不需要補(bǔ)充灌溉。

2.2 監(jiān)測指標(biāo)與獲取方法

本文所選擇的田間試驗對象為夏玉米(紀(jì)元16)。夏玉米為2016年6月16日播種，2016年9月27日收獲，播種密度約為4.54 萬株/hm2。試驗的氣象數(shù)據(jù)來自于試驗站的氣象站，試驗區(qū)的土壤主要為沙壤土。土壤主要性質(zhì)見表1。

研究涉及到的觀測數(shù)據(jù)主要為夏玉米的產(chǎn)量、干物質(zhì)量和葉面積指數(shù)，夏玉米的葉面積指數(shù)每7～10 d左右測一次。玉米的干物質(zhì)量每10 d左右測定一次。

3 結(jié)果與討論

3.1 全局敏感性分析

3.1.1 產(chǎn)量的參數(shù)敏感性分析

對于產(chǎn)量，一階敏感性指數(shù)Si>0.1的前6位的參數(shù)依次為二氧化碳同化率在12 ℃下的矯正因子(TMNFTB12)、最大光合速率在30 ℃下的校正因子(TMPFTB30)、開花到成熟之間的積溫(TSUM2)、貯存器官生長同化物轉(zhuǎn)化效率(CVO)、生育期為2(DVS=2)時的可見光的消散系數(shù)(KDIFTB2)、溫度為0 ℃時的單葉片初始量子效率(EFFTB0)，其余參數(shù)一階敏感性指數(shù)均小于0.1；全局敏感性指數(shù)(ST,i>0.10)前六位參數(shù)的分布與一階敏感性指數(shù)的分布一致，其值分別為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。其余參數(shù)全局敏感性指數(shù)均小于0.1。具體敏感性參數(shù)分布見圖1。

3.1.2 干物質(zhì)量隨時間變化的參數(shù)敏感性

干物質(zhì)量(TAGP)在整個生育期不斷的變化，分析參數(shù)的敏感性隨時間變化將有助于分析各個參數(shù)在不同生育期的作用，選取不同生育期一階敏感性指數(shù)Si>0.05、全局敏感性分ST,i>0.1的參數(shù)，在此基礎(chǔ)上綜合取全局敏感性指數(shù)較大的前10個作物參數(shù)作為分析對象，分別為：CO2同化率在12 ℃下的矯正因子(TMNFTB12)、在35 ℃時葉面積的生命周期(SPAN)、生育期為0時的比葉面積(SLATB0)、生育期為1.1時的根干物質(zhì)的分配系數(shù)(FRTB11)、30 ℃下積溫的增長量(DTSMTB30)、莖同化物轉(zhuǎn)換效率(CVS)、根的同化物轉(zhuǎn)換效率(CVR)、儲存器官的同化物轉(zhuǎn)換效率(CVO)、葉片的同化物轉(zhuǎn)換效率(CVL)和生育期為1時的單葉片CO2同化效率(AMAXTB100)，干物質(zhì)量隨時間變化的一階以及全局敏感性參數(shù)變化如圖2所示。

圖1 產(chǎn)量的全局和一階敏感性分析結(jié)果Fig.1 First order and total sensitivity analysis result of summer maize yield

圖2 干物質(zhì)量隨時間變化的全局和一階敏感性分析結(jié)果Fig.2 First order and total sensitivity analysis results of summer maize TAGP

可以看出，選取參數(shù)的一階敏感性指數(shù)和全局敏感系指數(shù)的變化趨勢基本趨于一致，但全局敏感性分析，不僅包含某一參數(shù)對模型結(jié)果的貢獻(xiàn)率，同時也包含不同參數(shù)之間的交互作用。CO2同化率在12℃下的矯正因子(TMNFTB12)在夏玉米的整個生長周期，都表現(xiàn)出很大的敏感性，在出苗后的73 d內(nèi)一階和全局敏感性指數(shù)都保持穩(wěn)定的增長，但在73 d后至整個生育期的結(jié)束，TMNFTB12的敏感性指數(shù)出現(xiàn)大幅度的下降。在夏玉米生長前期，SLATB0、CVR、CVS的一階、全局敏感系指數(shù)都有較高的數(shù)值，但在整個生育期都有不同程度的下降，全局敏感系指數(shù)在后期有所上升，但上升幅度較小，均未超過0.1。在夏玉米生長前期，CVO、SPAN的一階、全局敏感系指數(shù)數(shù)值較小，但在生育的中后期，CVO的一階以及全局敏感系指數(shù)逐漸上升至較高的數(shù)值，SPAN的全局敏感系指數(shù)也有較大幅度的上升。

3.1.3 葉面積指數(shù)隨時間變化的參數(shù)敏感性

選取不同生育期一階敏感性指數(shù)Si>0.05、全局敏感性分析ST,i>0.1的參數(shù)，在此基礎(chǔ)上取全局敏感性指數(shù)較大的前10個作物參數(shù)作為分析對象，分別為：根的同化物轉(zhuǎn)換效率(CVR)、生育期為0時的比葉面積(SLATB0)、CO2同化率在12 ℃下的矯正因子(TMNFTB12)、生育期為1時的單葉片CO2同化效率(AMAXTB100)、最大光合速率在0 ℃時的矯正因子(TMPFTB0)、生育期為1時的比葉面積(SLATB1)、莖同化物轉(zhuǎn)換效率(CVS)、水分限制時葉片的相對死亡率(PERDL)、出苗至開花有效積溫(TSUM1)和在35 ℃時葉面積的生命周期(SPAN)，葉面積指數(shù)隨時間變化的一階以及全局敏感性參數(shù)變化如圖3所示。

圖3 葉面積指數(shù)隨時間變化的全局和一階敏感性分析結(jié)果Fig.3 First order and total sensitivity analysis results of summer maize LAI

可以看出CVR、SLATB0、TMNFTB12的一階以及全局敏感性指數(shù)在前期有較高的數(shù)值，但隨著時間的推進(jìn)，一階以及全局敏感性指數(shù)降低至0.1以下；SPAN及SLATB1參數(shù)在生育前期，一階以及全局敏感性指數(shù)均較小，但在出苗72 d后，兩個參數(shù)的敏感性指數(shù)均出現(xiàn)較大幅度的增長。

3.2 參數(shù)優(yōu)化驗證結(jié)果

3.2.1 優(yōu)化參數(shù)的選擇

綜合選取CO2同化率在12 ℃下的矯正因子(TMNFTB12)、在35 ℃時葉面積的生命周期(SPAN)、生育期為0時的比葉面積(SLATB0)、生育期為1時的比葉面積(SLATB1)、莖同化物轉(zhuǎn)換效率(CVS)、根的同化物轉(zhuǎn)換效率(CVR)、儲存器官的同化物轉(zhuǎn)換效率(CVO)、最大光合速率在30 ℃下的校正因子(TMPFTB30)、生育期為2(DVS=2)時的可見光的消散系數(shù)(KDIFTB2)、溫度為0 ℃時的單葉片初始量子效率(EFFTB0)等10個對夏玉米產(chǎn)量、干物質(zhì)量、葉面積指數(shù)敏感的參數(shù)，運(yùn)用PEST軟件，對2016年的夏玉米數(shù)據(jù)進(jìn)行調(diào)參，并利用2015-2016年的數(shù)據(jù)進(jìn)行模擬驗證。

3.2.2 參數(shù)優(yōu)化驗證

夏玉米干物質(zhì)量、葉面積指數(shù)的優(yōu)化及驗證結(jié)果如圖4，由圖4可知，2016年模擬干物質(zhì)量(TAGP)一致性指數(shù)d=0.995，接近于1，標(biāo)準(zhǔn)均方根誤差nRMSE=12%，在10%～20%之間，模擬結(jié)果良好。利用2015年干物質(zhì)量(TAGP)的實測數(shù)據(jù)驗證調(diào)整參數(shù)，得到一致性指數(shù)d=0.993、標(biāo)準(zhǔn)均方根誤差nRMSE=12%，擬合結(jié)果良好。2016年模擬葉面積指數(shù)(LAI)一致性指數(shù)d=0.948，標(biāo)準(zhǔn)均方根誤差nRMSE=19%，2015年葉面積指數(shù)(LAI)擬合結(jié)果一致性指數(shù)d=0.919、標(biāo)準(zhǔn)均方根誤差nRMSE=18%。2016年模擬產(chǎn)量(WSO)的相對誤差RE=0.15%，2015年驗證產(chǎn)量的相對誤差RE=2.9%。

從模型驗證角度看，對2015、2016年的模擬結(jié)果與觀測值對比，相對誤差RE小于5%，一致性指數(shù)d均大于0.9，標(biāo)準(zhǔn)均方根誤差nRMSE均小于20%。綜合上述得出，WOFOST模型調(diào)參后對于北京大興地區(qū)的夏玉米的適應(yīng)性良好。

圖4 干物質(zhì)量、葉面積指數(shù)的優(yōu)化驗證結(jié)果Fig.4 Amount of dry matter、leaf area index optimization results

4 結(jié) 語

本研究運(yùn)用中國水利水電科學(xué)研究院大興試驗基地2016年夏玉米數(shù)據(jù)，利用擴(kuò)展的傅里葉幅度檢驗法(EFAST)對WOFOST模型的49個作物參數(shù)進(jìn)行敏感性分析，選出對產(chǎn)量、干物質(zhì)量、葉面積指數(shù)敏感性較高的10個參數(shù)進(jìn)行調(diào)參，并利用2015年的試驗數(shù)據(jù)進(jìn)行模擬驗證。得到的主要結(jié)論如下。

(1)對產(chǎn)量敏感的參數(shù)主要有：二氧化碳同化率在12 ℃下的矯正因子(TMNFTB12)、最大光合速率在30 ℃下的校正因子(TMPFTB30)、開花到成熟之間的積溫(TSUM2)、貯存器官生長同化物轉(zhuǎn)化效率(CVO)、生育期為2(DVS=2)時的可見光的消散系數(shù)(KDIFTB2)。

(2)對干物質(zhì)量敏感的參數(shù)主要有：CO2同化率在12 ℃下的矯正因子(TMNFTB12)、在35 ℃時葉面積的生命周期(SPAN)、生育期為0時的比葉面積(SLATB0)、生育期為1.1時的根干物質(zhì)的分配系數(shù)(FRTB11)、30 ℃下積溫的增長量(DTSMTB30)。

(3)對葉面積指數(shù)敏感的參數(shù)主要有：根的同化物轉(zhuǎn)換效率(CVR)、生育期為0時的比葉面積(SLATB0)、CO2同化率在12 ℃下的矯正因子(TMNFTB12)、生育期為1時的單葉片CO2同化效率(AMAXTB100)、最大光合速率在0 ℃時的矯正因子(TMPFTB0)。

(4)2016年調(diào)參后的夏玉米模擬值與試驗觀測值對比，產(chǎn)量的相對誤差RE=0.15%，干物質(zhì)量的一致性指數(shù)d=0.995，標(biāo)準(zhǔn)均方根誤差nRMSE=12%，葉面積指數(shù)的一致性指數(shù)d=0.948，標(biāo)準(zhǔn)均方根誤差nRMSE=19%。利用2015年數(shù)據(jù)對參數(shù)進(jìn)行驗證，夏玉米的模擬值與試驗觀測值對比，產(chǎn)量的相對誤差RE=2.9%，干物質(zhì)量的一致性指數(shù)d=0.993、標(biāo)準(zhǔn)均方根誤差nRMSE=12%，葉面積指數(shù)的一致性指數(shù)d=0.919、標(biāo)準(zhǔn)均方根誤差nRMSE=18%。表明優(yōu)化后的參數(shù)對于北京地區(qū)夏玉米模擬效果良好。

敏感性分析方法及調(diào)參軟件在作物生長模型中的應(yīng)用減少了人為主觀因素的影響，極大的縮短調(diào)參時間。

□