米榮娟，聶志剛,2

(1.甘肅農(nóng)業(yè)大學(xué)信息科學(xué)技術(shù)學(xué)院,甘肅 蘭州 730070；2.甘肅農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院，甘肅 蘭州 730070)

建立和引入作物模型，為農(nóng)作物生產(chǎn)調(diào)控提供預(yù)測依據(jù)[1].模型的有效應(yīng)用依賴于對參數(shù)進(jìn)行準(zhǔn)確估計[2]，通過合適的方法定量評估影響模型的敏感參數(shù)是模型校準(zhǔn)和優(yōu)化的前提，同時是模型本地化應(yīng)用的基礎(chǔ).近年來，利用敏感性分析評估作物模型敏感參數(shù)的研究屢見報道，Confalonieri[3]和Makowski等[4]分別通過敏感性分析對不同作物模型進(jìn)行評估，以確定參數(shù)對模型預(yù)測方差的貢獻(xiàn)，擬合觀測值與模擬值，證明了敏感性分析對簡化模型的重要性.局部敏感性分析法通過改變單一輸入評價單個變量對模型輸出的影響，操作簡單[5]，多用于線性模型分析，但未充分考慮參數(shù)之間互作效應(yīng)，分析不全面.全局敏感性分析法克服了局部敏感性分析法的局限性，被廣泛運用于復(fù)雜機理模型分析中.其中，Sobol法、Morris參數(shù)選擇法及擴展傅里葉幅度檢驗法[6]等均為常見的全局敏感性分析法.Morris參數(shù)選擇法和Sobol法相較于局部敏感性分析法更精確和全面，效率和精度有所改善，但不能量化分析結(jié)果.實踐證明，EFAST全局敏感性分析法針對非線性復(fù)雜作物模型敏感參數(shù)篩選優(yōu)勢明顯，如姜志偉等[7]以河南洛陽為試驗區(qū)，運用EFAST法篩選出CERES-Wheat模型“本地化”關(guān)鍵作物參數(shù)及田間管理參數(shù)；邢會敏等[8]以北京地區(qū)試驗為例，運用EFAST法分不同時期評估Aqua-Crop 模型中的42個作物參數(shù)，以選擇敏感參數(shù)優(yōu)化模型；Li等[9]運用EFAST法評估DSSAT-CERES模型輸出響應(yīng)對39種作物基因型參數(shù)和6種土壤參數(shù)的敏感性，均取得良好效果.目前，國內(nèi)基于EFAST法針對APSIM模型參數(shù)的敏感性分析的研究尚不多見.此外，模型存在不確定性，敏感性分析結(jié)果會因不同地域氣候條件不同產(chǎn)生差異.定西是典型的半干旱黃土丘陵區(qū)，小麥?zhǔn)窃摰貐^(qū)重要的糧食作物，運用EFAST法在該地區(qū)開展小麥產(chǎn)量形成參數(shù)敏感性分析的研究鮮見報道.因此，本研究選用具有較好模擬效果且運用廣泛的APSIM模型為研究對象，以旱地春小麥試驗為例，依托試驗區(qū)實測數(shù)據(jù)資料，采用EFAST全局敏感性分析法對APSIM模型中影響春小麥產(chǎn)量形成的敏感參數(shù)進(jìn)行篩選，以簡化模型，并為APSIM模型優(yōu)化及本地化應(yīng)用提供依據(jù).

1 材料與方法

1.1 試驗區(qū)概況與數(shù)據(jù)來源

試驗區(qū)位于典型黃土高原丘陵區(qū)甘肅省定西市安定區(qū)鳳翔鎮(zhèn)安家溝村，該地區(qū)屬中溫帶旱作雨養(yǎng)農(nóng)業(yè)區(qū)，平均海拔約2 000 m，晝夜溫度變化大，光熱資源較為豐富，但水分嚴(yán)重不足，且基本來源于自然降水，無灌溉條件.表1列出了試驗區(qū)基本土壤特征參數(shù).

春小麥定西35為試驗作物，播種采用免耕播種機，設(shè)置播種量187.5 kg/hm2，一般播深7 cm，行距 0.25 m.播種試驗小區(qū)共計15個，面積6 m×4 m，保護行 0.5 m，完全隨機區(qū)組設(shè)置，施肥及田間管理等參考當(dāng)?shù)貙嶋H情況.當(dāng)?shù)剡m宜正常播期3月19日，7月下旬春小麥成熟后收獲.小麥產(chǎn)量以收獲時各小區(qū)打碾產(chǎn)量折算公頃產(chǎn)量.

1971～2019年定西市安定區(qū)氣象資料均來源于甘肅省氣象局，所需土壤參數(shù)和田間管理參數(shù)為試驗點實測數(shù)據(jù)，模型檢驗所需2015～2017年產(chǎn)量數(shù)據(jù)經(jīng)田間實測獲得，2018～2019年產(chǎn)量數(shù)據(jù)為查閱年鑒獲得，APSIM 模型模擬獲得其他有關(guān)產(chǎn)量數(shù)據(jù).

表1 土壤屬性特征

1.2 APSIM模型

APSIM是一個大型農(nóng)業(yè)生產(chǎn)系統(tǒng)模擬模型[10].該模型基于組件設(shè)計，通過插件機制容納各個模塊(作物，土壤，氣候，管理等)，將模型與用戶界面分離；用戶可以將模塊不同組合插入在一起，以針對不同的仿真配置APSIM，具有靈活性.模型以土壤模塊為中心，以日步長為模擬跨度，受土壤、水分及管理措施等驅(qū)動進(jìn)行農(nóng)作物生長的動態(tài)模擬[11-14].其它模型說明和運用實例的詳盡信息參考官網(wǎng)：https://www.apsim.info/.

1.3 全局敏感性分析方法

Saltelli等[15]基于方差分解原理提出EFAST法，綜合了Sobol法和FAST法，并在兩者基礎(chǔ)上進(jìn)行改進(jìn).它對于樣本要求低，計算高效，穩(wěn)定性高[16]，適用于包含少量參數(shù)的系統(tǒng)，具有理論優(yōu)勢和應(yīng)用潛力.主要采用2個指標(biāo)：一階敏感性系數(shù)和全局敏感性系數(shù)來評價敏感性高低的標(biāo)準(zhǔn).其中，一階敏感性系數(shù)評價單一輸入?yún)?shù)對輸出的影響，全局敏感性系數(shù)評價參數(shù)之間交互作用對輸出的影響.基本原理如下：

假設(shè)模型y=f(X)，X=(x1，x2，x3，…，xn)為滿足一定概率分布的多維參數(shù)空間，通過適當(dāng)?shù)霓D(zhuǎn)換函數(shù)Gi將模型轉(zhuǎn)換為y=f(s)，轉(zhuǎn)換函數(shù)Gi與參數(shù)xi的概率密度分布函數(shù)相關(guān):

Xi(s)=Gi(sinωis)，?i∈{1,2,3,…,n}

(1)

式中：s為標(biāo)量且s∈(-∞，+∞)；{ωi}為參數(shù)xi所定義的整數(shù)頻率；

對y=f(s)進(jìn)行傅里葉變化：

(2)

式中：通過對s在[-π，π]中等間距取值，對取值結(jié)果輸入模型，多次運行，得到Ak與Bk的近似值如下：

Λk=Ak2+Bk2

(5)

式中：k為傅里葉變換參數(shù)且k∈Z，A-k=Ak,B-k=Bk,Λ-k=Λk.參數(shù)xi變化引起模型方差為：

(6)

式中：Z0為非零整數(shù)，ωi是參數(shù)xi對應(yīng)的頻率.模型總方差為：

(7)

模型總方差可分解為：

(8)

式中，Vi是參數(shù)xi變化引起的模型方差，Vij是參數(shù)xi通過xj作用引起的模型方差，V1,2,…,k為參數(shù)xi通過x1、x2、…、xk引起的模型方差.歸一化處理后，參數(shù)xi的一階敏感性系數(shù)為：

(9)

全局敏感性系數(shù)為：

(10)

Dejonge[17]曾界定Si>0.05，STi>0.10(Si+STi>0.15)作為判斷敏感性的標(biāo)準(zhǔn).

1.4 模擬試驗設(shè)計和待分析參數(shù)選擇

表2是待分析參數(shù)說明及范圍，本研究主要分析小麥品種參數(shù).品種參數(shù)選擇來源于APSIM官網(wǎng)小麥模塊文檔說明(網(wǎng)址：https://www.apsim.info/)及Zhao等[13]與何亮等[18]文獻(xiàn)參考；品種參數(shù)上下限是根據(jù)模型初始值上下浮動±50%.基礎(chǔ)氣象數(shù)據(jù)來源于甘肅省氣象局，土壤數(shù)據(jù)、管理數(shù)據(jù)來源于田間實測，這些參數(shù)在模型輸入中均采用固定值.

本研究首先通過5 a的數(shù)據(jù)驗證APSIM模型的適用性；其次結(jié)合模型，借助Simlab軟件并采用其內(nèi)置EFAST方法，對2016～2017年試驗區(qū)影響小麥產(chǎn)量的部分參數(shù)進(jìn)行全局敏感性分析，取2 a敏感性系數(shù)平均值作為結(jié)果；然后對敏感性較高參數(shù)設(shè)置梯度模擬，確定參數(shù)變化范圍，檢驗單因素變化效應(yīng).

表2 選擇參數(shù)的上下限及分布和模型輸出

1.5 模型檢驗方法

基于已有研究，利用已校準(zhǔn)的APSIM模型對2015～2017年(田間實測)、2018～2019年(查閱年鑒)定西35號春小麥產(chǎn)量實測值及模型模擬值進(jìn)行擬合(圖1).歸一化均方根誤差(NRMSE)用以檢驗誤差，模型有效性指數(shù)(ME)用以檢驗?zāi)M精度.其中，NRMSE值越小，誤差越?。籑E>0.5，則模擬效果良好[19].相關(guān)公式為：

(11)

(12)

式中：Yobs，Ysim和Ymean分別為實測值、模擬值和實測平均值.

1.6 數(shù)據(jù)處理

采用 Microsoft Excel 2010 軟件整理匯總產(chǎn)量數(shù)據(jù)，利用 Origin 2017 軟件擬合實測值和模擬值，并繪制旱地春小麥產(chǎn)量形成單因素效應(yīng)曲線和不同比例下參數(shù)引起產(chǎn)量波動圖.

2 結(jié)果與分析

2.1 模型驗證

經(jīng)檢驗(圖1)，實測值與模擬值擬合較好，數(shù)據(jù)點均分布在限定誤差范圍線內(nèi).通過線性擬合發(fā)現(xiàn)，線性回歸決定系數(shù)R2=0.998 5，歸一化均方根誤差NRMSE=5.314%，模型有效性指數(shù)ME=0.976，模型模擬產(chǎn)量較為準(zhǔn)確.

圖1 小麥產(chǎn)量模擬值與實測值線性擬合Figure 1 Linear fit between simulated and observed values of wheat yield

2.2 敏感性分析結(jié)果

APSIM小麥產(chǎn)量形成參數(shù)敏感性分析借助SimLab[20]專業(yè)軟件完成，該軟件由預(yù)處理模塊、模型執(zhí)行模塊和后處理模塊構(gòu)成.具體操作步驟：輸入?yún)?shù)因子，選擇參數(shù)范圍及分布概率函數(shù)；選擇抽樣方法，生成.sam格式的樣本數(shù)據(jù)；將樣本數(shù)據(jù)輸入模型，獲得模型輸出結(jié)果；對參數(shù)進(jìn)行敏感性分析.

根據(jù)圖2可知，一階敏感性系數(shù)和全局敏感性系數(shù)的分析結(jié)果總體走勢基本呈現(xiàn)一致性.在參數(shù)選擇范圍內(nèi)，按照前文所述判斷依據(jù)，得到敏感參數(shù)及排序：每克莖籽粒數(shù)量>灌漿期籽粒日潛在灌漿速率>單株最大籽粒質(zhì)量>作物春化敏感性指數(shù)>灌漿到成熟積溫>出苗到拔節(jié)積溫.其中，每克莖籽粒數(shù)量、灌漿期籽粒日潛在灌漿速率和單株最大籽粒質(zhì)量與產(chǎn)量形成密切相關(guān)，對模擬產(chǎn)量貢獻(xiàn)率達(dá)50.6%、31.9%和27.5%.敏感參數(shù)主要包括直接關(guān)系產(chǎn)量形成的參數(shù)及對生育期構(gòu)成影響的物候參數(shù)(不同生育階段的積溫等)等，均與實際生產(chǎn)情況相符合.將對產(chǎn)量形成較敏感的6個參數(shù)，作為模型進(jìn)一步調(diào)參優(yōu)化基礎(chǔ).

圖2 APSIM模型參數(shù)敏感性系數(shù)Figure 2 Parameter sensitivity coefficient of APSIM model

2.3 敏感系數(shù)對小麥產(chǎn)量的影響

對2016～2017年影響旱地春小麥產(chǎn)量形成的敏感參數(shù)設(shè)置梯度檢驗，觀察敏感參數(shù)影響產(chǎn)量變化的趨勢.梯度設(shè)置以模型默認(rèn)值(base)為基準(zhǔn)，增加10%、20%、30%、40%、50%及減少10%、20%、30%、40%、50%.

2.3.1 單因素總體效應(yīng) 圖3表明，每克莖籽粒數(shù)量在默認(rèn)值基礎(chǔ)上減少比例時，模擬產(chǎn)量出現(xiàn)持續(xù)下降趨勢；而在默認(rèn)值基礎(chǔ)上增加比例，模擬產(chǎn)量趨于平穩(wěn).灌漿期籽粒日潛在灌漿速率在默認(rèn)值基礎(chǔ)上減少10%時模擬產(chǎn)量達(dá)到峰值，繼續(xù)增長后模擬產(chǎn)量基本趨于平穩(wěn)；而減少比例時，模擬產(chǎn)量整體呈下降趨勢.單株最大籽粒質(zhì)量與每克莖籽粒數(shù)量的變化趨勢基本呈現(xiàn)一致性.出苗到拔節(jié)積溫在模擬上出現(xiàn)了年份差異，兩者在默認(rèn)值基礎(chǔ)上增加比例時，產(chǎn)量基本處于遞減模式；而在減少比例時，2016年模擬產(chǎn)量遞減，2017年遞增，可能是由于不同播期導(dǎo)致外部自然環(huán)境不同產(chǎn)生的差異.灌漿到成熟積溫在默認(rèn)值基礎(chǔ)上減少比例時，產(chǎn)量基本是下降狀態(tài)，僅僅在2017年出現(xiàn)部分偏差；但當(dāng)其在默認(rèn)值基礎(chǔ)上增加時，卻出現(xiàn)了不同年份產(chǎn)量增加和減少的差異性；作物春化敏感性指數(shù)隨數(shù)值增加時模擬產(chǎn)量基本呈下降趨勢.總體來看，當(dāng)變化比例處于-10%～+10%時，產(chǎn)量變化幅度極小，趨于穩(wěn)定；當(dāng)變化比例處于-50%～-10%時，產(chǎn)量變化幅度較大.

圖3 旱地小麥產(chǎn)量形成因素單因素效應(yīng)Figure 3 Single factor effect of wheat yield formation

2.3.2 等比例產(chǎn)量變化 圖4表示敏感參數(shù)在模型默認(rèn)值基礎(chǔ)上增加和減少比例時2016和2017年模擬產(chǎn)量平均值相較于默認(rèn)值狀態(tài)下模擬產(chǎn)量平均值的具體變化量.總體來看，產(chǎn)量變化絕對值介于0～500 kg/hm2之間，且只有降低比例時，春化敏感性指數(shù)引起模擬產(chǎn)量變化是在默認(rèn)值基礎(chǔ)上增加，其余均是在默認(rèn)值基礎(chǔ)上減少.分開來看，在減少同等比例時，單株最大籽粒質(zhì)量引起產(chǎn)量變化幅度均最高，在減少50%比例時，產(chǎn)量差值絕對值達(dá)到最高493.55 kg/hm2；而在增加同等比例時，出苗到拔節(jié)積溫引起產(chǎn)量變化幅度最高，在增加50%比例時，產(chǎn)量差值絕對值達(dá)到最高249.85 kg/hm2.

3 討論

在影響小麥產(chǎn)量形成參數(shù)中，每克莖籽粒數(shù)量、灌漿期籽粒日潛在灌漿速率和單株最大籽粒質(zhì)量是敏感性系數(shù)最高的3個參數(shù)，Zhao等[13]研究中有相似結(jié)論.其中，每克莖籽粒數(shù)量和單株最大籽粒質(zhì)量2個參數(shù)是構(gòu)成小麥產(chǎn)量要素的基本參數(shù)，直接描述了小麥品種的產(chǎn)量特性.籽粒數(shù)量和籽粒大小的乘積用以計算產(chǎn)量，籽粒數(shù)是在開花期根據(jù)每克莖籽粒數(shù)量及抽穗期和籽粒灌漿開始之間干物質(zhì)積累的函數(shù)計算，小麥籽粒大小是通過籽粒的生長速率(籽粒潛在灌漿速率)和開花后生殖階段生物量向籽粒的分配來計算的[21].灌漿期籽粒日潛在灌漿速率是調(diào)節(jié)小麥生育期籽粒灌漿速率的參數(shù).小麥灌漿期灌漿速率一般是“慢-快-慢”，春小麥灌漿時間長短和灌漿速率快慢，直接影響粒重高低.加快灌漿速率能有效增加千粒質(zhì)量，千粒質(zhì)量作為產(chǎn)量構(gòu)成重要因素，增加可使產(chǎn)量提高.積溫影響作物生長發(fā)育，積溫積累快，能夠加速籽粒形成.研究表明，春小麥生長中后期，積溫影響生育期長短；生育期長短決定干物質(zhì)積累時間，干物質(zhì)積累與籽粒形成有關(guān)，籽粒形成影響最終產(chǎn)量[18].根據(jù)前人研究發(fā)現(xiàn)，從生長發(fā)育和小麥生產(chǎn)2個方面看，出苗到拔節(jié)積溫都是一項可利用的生態(tài)指標(biāo)[22]，期間積溫增加能相應(yīng)地縮短出苗天數(shù).此外，溫度適當(dāng)升高有利于壯苗形成，利于后期生長發(fā)育，促進(jìn)產(chǎn)量形成[23].灌漿期籽粒日潛在灌漿速率是籽粒灌漿過程的重要因素.近期研究指出，與出苗到拔節(jié)期間溫度變化相比，灌漿到成熟階段溫度變化對小麥產(chǎn)量形成的影響更顯著[24].在適宜溫度范圍內(nèi)，積溫與灌漿速率存在正相關(guān)性，積溫增加會促使灌漿速率上升[25].但在過高溫影響下，易導(dǎo)致灌漿過程中止，成熟期提前，干物質(zhì)積累時間縮短，籽粒干癟，影響千粒質(zhì)量，從而導(dǎo)致減產(chǎn).灌漿期是關(guān)乎小麥籽粒形成關(guān)鍵期，易受自然災(zāi)害影響.因此，做好管理，適時適當(dāng)澆水施肥，能夠有效延長灌漿時間，更好促進(jìn)小麥產(chǎn)量形成[26].春化階段是小麥必經(jīng)發(fā)育期，春化對于幼苗階段影響尤其顯著，此階段的生長發(fā)育直接決定了后續(xù)發(fā)育進(jìn)程，影響最終產(chǎn)量形成.在APSIM中，春化影響從出苗到開花期.春化敏感性指數(shù)，在小麥文件中為品種特定參數(shù)，影響春化因子，造成累計春化的差異性，具有較強敏感性，需要進(jìn)一步校準(zhǔn).結(jié)果表明，小麥產(chǎn)量形成與其自身品種遺傳參數(shù)相關(guān)外，也受到物候期影響.

圖4 不同比例下參數(shù)引起產(chǎn)量波動Figure 4 Parameters at different ratios cause output fluctuations

敏感性分析得到了 APSIM 模型中對小麥產(chǎn)量模擬結(jié)果影響顯著的品種參數(shù).為了進(jìn)一步實現(xiàn)參數(shù)“本地化”，使參數(shù)適用于特定研究區(qū)域[27]—定西市安定區(qū)，須對選定參數(shù)以優(yōu)化調(diào)整，使模型具有高可適用性.對敏感性系數(shù)未達(dá)到判斷指標(biāo)的模型參數(shù)，采取經(jīng)驗值或?qū)崪y值及查閱文獻(xiàn)等方法獲得.

前人研究證實[28-29]EFAST方法為評估整個可行參數(shù)空間中參數(shù)敏感性提供了一種簡單、快速且全局的方法.但在實施EFAST方法時，模型參數(shù)變化范圍及樣本數(shù)量選擇對敏感性分析結(jié)果會產(chǎn)生影響.一般而言，參數(shù)范圍選擇基于不同作物種類及模型文檔說明，參考實際應(yīng)用情況界定.本研究是為了客觀選擇敏感性較強的參數(shù)，暫時忽略參數(shù)范圍變化影響，這在非特殊情況下基本可行.本研究中對管理參數(shù)(如播期、行間距等)、土壤參數(shù)(田間持水率、凋萎系數(shù)等)和日最高氣溫、最低氣溫、日照時長等氣象參數(shù)采用固定值一次性輸入，未作具體討論.此外，參數(shù)取值范圍差異可能會造成分析結(jié)果的不確定性.因此，繼續(xù)補充其他參數(shù)的敏感性分析，合理地確定參數(shù)分布范圍以及對敏感參數(shù)進(jìn)行優(yōu)化則是下一步需要研究的方向.

4 結(jié)論

本研究利用擴展傅里葉幅度檢驗法對定西春小麥APSIM產(chǎn)量形成參數(shù)進(jìn)行敏感性分析，得到6個敏感參數(shù)：每克莖籽粒數(shù)量、灌漿期籽粒日潛在灌漿速率、單株最大籽粒質(zhì)量、出苗到拔節(jié)積溫、灌漿到成熟積溫和作物春化敏感性指數(shù)，可優(yōu)先校準(zhǔn).通過對敏感參數(shù)設(shè)置單因素梯度檢驗表明，當(dāng)各個參數(shù)變化比例處于-10%～+10%時，產(chǎn)量變化幅度極小，趨于穩(wěn)定；當(dāng)變化比例處于-50%～-10%時，產(chǎn)量變化幅度較大.同時，產(chǎn)量變化絕對值介于0～500 kg/hm2之間，在同等增加和減少比例時分別是出苗到拔節(jié)積溫和單株最大籽粒質(zhì)量引起產(chǎn)量變化幅度最大.