不同參數(shù)取值范圍下CROPGRO-cotton模型全局敏感性和不確定性分析

高浚 周保平 王昱 于晗 王君

摘要：在對新作物模擬的過程中，需要對品種參數(shù)重新進行率定，但是很多參數(shù)難以測量和校正，所以需要確定參數(shù)的取值范圍。目前，參數(shù)的取值范圍大部分通過2種方式得到：（1）實測數(shù)據(jù)和前人文獻;（2）參數(shù)初始默認(rèn)值的某個百分比上下變動。前者因為很多實測數(shù)據(jù)和率定后的參數(shù)未公布，所以選擇后者確定參數(shù)取值范圍。在5個參數(shù)取值范圍（10%、20%、30%、40%、50%）下，通過Simlab軟件，采取擴展傅里葉幅度敏感性檢驗（EFAST）法，對不同參數(shù)取值范圍的輸出變量進行敏感性和不確定性分析可知，30%的參數(shù)取值范圍能很好地模擬作物生長。敏感性分析與不確定性分析標(biāo)明其輸出變量成熟期（MDAP）、產(chǎn)量（HWAM）、生物量（CWAM）的一致性檢驗系數(shù)（TDCC）分別為0.86、0.9、0.88，P值（顯著性分析）均小于0.05，具有顯著一致性，開花期（ADAP）的TDCC較小，為0.66，這因為只有參數(shù)EM-FL為敏感性參數(shù)，其余參數(shù)均不敏感。模型的輸出結(jié)果中ADAP、MDAP、HWAM的實測值均處于95%置信區(qū)間內(nèi)，CWAM略差。本研究可幫助CROPGRO-cotton模型在南疆地區(qū)的應(yīng)用提升模擬精度和模擬效果。

關(guān)鍵詞：CROPGRO-cotton模型;敏感性分析;不確定性分析;參數(shù)取值范圍;EFAST法

中圖分類號： S126? 文獻標(biāo)志碼： A

文章編號：1002-1302（2022）09-0195-07

作物模型運用數(shù)學(xué)模型把天氣、土壤、田間管理等體系組合成一個整體，運用計算機定量計算實現(xiàn)動態(tài)模擬[1]，了解作物的生長情況。因為生長過程采用的是數(shù)學(xué)公式表達(dá)，所以模型輸出的結(jié)果具有一定的不確定性。不確定性分析（UA）用來幫助量化試驗數(shù)據(jù)和參數(shù)估值的可信度，當(dāng)前運用概率分布函數(shù)來說明不確定性。敏感性分析（SA）用來量化某個不確定參數(shù)對模型輸出的作用。在模型應(yīng)用前，都要對模型參數(shù)進行率定，實現(xiàn)模型的本地化，參數(shù)的率定結(jié)果能很大程度上影響模型模擬的結(jié)果，所以，不確定性分析和敏感性分析非常重要。

當(dāng)今，作物模型探索的一個關(guān)鍵點就是對模型品種參數(shù)、土壤參數(shù)、田間管理參數(shù)等進行敏感性分析，國內(nèi)外有不少研究專家運用了各種試驗方案去進行探索。例如，不同灌溉、不同氣候等方法對ALMANAC模型[2-3]、WOFOST模型[4-6]、COTTON2K模型[7-9]、農(nóng)業(yè)技術(shù)轉(zhuǎn)移決策支持系統(tǒng)（DSSAT）模型[10-12]等模型實現(xiàn)不確定性分析和敏感性分析。但是，目前對于不同參數(shù)取值范圍下CROPGRO-cotton模型[13-15]敏感性分析和不確定性分析還鮮有人研究。

本研究主要研究南疆地區(qū)的棉花，實測數(shù)據(jù)來自于南疆阿克蘇地區(qū)的天氣、土壤和大田數(shù)據(jù)。運用擴展傅里葉幅度敏感性檢驗（EFAST）法[16-17]對模型的品種參數(shù)進行敏感性分析，然后用EFAST法的輸出結(jié)果進行不確定性分析，比較不同參數(shù)取值范圍下模擬精度的優(yōu)劣性，這為模型本地化參數(shù)的率定取值提供一定的理論和技術(shù)支持。

1 材料與方法

1.1 CROPGRO-cotton模型

DSSAT模型中包含CROPGRO模型，此模型囊括了許多作物模型，如番茄、玉米、水稻、棉花模型等。這些作物的子模型構(gòu)造功能基本一樣，主要模塊都是田間管理模塊、氣象參數(shù)、土壤參數(shù)和作物品種。田間管理模塊主要控制各種管理方式，如種植時間、種植方式、灌溉、施肥、收獲等。氣象參數(shù)主要包括每天太陽輻射量、每天最高氣溫和最低氣溫、風(fēng)速、降水量、濕度等。土壤參數(shù)主要包括土壤成分、土壤有機質(zhì)含量、土壤水分等。CROPGRO模型作物的模擬涉及到三大類參數(shù)：品種參數(shù)、物種參數(shù)、生態(tài)型參數(shù)。其中，品種參數(shù)的敏感性分析是重中之重。

1.2 田間試驗及數(shù)據(jù)

CROPGRO-cotton模型模擬所需要的田間管理及天氣數(shù)據(jù)來自于新疆生產(chǎn)建設(shè)兵團第一師阿拉爾水利局十團灌溉試驗站棉花田，位于81.31°E、40.13°N，水源充足，年均日照2 680 h，最高氣溫 42 ℃，最低氣溫-29 ℃，年平均氣溫14 ℃，年平均降水量為51 mm，年蒸發(fā)量2 684 mm。試驗區(qū)的天氣數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)，田間管理數(shù)據(jù)來源于大田觀測，土壤剖面數(shù)據(jù)來源于實際勘測，土壤的基本理化性質(zhì)見表1。

本研究的灌溉數(shù)據(jù)來源于大田試驗，試驗地根據(jù)棉花生育期內(nèi)所需的水分不同，通過輪灌工作制度，共設(shè)計3個灌水處理，具體的灌水時間及灌水定額見表2。

1.3 模型輸入?yún)?shù)與輸出變量

本研究選用了作物棉花的18個品種參數(shù)作為棉花的輸入?yún)?shù)，參數(shù)值為DSSAT模型Gentype文件夾下棉花作物的默認(rèn)初始值，取值范圍采用DSSAT模型的參數(shù)取值范圍，具體的參數(shù)情況見表3。

1.4 擴展傅里葉檢驗法

EFAST法是一種基于方差分解的敏感性分析方法。它通過傅里葉變換得到傅里葉級數(shù)的頻譜曲線，然后用曲線獲得所有參數(shù)以及參數(shù)相互作用所引起的模型方差，其計算公式為：

V（y）=∑ni=1V i+∑nij≤1V ij+…+V（1，2，…，n）;（1）

V ij=V[E（y|x ix j）]-V i-V j;（2）

S i=V iV（y）;（3）

S τi=V（y）-V -iV（y）。（4）

式中：V（y）為所有變量總方差;V i為第i個變量的方差;V j為第j個變量的方差;E（y|x ix j）為對于確定變量x i所對應(yīng)的y的期望;V ij為變量間相互作用的方差;V -i為去除變量后的方差之和;S i為變量x i的一階敏感性指數(shù);S τi為變量x i的全局敏感性指數(shù)。

1.5 敏感性與不確定性分析方法

采用Simlab軟件[18]輸入?yún)?shù)名稱、參數(shù)分布、參數(shù)取值范圍、敏感性分析方法和抽樣次數(shù)。首先，本研究設(shè)置了5種處理，EFAST法共抽樣56 000次，參數(shù)取值范圍分別為參數(shù)默認(rèn)值（1±10%）、參數(shù)默認(rèn)值（1±20%）、參數(shù)默認(rèn)值（1±30%）、參數(shù)默認(rèn)值（1±40%）、參數(shù)默認(rèn)值（1±50%），分別記為10%、20%、30%、40%、50%，對應(yīng)生成的參數(shù)樣本數(shù)量分別為3 600、7 200、10 800、14 400、18 000個。其次，用R語言編寫的批處理程序，進行DSSAT參數(shù)的修改、調(diào)用和輸出。模型的參數(shù)取值不同時，模型的輸出結(jié)果也會變化，具備了一定的不確定性，運用均值、標(biāo)準(zhǔn)差、變異系數(shù)對模型的輸出結(jié)果進行評價。最后，采用一致性檢驗系數(shù)（TDCC）[19-20]來衡量EFAST法參數(shù)敏感性排序。

2 結(jié)果與分析

2.1 參數(shù)取值范圍對輸出結(jié)果影響

按照對應(yīng)的參數(shù)取值范圍，選擇輸出結(jié)果中開花期（ADAP）、成熟期（MDAP）、產(chǎn)量（HWAM）、生物量（CWAM）這4個輸出變量，隨后分別計算它們的變異系數(shù)（CV），其變化趨勢見圖1。

由圖1可知，參數(shù)取值范圍為10%、20%、30%、40%、50%時，ADAP的CV變化范圍為4.0%～20.1%，MDAP的CV變化范圍為5.0%～14.2%，HWAM的CV變動范圍為10.9%～56.6%，CWAM的CV變化范圍為8.7%～27.0%。當(dāng)參數(shù)取值范圍變大時，輸出變量ADAP、MDAP、HWAM、CWAM的CV變大，呈線性相關(guān)關(guān)系。

2.2 參數(shù)取值范圍對EFAST法敏感性分析的影響

用EFAST法對不同參數(shù)取值范圍下的輸出變量進行敏感性分析，然后通過一致性檢驗系數(shù)（TDCC）對不同參數(shù)取值范圍下參數(shù)敏感性排序進行評價，結(jié)果見表4。由表4可知，不同參數(shù)取值范圍的P值都小于0.05，MDAP、HWAM、CWAM的TDCC都較大，但是ADAP的TDCC較小，因為只有參數(shù)EM-FL為敏感性參數(shù)，其余參數(shù)均不敏感，所以導(dǎo)致TDCC計算有一定的誤差。說明不同參數(shù)取值范圍下參數(shù)敏感性排序無明顯差異。

P值<0.05、TDCC較大， 只能表明參數(shù)在敏感性排序的統(tǒng)計意義方面有一致性，各個參數(shù)的重要程度略有不同。為了解釋該問題，用EFAST法把所有敏感性指數(shù)>0.1的作為敏感性參數(shù)進行排序，具體結(jié)果見表5。

由表5可知，不管參數(shù)取值范圍如何變化，在輸出變量ADAP下，只有EM-FL對開花期最敏感，主要是因為現(xiàn)在有很多棉花品種已發(fā)展成光周期鈍化品種，只有出苗到初光期的光熱時間的變化會對開花產(chǎn)生影響。在輸出變量MDAP下，SD-PM、EM-FL、FL-SD屬于極敏感參數(shù)，不同參數(shù)取值范圍下，排名基本無變化。在輸出變量HWAM下，THRSH、XFRT、EM-FL屬于極敏感參數(shù)，不同參數(shù)取值范圍下，排名變化幅度很小。在輸出變量CWAM下，SD-PM、EM-FL、LFMAX屬于極敏感參數(shù)，不同參數(shù)取值范圍下，排名稍微有些變化。CSDL和PPSEN排名變化幅度較小，不管在哪個參數(shù)范圍下都屬于極不敏感參數(shù)。

為方便研究參數(shù)敏感性在不同參數(shù)取值范圍下的變化情況，把輸出變量ADAP、MDAP、HWAM、CWAM中的極敏感參數(shù)提取出來，探討它們與其取值范圍之間的變動趨勢，詳細(xì)趨勢見圖2。

對ADAP來說，圖2-j表明，參數(shù)EM-FL的敏感性指數(shù)在整個參數(shù)取值范圍內(nèi)均無明顯變化。

對MDAP來說，圖2-a表明，參數(shù)SDPM取值范圍在50%時，敏感性指數(shù)最大，在10%～30%時，敏感性指數(shù)先增大再減小，在30%～50%時，敏感性指數(shù)逐漸變大;圖2-b表明，參數(shù)EM-FL隨取值范圍的增大而緩慢減小;圖2-c表明，參數(shù)FL-SD在10%時，敏感性指數(shù)最大，在20%時敏感性指數(shù)最小。

對HWAM來說，圖2-d表明，參數(shù)THRSH在10%～20%時基本無變化，在20%～50%時，參數(shù)THRSH敏感性指數(shù)隨取值范圍的增大而減小;圖2-e 表明，參數(shù)XFRT在30%時，敏感性指數(shù)最大，在40%時敏感性指數(shù)最小;圖2-f表明，參數(shù)EM-FL在20%時，敏感性指數(shù)最小，在40%時，敏感性指數(shù)最大。

對CWAM來說，圖2-g表明，敏感性指數(shù)隨參數(shù)SDPM隨取值范圍的增大而減小，分別在10%、50%時最大、最小;圖2-h表明，參數(shù)EM-FL在10%時，敏感性指數(shù)最大，在20%～50%時，變化趨勢比較穩(wěn)定;圖2-i表明，參數(shù)LFMAX在30%時，敏感性指數(shù)會略微下降， 但在參數(shù)LFMAX的整個取值范圍內(nèi)總體趨勢是上升的，在10%、50%時敏感性指數(shù)最小、最大。

2.3 參數(shù)取值范圍對EFAST法不確定性分析影響

不同參數(shù)取值范圍會對輸出結(jié)果有不同的影響。當(dāng)參數(shù)取值范圍過小時，輸出變量就不容易呈現(xiàn)出實測值的變化范圍;當(dāng)參數(shù)取值范圍過大時，輸出變量就容易和實測值不符。所以，就要通過試驗篩選出一個合適的參數(shù)取值范圍。本研究通過大量試驗得出，當(dāng)參數(shù)取值為默認(rèn)值的30%時，模型的輸出結(jié)果基本都分布在95%置信區(qū)間[21-22]內(nèi)，模型整體模擬結(jié)果擁有比較大的可信度。模型的不確定性分析通過累計分布函數(shù)[23-24]來進行量化，具體見圖3和表6。

由表6可知，不同參數(shù)取值范圍下各灌水處理與不確定分析結(jié)果無關(guān)。當(dāng)參數(shù)取值范圍為10%時，模型的輸出結(jié)果中ADAP、CWAM的實測值均不處于95%置信區(qū)間內(nèi);當(dāng)參數(shù)取值范圍為20%時，模型的輸出結(jié)果中ADAP、CWAM的實測值均不處于95%置信區(qū)間內(nèi);當(dāng)參數(shù)取值范圍為30%時，模型的輸出結(jié)果中CWAM的實測值不處于95%置信區(qū)間內(nèi);當(dāng)參數(shù)取值范圍為40%時，模型的輸出結(jié)果中MDAP、CWAM的實測值均不處于95%置信區(qū)間內(nèi);當(dāng)參數(shù)取值范圍為50%時，模型的輸出結(jié)果中MDAP、CWAM的實測值均不處于95%置信區(qū)間內(nèi)。

通過不確定性分析輸出結(jié)果可得，在10%時，ADAP的95%置信區(qū)間為64～73 d，出現(xiàn)頻率最高的時間是69 d，出現(xiàn)了416次，跟實測值相差較大;在20%時ADAP的95%置信區(qū)間為59～77 d，出現(xiàn)頻率最高的時間是76 d，出現(xiàn)了420次，跟實測值比較接近;在30%時，ADAP的95%置信區(qū)間為54～82 d，出現(xiàn)頻率最高的時間是76 d，出現(xiàn)了416次，跟實測值比較接近;在40%時，ADAP的95%置信區(qū)間為49～86 d，出現(xiàn)頻率最高的時間是50 d，出現(xiàn)了460次，跟實測值相差較大;在50%時，ADAP的95%置信區(qū)間為45～90 d，出現(xiàn)頻率最高的時間是50 d，出現(xiàn)了455次，跟實測值相差較大。MDAP的95%置信區(qū)間下限隨參數(shù)取值范圍的變大而減小，95%置信區(qū)間上限隨參數(shù)取值范圍的變大基本無變化。HWAM的95%置信區(qū)間上限隨參數(shù)取值范圍的變大而總體呈變大趨勢，取值范圍為10%～40%時95%置信區(qū)間上限隨參數(shù)取值范圍的變大而變大， 在50%時下降了一些。CWAM的95%置信區(qū)間上限隨參數(shù)取值范圍的變大而總體呈變大趨勢，取值范圍為10%～40%時95%置信區(qū)間上限隨參數(shù)取值范圍的變大而變大，在50%時同樣下降了一些，說明參數(shù)取值范圍在50%時，模擬情況會有比較大的誤差。另外，當(dāng)參數(shù)取值范圍為10%時，模型模擬的MDAP有689次為0的情況;當(dāng)參數(shù)取值范圍為20%時，模型模擬的MDAP有2 340次為0的情況;當(dāng)參數(shù)取值范圍為30%時，模型模擬的MDAP有4 212次為0的情況（圖3）;當(dāng)參數(shù)取值范圍為40%時，模型模擬的MDAP有6 157次為0的情況，HWAM有451次為0的情況;當(dāng)參數(shù)取值范圍為50%時，模型模擬的MDAP有8 099次為0的情況，HWAM有 1 540 次為0的情況。由此可得，取值范圍為30%時，模擬的效果最好，但DSSAT模型對HWAM的模擬效果略差。

3 討論與結(jié)論

由圖1可知，輸出變量的CV會受到參數(shù)取值范圍的影響，并且同一輸出變量的CV在不同參數(shù)取值范圍內(nèi)呈線性增加。但是HWAM的CV過大，這可能跟參數(shù)取值范圍過大有關(guān)，某些參數(shù)的取值超過了它的極限。

由表3可知，P值和TDCC顯示出不同參數(shù)的取值范圍對敏感性參數(shù)的排名并無顯著差異，但是排名并不一樣，并且最敏感參數(shù)基本完全一致。可是在敏感性分析時，當(dāng)參數(shù)取值范圍在10%或50%時，可能會使部分敏感性參數(shù)變成不敏感參數(shù)，由此得到一個錯誤的敏感性分析結(jié)果。

本研究由于低估了灌溉對棉花生物量的影響，導(dǎo)致CWAM的模擬效果比較差。在敏感性和不確定性分析中，應(yīng)統(tǒng)計實測值的大概變化范圍，因為當(dāng)參數(shù)取值范圍過小時，就不能反映出實測值的變化程度，而當(dāng)參數(shù)取值范圍過大時，又會大大超出實測值的變化程度，所以確定好實測值的大概變化范圍，才能選擇出恰當(dāng)?shù)膮?shù)取值范圍。

本研究結(jié)論：（1）參數(shù)取值范圍變大，輸出結(jié)果的CV也會隨著增大，但是它們的平均值會表現(xiàn)出下降的趨向;（2）不同的參數(shù)取值范圍會對參數(shù)的敏感性分析結(jié)果產(chǎn)生一定程度的作用，即使它不會影響判別最敏感參數(shù)的個數(shù)，但是當(dāng)參數(shù)的取值范圍處于極端時（10%或50%），就會致使某些敏感性參數(shù)成為不敏感參數(shù);（3）當(dāng)參數(shù)取值范圍太小時，輸出變量就不容易呈現(xiàn)出實測值的變化范圍，當(dāng)參數(shù)取值范圍過大時，輸出變量就容易和實測值不符;（4）根據(jù)試驗得出參數(shù)默認(rèn)值±30%是CROPGRO-cotton模型敏感性分析和不確定性分析最佳參數(shù)取值范圍，本研究方法可以幫助其他模型在進行敏感性分析和不確定性分析時提供對照。

參考文獻：

[1]Li N，Yao N，Li Y，et al. A meta-analysis of the possible impact of climate change on global cotton yield based on crop simulation approaches[J]. Agricultural Systems，2021，193：103221.

[2]Jacot J，Kiniry J R，Williams A S，et al. Use of PhenoCam measurements and image analysis to inform the ALMANAC process-based simulation model[J]. Journal of Experimental Agriculture International，2021：120-144.

[3]劉 剛，謝 云，高曉飛，等. ALMANAC作物模型參數(shù)的敏感性分析[J]. 中國農(nóng)業(yè)氣象，2008，29（3）：259-263.

[4]李 穎，趙國強，陳懷亮，等. 基于冬小麥農(nóng)業(yè)氣候分區(qū)的WOFOST模型參數(shù)標(biāo)定[J]. 應(yīng)用氣象學(xué)報，2021，32（1）：38-51.

[5]蔡 福，米 娜，明惠青，等. WOFOST模型蒸散過程改進對玉米干旱模擬影響[J]. 應(yīng)用氣象學(xué)報，2021，32（1）：52-64.

[6]李 明. 基于WOFOST模型的春小麥產(chǎn)量預(yù)測方法研究[D]. 泰安：山東農(nóng)業(yè)大學(xué)，2020.

[7]楊艷敏，歐陽竹，王淑芬. 基于COTTON2K的華北平原和新疆2個棉區(qū)棉花耗水特征比較[J]. 華北農(nóng)學(xué)報，2012，27（增刊1）：229-233.

[8]Attia A，Rajan N，Nair S S，et al. Modeling cotton lint yield and water use efficiency responses to irrigation scheduling using Cotton2K[J]. Agronomy Journal，2016，108（4）：1614-1623.

[9]楊艷敏，劉小京，歐陽竹. COTTON2K在新疆棉花精準(zhǔn)種植管理中的應(yīng)用[J]. 中國生態(tài)農(nóng)業(yè)學(xué)報，2005，13（3）：123-126.

[10]徐瑞陽. 基于DSSAT-SUBSTOR模型的東北地區(qū)馬鈴薯單產(chǎn)模擬和種植成本效益分析[D]. 北京：中國農(nóng)業(yè)科學(xué)院，2021.

[11]徐春萌，田芷源，陳 威，等. 基于DSSAT作物模型的中美大豆主產(chǎn)區(qū)單產(chǎn)模擬與驗證[J]. 農(nóng)業(yè)工程學(xué)報，2021，37（3）：132-139.

[12]李翊塵，史培軍. 基于DSSAT模型的河套平原春小麥灌溉方案優(yōu)化[J]. 北京師范大學(xué)學(xué)報（自然科學(xué)版），2020，56（5）：728-739.

[13]Mack L，Boote K J，Munz S，et al. Adapting the CROPGRO model to simulate chia growth and yield[J]. Agronomy Journal，2020，112（5）：3859-3877.

[14]李 波，孫翔龍，姚名澤，等. 溫室秸稈不同還田量條件下DSSAT-CROPGRO-Tomato模型的調(diào)參與驗證[J]. 生態(tài)學(xué)雜志，2021，40（3）：908-918.

[15]趙子龍. 基于CROPGRO-Tomato模型東北寒區(qū)溫室番茄水分生產(chǎn)函數(shù)研究[D]. 沈陽：沈陽農(nóng)業(yè)大學(xué)，2018.

[16]王連喜，張 陽，李 琪，等. 作物模型參數(shù)敏感性分析現(xiàn)狀與展望[J]. 氣象科技，2018，46（2）：382-389.

[17]崔金濤，丁繼輝，Yesilekin N，等. 基于EFAST的CERES-wheat模型土壤參數(shù)敏感性分析[J]. 農(nóng)業(yè)機械學(xué)報，2020，51（12）：276-283.

[18]Foscarini F，Bellocchi G，Confalonieri R，et al. Sensitivity analysis in fuzzy systems：integration of SimLab and DANA[J]. Environmental Modelling & Software，2010，25（10）：1256-1260.

[19]宋明丹，馮 浩，李正鵬，等. 基于Morris和EFAST的CERES-wheat模型敏感性分析[J]. 農(nóng)業(yè)機械學(xué)報，2014，45（10）：124-131，166.

[20]李 波，李長信，姚名澤，等. 不同灌水處理下CROPGRO-Tomato模型全局敏感性和不確定性分析[J]. 沈陽農(nóng)業(yè)大學(xué)學(xué)報，2020，51（2）：153-161.

[21]Scosyrev E，Pethe A. Confidence intervals for exposure-adjusted rate differences in randomized trials[J]. Pharmaceutical Statistics，2021.

[22]Hemming K，Taljaard M. Why proper understanding of confidence intervals and statistical significance is important[J]. The Medical Journal of Australia，2021，214（3）：116-118.

[23]丁潤霞，黃韓亮. 正態(tài)分布概率區(qū)間模糊數(shù)及測度問題[J]. 佳木斯大學(xué)學(xué)報（自然科學(xué)版），2021，39（4）：147-151.

[24]梁 駿. 水文設(shè)計值置信區(qū)間計算方法研究[D]. 楊凌：西北農(nóng)林科技大學(xué)，2016.