基于EFAST的不同生產水平下WOFOST模型參數敏感性分析

興 安 卓志清 趙云澤 李 勇 黃元仿,2

(1.中國農業大學土地科學與技術學院, 北京 100193; 2.自然資源部農用地質量與監控重點實驗室, 北京 100035)

0 引言

目前，作物模型在生產潛力評價、氣候變化影響評估、作物生長監測、作物產量預報、農業風險評估、農業管理決策等諸多領域的應用越來越廣泛[1-6]，其中模型參數校正和優化仍然是模型本地化應用的關鍵過程，而敏感性分析是進行參數校正和優化的重要前提[7-9]。敏感性分析能夠將不確定性結果分配至不同的模型參數，量化模型參數對輸出結果的影響，識別并篩選出模型中的關鍵參數，從而降低參數優化過程中的工作量和輸出結果的不確定性[10]。敏感性分析通常分為局部敏感性分析方法和全局敏感性分析方法[11]，其中，全局敏感性分析方法有EFAST方法[12-13]、Morris方法[14]、Sobol’方法[15]等，均在作物模型敏感性分析中得到了廣泛的應用[16-20]。

作物模型中參數的重要性不僅與模型結構有關，還受到其他參數和輸入數據的影響。特別是在區域尺度上，氣候、土壤以及管理條件的復雜性導致區域尺度上參數敏感性分析結果存在一定的不確定性[21-22]。因此，許多研究針對不同區域年份的氣候、土壤、管理條件的差異和作物參數取值范圍進行了參數敏感性的影響分析[22-24]。如文獻[25]選取不同氣候類型下的4個冬小麥代表性站點，對APSIM-Wheat模型品種參數、土壤參數在不同氣候類型及不同產量水平下的敏感性差異進行分析，結果表明，不同氣候區間產量和生育期蒸散對參數的敏感性分析結果存在差異，不同產量水平間差異不大；文獻[26]采用Morris和Sobol’方法，分別計算了歐洲5個國家的主要水稻種植區及代表不同大陸性氣候條件的3個年份下的WARM模型參數敏感性，結果表明，模型對作物參數的敏感性在不同氣候條件下并不總是一致的；文獻[27]采用Morris和EFAST方法對歐洲不同氣候條件(包括未來氣候變化)下WOFOST模型不同作物參數的敏感性進行了分析，結果表明，敏感性分析受氣候條件，尤其是極端氣候條件的影響，從而表現出不同的結果，但整體上最敏感參數趨于一致。以上研究均強調了在區域尺度和長時間序列上進行參數敏感性分析的必要性，從結果上看，在較大時空尺度上敏感性分析結果因氣候類型的多樣性及復雜性而存在較大差異，而在較小尺度相同氣候類型下，不同氣候條件對參數敏感性的影響還有待進一步探討。此外，WOFOST模型作為應用最廣泛的機理作物模型之一，其參數敏感性分析也是目前的研究熱點。然而，大部分研究在進行參數敏感性分析[24,27-29]和參數校驗[3-4,30-32]時未考慮蒸散速率修正因子(CFET)和土壤水分消耗作物群數量(DEPNR)等作物水分利用相關品種參數，而這些參數能夠表征作物對水分脅迫的響應，這會影響模型模擬精度，尤其是水分限制生產水平下的作物生長模擬。同時WOFOST模型對所有作物生長的描述相同，在同一種植制度、不同作物間的參數敏感性差異分析有助于進一步理解模型運作過程。

本文選取溫帶季風氣候類型的黃淮海平原旱作區不同積溫區內的3個站點，基于農業氣象試驗站點數據以及實地采樣數據，采用EFAST方法分析WOFOST模型冬小麥和夏玉米產量在不同生產水平、不同氣候條件下的參數敏感性及其排序一致性，并探討敏感性差異形成的原因，為后續區域尺度模型參數標定和驗證提供基礎。

1 材料與方法

1.1 研究區概況

黃淮海平原旱作區[33-34]是我國重要旱作作物種植區及商品糧基地，其大部分屬溫帶季風氣候，年均溫度14～16℃，年均降雨量500～1 200 mm，土壤類型主要為潮土、褐土和砂漿黑土，主要種植糧食作物為冬小麥和夏玉米，冬小麥播種時間一般為10月，次年5月下旬或6月上、中旬成熟，夏玉米常年播種在6月，9月上、中旬成熟。

本研究基于黃淮海旱作區2000—2015年日平均溫度大于等于10℃的年均積溫數據，以200℃為間隔劃分不同積溫區[35]，在積溫為4 600～4 800℃、4 800～5 000℃、5 200℃以上的積溫區中選取黃驊、商丘和駐馬店3個農業氣象試驗站點作為研究地區(圖1)，其基本信息和氣候特征見表1。

圖1 黃淮海平原旱作區分布圖

表1 3個代表站點的氣候基本特征

1.2 WOFOST模型

WOFOST模型由荷蘭瓦赫寧根大學和世界糧食研究中心(CWFS)共同開發研制，以逐日氣象數據為驅動，通過土壤、管理和作物參數限制與調整作物生長過程的動態解釋性模型，可以模擬潛在、水分限制和養分限制3種水平下的作物生長。WOFOST模型中物候期用無量綱的狀態變量——生長階段(Development stage, DVS)來描述，DVS為0表示出苗期、1表示開花期、2表示成熟期[36]。

模型運行所需的氣象數據(逐日最高氣溫、最低氣溫、降雨量)、管理數據(播種日期、出苗期、成熟期)從中國氣象數據網(http:∥data.cma.cn/site/index.html)獲取，逐日輻射、早晨水汽壓、地上2 m處平均風速采用FAO提出的Angstrom公式將從中國氣象數據網獲取的日照時數、平均相對濕度以及平均風速進行轉換獲取。土壤數據(土壤含水率、田間持水量、萎蔫系數、水分特征曲線、飽和導水率)由旱作區取樣數據和農業氣象站數據獲取，土樣取于2017年5—6月，采用重鉻酸鉀外加熱法、環刀法測定土壤有機質含量和容重，采用激光粒度儀測定土壤顆粒組成，并采用Soil Water Characteristics軟件[37]計算轉換為所需土壤參數。

1.3 研究方法

1.3.1參數選取及取值范圍

根據冬小麥和夏玉米不同生長發育特點以及前人關于WOFOST模型作物參數敏感性分析結果[24,27,29]，選取初始、綠葉、同化、同化物轉化、呼吸作用、干物質分配、死亡、根系以及水分利用相關參數進行敏感性分析，2種作物參數各50個，其中生育期積溫參數未考慮，運行模型時各生育期參數根據農業氣象站點數據中的實際生育期計算獲取。模型作物參數缺省值基礎上浮動±10%并考慮模型參數默認范圍作為參數的取值范圍，敏感性分析中將所有參數的分布假定為均勻分布，詳見表 2。表2中TDWI為初始總干物質質量；LAIEM為出苗時葉面積指數；RGRLAI為葉面積指數最大增長率；SLATB為不同生育期的比葉面積；SPAN為35℃時葉片生命周期；TBASE為葉齡的低溫閾值；KDIFTB為不同生育期散射光消光系數；EFFTB為日平均氣溫為0℃或40℃時單葉同化CO2的光能利用效率；AMAXTB為單葉最大CO2同化速率；TMPFTB為最大同化速率在不同均溫條件下的校正因子；TMNFTB為總同化速率低溫校正因子；CVL為葉同化物轉化效率；CVO為儲存器官同化物轉化效率；CVR為根同化物轉化效率；CVS為莖同化物轉化效率；Q10為溫度變化10℃時呼吸作用變化率；RML為葉維持呼吸速率；RMO為儲存器官的維持呼吸速率；RMR為根的維持呼吸速率；RMS為莖的維持呼吸速率；RFSETB為不同生育期衰老校正因子；FRTB為不同生育期莖干物質的分配系數；FLTB為不同生育期葉干物質的分配系數；FOTB為不同生育期儲存器官干物質的分配系數；RDRRTB為不同生育期根相對死亡速率；RDRSTB為不同生育期莖相對死亡速率；CFET為蒸散速率修正因子；DEPNR為土壤水分消耗作物群數量；RDI為初始根深；RRI為根深最大日增加量；RDMCR為最大根深。

表2 WOFOST模型冬小麥-夏玉米參數取值范圍

1.3.2敏感性分析方法

EFAST法是SALTELLI等[13]在傅里葉幅度敏感性檢驗法(Fourier amplitude sensitivity test, FAST)的基礎上，結合Sobol’法的優點所提出的一種基于方差分解的全局敏感性分析方法，具有穩定、計算快速等特點。EFAST法原理是選取合適的搜索曲線在多維參數空間中運行，將一組非線性相關的整數頻率分配給模型的所有輸入參數，并對模型中選取的參數引入一個具有共同獨立參數的函數，使模型作為獨立參數的周期函數，將多維積分降為一維積分。通過把目標函數轉換成傅里葉級數，得到各頻率的傅里葉頻譜曲線，由頻譜曲線可計算由參數xi引起的模型輸出方差和總方差，其比值即為該參數的敏感度，具體算法詳見文獻[13,38-39]，其簡單算法公式為

(1)

(2)

(3)

STi=Si+Sij+…+Sij…k

(4)

式中V——模型總方差

Vi——xi輸入變化單獨引起的模型方差

Vij——xi通過xj作用貢獻的耦合方差

Vij…k——xi通過與xj、…、xk相互作用貢獻的耦合方差

Si——xi的一階敏感性指數

Sij——xi的二階敏感性指數

Sij…k——xi的多階敏感性指數

STi——xi的總敏感性指數

1.3.3敏感性分析方案

基于農業氣象試驗站點多年降雨量數據，選取2000年(豐水)、2007年(平水)、2013年(枯水)為研究年份，采用EFAST方法，分別對黃驊、商丘、駐馬店等站點冬小麥和夏玉米不同生產水平下的產量進行參數敏感性分析，旨在分析不同氣候條件下和不同生產水平下參數敏感性的差異。模型產量的輸出結果為總儲存器官干物質量(Total dry weight of living storage organs，TWSO)，不同氣候條件由不同站點(不同積溫區)和不同年份(不同降水年型)下的氣候條件來表征，即本研究中共考慮9個不同氣候條件下的參數敏感性，不同生產水平考慮潛在和水分限制生產水平。全局敏感性分析試驗借助軟件Simlab2.2的EFAST模塊實現。EFAST方法認為采樣個數大于等于65倍的參數個數為有效，本研究采樣數取73，總采樣數為3 650，能夠滿足EFAST方法基本要求。模型3年3站點2種作物共需運行3 650×3×3×2=65 700次，本研究全部的模型運算通過編寫Python程序實現。

通過對不同氣候條件下各參數敏感性指數進行平均，計算出各參數總敏感指數，采用參數敏感性指數大于或等于0.10作為界定主要敏感參數的標準[40]。

1.3.4一致性檢驗方法

不同生產水平、氣候條件的敏感性排序的一致性由TDCC(Top-down concordance coefficient)系數來衡量[41]。

假設有排序矩陣Rij，其中有m個因子，n個觀測值。rnm為變量m在所有n個觀測值的排序序號。在本研究中，m為不同生產水平或不同氣候條件，n為待分析的參數個數。

(5)

計算排序矩陣Rij的savage得分

(6)

并代入TDCC系數計算公式

(7)

式中CTDCC——TDCC系數

CTDCC的顯著性P值由統計值T計算，計算式為

T=m(n-1)CTDCC

(8)

T為自由度n-1的χ2分布。一般認為，CTDCC越接近1，且P值小于0.05，表明各變量排序具有顯著的一致性。

2 結果分析

2.1 冬小麥-夏玉米參數敏感性

2.1.1冬小麥參數敏感性

WOFOST模型冬小麥產量在2種不同生產水平下的參數敏感性分析結果如圖2所示。潛在生產水平下產量主要敏感參數有TBASE、CVO、TMNFTB3、AMAXTB1.50、FLTB0.25、FOTB1.00等，其中TBASE和CVO的全局敏感性指數最高，為0.21和0.17。水分限制生產水平下對產量敏感的冬小麥參數有CFET、CVO、TBASE、KDIFTB0.00、TMNFTB3、SLATB0.00、DEPNR等，其中CFET為最敏感的參數，其全局敏感性指數為0.56。由此可發現，2種生產水平下的主要敏感參數有所區別，潛在生產水平下，冬小麥的生長發育主要受溫度和太陽輻射影響，因此溫度相關的參數如葉齡低溫閾值與總同化速率低溫校正因子有較高的敏感性；而水分限制生產水平下，冬小麥所需水分僅由降水提供，當降水不足導致水分虧缺會影響冬小麥的正常生長發育，因此產量對蒸散速率修正因子等作物水分利用相關參數高度敏感。

圖2 冬小麥參數敏感性分析結果

2.1.2夏玉米參數敏感性

圖3 夏玉米參數敏感性分析結果

夏玉米產量在不同生產水平下的參數敏感性分析結果如圖3所示。潛在生產水平下，夏玉米產量主要敏感參數有TMNFTB10、EFFTB40、FLTB0.33、SPAN、CVL、SLATB0.00、AMAXTB1.75以及FLTB0.00等，其中最敏感的參數為TMNFTB10和EFFTB40，其全局敏感性指數分別為0.88和0.54。水分限制生產水平下，夏玉米作物參數敏感性分析結果與潛在生產水平差異不大，除了上述潛在生產水平主要敏感參數以外，RDI和SLATB0.78也表現較高的敏感性，其中TMNFTB10依然是最敏感的參數，其全局敏感性指數為0.69。該結果表明，在潛在和水分限制生產水平下，受溫度影響的同化速率相關參數是影響夏玉米產量的主要敏感參數，而水分則沒有對夏玉米生長發育起到太大的限制作用。

2.2 冬小麥-夏玉米參數敏感性排序一致性

2.2.1不同生產水平

采用TDCC系數量化冬小麥和夏玉米產量在不同生產水平作物參數敏感性的排序一致性，結果表明冬小麥和夏玉米不同生產水平TDCC系數分別為0.82和0.98，P均小于0.01。可見冬小麥和夏玉米在不同生產水平下參數敏感性排序均有較高的顯著一致性，但相較于夏玉米，冬小麥不同生產水平TDCC系數較低，其不同生產水平間敏感性差異較夏玉米大，這與上述冬小麥和夏玉米參數敏感性分析結果相一致。

2.2.2不同氣候條件

不同氣候條件冬小麥作物參數敏感排序一致性檢驗結果如表3所示，潛在生產水平下由不同站點和年份組成的不同氣候條件的總排序一致性較好，TDCC系數為0.92，不同站點和不同年份上也均表現出較高的一致性，TDCC系數最低達0.91，最高達0.97，P均小于0.01，說明不同氣候條件下潛在生產水平參數敏感性排序具有顯著的一致性，不同氣候條件對參數敏感性的影響較低。而水分限制生產水平下不同氣候條件排序一致性較差，TDCC系數僅為0.61，同時在不同站點和不同年份參數敏感排序一致性檢驗結果中可發現，2000、2007、2013年各年份不同站點間TDCC系數分別為0.84、0.94和0.77，黃驊、商丘、駐馬店各站點不同年份間TDCC系數為0.60、0.70和0.78，表明各站點在不同年份間的一致性較差，而同一年份在不同站點間的一致性較好，可見不同氣候條件對水分限制生產水平下的參數敏感性有較大的影響，且不同降水年型氣候條件差異造成的敏感性差異較大。

表3 不同氣候條件下冬小麥參數敏感性的排序一致性

不同氣候條件夏玉米參數敏感性的排序一致性檢驗結果如表 4所示，潛在生產水平下夏玉米參數敏感性TDCC系數為0.98，在不同年份和不同站點間TDCC系數均在0.97以上，P均小于0.01，表明潛在生產水平下夏玉米參數敏感性排序在不同氣候條件下具有顯著的一致性。而在水分限制生產水平下，其TDCC系數為0.86，與潛在生產水平相比，不同氣候條件下夏玉米水分限制生產水平參數敏感性排序一致性相對較差。不同年份和不同站點的一致性檢驗表明商丘和駐馬店站點在3個不同年份和2013年3個不同站點TDCC系數較低，參數敏感性排序一致性較差，可見不同氣候條件對水分限制生產水平下的參數敏感性有影響，其中不同降水年型間氣候條件差異是造成總體一致性較差的主要原因。

表4 不同氣候條件下夏玉米參數敏感性的排序一致性

2.3 不同氣候條件下水分限制生產水平參數敏感性差異

對于冬小麥產量而言，2000、2013年所有站點產量對CFET、DEPNR等水分利用相關參數高度敏感，而2007年所有站點產量對CFET、DEPER等參數不敏感，相反TBASE等受溫度影響的同化速率相關參數是其主要敏感參數。夏玉米相較于冬小麥，在不同氣候條件下的參數敏感性差異不大，在2013年商丘和駐馬店站點夏玉米主要敏感參數為SPAN、CFET、DEPNR、RDI等葉片和水分利用相關參數，而其余站點、年份主要敏感的參數均為TMNFTB10、EFFTB40等(圖4，圖中站點年份編號中，HH表示黃驊，SQ表示商丘，ZMD表示駐馬店，后面數字表示年份)?？傮w而言，水分限制生產水平下不同氣候條件對參數敏感性的影響有以下2種表現形式：降水不足、水分脅迫程度較大時(枯水年)，作物產量對CFET、DEPER、RDI等水分利用和根系相關參數敏感；降水較充足、水分脅迫程度較低時(豐水年和平水年)，作物產量表現出對TBASE、TMNFTB10及EFFTB40等受溫度影響的同化速率和光能利用效率相關參數敏感，這與潛在和水分限制生產水平下的參數敏感性分析結果相似。此外，2000年雖為豐水年，但因1999年(枯水年)降雨量僅為517.3 mm，冬小麥生育期前期水分脅迫程度較大，因此也會對CFET和DEPER等參數有較大的敏感性。

3 討論

3.1 參數全局敏感性

上述研究表明冬小麥潛在和水分限制生產水平下產量主要敏感參數有所區別。潛在生產水平下主要對TBASE、TMNFTB3、CVO等參數敏感，其中前2個參數表征冬小麥低溫條件下進入越冬期同化速率減慢或停止的過程，CVO表示的是同化物轉化為儲存器官干物質的效率，這些過程均對作物生長有著重要的影響，文獻[27,30,42]等研究中也表明了這些參數的敏感性和重要性。冬小麥在水分限制生產水平下主要敏感參數為CFET、CVO、TBASE、TMNFTB3等，其中敏感性指數最高的CFET參數表征作物面對水分脅迫時的響應[43]。目前大部分研究在進行敏感性分析時對參數CFET的考慮較少，但該參數在本研究中表現出對水分限制生產水平下產量有較大的影響，該參數準確校驗可能對黃淮海旱作區作物產量模擬有重要意義。

夏玉米潛在和水分限制生產水平下的產量主要敏感參數相一致，主要為TMNFTB10、EFFTB40、SPAN，當水分充足或水分脅迫較小時，作物產量形成主要由光溫影響，因此在低溫或高溫下產量對光能利用效率和最大同化速率有較高的敏感性。但文獻[24,27-28]關于夏玉米參數敏感性的研究表明，除EFFTB40、SPAN 等參數外，夏玉米產量還對CVO、TBASE、FOTB1.10、FLTB0.88、SLATB0.00和SLATB0.5等參數有著較高的敏感性，其中部分參數在本研究中并沒有表現出很強的敏感性，可能的原因是參數的取值范圍、區域尺度、模擬的氣象條件、環境或田間管理措施等因素不同[44]，這也說明了在特定操作環境下運用模型之前進行敏感性分析的重要性[27]。

3.2 參數敏感性一致性檢驗

冬小麥和夏玉米不同生產水平參數敏感性分析結果表明冬小麥不同生產水平主要敏感參數有所差別，而夏玉米則差異不大，冬小麥不同生產水平一致性相對于夏玉米較差。這主要是由2種作物生育期內降雨量的不同所導致，冬小麥生長期主要為10月至次年6月，在黃淮海旱作地區正是降水較少時期，水分缺乏情況下會提高作物水分利用相關參數的敏感性，而夏玉米生育期為6月至10月，這期間雨水較充足，基本能夠滿足作物生長發育需求，夏玉米生長發育主要受光溫條件限制，因此其不同生產水平的參數敏感排序一致性較高，主要敏感參數間差異不大。

冬小麥和夏玉米參數敏感在不同氣候條件下的排序一致性檢驗結果表明：在不同氣候條件下WOFOST模型潛在生產水平下的產量對作物參數的敏感性有較高的一致性和穩定性，而在水分限制生產水平下一致性較差。這可能是由于在黃淮海平原旱作區處于同一氣候類型區，光溫條件差異并不大，模型對外界光溫條件的響應較為一致，因此會導致不同氣候條件下潛在生產水平敏感性參數有較高的一致性。相反，研究區內無論年際間或是不同地區降水差異均較大，從而增加了不同氣候條件下水分限制生產水平參數敏感性差異。本研究的TDCC系數相對低于文獻[17]的研究結果，其主要原因可能在于本研究計算一致性時考慮了全部的參數，而文獻[17]選取前9個敏感參數進行一致性檢驗；而不同氣候條件參數敏感性TDCC系數高于文獻[25]的研究結果，這主要是由于文獻[25]的研究主要是針對大區域尺度，選取的4個站點分屬不同的氣候類型，氣候條件相差較大，因此各氣候區間敏感參數TDCC系數較低。同時，本研究中不同生產水平間參數敏感性TDCC系數高于前人的研究結果，其原因可能是由于本研究只考慮了潛在和水分限制生產水平，如果考慮養分限制甚至是實際產量，其影響因素更加復雜，相應的敏感性分析結果可能差異更大。由此可進一步明確作物模型是通過作物參數來表達作物對不同外界條件的不同響應，不同外界條件會導致不同的參數敏感性，考慮的外界條件越復雜，參數敏感性的差異會越大。因此，在特定環境下運用模型時進行敏感性分析是不可缺少的過程，同時，多種生產水平及更為復雜的影響因素下的模型敏感性分析可能是進一步了解模型并提高區域上應用的有效措施之一。

4 結論

(1)WOFOST模型冬小麥2種生產水平的TDCC系數為0.82，主要敏感參數間有所差別，潛在生產水平下的主要敏感參數有葉齡的低溫閾值(TBASE)、儲存器官同化物轉化效率(CVO)、總同化速率在低溫3℃時的校正因子(TMNFTB3)等，而冬小麥水分限制生產水平下蒸散速率修正因子(CFET)和土壤水分消耗作物群數量(DEPNR)等水分利用相關參數以及儲存器官同化物轉化效率(CVO)等參數表現出較大的敏感性。

(2)WOFOST模型夏玉米2種生產水平TDCC系數為0.98，主要敏感參數間差異不大，主要為總同化速率低溫10℃時的校正因子(TMNFTB10)、每日溫度為40℃時單葉片同化CO2的初始光能利用效率(EFFTB40)、35℃時葉片生命周期(SPAN)等。

(3)相同氣候類型下不同氣候條件對不同生產水平下的參數敏感性的影響不同。潛在生產水平下冬小麥和夏玉米不同氣候條件TDCC系數分別為0.92和0.98，不同氣候條件對潛在生產水平下影響較小，而水分限制生產水平下2種作物不同氣候條件TDCC系數分別為0.61和0.86，不同氣候條件對水分限制生產水平影響較大，其中冬小麥水分限制生產水平下由不同氣候條件引起的參數敏感性差異相對較大，這均與氣候條件在不同時空上的差異有關。

(4)作物模型是通過作物參數來表達作物對不同外界條件的響應，不同外界條件會導致不同的參數敏感性，外界條件越復雜，參數敏感性差異越大。因此，在特定環境下運用模型時進行敏感性分析是不可缺少的過程。