基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數(shù)優(yōu)化

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

?

基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數(shù)優(yōu)化

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

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

何亮，侯英雨，趙剛，鄔定榮，于強(qiáng). 基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數(shù)優(yōu)化[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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　引　言

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

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

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

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

1　材料與方法

1.1試驗(yàn)站點(diǎn)和數(shù)據(jù)

本文所用到的試驗(yàn)觀測數(shù)據(jù)在欒城站進(jìn)行。欒城站位于河北省石家莊市欒城縣（114°40′E，37°50′N，海拔50.1 m），為太行山前平原農(nóng)業(yè)高產(chǎn)區(qū)的典型代表。本地主要?dú)夂蝾愋褪桥瘻貛О霛駶櫦撅L(fēng)氣候，多年平均太陽輻射為5 433 MJ/m2，平均氣溫為12.3℃，大于0℃積溫為4 710℃·d，大于10℃積溫為4 232℃·d，無霜期為200 d。豐富的光熱資源滿足一年二熟或兩年三熟的作物生長。年降雨量大致為480 mm，降水分配不均，超過65%的降水集中在夏季。冬小麥生長季降水量僅僅120 mm左右，而整個(gè)生長季的需水量達(dá)460 mm以上，降水量難以滿足冬小麥生長需要[18-21]。

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

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

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

1.2WOFOST模型

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

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

1.3全局敏感性分析方法

1.3.1擴(kuò)展傅里葉幅度檢驗(yàn)法

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

模型y＝f（x1，x2，…，xk）可用合適的轉(zhuǎn)換函數(shù)轉(zhuǎn)換為y＝f（s），對f（s）進(jìn)行傅里葉變換，

式中Ns為取樣數(shù)，。

模型的總方差分解為

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

參數(shù)xi的總敏感性為

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

1.3.2敏感性參數(shù)選擇和模擬設(shè)計(jì)

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

1.4貝葉斯方法

1.4.1馬爾科夫蒙特卡洛

對于作物模型的驅(qū)動數(shù)據(jù)、模型參數(shù)和模型輸出，公式描述為

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

參數(shù)估計(jì)的問題就是已知xi和yi，來估計(jì)模型的參數(shù)θ。統(tǒng)計(jì)學(xué)中核心概念似然函數(shù)（Likelihood function），它是觀測量的聯(lián)合概率分布，滿足

根據(jù)貝葉斯統(tǒng)計(jì)觀點(diǎn)，參數(shù)也滿足一定的概率分布，假設(shè)參數(shù)的先驗(yàn)分布為p（θ）。根據(jù)貝葉斯公式，得到參數(shù)的后驗(yàn)分布為

要得到后驗(yàn)概率密度，關(guān)鍵是要解出高維的聯(lián)合概率分布函數(shù)p（yi|θ, xi），貝葉斯統(tǒng)計(jì)方法求這個(gè)高維的概率密度函數(shù)的思路是求出每一個(gè)參數(shù)的邊緣分布，如

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

馬爾科夫蒙特卡洛（MCMC）基本思想是產(chǎn)生一個(gè)馬爾科夫鏈，以目標(biāo)分布為平穩(wěn)分布，根據(jù)馬爾科夫理論，一個(gè)馬爾科夫鏈從任意初值出發(fā)，都會收斂到平穩(wěn)分布。

馬爾科夫的基本原理是：時(shí)間序列上，某一時(shí)刻t的狀態(tài)，只與前一個(gè)時(shí)刻t-1有關(guān)，與其他時(shí)刻無關(guān)，這樣的隨機(jī)時(shí)間序列θ為馬爾科夫鏈。

設(shè){θt}t＞0為空間Θ上的齊次馬爾科夫鏈，即：轉(zhuǎn)移概率函數(shù)p（.，.）與時(shí)間t無關(guān)，轉(zhuǎn)移概率函數(shù)為

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

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

1.4.2MCMC采樣方法

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

1.4.3MCMC優(yōu)化參數(shù)選擇、模型優(yōu)化和驗(yàn)證方案

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

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

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

2　結(jié)果與分析

2.1全局敏感性結(jié)果

2.1.1潛在產(chǎn)量水平下的參數(shù)敏感性

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

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

2.1.2雨養(yǎng)產(chǎn)量水平下的參數(shù)敏感性

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

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

圖1　潛在產(chǎn)量水平下的擴(kuò)展傅里葉幅度檢驗(yàn)法的敏感性分析結(jié)果Fig.1　Sensitivity results of extended Fourier amplitude sensitivity test in potential yield level

2.2馬爾科夫蒙特卡洛（MCMC）優(yōu)化

2.2.1馬爾科夫優(yōu)化參數(shù)選擇、初值和先驗(yàn)分布

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

2.2.2不同校正優(yōu)化方案得到的參數(shù)比較

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

圖2　雨養(yǎng)產(chǎn)量水平下的擴(kuò)展傅里葉幅度檢驗(yàn)法的敏感性分析結(jié)果Fig.2　Sensitivity results of extended Fourier amplitude sensitivity test in rain-fed yield level

表5　WOFOST中11個(gè)待優(yōu)化的品種參數(shù)初值和先驗(yàn)分布Table 5　Initial value and prior distribution of 11 parameters for optimization

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

2.2.3不同優(yōu)化方案的模型校正和驗(yàn)證結(jié)果

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

圖3　方案1（1998-1999季節(jié)數(shù)據(jù)作為校正）模型優(yōu)化結(jié)果Fig.3　Model calibration of plan one （data of 1998-1999 as calibration）

表7　方案1的優(yōu)化參數(shù)在其他2季的模型驗(yàn)證Table 7　Model verification in other two seasons in plan one

表8　3個(gè)優(yōu)化方案葉面積指數(shù)和干物質(zhì)模擬校正和驗(yàn)證結(jié)果精度比較Table 8　Comparison of LAI and biomass simulation accuracies of different plans of model calibration and verification

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

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

3　討　論

3.1參數(shù)的全局敏感性

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

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

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

3.2馬爾科夫蒙特卡洛參數(shù)優(yōu)化

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

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

相比傳統(tǒng)的試錯(cuò)法，本研究中的馬爾科夫蒙特卡洛方法對WOFOST參數(shù)的自動優(yōu)化方法，更具客觀性。與之其他的調(diào)參研究[1-2]，本研究在調(diào)參之前進(jìn)行的參數(shù)敏感性分析，可以有的放矢地選擇最敏感的參數(shù)，一來可以減少優(yōu)化參數(shù)的個(gè)數(shù)，二來也可以正確地選擇優(yōu)化對象。

4　結(jié)　論

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

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

2）潛在和雨養(yǎng)產(chǎn)量水平下，最大葉面積和地上生物量對參數(shù)的敏感性差異不大。

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

4）優(yōu)化后的參數(shù)，相比有水分脅迫的雨養(yǎng)的生長情況，WOFOST在潛在的生長水平下模擬較好。

[參考文獻(xiàn)]

[1] 房全孝. 根系水質(zhì)模型中土壤與作物參數(shù)優(yōu)化及其不確定性評價(jià)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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模型的調(diào)參與驗(yàn)證[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 何亮，趙剛，靳寧，等. 不同氣候區(qū)和不同產(chǎn)量水平下APSIM-Wheat模型的參數(shù)全局敏感性分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 黃敬峰，陳拉，王秀珍. 水稻生長模型參數(shù)的敏感性及其對產(chǎn)量遙感估測的不確定性[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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 作物模型參數(shù)全局敏感性分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，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 模型的冬小麥生長模擬參數(shù)全局敏感性分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 王會肖，劉昌明. 作物水分利用效率內(nèi)涵及研究進(jìn)展[J].水科學(xué)進(jìn)展，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ìn)[J]. 中國農(nóng)業(yè)氣象，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]. 植物生態(tài)學(xué)報(bào)，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方法的概念性流域水文模型參數(shù)優(yōu)選及不確定性研究[D]. 杭州：浙江大學(xué)，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] 茆師松，湯銀才. 貝葉斯統(tǒng)計(jì)[M]. 北京：中國統(tǒng)計(jì)出版社，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模型的區(qū)域小麥生育期預(yù)測與評價(jià)[J]. 中國農(nóng)業(yè)科學(xué)，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

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

基金項(xiàng)目：公益性行業(yè)（氣象）科研專項(xiàng)（GYHY201306052，GYHY201506001）

收稿日期：2015-11-13

修訂日期：2015-12-21

中圖分類號：S512.1

文獻(xiàn)標(biāo)志碼：A

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

doi：10.11975/j.issn.1002-6819.2016.02.025