黃健熙，黃 海，馬鴻元，李 穎，侯英雨，何 亮，朱德海

(1. 中國農(nóng)業(yè)大學土地科學與技術學院，北京 100083；2. 農(nóng)業(yè)部農(nóng)業(yè)災害遙感重點實驗室，北京 100083；3. 中國氣象局河南省農(nóng)業(yè)氣象保障與應用技術重點開放實驗室，鄭州 450003；4. 河南省氣象科學研究所，鄭州 450003；5. 中國氣象局國家氣象中心，北京 100081)

0 引 言

氣候變化、人口膨脹、資源短缺等一系列問題對農(nóng)業(yè)的高效持續(xù)發(fā)展提出了更高要求，如何在宏觀上準確地進行作物估產(chǎn)，獲取作物生長關鍵性參數(shù)，以幫助確保糧食安全，已經(jīng)越來越成為研究的熱點?；跈C理過程模擬的作物生長模型普遍用于農(nóng)作物產(chǎn)量預報、參數(shù)估計、災害影響評估[1-7]，還有學者將其成功應用于作物重金屬、營養(yǎng)品質(zhì)監(jiān)測等[8-9]，在不同區(qū)域環(huán)境、不同作物研究中均得到廣泛驗證。

作物模型機理性模擬的準確性依賴于眾多模型參數(shù)、氣象驅(qū)動數(shù)據(jù)、土壤參數(shù)的準確獲取，模型應用的區(qū)域越大，數(shù)據(jù)的直接、準確獲取越困難，模擬的精度也越難以保證。實際應用中，往往將一定的地理區(qū)域進行人為的均質(zhì)化，在此區(qū)域內(nèi)共用全部或者部分模型參數(shù)，由此不可避免的引起代表性誤差。為了更準確反映出區(qū)域模型模擬的空間變異性，降低以點代面導致的誤差，近年來眾多學者將遙感觀測數(shù)據(jù)對作物模型進行約束，通過數(shù)據(jù)同化在空間上對模型參數(shù)進行優(yōu)化，取得了大量的成果[3-5,10-12]。然而，諸多研究中鮮有對模型參數(shù)標定的詳盡描述，也缺乏模型不確定性的定量表達，使得模型在更大區(qū)域上的應用可靠性受到質(zhì)疑。

所謂模型參數(shù)標定，就是調(diào)整模型參數(shù)使其適用于特定的研究區(qū)域。傳統(tǒng)方法是采用一定數(shù)量的實測值，通過試錯法標定出一套“最適”參數(shù)，這種標定方法很難準確代表空間變異性的結果的綜合。并且由于數(shù)據(jù)代表性、測量精度、標定方法等限制，常?；谝唤M觀測值會得到多組相異參數(shù)，即模型的“異參同效”現(xiàn)象。因此，無論是“以點代面”的區(qū)域模擬，還是更大范圍的模型同化應用，都迫切需要更為科學精確的參數(shù)標定方法以及標定后參數(shù)不確定性的定量評價方法。

基于貝葉斯理論的 MCMC（Markov chain Monte Carlo，馬爾科夫鏈-蒙特卡洛）方法能夠結合觀測數(shù)據(jù)和參數(shù)的先驗知識得到模型參數(shù)的后驗分布，從而利用參數(shù)后驗樣本的均值或中值作為參數(shù)標定結果，并且可以通過后驗樣本的均方根偏差定量的表達模型參數(shù)的不確定性。MCMC方法在水文模型[13-15]、森林管理模型[16]、林冠蒸騰模型[17]等模型參數(shù)標定研究中已得到較多應用，并且顯示出較強的適用性，其在作物模型參數(shù)標定與不確定性評價中也逐漸得到應用[18-20]。中國的作物模型應用面臨地形復雜、氣候多樣空間異質(zhì)性挑戰(zhàn)，因此模型參數(shù)也具有更大的不確定性。利用MCMC方法對中國小麥種植區(qū)作物模型參數(shù)較大范圍的標定、不確定性評價以及比較目前還少有報道，因此需要進一步研究探索。

本研究實現(xiàn)利用融入 snooker更新算法的 DE-MC（differential evolution Markov chain ，差分進化馬爾科夫鏈）方法在Python環(huán)境下對WOFOST作物生長模型參數(shù)進行標定和不確定性分析。標準的 DE-MC算法由Braak等[21]于2006年提出，并通過引入snooker更新算法進一步改進，研究表明融入snooker更新的DE-MC算法能夠通過少至3條的并行鏈實現(xiàn)對多達50～100維參數(shù)空間的后驗估計[22]。在河南鄭州農(nóng)業(yè)氣象站點上，結合多年實測LAI（leaf area index，葉面積指數(shù)）和產(chǎn)量數(shù)據(jù)得到參數(shù)后驗分布，實現(xiàn)對模型參數(shù)的估計和不確定性評價。

1 數(shù)據(jù)及模型

1.1 研究數(shù)據(jù)

本文所用的試驗觀測數(shù)據(jù)來自河南鄭州農(nóng)業(yè)氣象試驗站2008—2012年的觀測數(shù)據(jù)，其中2008—2011年數(shù)據(jù)用于模型參數(shù)標定，2012年數(shù)據(jù)用于標定結果驗證。研究區(qū)均位于種植區(qū)劃中的黃淮冬麥區(qū)，耕作模式主要為一年兩熟，小麥生育期降水約 120～140 mm，播種期至成熟期大于0 ℃積溫2 100～2 300 ℃，無霜期180～200 d。冬小麥播種期一般為10月份上中旬，成熟期一般為次年5月下旬到6月上旬。

作物數(shù)據(jù)主要包括從播種到成熟各個生育期的時間、生長狀況、管理措施、土壤情況、產(chǎn)量及產(chǎn)量結構等。其中，LAI觀測值分別對應三葉期、分蘗期、越冬始期、返青期、拔節(jié)期、抽穗期和乳熟期 7個生育時期；模型運行所需氣象數(shù)據(jù)來自中國科學院青藏高原研究所制作中國區(qū)域地面氣象要素數(shù)據(jù)集[23]，該數(shù)據(jù)是用于驅(qū)動陸面模型模擬的柵格數(shù)據(jù)集，包括近地面氣溫、近地面氣壓、近地面空氣比濕、近地面全風速、地面向下短波輻射、地面向下長波輻射與降水率共7個要素。

1.2 WOFOST模型

WOFOST模型是荷蘭瓦赫寧根大學開發(fā)的一個機理性過程模型，基于光截獲和二氧化碳同化作為驅(qū)動過程來描述作物生長，并將物候發(fā)育作為生長控制過程，分別用 0、1、2代表作物出苗、開花、成熟 3個生育期（development stage，DVS）。

本研究采用 Python下的 PCSE 5.3.1（python crop simulation environment，PCSE）軟件包進行WOFOST的運行。PCSE提供了作物模型模擬的環(huán)境、讀取數(shù)據(jù)（如氣象、土壤、農(nóng)藝管理等）的工具、模擬生物物理過程（如物候、呼吸、蒸散發(fā)等）的模塊。相比于傳統(tǒng)FORTRAN版本的作物模型，PCSE是用純Python代碼編寫，使其具有更加靈活、易于修改和拓展的特性，并可輕松連接數(shù)據(jù)庫、圖形用戶界面、可視化工具和數(shù)值統(tǒng)計軟件包。有關PCSE更詳細的介紹可以參考網(wǎng)站：https://pcse.readthedocs.org或者http://github.com/ajwdewit/pcse。

1.3 標定參數(shù)選擇

本研究利用LAI和產(chǎn)量觀測值對模型參數(shù)進行標定，與積溫有關的作物參數(shù)直接依據(jù)作物物候發(fā)育日期與近地面氣溫計算得到。其余作物參數(shù)通過 Sobol方法[24-25]進行全局敏感性分析，以Sobol全局敏感性指數(shù)大于0.05作為參數(shù)敏感的標準，選取對不同生育期LAI以及成熟期產(chǎn)量敏感的參數(shù)作為待標定參數(shù)。

2 研究方法

2.1 MCMC方法

貝葉斯理論能夠結合模型參數(shù)的先驗知識和模型輸出對應的觀測值，實現(xiàn)對模型參數(shù)的后驗估計。MCMC方法即在貝葉斯理論框架下，構造以參數(shù)后驗分布為平穩(wěn)分布的馬爾科夫鏈，從而得到參數(shù)的后驗樣本，基于這些樣本實現(xiàn)對參數(shù)數(shù)字特征的推斷。貝葉斯理論用公式表達如下

式中θ和y分別代表WOFOST模型的參數(shù)和模擬輸出值（如LAI、產(chǎn)量等）；p(θ/y)為參數(shù)的后驗概率密度函數(shù)；f(y/θ)為觀測值似然函數(shù)；g(θ)為參數(shù)的先驗分布。研究中，將模型的驅(qū)動數(shù)據(jù)，如氣象數(shù)據(jù)、降水等，記為x，則式（1）可以改為

引入MCMC方法的目的就是通過結合實測數(shù)據(jù)得到模型參數(shù)后驗樣本和并定量評價其不確定性，研究中，模型參數(shù)的不確定性用MCMC方法標定得到的參數(shù)后驗樣本的均方根偏差RMSDq計算

式中N表示參數(shù)θ后驗樣本個數(shù)；θi表示第i個后驗參數(shù)；表示參數(shù)θ后驗樣本均值，當θ服從高斯分布時，RMSDθ可以作為θ總體標準差的無偏估計。

為便于比較不同參數(shù)樣本內(nèi)部的離散程度，對每個參數(shù)后驗分別再計算相對均方根偏差RRMSDq，定義為

式中RMSDθ表示參數(shù)θ的均方根偏差，`表示參數(shù)θ后驗樣本均值。

2.2 DE-MC算法及snooker更新

差分進化馬爾科夫鏈(differential evolution Markov chain, DE-MC)算法[21]通過多條平行鏈的運行，以實現(xiàn)更好的探索參數(shù)空間。但標準DE-MC算法中并行鏈數(shù)必須大于參數(shù)空間的維數(shù)。為了提高采樣效率，Braak等[22]在提出DE-MC算法的基礎上，引入snooker更新算法[26-27]以部分取代其中的平行方向更新（parallel direction update），克服了標準DE-MC中并行鏈數(shù)必須大于維數(shù)的限制，發(fā)展成為DE-MC(Z)S算法。本文基于這種思想實現(xiàn)算法的構建，主要步驟為：

1）初始化M0′d維的矩陣Z，即依據(jù)先驗分布隨機生成M0組參數(shù)，其中M0>max(d,N),然后復制矩陣Z中前N行，記為種群（population）X，并設置M←M0；

2）設置一定的次數(shù)K，每逢K倍對X進行更新（update）。

3）將更新后的X添加到Z中，即M←M+N；

4）如果X已經(jīng)達到收斂（本研究通過比方差法（variance ratio）[28]進行判斷），或者總的更新次數(shù)超過了程序設定值，則進行到第5步，否則回到步驟2；

5）丟棄Z中一部分較早狀態(tài)的樣本（burn in），并進行一定比例的細化（thin），最終得到參數(shù)的后驗樣本。

其中每次對X的更新形成一個世代周期（generation cycle），并依次對鏈x1,x2,…,xN進行更新，首先隨機生成0到1之間的隨機數(shù)，如果小于0.1，則進行snooker更新，否則按照標準的 DE-MC算法進行平行方向更新[22]。snooker更新的步驟為：

1）選擇另外一條處于狀態(tài)z的鏈；

2）沿著xiz方向按如下步驟依概率密度采樣，

① 選擇另外2條分別處于狀態(tài)zR1和zR2的隨機鏈R1和R2，

② 將zR1和zR2正射投影到xiz方向，并分別生成zP1和zP2，

③ 提議x*←xi+γs(zP1-zP1)，其中γs為縮放因子，取1.2～2.2之間的隨機數(shù)，

④計算Metropolis系數(shù)：

⑤ 以概率min(1,r)接受提議，否則保持xi狀態(tài)不變。

3 結果與分析

3.1 Sobol敏感性分析結果

Sobol全局敏感性分析結果表明，對三葉期LAI敏感的參數(shù)包括SLATB00（specific leaf area at development stage of 0，生育期為0時的比葉面積）、TBASEM（lower threshold temperature for emergence，出苗的低溫閾值）、TDWI（initial total crop dry weight，作物初始干物質(zhì)質(zhì)量）、SLATB020（specific leaf area at development stage of 0.2，生育期為0.2時的比葉面積）；對分蘗期LAI敏感的參數(shù)包括SLATB00、SLATB020、TBASEM、AMAXTB00（maximum CO2assimilation rate at development stage of 0，生育期為 0時的單葉最大二氧化碳同化速率）；對越冬期 LAI敏感的參數(shù)包括 SLATB00、SLATB020、TBASEM、AMAXTB00、TBASE（lower threshold temperature for aging of leaves，葉齡的低溫閾值）、TMNFTB_min（low temperature threshold when reduction factor of gross assimilation rate is 0，總同化率降低因子為0時的低溫閾值）；對返青期 LAI敏感的參數(shù)包括SLATB00、SLATB020、TMNFTB_min、TBASEM、TBASE、AMAXTB00；對拔節(jié)期 LAI敏感的參數(shù)包括SLATB040（specific leaf area at development stage of 0.4，生育期為 0.4時的比葉面積）、SLATB00、TBASE、SLATB020、AMAXTB00、 SLATB070（specific leaf area at development stage of 0.7，生育期為0.7時的比葉面積）、TMNFTB_min、RGRLAI（maximum relative increase in LAI，葉面積指數(shù)最大增長率）；對抽穗期LAI敏感的參數(shù)包括 SLATB070、SLATB040、SLATB00、TBASE、KDIFTB100（extinction coefficient for diffuse visible light at development stage of 1，生育期為1時的散射消光系數(shù)）、AMAXTB00、SLATB020、AMAXTB100（maximum CO2assimilation rate at development stage of 1，生育期為1時的單葉最大二氧化碳同化速率）；對乳熟期LAI敏感的參數(shù)包括SPAN（life span of leaves growing at 35 Celsius，在35 ℃時葉面積的生命周期）、SLATB070、KDIFTB100；對成熟期產(chǎn)量敏感的參數(shù) SPAN、CVO（efficiency of conversion into storage organs，儲存器官的同化物轉(zhuǎn)換效率）、AMAXTB130（maximum CO2assimilation rate at development stage of 1.3，生育期為1.3時的單葉最大二氧化碳同化速率）、SLATB070、AMAXTB100、KDIFTB100、SLATB040。

上述敏感性分析是以2008年鄭州農(nóng)業(yè)氣象試驗站上數(shù)據(jù)運行分析得到的，由于TBASEM的改變會引起生育期在時間上的“漂移”，因而設置該參數(shù)為默認值，不進行標定。最終綜合敏感性分析結果，選取 SLATB00、SLATB020、AMAXTB00、SLATB070、TBASE、AMAXTB100、KDIFTB100、SLATB040、TMNFTB_min、SPAN、CVO、AMAXTB130、RGRLAI、TDWI共計 14個參數(shù)作為待標定參數(shù)。

3.2 MCMC標定

3.2.1 參數(shù)先驗分布

對14個參數(shù)設置為以模型默認值為初始值、在參數(shù)范圍內(nèi)的均勻分布。參數(shù)初始值及范圍如表1。

表1 待優(yōu)化的14個參數(shù)的初始值及范圍Table 1 Initial value and range of 14 parameters for optimization

3.2.2 似然函數(shù)的確定

本研究中認為LAI和產(chǎn)量服從以觀測值為期望的高斯分布，由于概率密度可能非常小，為了避免位數(shù)舍入引起的誤差，所有概率密度的計算均通過取對數(shù)形式計算，由此確立似然函數(shù)

式中LLAI和Lyield分別表示LAI和產(chǎn)量的似然函數(shù)；L表示似然函數(shù)（likelihood function）的返回值；矢量x和xobs分別表示對應不同生育期LAI的模型模擬值和觀測值；Σ表示LAI觀測值的協(xié)方差矩陣，本研究中不同生育期LAI觀測值相互獨立，且方差均取0.2；K表示空間維數(shù)，即LAI觀測值的個數(shù)；detΣ表示 Σ的行列式值；Y和Yobs分別表示成熟期產(chǎn)量的模型模擬值和觀測值；σ表示產(chǎn)量觀測值的標準差，本研究中取觀測值的10%。

圖1 14個參數(shù)的先驗及后驗分布Fig.1 Prior and posterior distribution of 14 parameters

3.2.3 參數(shù)的后驗分布

本研究中設定DE-MC種群（并行鏈）數(shù)目為3，每5次采樣進行一次種群進化（鏈的更新），每5 000次迭代通過比方差法進行一次收斂性判斷，計算診斷指標?R[28]，?R>1時，表明鏈沒有達到收斂，若?R≈1，表明鏈已經(jīng)處于靜止狀態(tài)。本研究中當連續(xù)超過3次診斷指標?R<1.05，則認為馬爾科夫鏈達到收斂。收斂后舍棄前2000次采樣（burn-in），為解決自相關問題，并每隔5次對鏈進行細化，最終得到參數(shù)的后驗樣本及其分布。圖1為14個參數(shù)先驗及后驗樣本的分布圖。表2列出了最終得到參數(shù)后驗樣本的均值、中值、95%置信區(qū)間、RMSD以及 RRMSD，結果均保留小數(shù)點后4位。

從圖 1可以看出，相對于各參數(shù)的先驗分布，參數(shù)的后驗分布均出現(xiàn)較大改變，表明標定參數(shù)選取的有效性。由于WOFOST模型無法準確模擬越冬期小麥實際生長情況，因而在標定參數(shù)包含越冬期前后LAI觀測值的情況下，部分參數(shù)的后驗分布表現(xiàn)出極端變化，表現(xiàn)為：SLATB020集中分布在參數(shù)區(qū)間的左邊界，這是因為模型運行時，比葉面積與LAI模擬值成正比例關系，生育期為0.2較為接近越冬期，由于WOFOST無法模擬低溫導致的葉片損傷減少，其LAI“模擬期望值”很大程度上高于實際LAI值，因而在實際LAI觀測值的標定下，導致該時期比葉面積極端減小以使模型模擬出較小的LAI值。同理，生育期為 1時大致處于抽穗期附近，此階段 LAI值基本處于整個生育期的峰值，因而在抽穗期LAI峰值的影響下，AMAXTB100和KDIFTB100集中分布在參數(shù)區(qū)間的右邊界。

除個別極端分布的參數(shù)，表 2中參數(shù)的后驗樣本均值和中值基本相一致。后驗參數(shù)的 95%置信區(qū)間一定程度上反映了參數(shù)的合理取值區(qū)間，與表 1中參數(shù)的初始范圍相比，其范圍縮小的越多，也一定程度表明該參數(shù)標定后的不確定越小。RMSD和RRMSD描述了參數(shù)樣本內(nèi)部相對于均值的離散程度，可作為參數(shù)不確定性的定量評價。在參數(shù)符合高斯分布時，參數(shù)樣本的 RMSD是參數(shù)總體高斯分布的標準差σ的無偏估計。由圖1分布，可近似認為 SLATB00、SLATB070、SLATB040、SPAN、AMAXTB130符合高斯分布，以RRMSD為指標，不確定性由小到大依次為 SPAN、SLATB070、SLATB040、AMAXTB130、SLATB00。

表2 參數(shù)后驗樣本的均值、中值、95%置信區(qū)間、RMSD、RRMSDTable 2 Mean, median, 95% confidence interval, RMSD and RRMSD of parameters posterior sample

3.3 標定運行結果

3.3.1 標定精度驗證

以2012年生長季站點上氣象數(shù)據(jù)為驅(qū)動，分別以參數(shù)初始值、后驗均值和中值帶入模型運行，與觀測值進行比較，結果如圖2所示。計算LAI模擬值和對應觀測值間的均方根誤差（RMSE）作為精度評價指標，公式為

式中LAIi表示第i個LAI模擬值，LAIiobs表示第i個對應的觀測值，n表示LAI觀測值個數(shù)。結果顯示：參數(shù)初始值（模型默認參數(shù)）、后驗均值和后驗中值帶入模型模擬的LAI值的RMSE分別為1.79、0.84、0.87，相比于初始值，標定后LAI模擬精度分別提高了53.07%（后驗均值）和51.40%（后驗中值）；產(chǎn)量模擬精度從初始值的75.15%提高至 83.40%（后驗均值）和 84.03%（后驗中值），模擬精度提高8.25%~8.88%。這表明利用MCMC標定得到的后驗樣本均值和中值在很大程度能夠?qū)崿F(xiàn)對觀測數(shù)據(jù)的擬合。此外，由圖 2可以發(fā)現(xiàn)，相比于初始狀態(tài)，標定后并非在每個觀測節(jié)點上都提高了精度，但是卻在整體上提高了與觀測數(shù)據(jù)的擬合精度。標定后的LAI曲線較好的捕捉到了LAI的峰值，但是在越冬期前后的LAI模擬值具有較大誤差，主要是由于WOFOST模型開發(fā)之初主要是為了模擬歐洲的作物生長，并未實現(xiàn)中國冬小麥越冬期情況模擬[29]。

標定后的產(chǎn)量模擬精度較初試狀態(tài)有所提高，但相對觀測值仍有一定誤差，這很大程度上是因為用于標定的2008—2011年的產(chǎn)量數(shù)據(jù)僅僅是4個標量值，并不能夠充分實現(xiàn)對產(chǎn)量形成相關參數(shù)的準確標定，進一步研究可以引入生育期內(nèi)的穗重數(shù)據(jù)進行標定，以提高產(chǎn)量標定精度。

圖2 參數(shù)標定前后模型模擬結果對比Fig.2 Comparison of model simulated results before and after parameters calibration

3.3.2 模型輸出的不確定性評價

將后驗樣本作為參數(shù)集合，分別帶入WOFOST模型并輸出得到LAI和產(chǎn)量時間序列上的模擬集合，結果如圖 3所示。由陰影部分的走勢可以直觀的反映出后驗參數(shù)集模擬結果所代表的不確定性，整體上LAI模擬的不確定性在三葉期至返青期之間以及LAI曲線峰值（拔節(jié)期至抽穗期之間）較大，其它時期相對穩(wěn)定。產(chǎn)量模擬的不確定性在開始模型模擬籽粒形成開始逐漸擴大，并在乳熟期達到最大，至成熟期基本趨于穩(wěn)定。

圖3 參數(shù)標定后模型模擬LAI和產(chǎn)量的不確定性Fig.3 Uncertainty of model simulated LAI and yield with calibrated parameters dataset

4 結 論

本研究以鄭州農(nóng)業(yè)氣象試驗站2008—2011年觀測數(shù)據(jù)作為標定數(shù)據(jù)，選取14個敏感參數(shù)，運用融入snooker更新的DE-MC方法實現(xiàn)參數(shù)的標定，并用2012年的觀測數(shù)據(jù)進行了驗證，得到結論如下：

1）以模型默認參數(shù)為初始值，通過DE-MC方法可以實現(xiàn)參數(shù)的自動標定，相比于模型默認參數(shù)模擬結果，參數(shù)標定后LAI模擬精度可提高51.40%~53.07%，產(chǎn)量模擬精度提高8.25%~8.88%。

2）以生育期內(nèi)LAI觀測值和成熟期產(chǎn)量為標定數(shù)據(jù)，SPAN、SLATB070、SLATB040、AMAXTB130 和 SLATB00的后驗分布可近似為高斯分布，其中 SPAN的不確定性最低。

3）帶入后驗參數(shù)集合進行模型，LAI在三葉期至返青期之間以及拔節(jié)期至抽穗期之間模擬的不確定性較大；產(chǎn)量模擬的不確定性隨時間不斷增大，至乳熟期前后達到穩(wěn)定。

本研究使用的標定數(shù)據(jù)為農(nóng)業(yè)氣象站點上的觀測數(shù)據(jù)，雖然作物品種一致，但年際間的栽培管理措施并未嚴格進行控制（如播種密度等），導致標定結果的精度可能受到一定影響。后續(xù)可在更多地區(qū)，采用具有定量控制的試驗數(shù)據(jù)開展研究，以進一步驗證本方法的有效性。

[1]Wang H, Zhu Y, Li W, et al. Integrating remotely sensed leaf area index and leaf nitrogen accumulation with RiceGrow model based on particle swarm optimization algorithm for rice grain yield assessment[J]. Journal of Applied Remote Sensing, 2014, 8(1): 83674.

[2]Huang J, Tian L, Liang S, et al. Improving winter wheat yield estimation by assimilation of the leaf area index from Landsat TM and MODIS data into the WOFOST model[J].Agricultural & Forest Meteorology, 2015, 204: 106－121.

[3]He B, Li X, Quan X, et al. Estimating the aboveground dry biomass of grass by assimilation of retrieved LAI into a crop growth model[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2015, 8(2): 550－561.

[4]王東偉，王錦地，梁順林. 作物生長模型同化 MODIS反射率方法提取作物葉面積指數(shù)[J]. 中國科學：地球科學，2010(1)：73－83.Wang Dongwei, Wang Jindi, Liang Shunlin. Assimilating MODIS reflectance into crop growth model to extract crop leaf area index [J]. Scientia Sinica Terrae, 2010(1): 73－83.(in Chinese with English abstract)

[5]Fang H, Liang S, Hoogenboom G, et al. Corn-yield estimation through assimilation of remotely sensed data into the CSM-CERES-Maize model[J]. International Journal of Remote Sensing, 2008, 29(10): 3011－3032.

[6]張建平，趙艷霞，王春乙，等. 基于 WOFOST作物生長模型的冬小麥干旱影響評估技術[J]. 生態(tài)學報，2013，33(6)：1762－1769.Zhang Jianping, Zhao Yanxia, Wang Chunyi, et al. Evaluation technology on drought disaster to yields of winter wheat based on WOFOST crop growth model[J]. Acta Ecologica Sinica, 2013, 33(6): 1762－1769. (in Chinese with English abstract)

[7]馬玉平，王石立，李維京. 基于作物生長模型的東北玉米冷害監(jiān)測預測[J]. 作物學報，2011，37(10)：1868－1878.Ma Yuping, Wang Shili, Li Weijing. Monitoring and predicting of maize chilling damage based on crop growth model in Northeast China[J]. Acta Agronomica Sinica, 2011,37(10): 1868－1878. (in Chinese with English abstract)

[8]Jin M, Liu X, Wu L, et al. An improved assimilation method with stress factors incorporated in the WOFOST model for the efficient assessment of heavy metal stress levels in rice[J].International Journal of Applied Earth Observations &Geoinformation, 2015, 41: 118－129.

[9]李振海. 基于遙感數(shù)據(jù)和氣象預報數(shù)據(jù)的 DSSAT模型冬小麥產(chǎn)量和品質(zhì)預報[D]. 杭州：浙江大學，2016.Li Zhenhai. Predicting Winter Wheat Yield and Quality by Integrating of Remote-sensing Data and the Weather Forecast Data Into the DSSAT Model[D]. Hangzhou: Zhejiang University, 2016. (in Chinese with English abstract)

[10]Huang J, Sedano F, Huang Y, et al. Assimilating a synthetic Kalman filter leaf area index series into the WOFOST model to improve regional winter wheat yield estimation[J].Agricultural & Forest Meteorology, 2016, 216: 188－202.

[11]Ma H, Huang J, Zhu D, et al. Estimating regional winter wheat yield by assimilation of time series of HJ-1 CCD NDVI into WOFOST–ACRM model with Ensemble Kalman Filter[J].Mathematical & Computer Modelling, 2013, 58(3-4):753－764.

[12]Morel J, Todoroff P, Bégué A, et al. Toward a satellite-based system of sugarcane yield estimation and forecasting in smallholder farming conditions: A case study on Reunion Island[J]. Remote Sensing, 2014, 6(7): 793708.

[13]Rwasoka D T, Madamombe C E, Gumindoga W, et al.Calibration, validation, parameter indentifiability and uncertainty analysis of a 2 – parameter parsimonious monthly rainfall-runoff model in two catchments in Zimbabwe[J].Physics & Chemistry of the Earth Parts A/b/c, 2014,67–69(12): 36－46.

[14]Hassan A E, Bekhit H M, Chapman J B. Using Markov Chain Monte Carlo to Quantify Parameter Uncertainty and Its Effect on Predictions of a Groundwater Flow Model.[M].Elsevier Science Publishers B. V., 2009.

[15]Gallagher M, Doherty J. Parameter estimation and uncertainty analysis for a watershed model[J]. Environmental Modelling & Software, 2007, 22(7): 1000－1020.

[16]Van O 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 &Management, 2013, 289(2): 255－268.

[17]Li X, Yang P, Ren S, et al. An improved canopy transpiration model and parameter uncertainty analysis by Bayesian approach[J]. Mathematical & Computer Modelling An International Journal, 2010, 51(11, 12): 1368－1374.

[18]Toshichika I, Masayuki Y, Motoki N. Parameter estimation and uncertainty analysis of a large-scale crop model for paddy rice: Application of a Bayesian approach[J].Agricultural & Forest Meteorology, 2009, 149(2): 333－348.

[19]Marin F, Jones J W, Boote K J. A stochastic method for crop models: including uncertainty in a sugarcane model[J].Agronomy Journal, 2017, 109(2).

[20]何亮，侯英雨，趙剛，等. 基于全局敏感性分析和貝葉斯方法的 WOFOST作物模型參數(shù)優(yōu)化[J]. 農(nóng)業(yè)工程學報，2016，32(2)：169－179.He Liang, Hou Yingyu, Zhao Gang, et al. 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)

[21]Braak C J F T. A Markov chain Monte Carlo version of the genetic algorithm differential evolution: Easy Bayesian computing for real parameter spaces[J]. Statistics &Computing, 2006, 16(3): 239－249.

[22]Braak C J T, Vrugt J A. Differential evolution Markov chain with snooker updater and fewer chains[J]. Statistics &Computing, 2008, 18(4): 435－446.

[23]何杰，陽坤. 中國區(qū)域高時空分辨率地面氣象要素驅(qū)動數(shù)據(jù)集[Z]. 寒區(qū)旱區(qū)科學數(shù)據(jù)中心，2011.

[24]Sobol, I. M. Global Sensitivity Indices for Nonlinear Mathematical Models and Their Monte Carlo Estimates[M].Elsevier Science Publishers B. V., 2001.

[25]Zhang X Y, Trame M N, Lesko L J, et al. Sobol sensitivity analysis: A tool to guide the development and evaluation of systems pharmacology models[J]. Cpt Pharmacometrics &Systems Pharmacology, 2015, 4(2): 69.

[26]Gilks W R, Roberts G O, George E I. Adaptive direction sampling[J]. Journal of the Royal Statistical Society, 1994,43(1): 179－189.

[27]Liang F, Wong W H. Real-parameter evolutionary Monte Carlo with applications to Bayesian mixture models[J].Publications of the American Statistical Association, 2001,96(454): 653－666.

[28]Brooks S P, Roberts G O. Convergence assessment techniques for Markov chain Monte Carlo[J]. Statistics &Computing, 1998, 8(4): 319－335.

[29]馬玉平，王石立，張黎. 針對華北小麥越冬的 WOFOST模型改進[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)