基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數優化

何　亮，侯英雨，趙　剛，鄔定榮，于　強（. 國家氣象中心，北京 0008；　. 德國波恩大學作物科學研究組，波恩 D-55；. 中國氣象科學研究院，北京 0008；　4. 中國科學院地理科學與資源研究所陸地水循環與地表過程重點實驗室，北京 000； 5. 澳大利亞悉尼科技大學生命科學學院，悉尼 007）

?

基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數優化

何亮1，侯英雨1，趙剛2，鄔定榮3，于強4，5
（1. 國家氣象中心，北京 100081；2. 德國波恩大學作物科學研究組，波恩 D-53115；3. 中國氣象科學研究院，北京 100081；4. 中國科學院地理科學與資源研究所陸地水循環與地表過程重點實驗室，北京 100101； 5. 澳大利亞悉尼科技大學生命科學學院，悉尼 2007）

摘要：作物模型參數的敏感性分析、標定和驗證可以提高模型的效率和精準度，進而為模型應用做好準備工作。該研究結合參數全局敏感性分析方法以及貝葉斯后驗估計理論的馬爾科夫蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法，以華北欒城站三年的冬小麥觀測數據（葉面積和地上生物量）為參照，對WOFOST模型的55個品種參數進行了敏感性分析、篩選和優化。發現：1）對葉面積影響較大的參數為：生育期為0、0.5、0.6和0.75時的比葉面積、生育期為1.5時的最大光合速率、葉面積指數最大增長率；對地上干物質影響較大的參數為：生育期為1.5時的最大光合速率、生育期為0時的比葉面積、35℃時葉面積的生命周期、生育期為0時的散射消光系數、生育期為1.8時的最大光合速率、儲存器官的同化物轉換效率。2）潛在和雨養產量水平下，最大葉面積和地上生物量對參數的敏感性差異不大。3）馬爾科夫蒙特卡洛方法（MCMC）可以對WOFOST模型品種參數較好地優化；設計的3種校正-驗證方案中，第1種方案（用1998－1999年作為校正年份，1999－2000年，2000－2001年作為驗證年份）模擬效果最好。4）優化后的參數，模型對潛在產量水平模擬較好，一致性指數均大于0.9，相對均方根誤差小于20%；而對有水分脅迫的雨養情況下比潛在產量水平的模擬結果差，表明模型對水分脅迫的模擬不足。該研究為WOFOST模型區域應用和模型調整優化提供科學理論依據。關鍵詞：模型；作物；優化；WOFOST；全局敏感性分析；MCMC；模型參數優化

何亮，侯英雨，趙剛，鄔定榮，于強. 基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數優化[J]. 農業工程學報，2016，32（2）：169－179.doi：10.11975/j.issn.1002-6819.2016.02.025http://www.tcsae.org

He Liang， Hou Yingyu， Zhao Gang， Wu Dingrong， Yu Qiang. Parameters optimization of WOFOST model by integration of global sensitivity analysis and Bayesian calibration method[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE）， 2016， 32（2）: 169－179. （in Chinese with English abstract）doi：10.11975/j.issn.1002-6819.2016.02.025 http://www.tcsae.org

Email：heliang_hello@163.com

0　引　言

模型的標定（calibration）和驗證（validation）等準備工作是基于機理作物模型應用的前提。對于基于機理的作物模型，參數眾多。作物模型的很多參數（如比葉面積）不能通過田間觀測直接獲得，而是需要通過可觀測的變量（如葉面積指數）進行反演。通過準確的方法獲得模型參數是進行模型應用、提高模型可預報性的前提[1]。

根據觀測值反推參數屬于參數估計問題。對于線性方程或者簡單的非線性方程可以用最小二乘法來解決。基于過程的作物模型，刻畫了光合、干物質分配、土壤水分運移和蒸發等眾多生物物理過程，往往包含的方程較多，模型的非線性效應很明顯。對于一般的最小二乘法和非線性參數估計方法難以獲得全局的最優解。一種常用的模型參數標定是基于蒙特卡洛試錯法。這種方法根據以往的經驗或者參數文獻參考值重復和隨機地選擇參數，使得觀測值和模擬值擬合程度指標：如決定系數（R2）、相關系數、均根方差等指標達到預期的要求，便認為此組參數為最優參數。這種方法帶有很大的主觀性并且工作繁雜，而且很難獲得最可靠、最優的參數[2]。并且，這種參數估計方法計算量巨大，計算量隨著參數數量呈幾何級數增加。

為了避免上述缺點，基于過程的模型開始采用非線性參數優化方法，例如遺傳算法[3]，普適似然不確定估計法（general likelihood uncertainty estimation，GLUE）[4]等，也有類似的參數估計軟件出現，例如PEST（parameter estimation software）[1]。He等[5-6]利用GLUE方法對CERES-Maize的品種和土壤參數進行了估計，取得了很好的效果；房全孝[1]利用PEST對根系水質模型（root zone water quality model，RZWQM）的土壤和根系生長參數和作物遺傳參數進行了優化，相比傳統試錯法的校正結果，有明顯的參數優化效率。馬爾科夫蒙特卡洛方法（Markov Chain Monte Carlo，MCMC）是基于貝葉斯統計理論的一種參數估計方法，其在模型參數估計已有應用。近年來，運用該方法來對基于過程的生物模型也有若干運用[7-10]。Makowski等[11]對比了GLUE和MCMC方法的參數優化效率，發現對一個非線性農業模型的22個參數估計時MCMC方法估計的誤差比GLUE方法要小。對于這些非線性優化方法的應用過程當中，一個重要的問題是優化參數的選擇和計算量。對非敏感參數進行優化不會提高模型的精度反而會成倍加大計算量。模型參數敏感性分析可以有效地區分和界定參數的敏感性和重要性，為優化參數選擇提供有效篩選依據。參數的敏感性分析可以分為局部敏感性分析和全局敏感性分析。近年來，在過程模型的參數敏感性分析中，全局敏感性分析方法受到了學者們的親睞[12-17]。這主要是因為全局敏感性分析方法不僅僅反應了單個參數對模擬結果的影響，而且可以定量參數與參數之間的相互作用對結果的影響。鑒于此，在對參數進行優化之前，先進行參數的敏感性分析，篩選出對模擬結果影響大的參數，然后再來優化。這可以使優化算法做到有的放矢、提高效率并且減少計算時間。

本研究選擇全球得到廣泛應用的WOFOST作物模型為研究對象。以冬小麥為實例，首先用全局敏感性分析方法對模型的參數進行參數篩選。然后采用馬爾科夫蒙特卡洛（MCMC）的參數優化方法，對敏感的作物品種參數進行優化。以期探索一種普適的作物模型參數優化方法，為作物模型參數校正和模型應用提供一定的指導。

1　材料與方法

1.1試驗站點和數據

本文所用到的試驗觀測數據在欒城站進行。欒城站位于河北省石家莊市欒城縣（114°40′E，37°50′N，海拔50.1 m），為太行山前平原農業高產區的典型代表。本地主要氣候類型是暖溫帶半濕潤季風氣候，多年平均太陽輻射為5 433 MJ/m2，平均氣溫為12.3℃，大于0℃積溫為4 710℃·d，大于10℃積溫為4 232℃·d，無霜期為200 d。豐富的光熱資源滿足一年二熟或兩年三熟的作物生長。年降雨量大致為480 mm，降水分配不均，超過65%的降水集中在夏季。冬小麥生長季降水量僅僅120 mm左右，而整個生長季的需水量達460 mm以上，降水量難以滿足冬小麥生長需要[18-21]。

在欒城站進行了冬小麥的水分池試驗，試驗的時間是1998年10月到2001年6月，一共設計了16個水分試驗池，每個試驗池的面積為50 m2（5 m×10 m），四周由厚24.5 cm和深1.5 m的水泥墻隔離，防止不同水分處理的水分側向運動。一共3個冬小麥的生長年度，冬小麥的品種為高優503，小麥播種量為135 kg/hm2，采用人工播種，播種前深翻土壤到15 cm，播種時施肥磷酸二氨480 kg/hm2，復合肥1 600 kg/hm2。冬小麥生育期共設計5種不同水分處理，每個水分處理重復3次，處理包括：充分灌溉，雨養，返青期水分脅迫，拔節期水分脅迫，灌漿期水分脅迫。模型標定一般選擇充分灌溉和雨養的兩種情況，結合WOFOST模型的潛在生長和水分限制生長的兩種產量水平，本研究選擇了充分灌溉和雨養兩種情況的數據。太陽輻射、最高最低溫度、降水、氣壓、濕度等氣象因子由試驗站的自動氣象觀測站測定。干物質和葉面積每隔5～7 d測定一次。土壤水分測定由中子儀在田間每5 d測定一次，測定深度為20～160 cm，間隔20 cm。1998年到2001年3個冬小麥生育期的氣候基本特征如表1。欒城站的WOFOST土壤參數如表2。其他有關本次試驗的更詳細的描述可見Zhang等[19]和Chen等文獻[20-22]。

表1　1998－2001年3個冬小麥生育期的基本氣候特征Table 1　Climate characteristics in three winter wheat seasons between 1998 and 2001

表2　欒城站的土壤屬性特征Table 2　Soil characteristic parameters in Luancheng

1.2WOFOST模型

WOFOST模型是由荷蘭瓦赫寧根大學開發的一個根據氣象、土壤條件和管理措施模擬作物根、莖、穗生物量和土壤水分動態的模型。幾十年來，它已經在學術和工業界得到了廣泛應用。以下對模型機理進行簡略介紹，更加詳細的模型說明可以參考官方網站：www.supit.net。它模擬的主要過程包括作物發育、二氧化碳同化、呼吸作用、作物蒸騰、干物質分配、葉面積增長、干物質和葉片衰老、死亡、土壤水分平衡等過程[23-24]。WOFOST模型是一個抽象化通用的作物模型，也就是對各種作物的生長發育過程描述是一致的。它通過改變作物干物質分配和植物結構有關遺傳參數來實現對不同作物的區分和模擬。它可以模擬3種產量水平即：光溫限制的潛在產量、光溫水限制的雨養產量以及光溫水肥限制可獲得產量。它通過逐日氣象數據進行驅動，通過土壤、管理和作物參數數據限制和調整作物的生長過程。氣象數據包括太陽輻射、最高、最低氣溫、早晨的水汽壓、2米高度的平均風速、降雨量。土壤數據主要為田間持水率、飽和含水率、凋萎系數和導水率，如果需要模擬地下水的影響，則還需提供土壤水分特征曲線和導水率曲線參數。最重要的就是作物參數，其中包括不同發育階段所需要的積溫、光周期影響因子、不同生育期的最大光合速率、不同生育期的比葉面積、干物質分配系數、干物質和葉片的死亡率等。其中，不確定性最大的參數也就是作物品種參數。對于作物品種參數的詳細描述可見表3。

表3　WOFOST品種參數上下限和分布Table 3　Upper and down bound of cultivar parameters in WOFOST

1.3全局敏感性分析方法

1.3.1擴展傅里葉幅度檢驗法

全局敏感性分析方法很多，包括：Morris參數篩選法、基于方差的Sobol法、傅里葉幅度檢驗法（Fourier amplitude sensitivity test，FAST）、擴展傅里葉幅度檢驗法（extended Fourier amplitude sensitivity test，EFAST）等[25]。本研究采用樣本數要求低，計算高效的EFAST法，它是Saltelli等[25]結合了Sobol法和傅里葉幅度法兩種方法的優點提出了一種新的全局敏感性分析方法。其算法的基本思想是分解參數對模型結果的方差，把參數敏感性分為兩種類型：單個參數對結果的影響和參數之間耦合對模型結果的影響。其中單個參數獨立作用是用主敏感性度指數（main effect）衡量，而參數相互作用則用總敏感度（total effect）和主敏感度的差別衡量。算法簡單介紹如下：

模型y＝f（x1，x2，…，xk）可用合適的轉換函數轉換為y＝f（s），對f（s）進行傅里葉變換，

式中Ns為取樣數，。

模型的總方差分解為

式中Vi為參數xi輸入變化單獨引起的模型方差，Vij為參數xi通過參數xj作用貢獻的耦合方差，Vijm為參數xi通過參數xj，xm作用貢獻的方差，則依次類推，V1，2，…，k為參數xi通過x1，2，…，k貢獻的方差。通過歸一化處理參數xi的一階敏感性指數Si定義為

參數xi的總敏感性為

式中V-i為不包括參數xi的所有參數方差之和。算法詳細的介紹參考Saltelli等文獻[25]。

1.3.2敏感性參數選擇和模擬設計

主要的參數說明和參數范圍見表3。其中品種參數的范圍選擇來源于模型文檔提供的合理的范圍。以上所有的變量都服從均一分布，見表3。其中關于WOFOST模型中控制生育期的參數：出苗到開花的積溫（TSUM1）和開花到成熟的積溫（TSUM2）是根據觀測實際算出。模型的其他參數：氣象數據、播種日期、土壤參數都為實際的觀測。模型考慮兩個輸出，最大的葉面積指數（MAXLAI）和地上生物量（TAGP）。考慮兩種產量水平即潛在和水分限制。在EFAST方法中每次的敏感性計算需要運行n×p次，其中n為采樣數量，p為參數個數。在EFAST方法中認為參數采樣個數大于65倍的參數個數為有效（即采樣個數≥參數個數×65），因此，為減少模擬次數，先計算敏感性結果隨n的收斂性，大致在n 取80左右結果收斂，本研究采樣大小n取150，采用EFAST方法，一共需要模擬3（3年）×150×55（參數）＝24 750次。敏感性計算框架結合R語言（http://www.r-project.org/）的Sensitivity包。用R語言編寫程序調用WOFOST的函數，用Sensitivity包生成參數樣本，然后交給WOFOST計算，最后利用Sensitivity包分析和計算模擬結果的敏感性。更為詳細說明見何亮等的文獻[13]。

1.4貝葉斯方法

1.4.1馬爾科夫蒙特卡洛

對于作物模型的驅動數據、模型參數和模型輸出，公式描述為

式中yi代表模型n個輸出，例如生物量，葉面積等。xi為模型驅動因子，如溫度，太陽輻射和降水等。θ為模型參數，例如作物品種參數，土壤參數等，εi為模型的隨機誤差，均值為0，方差為σ2，未知。

參數估計的問題就是已知xi和yi，來估計模型的參數θ。統計學中核心概念似然函數（Likelihood function），它是觀測量的聯合概率分布，滿足

根據貝葉斯統計觀點，參數也滿足一定的概率分布，假設參數的先驗分布為p（θ）。根據貝葉斯公式，得到參數的后驗分布為

要得到后驗概率密度，關鍵是要解出高維的聯合概率分布函數p（yi|θ, xi），貝葉斯統計方法求這個高維的概率密度函數的思路是求出每一個參數的邊緣分布，如

計算該邊緣分布有兩個難點：第一是歸一化常數未知，第二是高維數值積分的困難。鑒于此，為解決此問題，出現了一系列的馬爾科夫鏈蒙特卡羅（MCMC）算法。

馬爾科夫蒙特卡洛（MCMC）基本思想是產生一個馬爾科夫鏈，以目標分布為平穩分布，根據馬爾科夫理論，一個馬爾科夫鏈從任意初值出發，都會收斂到平穩分布。

馬爾科夫的基本原理是：時間序列上，某一時刻t的狀態，只與前一個時刻t-1有關，與其他時刻無關，這樣的隨機時間序列θ為馬爾科夫鏈。

設{θt}t＞0為空間Θ上的齊次馬爾科夫鏈，即：轉移概率函數p（.，.）與時間t無關，轉移概率函數為

p（θ,θ*）為馬爾科夫鏈的轉移核，即為推薦分布。對于分布π（θ），滿足

如果?θ*∈Θ，則π（θ）為轉移核p（θ,θ*）的平穩分布。當某θt～π（ θ），則θj～π（θ）， j＝t＋1，…理論上講，無論θ0取何種分布，經過長時間搜索后，θt的邊緣分布即為平穩分布π（θ），則稱馬爾科夫收斂[26]。

1.4.2MCMC采樣方法

常用的MCMC采樣方法有Metropolis 算法、M-H （metropolis-hastings）算法、Gibbs采樣和自適應的Metropolis算法[27]。本研究用M-H抽樣方法，步驟為：1）構造合適的建議分布q（?| θt）；2）根據先驗分布g產生θ0；3）從（|t）q? θ中產生候選點θ；4）從均勻分布U（0，1）中產生參數矩陣U；5）判斷：若，則接受θ′，并令θt＋1＝ θ'，否則令θt＋1＝ θt；6）增加t，返回到第3步。

1.4.3MCMC優化參數選擇、模型優化和驗證方案

本研究中，共涉及3個年度的試驗數據，如表1。用來優化的觀測數據為生育期內不同時間段觀測的地上生物量和葉面積指數，每年包括雨養和充分灌溉2套地上生物量和葉面積觀測數據。參數估算方法中似然函數（式9）的構建是同時考慮了葉面積和生物量。待優化的作物品種參數為參數敏感性分析后，選取的較敏感的參數。由于3個季節是同一個品種，模型優化和驗證設計3個情景，如表4：分別用其中一年的所有全灌溉和雨養的數據作為模型的校正優化，其他剩下兩年作為模型的驗證。模型待優化的參數選擇依賴于全局敏感性分析的結果，這將在2.2.1節中詳細介紹。模型校正和驗證的過程用觀測和模擬的一致性指數（index of agreement，d，見式14）和相對均方根誤差（relative root mean square error，RRMSE，見式16）來評價，d越接近1，RRMSE越小，說明模擬精度越高。

式中d為一致性指數，RMSE（root mean square error）為均方根誤差，RRMSE為相對均方根誤差，Si為第i個模擬值，Oi為第i個觀測值，O為觀測值的均值，n為樣本數。

表4　WOFOST模型的不同校正-驗證方案Table 4　Different plans for calibration and verification of WOFOST model

2　結果與分析

2.1全局敏感性結果

2.1.1潛在產量水平下的參數敏感性

在潛在產量水平下，如圖1a，對于最大葉面積指數（MAXLAI），一階敏感性指數最大的前5個因子依次為：生育期為0時的比葉面積（SLATB00）、生育期為0.5時的比葉面積（SLATB050）、生育期為0.6時的比葉面積（SLATB060）、生育期為0.75時的比葉面積（SLATB075）、生育期為1.5時的最大光合速率（AMAXTB150），其敏感性值分別為12.8%、2.7%、2.7%、2.5%、2.5%，其他參數的值分別小于1%；全局敏感性指數最大的前5個參數與一階敏感性指數一樣，其值分別為17.4%、4.4%、3.6%、3.4%、3.3%，但是對于參數SPAN（35℃時葉面積的生長周期）也較大，排在第六，值為1.9%，其余的參數分別小于1%。

如圖1b，對于地上生物量（TAGP），一階敏感性指數最大的前5個因子依次為：生育期為0時的比葉面積（SLATB00）、35℃時葉面積的生命周期（SPAN）、生育期為1.5時的最大光合速率（AMAXTB150）、儲存器官的同化物轉換效率（CVO）、生育期為1.8時的最大光合速率（AMAXTB180），其值分別為6.1%、4.1%、2.6%、2.2%、1.8%，其他參數的值分別小于1%；全局敏感性指數最大的前5個參數與一階敏感性指數一樣，其值分別為8.2%、5.3%、3.1%、2.9%、2.3%，其余的參數分別小于2%。

2.1.2雨養產量水平下的參數敏感性

在雨養產量水平下，如圖2a，對于最大葉面積指數（MAXLAI），一階敏感性指數最大的前5個因子依次為：生育期為0時的比葉面積（SLATB00）、生育期為0.5時的比葉面積（SLATB050）、生育期為0.6時的比葉面積（SLATB060）、生育期為1. 5時的最大光合速率（AMAXTB150）、生育期為0.75時的比葉面積（SLATB075）、其值分別為12.1%、4.8%、2.1%、2.1%、1.2%，其他參數的值分別小于1%；全局敏感性指數最大的前4個參數與一階敏感性指數一樣，其值分別為15.6%、6.3%、3.5%、2.8%，參數RGRLAI（葉面積指數最大增長率）也較大，排在第5，值為2.7%，其余的參數分別小于2%。

如圖2b，對于地上生物量（TAGP），一階敏感性指數最大的前5個因子依次為：生育期為1.5時的最大光合速率（AMAXTB150）、生育期為0時的比葉面積（SLATB00）、35℃時葉面積的生命周期（SPAN）、生育期為0時的散射消光系數（KDIFFTB00）、生育期為1. 8時的最大光合速率（AMAXTB180），其值分別為4.8%、4.4%、2.6%、2.0%、1.3%，其他參數的值分別小于1%；全局敏感性指數最大的前5個參數分別是生育期為0時的比葉面積（SLATB00）、生育期為0.5時的比葉面積（SLATB050）、35℃時葉面積的生命周期（SPAN）、生育期為1.8時的最大光合速率（AMAXTB180）、葉齡的低溫閾值（TBASE），其值分別為23.0%、7.1%、6.7%、6.6%、4.2%，其余的參數分別小于4%。

圖1　潛在產量水平下的擴展傅里葉幅度檢驗法的敏感性分析結果Fig.1　Sensitivity results of extended Fourier amplitude sensitivity test in potential yield level

2.2馬爾科夫蒙特卡洛（MCMC）優化

2.2.1馬爾科夫優化參數選擇、初值和先驗分布

根據上節的敏感性結果，選擇葉面積和干物質敏感較大的參數作為優化對象。其中主要包括以下11個參數：生育期為0時的比葉面積（SLATB00）、生育期為0.5時的比葉面積（SLATB050）、生育期為0.6時的比葉面積（SLATB060）、生育期為0.75時的比葉面積（SLATB075）、生育期為1.5時的最大光合速率（AMAXTB150）、生育期為1.8時的最大光合速率（AMAXTB180）、RGRLAI（葉面積指數最大增長率）、35℃時葉面積的生命周期（SPAN）、生育期為0時的散射消光系數（KDIFFTB00）、葉齡的低溫閾值（TBASE）儲存器官的同化物轉換效率（CVO）；這11個參數的初值是根據觀測值，先用試錯法手動大致調整參數，使觀測葉面積和干物質趨勢大致一樣時的參數值。這11個參數的初值如表5，每個參數的先驗分布都為均一分布。在優化的過程中，對于11個參數給予±10%的上下擾動。其他敏感性不大的參數采樣模型中默認的參數。對于控制生育期的參數：出苗到開花的積溫（TSUM1）和開花到成熟的積溫（TSUM2）是根據觀測實際算出。

2.2.2不同校正優化方案得到的參數比較

模型優化采樣3000次，收斂的判斷是做出跡圖（trace plot），即將所產生的樣本對迭代次數作圖，生成馬氏鏈的一條樣本路徑。如果當采樣次數足夠大時，路徑表現出穩定性沒有明顯的周期和趨勢，就可以認為是收斂（由于篇幅，省略了跡圖）。MCMC優化后，待參數趨于收斂后，11個參數得后驗分布情況（由于篇幅，省略了11個參數的后驗概率分布圖）。根際參數的后驗分布，把后驗分布的均值作為參數的優化值。3種優化方案得到的參數如表6。總體來看，三種優化方案的各個參數變異系數相差較小，最大的葉齡的低溫閾值（TBASE），為3%。不同的優化方案下，優化后的參數變異性小，說明馬爾科夫蒙特卡洛（MCMC）方法在不同環境下的優化可靠性較好。

圖2　雨養產量水平下的擴展傅里葉幅度檢驗法的敏感性分析結果Fig.2　Sensitivity results of extended Fourier amplitude sensitivity test in rain-fed yield level

表5　WOFOST中11個待優化的品種參數初值和先驗分布Table 5　Initial value and prior distribution of 11 parameters for optimization

表6　不同參數優化策略下的WOFOST作物品種參數取值Table 6　Estimation results of crop parameters of winter wheat for WOFOST in different parameter estimation plans

2.2.3不同優化方案的模型校正和驗證結果

本研究設計了3種模型校正-驗證方案。對觀測的葉面積指數、地上干物質進行了驗證。通過對比分析不同的校正-驗證方案，選擇較優的方案。圖3和表7是方案1的詳細校正和驗證結果。其他2種方案，處于簡潔需要，給出總體的驗證結果，如表8。

圖3　方案1（1998-1999季節數據作為校正）模型優化結果Fig.3　Model calibration of plan one （data of 1998-1999 as calibration）

表7　方案1的優化參數在其他2季的模型驗證Table 7　Model verification in other two seasons in plan one

表8　3個優化方案葉面積指數和干物質模擬校正和驗證結果精度比較Table 8　Comparison of LAI and biomass simulation accuracies of different plans of model calibration and verification

從方案1的校正結果看（圖3），模型對潛在的葉面積和干物質模擬的很好，一致性指數d分別為0.92，0.95，相應的RRMSE分別為15%，22%；對于雨養情況下，葉面積指數模擬較潛在情況下要差，一致性指數小于0.8，但模型對雨養情況下的干物質模擬也較好，一致性指數達到了0.94，RRMSE為22%；從方案1的模型驗證來看，如表7：潛在產量水平下，模型對2a的葉面積和干物質都模擬較好，一致性指數d都大于0.80，RRMSE小于21%。雨養產量水平下，模型在2000-2001年的葉面積指數模擬稍差，一致性指數僅僅為0.45，RRMSE達到了136%；

比較3個校正-驗證方案，如表8，從模型驗證的角度看，在潛在產量水平下三種方案都很好，一致性指數都不小于0.90，且RRMSE都小于20%。從雨養產量水平下一致性指數在0.74到0.77之間，雨養下的RRMSE都較大。綜合潛在和雨養兩種情況而言，方案1較其他兩個方案略好。

3　討　論

3.1參數的全局敏感性

本研究首先用擴展傅里葉幅度檢驗法（EFAST），分析了潛在和雨養兩種生產水平下，55個WOFOST品種參數的敏感性的大小。從結果看，最大葉面積（MAXLAI）、地上生物量（TAGP）對于光合速率相關的參數，如生育期為1.5時的最大光合速率（AMAXTB150）、生育期為1.8時的最大光合速率（AMAXTB180），以及有關比葉面積的參數，如，生育期為0時的比葉面積，SLATB00等反應較為敏感。這是因為不同生育階段的最大光合速率是限制光合作用大小的主要原因，是控制物質源的根本參數。而對于葉面積，比葉面積參數、最大葉面積增長率是控制葉片生長的關鍵參數，比葉面積控制了單位物質轉換成葉面積的大小，因而葉面積對這些參數敏感性較強。從其他主流的作物模型參數敏感性分析來看，例如APSIM模型中的谷粒最大灌漿速率和輻射利用效率[13]，CERES模型中的光能利用率[16]都是制約同化物質量大小的基本參數。與Wang 等[12]對WOFOST玉米的參數分析比較來看，本研究有些差別，Wang等的研究中葉片在35℃生命周期（SPAN）為對干物質積累最敏感的參數，差別的原因可能是作物類型的不一致，因為作物類型的不一樣，不同參數的取值范圍就有差異。與Ceglar等[28]對WOFOST玉米的分析比較，比葉面積都是比較敏感的參數。

從不同產量水平下的參數敏感性結果看，最敏感的參數基本上是一致的，如圖1和2。這也說明不同的產量水平下，模型的最大葉面積、地上生物量對模型參數的敏感性具有一致的表現，WOFOST模型的結果和作者對APSIM模型[13，29]分析結果基本是一致的。因而可以斷定，產量水平不是造成模型輸出對品種參數敏感性差異的原因。

Wang等[12]證實模型參數范圍、采樣次數對敏感性的結果影響很大。對于參數取值范圍而言，不同的作物類型有不同的取值。本研究的參數范圍是在前人的基礎上結合模型說明里面的取值給定，沒有考慮參數范圍變化對敏感性結果的影響，這主要是基于本研究中參數敏感性的主要目的是為了客觀地選擇出敏感性的參數，為下一步優化做準備。而對于采樣次數的影響，本研究認為，需要先對不同的采樣次數做計算，一直到敏感性結果趨于收斂，這才是最終的敏感性結果，這就消除了采樣次數對敏感性結果的影響。

3.2馬爾科夫蒙特卡洛參數優化

通過敏感性篩選結果，選擇了WOFOST模型中11個對葉面積和干物質影響大的參數作為優化對象，設計3種校正-驗證方案。從結果看，3種方案得到的參數輸入模型，模型都較好地模擬了觀測值。從3種優化方案的參數結果看，3套參數都趨于穩定，這也說明馬爾科夫蒙特卡洛方法可以較為穩定地反演出最優的參數，證實這套方法用在WOFOST模型調參的可行性。從模型驗證來看，模型對潛在水平的模擬比雨養下的要好，這可能和WOFOST的模型結構有關，雨養水平下，WOFOST模型把土壤水分循環模塊添加了進來，而對于潛在產量水平，模型默認是把土壤一直是處于田間持水量的狀態。雨養情況下模型的復雜度比潛在情況下要復雜，且由于降水少，作物生長處于水分脅迫的時候多，由此可見模型對于水分脅迫的模擬是不足的，類似的研究在CERES模型中也見有報道[2]。

對于非線性優化的過程中，優化收斂性對初值非常敏感。如果初始值選擇不好，有時候可能很難達到收斂。為了減少初值給優化帶來的收斂性難的問題，本研究采取的策略是先手動粗略地調整參數，使模擬趨勢大致差不多時的參數作為初始值。這可以減少收斂性帶來的不確定性。參數的確定，依賴于觀測數據，當觀測數據樣本能夠較好地反應總體時，估計得到的參數就越準確，也能夠避免“異參同效”的現象[30]。在反演品種參數的過程中，能有較多年份的觀測數據對反演品種參數是有幫助的。MCMC算法獲得參數的后驗分布而不是一組唯一的最優參數解，一般認為最優那一組參數，是對應的似然函數值最大的那一套參數。各個參數的后驗平均值作為一套品種參數，這套品種參數并不一定是最優的那一套品種參數。但是這個均值和MCMC獲得的參數最大后驗概率對應參數較為接近，說明參數分布是朝著高概率區進行收斂。

相比傳統的試錯法，本研究中的馬爾科夫蒙特卡洛方法對WOFOST參數的自動優化方法，更具客觀性。與之其他的調參研究[1-2]，本研究在調參之前進行的參數敏感性分析，可以有的放矢地選擇最敏感的參數，一來可以減少優化參數的個數，二來也可以正確地選擇優化對象。

4　結　論

本研究，運用欒城1998－2001年3個生育期的冬小麥在全灌溉和雨養的田間數據，對WOFOST模型55個參數進行了潛在生長、雨養生長水平下的參數全局敏感性分析，篩選出11個對葉面積和干物質影響最大的參數，運用馬爾科夫蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法，設計了3種校正-驗證方案，得到的結論如下：

1）在WOFOST冬小麥的55個作物參數中，對葉面積影響最大的參數為：生育期為0時的比葉面積（SLATB00）、生育期為0.5時的比葉面積（SLATB050）、生育期為0.6時的比葉面積（SLATB060）、生育期為0.75時的比葉面積（SLATB075）、生育期為1.5時的最大光合速率（AMAXTB150）、葉面積指數最大增長率（RGRLAI）。對地上干物質影響較大的參數為：生育期為1.5時的最大光合速率（AMAXTB150）、生育期為0時的比葉面積（SLATB00）、35℃時葉面積的生命周期（SPAN）、生育期為0時的散射消光系數（KDIFFTB00）、生育期為1.8時的最大光合速率（AMAXTB180）、儲存器官的同化物轉換效率（CVO）。

2）潛在和雨養產量水平下，最大葉面積和地上生物量對參數的敏感性差異不大。

3）馬爾科夫蒙特卡洛方法（MCMC）可以對WOFOST模型品種參數較好地優化，設計的三種校正-驗證方案中，第一種方案（用1998－1999年的數據作為校正年份，1999－2000年，2000－2001年的數據作為驗證年份）模擬效果最好。其中對潛在產量水平的葉面積和生物量的模擬，一致性指數d都大于0.9；對雨養產量水平的葉面積和生物量的模擬，一致性指數d＞0.75。

4）優化后的參數，相比有水分脅迫的雨養的生長情況，WOFOST在潛在的生長水平下模擬較好。

[參考文獻]

[1] 房全孝. 根系水質模型中土壤與作物參數優化及其不確定性評價[J]. 農業工程學報，2012，28（10）：118－23. Fang Quanxiao. Optimizing and uncertainty evaluation of soil and crop parameters in root zone water quality model[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE）， 2012， 28（10）: 118－123. （in Chinese with English abstract）

[2] 姚寧，周元剛，宋利兵，等. 不同水分脅迫條件下DSSAT-CERES-Wheat模型的調參與驗證[J]. 農業工程學報，2015，31（12）：138－150. Yao Ning， Zhou Yuangang， Song Libing， et al. Parameter estimation and verification of DSSAT-CERES-Wheat model for simulation of growth and development of winter wheat under water stresses at different growth stages[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE）， 2015， 31（12）: 138－150. （in Chinese with English abstract）

[3] Wang Q J. The genetic algorithm and its application to calibrating conceptual rainfall-runoff models[J]. Water Resourses Research， 1991， 27（9）: 2467－71.

[4] Beven K， Binley A. The future of distributed models: Model calibration and uncertainty prediction[J]. Hydrological Processes， 1992， 6（3）: 279－298.

[5] He Jiangqiang， Dukes， M D， Jones， J W， et al. Applying GLUE for estimating CERES-Maize genetic and soil parameters for sweet corn production[J]. Transactions of the ASABE， 2009， 52（6）: 1907－1921.

[6] He Jiangqiang， Jones， J.W.， Graham W.D.， et al. Influence of likelihood function choice for estimating crop model parameters using the generalized likelihood uncertaintyestimation method[J]. Agricultural Systems， 2010， 103（5）: 256－264.

[7] Iizumi T， Yokozawa M， Nishimori M. Parameter estimation and uncertainty analysis of a large-scale crop model for paddy rice: Application of a Bayesian approach[J]. Agricultural and Forest Meteorology， 2009， 149（2）: 333－348.

[8] van Oijen M， Cameron D R， Butterbach-Bahl K， et al. A Bayesian framework for model calibration， comparison and analysis: Application to four models for the biogeochemistry of a Norway spruce forest[J]. Agricultural and Forest Meteorology， 2011， 151（12）: 1609－1621.

[9] van Oijen M， Reyer C， Bohn F J， et al. Bayesian calibration，comparison and averaging of six forest models， using data from Scots pine stands across Europe [J]. Forest Ecology and Management， 2013， 289: 255－268.

[10] Tao Fulu， Yokozawa M， Zhang Zhao. Modelling the impacts of weather and climate variability on crop productivity over a large area: a new process-based model development，optimization， and uncertainties analysis[J]. Agricultural and Forest Meteorology， 2009， 149（5）: 831－850.

[11] Makowski D， Wallach D， Tremblay M. Using a Bayesian approach to parameter estimation: comparison of the GLUE and MCMC methods[J]. Agronomie， 2002， 22（2）: 191－203.

[12] Wang Jing， Li Xin， Lu Ling， et al. Parameter sensitivity analysis of crop growth models based on the extended Fourier Amplitude Sensitivity Test method[J]. Environmental Modelling ＆ Software， 2013， 48: 171－182.

[13] 何亮，趙剛，靳寧，等. 不同氣候區和不同產量水平下APSIM-Wheat模型的參數全局敏感性分析[J]. 農業工程學報，2015，31（14）：148－157. He Liang， Zhao Gang， Jin Ning， et al. Global sensitivity analysis of APSIM-Wheat parameters in different climate zones and yield levels[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE），2015， 31（14）: 148－157. （in Chinese with English abstract）

[14] 黃敬峰，陳拉，王秀珍. 水稻生長模型參數的敏感性及其對產量遙感估測的不確定性[J]. 農業工程學報，2012，28（19）：119－129. Huang Jingfeng， Chen La， Wang Xiuzhen. Sensitivity of rice growth model parameters and their uncertainties in yield estimation using remote sensing date[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE）， 2012， 28（19）: 119－129. （in Chinese with English abstract）

[15] 姜志偉，陳仲新，周清波，等. CERES-Wheat 作物模型參數全局敏感性分析[J]. 農業工程學報，2011，27（1）：236－242. Jiang Zhiwei， Chen Zhongxin， Zhou Qingbo， etal. Global sensitivity analysis of CERES-Wheat model parameters[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE）， 2011， 27（1）: 236－242. （in Chinese with English abstract）.

[16] 宋明丹，馮浩，李正鵬，等. 基于Morris 和EFAST 的CERES-Wheat 模型敏感性分析[J]. 農業機械學報，2014，45（10）：124－131. Song Mingdan， Feng Hao， Li Zhengpeng， et al. Global sensitivity analyses of DSSAT-CERES-wheat model using morris and EFAST Methods[J]. Transactions of the Chinese Society for Agricultural Machinery （Transactions of the CSAM）， 2014， 45（10）: 124－131. （in Chinese with English abstract）.

[17] 吳錦，余福水，陳仲新，等. 基于EPIC 模型的冬小麥生長模擬參數全局敏感性分析[J]. 農業工程學報，2009，25（7）：136－142. Wu Jin， Yu Fushui， Chen Zhongxin， et al. Global sensitivity analysis of growth simulation parameters of winter wheat based on EPIC model[J]. Transactions of the Chinese Society of Agricultural Engineering （Transactions of the CSAE），2009， 25（7）: 136－142. （in Chinese with English abstract）.

[18] 王會肖，劉昌明. 作物水分利用效率內涵及研究進展[J].水科學進展，2000，11（1）：99－104. Wang Huixiao， Liu Changming. Advances in crop water use eff ic iency research[J]. Advances in Water Sceience， 2000，11（1）: 99－104. （in Chinese with English abstract）.

[19] Zhang Yongqiang， Kendy E， Yu Qiang， et al. Effect of soil water deficit on evapotranspiration， crop yield， and water use efficiency in the North China Plain[J]. Agricultural Water Management， 2004， 64（2）: 107－122.

[20] Chen Chao， Wang Enli， Yu Qiang. Modeling wheat and maize productivity as affected by climate variation and irrigation supply in north China Plain[J]. Agronomy Jounal，2010， 102（3）: 1037－1049.

[21] Chen Chao， Wang Enli， Yu Qiang， et al. Quantifying the effects of climate trends in the past 43 years （1961－2003） on crop growth and water demand in the North China Plain[J]. Climatic Change， 2010， 100（3/4）: 559－578.

[22] Chen Chao， Wang Enli， Yu Qiang. Modelling the effects of climate variability and water management on crop water productivity and water balance in the North China Plain[J]. Agricultural Water Management， 2010， 97（8）: 1175－1184.

[23] 馬玉平，王石立，張黎. 針對華北小麥越冬的WOFOST模型改進[J]. 中國農業氣象，2005，26（3）：145－149. Ma Yuping， Wang Shili， Zhang Li. Study on improvement of WOFOST against overwinter of wheat in north China[J]. Chinese Journal of Agrometeorology， 2005， 26（3）: 145－149. （in Chinese with English abstract）

[24] 鄔定榮，歐陽竹，趙小敏，等. 作物生長模型WOFOST在華北平原的適用性研究[J]. 植物生態學報，2003，27（5）：594－602. Wu Dingrong， Ou Yangzhu， Zhao Xiaoming， et al. The applicability reseracher of WOFOST model in North China Plain[J]. Acta Phytoecologica Sinia， 2003， 27（5）: 594－602. （in Chinese with English abstract）.

[25] Saltelli A， Tarantola S， Campolongo F， et al. Sensitivity analysis in practice: a guide to assessing scientific models[M]. John Wiley and Sons， 2004: 20－78.

[26] 曹飛鳳. 基于MCMC方法的概念性流域水文模型參數優選及不確定性研究[D]. 杭州：浙江大學，2010：16－25. Cao Feifeng. Study on Paramater Optimization and Uncertainty Analysis for Concenptual Hydrological ModelBase on MCMC Method[D]. Hangzhou: Zhejiang University，2010: 16－25. （in Chinese with English abstract）

[27] 茆師松，湯銀才. 貝葉斯統計[M]. 北京：中國統計出版社，2012.

[28] Ceglar A， ?repin?ek Z， Kajfe?-Bogataj L， et al. The simulation of phenological development in dynamic crop model: The Bayesian comparison of different methods[J]. Agricultural and Forest Meteorology， 2011， 151（1）: 101－115.

[29] Zhao Gang， Bryan B A， Song Xiaodong. Sensitivity and uncertainty analysis of the APSIM-wheat model: Interactions between cultivar， environmental， and management parameters[J]. Ecological Modelling， 2014， 279: 1－11.

[30] 呂尊富，劉小軍，湯亮，等. 基于WheatGrow和CERES模型的區域小麥生育期預測與評價[J]. 中國農業科學，2013，46（6）：1136－1148.

Lu Zunfu， Liu Xiaojun， Tang Liang， et al. Regional prediction and evaluation of wheat phenology based on the WheatGrow and CERES models[J]. Scientia Agricultura Sinica， 2013， 46（6）: 1136－1148. （in Chinese with English abstract）

Parameters optimization of WOFOST model by integration of global sensitivity analysis and Bayesian calibration method

He Liang1, Hou Yingyu1, Zhao Gang2, Wu Dingrong3, Yu Qiang4,5
（1. Nɑtionɑl Meteorologicɑl Center， Beijing 100081， Chinɑ;2. Crop Science Group， Institute of Crop Science ɑnd Resource Conservɑtion （INRES）， University of Bonn， Kɑtzenburgweg 5， D-53115 Bonn， Germɑny;3. Chinese Acɑdemy of Meteorologicɑl Sciences， Beijing 100081，Chinɑ;4. Key Lɑborɑtory of Wɑter Cycle ＆ Relɑted Lɑnd Surfɑce Processes， Institute of Geogrɑphic Sciences ɑnd Nɑturɑl Resources Reseɑrch， University of Chinese Acɑdemy of Science， Beijing 100101， Chinɑ;5. School of Life Sciences， University of Technology Sydney，PO Box 123， Broɑdwɑy， NSW， 2007， Austrɑliɑ）

Abstract:Crop model calibration and validation are essential for model evaluation and application. It is important for model application to accurately estimate the values of crop model parameters and further improve the capacity of model prediction. In the previous researches， trial-and-error method was widely used in model calibration and validation. The deficiency of this method was subjective selection of parameter values and time-consuming processes. To overcome these issues， the optimization methods such as general likelihood uncertainty estimation （GLUE）， genetic algorithm （GA） and shuffled complex evolution （SCE-UA） algorithm were alternative method for model calibration and validation. However， it is a problem to decide which parameters for optimization. It is essential to select the most sensitive parameters among hundreds of parameters in the crop model for optimization. To avoid subjective selection of parameters for calibration and validation， we used the global sensitivity analysis method of model parameters and the Markov Chain Monte Carlo （MCMC） method based on Bayesian theory to optimize the crop genetic parameters in the WOFOST （world food studies）， and the data of three-year winter wheat field experiment in Luancheng in the North China Plain were adopted. The main objectives were: 1） to analyze the sensitivity and uncertainty of WOFOST brought by 55 crop genetic parameters using the extended Fourier amplitude sensitivity test; 2） to calibrate and validate the WOFOST using the MCMC method after sensitivity analysis. We found that: 1）The most sensitive parameters for maximum leaf area index （MAXLAI） in the crop growth period were successively: specific leaf area at development stage of 0， 0.5， 0.6， and 0.75， maximum CO2assimilation rate at development stage of 1.5， and maximum relative increase in LAI （RGTLAI）; 2） The most sensitive parameters for total above ground production （TAGP） in the crop growth period were successively: maximum CO2assimilation rate at development stage of 1.5 （AMAXTB150），specific leaf area at development stage of 0 （SLATB00）， life span of leaves growing at 35oC， extinction coefficient for diffuse visible light at development stage of 0 （KDIFFTB00）， maximum CO2assimilation rate at development stage of 1.8 （AMAXTB180）， efficiency of conversion into storage organs （CVO）; 3） The parameter sensitivity for MAXLAI and TAGP in potential and rain-fed production level was almost coincident， which indicated that yield level didn’t influence the parameter sensitivity results; 4） Eleven sensitive parameters were selected for optimization by using the MCMC method. The first calibration and validation strategy （i.e. the data in 1998-1999 for calibration and those in 1999-2000 and 2000-2001 for validation）， was better than other 2 strategies. 5） WOFOST simulation was much improved if the optimized parameters by the MCMC method were adopted. The index of agreement was higher than 0.9 and the relative root mean square error was less than 20%. However， WOFOST performed worse in rain-fed case because water stress factor was added to limit crop growth. The results indicate that more sensitive parameters should have priority in adjusting values for model calibration and validation. In addition， the MCMC method is a feasible optimization method for WOFOST calibration and validation.

Keywords:models; crops; optimization; WOFOST; global sensitivity; MCMC; model parameters optimization

作者簡介：何亮，男（漢族），湖南永興縣人，博士，主要從事作物模型、農業氣象和全球變化研究。北京國家氣象中心，100081。

基金項目：公益性行業（氣象）科研專項（GYHY201306052，GYHY201506001）

收稿日期：2015-11-13

修訂日期：2015-12-21

中圖分類號：S512.1

文獻標志碼：A

文章編號：1002-6819（2016）-02-0169-11

doi：10.11975/j.issn.1002-6819.2016.02.025