何力鴻，莫興國，劉蘇峽

（1. 中國科學(xué)院地理科學(xué)與資源研究所陸地水循環(huán)及地表過程重點實驗室，北京 100101；2. 中國科學(xué)院大學(xué)，北京 100049；3. 中國科學(xué)院大學(xué)資環(huán)學(xué)院/中丹學(xué)院，北京 100049）

0 引 言

基于動力學(xué)過程的生態(tài)系統(tǒng)模型具有結(jié)構(gòu)復(fù)雜、輸入?yún)?shù)多和參數(shù)空間變異性強等特點，使得參數(shù)校驗成為模型研究和應(yīng)用的難點[1]。而參數(shù)敏感性分析可有效識別模型敏感參數(shù)，通過固定不敏感參數(shù)，可有效降低運算成本[2]。

參數(shù)敏感性分析可分為局部和全局敏感性分析。其中全局敏感性分析方法充分考慮不同參數(shù)及其相互作用對模型模擬結(jié)果的影響，在非線性過程模型的參數(shù)校準(zhǔn)中得到廣泛應(yīng)用[3-5]。全局敏感性分析方法可分為定性和定量分析，其中定性敏感性方法研究不同參數(shù)對目標(biāo)函數(shù)影響的相對大小，如 Morris 法[6]、多元自適應(yīng)回歸樣條法（Multivariate Adaptive Regression Splines，MARS)[7]、Delta Test 法[8]和 Sum of Trees 法[9]等；而定量敏感性方法可獲得不同參數(shù)對目標(biāo)函數(shù)的貢獻率，如 Sobol'法[10]、McKay 法[11]和擴展傅里葉幅度檢驗法 （Extended Fourier Amplitude Sensitivity Test，EFAST)[12]等。近年來，國內(nèi)外學(xué)者對參數(shù)敏感性做了諸多分析，如梁浩等[13]通過Morris和Sobol'法分析WHCNS 模型中與土壤硝態(tài)氮和含水率相關(guān)的敏感參數(shù)；邢會敏等[14]采 EFAST 法對AquaCrop 模型中的作物參數(shù)進行分析。DeJonge 等[15]使用 Morris 和 Sobol'方法分析了不同灌溉處理下CERES-Maize 模型中作物參數(shù)對葉面積指數(shù)、蒸散和產(chǎn)量的影響，并基于敏感度分析結(jié)果提出了參數(shù)校準(zhǔn)方案。

在參數(shù)分析過程中，合適的試驗設(shè)計可減少抽樣次數(shù)，降低運算成本。目前試驗設(shè)計方法主要有蒙特卡洛（Monte Carlo，MC）法[16]、正交試驗設(shè)計（Orthogonal Experimental Design，OED）[17]和拉丁超立方（Latin Hypercube，LH）抽樣[18]等。其中，MC 通過參數(shù)空間上的概率密度函數(shù)隨機生成樣本點，其效率取決于隨機數(shù)的屬性。而隨機數(shù)抽樣容易產(chǎn)生聚類，出現(xiàn)空白抽樣空間，隨機添加新抽樣點不一定能填補空白抽樣空間[2]。有限數(shù)量的抽樣點的正交質(zhì)量取決于點分布的均勻性，為了保證所有的抽樣區(qū)域都能被抽樣點覆蓋，McKay 等[11]提出了LH 抽樣，該方法是分層抽樣的一種形式，可減少通過蒙特卡洛法估算時產(chǎn)生的誤差。后來研究對LH 抽樣進行了改進，如Tang[19]通過正交陣列（Orthogonal Array，OA）構(gòu)建了OALH（Orthogonal Array-Based Latin Hypercubes），與LH 相比，OALH 更適用于計算機試驗和數(shù)值積分。除上述抽樣方法外，還有針對特定敏感性方法設(shè)計的抽樣方式，如 Morris 抽樣[6]和Sobol'抽樣[10]等。

以往研究多基于一種或幾種敏感性分析方法篩選敏感性參數(shù)，對比分析不同敏感性方法有效性及效率的研究不夠深入。本研究基于VIP（Vegetation Interface Processes）模型探討不同敏感性方法和抽樣方式在參數(shù)敏感性分析過程中的有效性及其效率，VIP 模型綜合考慮了作物生長、碳氮循環(huán)過程、土壤水分運動等過程，經(jīng)過多年發(fā)展，模型已用于華北平原單點及區(qū)域作物產(chǎn)量、蒸散量及水分利用效率的模擬[20]。本研究以華北平原土壤硝態(tài)氮模擬為例，對比分析了PSUADE（Problem Solving environment for Uncertainty Analysis and Design Exploration）提供的8 種敏感性分析方法（包括6 種定性和2 種定量敏感性分析方法）在篩選農(nóng)田土壤硝態(tài)氮相關(guān)敏感參數(shù)時的有效性及其效率。

1 材料與方法

1.1 VIP 模型及參數(shù)

VIP 模型綜合考慮了植物生理過程、碳氮循環(huán)過程、土壤水分運動等過程，能夠合理模擬單點和區(qū)域作物生長和水量平衡過程。植被動態(tài)變化采用干物質(zhì)累積與分配方程來描述、土壤熱量傳輸采用熱擴散方程、土壤水量平衡方程遵循達西定律[20]。其中，VIP 模型土壤氮循環(huán)模塊考慮了土壤有機氮庫之間的周轉(zhuǎn)、礦化過程、無機氮庫之間的轉(zhuǎn)化過程。本研究選取了與土壤有機質(zhì)周轉(zhuǎn)、礦化、硝化過程、反硝化過程和氨揮發(fā)過程相關(guān)的參數(shù)（表1）來評估不同敏感性分析方法。

表1 VIP 模型中氮循環(huán)過程相關(guān)參數(shù)Table 1 Parameters related to nitrogen cycle process in VIP(Vegetation Interface Processes) model

1.2 試驗站點及數(shù)據(jù)

研究樣地位于華北平原 (31°14'～40°25'N，112°48'～122°45'E)，面積約33 萬km2，其中農(nóng)業(yè)用地占70%。屬半干旱半濕潤暖溫帶大陸性季風(fēng)氣候，夏季降雨量最多[21]。中部地區(qū)的年平均氣溫在10～15 ℃之間，日照時間約為2 800 h。土壤為黃河沉積物發(fā)育而成的輕壤質(zhì)潮土，耕層容重為1.48 g/m3，有機質(zhì)質(zhì)量分數(shù)為11 g/kg。該地區(qū)玉米于每年6 月初播種，9 月中下旬收獲。本研究選擇了3 個站點（封丘，禹城和欒城農(nóng)業(yè)生態(tài)站）的數(shù)據(jù)對VIP模型進行參數(shù)敏感性分析。

1.2.1 驅(qū)動數(shù)據(jù)

氣象數(shù)據(jù)、土壤及作物參數(shù)等數(shù)據(jù)用于驅(qū)動 VIP 模型，其中氣象數(shù)據(jù)包括最大溫度、最小溫度、平均溫度、降水量、日照時數(shù)、相對濕度、風(fēng)速。封丘，禹城和欒城站日尺度氣象數(shù)據(jù)從國家生態(tài)科學(xué)數(shù)據(jù)中心[21]獲得。另一類驅(qū)動數(shù)據(jù)是作物參數(shù)與土壤參數(shù)，其中涉及的作物參數(shù)主要有：有效積溫參數(shù)、作物不同時期含氮量；土壤參數(shù)主要有不同有機質(zhì)分解速率、碳氮比、氨揮發(fā)速率、反硝化速率的相關(guān)參數(shù)，由文獻資料獲得[22-23]。封丘、禹城和欒城農(nóng)業(yè)生態(tài)站田間管理措施數(shù)據(jù)從文獻資料[24-25]和國家生態(tài)科學(xué)數(shù)據(jù)中心[21]獲得。表 2 中列出了播種日期，收獲日期，氮肥和灌溉的詳細信息?？紤]到研究目的、數(shù)據(jù)的質(zhì)量等因素，本研究封丘站、禹城站和欒城站模擬時間段分別為2008年6月—2009年9月，2001 年 6 月—2002 年 9 月和 2004 年 6 月—2005 年 9 月，時間步長定為1 d。

表2 3 個農(nóng)業(yè)生態(tài)站玉米生長期及田間管理數(shù)據(jù)Table 2 Growing season and field management data of maize at three agro-ecological stations

1.2.2 驗證數(shù)據(jù)

本研究選定 3 個農(nóng)業(yè)生態(tài)站（封丘、禹城和欒城）樣地實測土壤硝態(tài)氮數(shù)據(jù)用于驗證 VIP 模型模擬結(jié)果。在玉米生長季，以20 cm 的距離分層收集土壤樣品，并在分層中的 2 個點隨機采集每個樣區(qū)。在每季作物收獲后采集土樣，用50 mL 2 mol/L KCL 提取10 g 每個土壤樣品，采用雙波長（220 nm，275 nm）法測定[26]。

1.3 敏感性分析方法

在以往農(nóng)業(yè)生態(tài)系統(tǒng)模型參數(shù)敏感性分析研究中，多數(shù)研究基于一種或幾種敏感性分析方法篩選敏感性參數(shù)，而對比分析不同敏感性方法有效性及效率研究較少。本研究選取了在生態(tài)模型中廣泛應(yīng)用的 Spearman秩相關(guān)系數(shù)（Spearman’s correlation coefficient，SPEA）法、Morris 法、多元自適應(yīng)回歸樣條（Multivariate Adaptive Regression Splines，MARS)、Delta Test（DT)法、Sum of Trees（SOT）法、Gaussian Process（GP）法、McKay 法和Sobol'法，對比分析這8 種全局敏感性分析方法在參數(shù)敏感性分析中的有效性及效率。8 種全局敏感性分析方法由 PSUADE（Problem Solving environment for Uncertainty Analysis and Design Exploration）軟件提供，PSUADE 為復(fù)雜系統(tǒng)模型提供不確定性分析的集成設(shè)計和分析環(huán)境，可通過 https://computing.llnl.gov/projects/psuade-uncertainty-quantification獲得。在參數(shù)敏感性分析過程中，首先指定抽樣方式，樣本量大小以及參數(shù)取值范圍，然后生成樣本組合。本研究使用抽樣方式有蒙特卡羅（Monte Carlo，MC）、METIS、Morris、拉丁超立方（Latin Hypercube，LH）、LPTAU、正交陣列 （Orthogonal Array，OA)和基于拉丁超立方的正交陣列（Orthogonal Array-Based Latin Hypercubes，OALH）抽樣方式。通過腳本將 VIP 模型可執(zhí)行文件（vip.exe）與PSUADE 進行鏈接，并將運行結(jié)果讀寫到輸出文件，利用PSUADE 提供的統(tǒng)計方法分析不同參數(shù)組合與目標(biāo)函數(shù)之間關(guān)系。

在評估敏感性分析方法有效性時，將敏感性分析結(jié)果進行歸一化，通過對比不同敏感性分析結(jié)果及關(guān)鍵參數(shù)的篩選來判斷敏感性方法的有效性，如潛在硝化速率（knit）對土壤硝態(tài)氮含量影響較大，若篩選結(jié)果判定其為敏感參數(shù)，則證明該敏感性分析方法有效。在評估敏感性方法效率時，通過篩選敏感參數(shù)所需樣本量大小來衡量不同敏感性方法效率，樣本量越大，效率越低，反之效率越高。在評估Morris 敏感性方法效率時，樣本量是n+1（n是參數(shù)個數(shù)）的整數(shù)倍，為探究篩選敏感參數(shù)所需最佳樣本量，設(shè)置樣本量分別為85, 170,…, 3400，每次增加 85 直至最大樣本量；Sobol'敏感性分析方法篩選敏感性參數(shù)過程中，樣本量是n+2（n是參數(shù)個數(shù)）的整數(shù)倍，本研究設(shè)置樣本量分別為90, 180,…, 3420，每次增加 90 直至最大樣本量；在評估 DT、MARS、Mckay和SOT 敏感性方法效率時，通過設(shè)置不同抽樣方式（MC，LH，OA，OALH，LPTAU 和METIS）和不同樣本量梯度分析有效敏感性分析方法的效率，其中 MC，LH，LPTAU 和 METIS 抽樣方式的樣本量設(shè)置為 85, 170,…,3400，樣本量每次增加85 直至最大樣本量；OA 和OALH抽樣法的樣本量為m·pt，本研究給定m=1 和t= 2，樣本量由素數(shù)p（17, 19, 23,…, 59）確定，樣本量分別為289,361, 529,…, 3481。樣本量從小依次增加到最大值（最大樣本量分析結(jié)果作為基準(zhǔn)），對比不同樣本量下敏感性分析結(jié)果與最大樣本量的結(jié)果，確定篩選出敏感參數(shù)所需的最少樣本量。

2 結(jié)果與分析

2.1 敏感性分析方法有效性

2.1.1 定性敏感性分析方法有效性

將 3 個農(nóng)業(yè)生態(tài)站定性敏感性方法（SPEA、DT、GP、MARS、Morris 和SOT）的分析結(jié)果進行加權(quán)平均，為了方便統(tǒng)一對比，對結(jié)果再進行歸一化處理，本研究將敏感性分析歸一化結(jié)果大于0.2 的參數(shù)作為敏感參數(shù)。對VIP 模型中16 個氮循環(huán)參數(shù)進行排序（圖1），結(jié)果表明，DT、MARS、Morris 和SOT 法篩選出與土壤硝態(tài)氮含量相關(guān)的敏感參數(shù)有：knit、kur、mnit、rmb、rsh、kmb、ksl和kml。其中，knit和kur在上述4 種敏感性分析方法篩選結(jié)果中表現(xiàn)出高敏感性，而SPEA 和GP 法篩選結(jié)果表明knit和kur為不敏感性參數(shù)，與其他定性敏感性分析結(jié)果不一致。此外，SPEA 法也無法有效識別kmb，GP 法無法識別ksl，而參數(shù)kmb和ksl在 DT、MARS、Morris 和 SOT法的結(jié)果中都表現(xiàn)為敏感參數(shù)。

2.1.2 定量敏感性分析方法有效性

在 McKay 和 Sobol'定量敏感性分析方法的有效性（圖 2）驗證過程中，2 種方法樣本量都設(shè)置為3 400。為了方便統(tǒng)一對比，將 3 個農(nóng)業(yè)生態(tài)站定量敏感性分析結(jié)果進行歸一化處理。結(jié)果表明，McKay 法和 Sobol'法都能有效篩選出 3 個農(nóng)業(yè)生態(tài)站土壤硝態(tài)氮相關(guān)的敏感參數(shù)，敏感參數(shù)分別為knit、kur、mnit、ksl、kml、ksh、kmb、rmb和rsh，且篩選出的最敏感參數(shù)都是knit，這與定性敏感性分析結(jié)果一致。

圖1 不同定性敏感性分析方法的參數(shù)敏感性Fig.1 Parameter sensitivity of different qualitative sensitivity analysis methods

2.2 敏感性分析方法效率

由于SPEA 和GP 法敏感性分析結(jié)果與其他敏感性分析方法結(jié)果不一致，未能將最敏感參數(shù)knit有效篩選出來，未進一步分析SPEA 和GP 法效率。本研究通過6 種不同的抽樣方法（MC，LH，OA，OALH，LPTAU 和METIS）來明確適用于MARS、DT、SOT 和Mckay 敏感性分析法的抽樣方式及最小樣本量。在對比不同敏感性分析方法效率時，將每種抽樣方式最大樣本量篩選結(jié)果作為敏感性方法有效性基準(zhǔn)，當(dāng)連續(xù) 3 個樣本量的敏感參數(shù)與最大樣本量篩選結(jié)果相同或接近時，則認為連續(xù)3 個樣本量中最小樣本量足以用于抽樣分析。對不同抽樣方式及樣本量的分析結(jié)果進行統(tǒng)計（表3），結(jié)果發(fā)現(xiàn)，MARS敏感性分析過程中，通過 OA 和 OALH 抽樣時，所需樣本量僅為361，這表明OA 和OALH 抽樣方式最適合MARS 法；在利用DT 和SOT 敏感性分析中發(fā)現(xiàn)，MC最適合DT 和SOT 敏感性分析方法，所需最小樣本量分別為425 和510；與MARS 和DT 敏感性分析方法相比，SOT 在篩選敏感參數(shù)過程中所需樣本量更多。在McKay 篩選敏感參數(shù)的過程中發(fā)現(xiàn)，與其他抽樣方式相比，OALH 抽樣需要的樣本點最?。颖玖繛?29），這表明，OALH 抽樣方式最適應(yīng)于McKay 敏感性分析方法。

圖2 3 個農(nóng)業(yè)生態(tài)站McKay 和Sobol'敏感性分析方法的參數(shù)敏感性Fig.2 Parameter sensitivity of Mckay and Sobol' methods at three agro-ecological stations .

表3 3 個農(nóng)業(yè)生態(tài)站不同抽樣方式下MARS、DT、SOT 和Mckay 敏感性分析方法所需最小樣本量Table 3 Minimum sample size required for MARS, DT, SOT and Mckay sensitivity analysis methods under different sampling methods at three agro-ecological stations

在分析Morris 法和Sobol'法篩選敏感參數(shù)效率時，設(shè)置不同樣本量來確定篩選華北平原 3 個農(nóng)業(yè)生態(tài)站（封丘、禹城和欒城）土壤硝態(tài)氮敏感參數(shù)所需的最小樣本量。不同參數(shù)對應(yīng)均值μ在4 kg/ hm2處會被明顯區(qū)分開，故本研究認為均值μ≥ 4 kg/ hm2的參數(shù)為敏感參數(shù)。結(jié)果表明，樣本量為 85 時，封丘、禹城和欒城站篩選出的敏感參數(shù)為knit和kur；隨著樣本量增加，mnit、ksl、kml、ksh、kmb、rmb和rsh被篩選出來；樣本量大于 340篩選出的敏感參數(shù)基本一致，這表明樣本量取340 時即能有效篩選 VIP 模型中土壤氮循環(huán)敏感參數(shù)。在分析Sobol'法篩選敏感參數(shù)效率時，在分析Sobol'法篩選敏感參數(shù)效率時，當(dāng)參數(shù) Sobol'一階指數(shù)大于 0.4 時，認為該參數(shù)為敏感參數(shù)。通過對比不同樣本量下敏感性分析結(jié)果發(fā)現(xiàn)（表 4），除了禹城生態(tài)站，當(dāng)樣本量為 720時，封丘和欒城生態(tài)站的kmb、ksh都不能被有效篩選出來；而樣本量大于810 的敏感性分析結(jié)果基本一致，篩選出的敏感參數(shù)分別為knit、kur、mnit、ksl、kml、ksh、kmb、rmb和rsh，這表明該研究中樣本量選取 810 就能有效篩選敏感參數(shù)。

表4 3 個農(nóng)業(yè)生態(tài)站Morris 和Sobol'敏感性分析方法所需最小樣本量Table 4 Minimum sample size required for Morris and Sobol'sensitivity analysis methods at three agro-ecological stations

3 討 論

3.1 不同敏感性分析方法有效性

8 種全局敏感性方法運用于土壤氮循環(huán)參數(shù)進行了全局敏感性分析。從全局敏感性分析結(jié)果來看，SPEA 和GP 法敏感參數(shù)篩選結(jié)果與其他方法的結(jié)果不同，主要由于 SPEA 指數(shù)是一種秩變換的統(tǒng)計量，適用于參數(shù)和模型響應(yīng)之間的線性關(guān)系，而大多數(shù)生態(tài)系統(tǒng)模型中參數(shù)與模擬結(jié)果呈非線性關(guān)系，如Duan 等[27]分析了SIXPAR模型模擬結(jié)果對參數(shù)的響應(yīng)，盡管 SIXPAR 模型結(jié)構(gòu)簡單，但結(jié)果表明不同參數(shù)之間表現(xiàn)出強烈的非線性交互作用。Gan 等[2]通過SAC-SMA 模型對比不同敏感性分析方法時發(fā)現(xiàn)，SPEA 和GP 法無法有效篩選出敏感參數(shù)，這與本研究結(jié)果類似。此外，與SOT 法相比，通過MARS和DT 法篩選出的敏感參數(shù)結(jié)果接近，且部分邊緣敏感參數(shù)敏感性高于SOT 法得到的結(jié)果，由于邊緣敏感參數(shù)與不敏感參數(shù)之間存在較大差異，可以推斷 MARS 和 DT法可能更有效識別不敏感參數(shù)；與MARS 和DT 法相比，SOT 法篩查所采用的不同抽樣方式所需的樣本量差異更大。這與Gan 等[2]研究結(jié)果類似。無論采用何種抽樣方式，SOT 法篩選比MARS 和DT 法篩選需要更多的樣本量來識別不敏感的參數(shù)，這支持了MARS 和DT 法擅長識別不敏感的參數(shù)，而SOT 法擅長識別高度敏感的參數(shù)的推斷。此外，與上述定性敏感性分析法相比，McKay 和Sobol′定量敏感性分析法都能有效識別敏感參數(shù)，但這2 種方法在敏感性分析過程中運算時間較長，所需計算成本較高。

3.2 不同敏感性分析方法效率

本研究為了確定有效篩選的最小樣本量，在使用Morris 方法時，當(dāng)樣本量低于 340 時，部分敏感參數(shù)未能被篩選出來，這表明樣本量低于 340 時不能有效篩選VIP 模型中氮循環(huán)敏感參數(shù)，這與以往研究結(jié)果類似[28]。通過對比不同樣本量下Sobol'一階指數(shù)結(jié)果發(fā)現(xiàn)，當(dāng)樣本量小于810 時，rmb和rsh不能被有效篩選出來。與Sobol'法相比，Morris 方法所需樣本量較小，對于復(fù)雜碳氮模型，隨著樣本量的增加，所需運算成本也增加，而Morris法在樣本量較少的情況下就能篩選出較敏感參數(shù)，所以在初步篩選復(fù)雜模型的敏感參數(shù)時推薦使用 Morris 方法。McKay 和 Sobol′方法的敏感性指標(biāo)豐富，可以從多角度考慮，得到可靠的結(jié)果，但McKay 和Sobol′方法需要基于較大樣本量才能得到可靠的定量敏感性分析結(jié)果，不適用于參數(shù)較多的復(fù)雜動力系統(tǒng)模式的參數(shù)篩選。在實際的應(yīng)用中，要根據(jù)研究目的選擇最恰當(dāng)?shù)姆椒?。如需定量分析參?shù)對模擬結(jié)果影響，并不考慮模型運算成本問題，推薦使用McKay 和Sobol′方法。

4 結(jié) 論

本研究以VIP（Vegetation Interface Processes）模型模擬華北平原 3 個農(nóng)業(yè)生態(tài)站（封丘、禹城和欒城）土壤硝態(tài)氮為例，評估不同敏感性分析方法有效性及其效率。通過不同定性敏感性分析結(jié)果發(fā)現(xiàn)，SPEA（Spearman’s correlation coefficient）和 GP（Gaussian Process）法篩選結(jié)果與 DT（Delta Test)、MARS（Multivariate Adaptive Regression Splines）、Morris 和 SOT（Sum of Trees）法的不一致；與MARS、DT 和SOT 定性敏感性分析方法相比，Morris 在篩選敏感參數(shù)所需樣本量最少；蒙特卡羅（Monte Carlo，MC），正交陣列（Orthogonal Array，OA）和基于拉丁超立方體的正交陣列（Orthogonal Array-Based Latin Hypercubes，OALH）抽樣方式都適用于MARS、DT 和SOT 敏感性分析方法。OA 和OALH 抽樣方式最適合MARS 法；MC 最適合DT和SOT 敏感性分析方法；OALH 最適應(yīng)于McKay 敏感性分析方法。與MARS 和DT 敏感性分析方法相比，SOT在篩選敏感參數(shù)過程中所需樣本量更多。

在篩選敏感性參數(shù)過程中，定性敏感性分析方法比定量敏感性方法所需樣本量少，但無法對敏感參數(shù)進行定量描述。對參數(shù)較多的復(fù)雜系統(tǒng)模型，可首先使用定性敏感性分析方法進行初步敏感性參數(shù)篩選，以較低的計算成本剔除不敏感的參數(shù)，然后采用定量敏感性分析方法對系統(tǒng)模型中的敏感參數(shù)進行定量分析。