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

謝松涯，張寶忠

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

0 引 言

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

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

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

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

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

1.1 WOFOST模型

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

1.2 敏感性分析方法

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

(1)

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

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

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

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

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

(6)

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

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

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

(7)

(8)

(9)

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

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

2.1 試驗區(qū)概況

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

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

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

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

3 結果與討論

3.1 全局敏感性分析

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

對于產量，一階敏感性指數(shù)Si>0.1的前6位的參數(shù)依次為二氧化碳同化率在12 ℃下的矯正因子(TMNFTB12)、最大光合速率在30 ℃下的校正因子(TMPFTB30)、開花到成熟之間的積溫(TSUM2)、貯存器官生長同化物轉化效率(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 干物質量隨時間變化的參數(shù)敏感性

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

可以看出，選取參數(shù)的一階敏感性指數(shù)和全局敏感系指數(shù)的變化趨勢基本趨于一致，但全局敏感性分析，不僅包含某一參數(shù)對模型結果的貢獻率，同時也包含不同參數(shù)之間的交互作用。CO2同化率在12℃下的矯正因子(TMNFTB12)在夏玉米的整個生長周期，都表現(xiàn)出很大的敏感性，在出苗后的73 d內一階和全局敏感性指數(shù)都保持穩(wěn)定的增長，但在73 d后至整個生育期的結束，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ù)，在此基礎上取全局敏感性指數(shù)較大的前10個作物參數(shù)作為分析對象，分別為：根的同化物轉換效率(CVR)、生育期為0時的比葉面積(SLATB0)、CO2同化率在12 ℃下的矯正因子(TMNFTB12)、生育期為1時的單葉片CO2同化效率(AMAXTB100)、最大光合速率在0 ℃時的矯正因子(TMPFTB0)、生育期為1時的比葉面積(SLATB1)、莖同化物轉換效率(CVS)、水分限制時葉片的相對死亡率(PERDL)、出苗至開花有效積溫(TSUM1)和在35 ℃時葉面積的生命周期(SPAN)，葉面積指數(shù)隨時間變化的一階以及全局敏感性參數(shù)變化如圖3所示。

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

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

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

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

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

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

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

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

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

4 結 語

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

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

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

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

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

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

□