基于EFAST的CERES-Wheat模型土壤參數敏感性分析

崔金濤 丁繼輝 YESILEKIN Nebi 鄧 升 邵光成

(1.河海大學水文水資源與水利工程科學國家重點實驗室， 南京 210098； 2.河海大學水文水資源學院， 南京 210098；3.河海大學農業科學與工程學院， 南京 210098； 4.佛羅里達大學農業生物工程學院， 蓋恩斯維爾 32608；5.江西省水利科學研究院農村水利研究所， 南昌 330029)

0 引言

隨著農業信息技術的興起，作物模型逐漸成為數字農業和現代農業研究的熱點[1]。自提出植被冠層光能截獲作物生長模型[2]后，經過50多年的發展已經涌現了大量的作物模型，其中較為成熟的作物生長模型主要有DSSAT(Decision support system for agrotechnology transfer)[3-4]、APSIM (Agricultural production system simulator)[5-6]、AquaCrop[7]、WOFOST(Worldfood studies)[8]和STICS (Simulateur multidisciplinaire pour Ies cultures standard)[9]等。與傳統田間試驗相比，作物生長模型借助一系列數學公式，通過對農田生態系統中天氣變化、養分循環、作物生長等過程進行模擬，為作物栽培、水肥管理、生態環境調控等方面提供決策，具有省時、省力、易重現的優點。因此，作物生長模型的應用越來越廣泛。

農業生態系統運行過程十分復雜，涉及眾多的物理過程及生化反應，因此作物模型在擁有眾多優點的同時，也面臨著參數過多、難以獲取的問題[10-11]。以DSSAT模型為例，僅要求輸入的土壤參數就有10多個，如果再考慮土壤異質性進行分層，則需要輸入的土壤參數數量會成倍增加。雖然土壤參數具有明確物理含義，可以根據田間實測數據確定，但是在實際應用中，由于田間土壤的異質性，一些土壤參數的測量值可能變化范圍較大[12]，同時一些參數的測量極其復雜繁瑣，很難對模型所涉及的全部參數進行測量和標定[13]。通常只測量和標定對模型模擬結果影響較大的參數，而固定或簡單處理對模型模擬結果影響較小的參數[11,13-15]。因此，篩選出對模擬結果敏感的關鍵控制參數，從而減少田間實測參數數量、降低模型輸入參數獲取難度，對模型參數率定及后續模型本地化、區域化應用至關重要。

敏感性分析可以確定參數對模擬結果的影響，從而篩選出對模擬敏感度較高的參數。全局敏感性分析不僅可以同時檢驗多個參數變化對模擬結果的影響，還評估了參數對模擬結果的直接和間接影響，因此廣泛應用于地學、農學等領域的模型參數敏感性分析。近年來，不少學者已將全局敏感性分析法成功應用到作物生長模型的參數敏感性分析中[10,13,16-20]。然而，已有研究大多集中于品種參數和田間管理對模擬結果的敏感性分析，研究目的多服務于模型調參率定，而對減少模型土壤參數實測數量及土壤參數簡單處理的研究尚不多見。近年來，在促進農業生產的同時，開始注重對肥料的高效利用和對農業面源污染的防治[21-22]。氮肥作為一種在農業生產中最為常見、施用最多的肥料[21]，常通過氨揮發、淋溶、徑流、硝化-反硝化作用等途徑損失。土壤狀況是影響氮素運移、轉化過程及不同形態氮素分布的重要因素之一[23]，研究氮素在農業生產中的運移、轉化過程及土壤對該過程的影響機理，有助于指導人們科學合理施肥[21]，模型氮素分布敏感性研究可為模型在肥料高效利用及面源污染防治方面的應用提供理論參考。然而，已有研究大多側重于模擬作物生長及產量方面，鮮見涉及氮素運移、轉化過程的研究。

本文以江蘇省漣水縣水利科學研究站為研究區，根據冬小麥田間實測數據，應用EFAST法分析CERES-Wheat模型中的土壤參數對冬小麥生長及氮素分布等模擬的敏感性，旨在降低模型土壤參數輸入要求，為模型本地化應用提供技術支撐。

1 材料與方法

1.1 CERES-Wheat模型及參數

CERES-Wheat模型是DSSAT-CERES模型中的一個子模塊，以天為步長動態模擬冬小麥在不同的氣象、土壤及灌溉、施肥等田間管理下的生長發育[24]。模型通過計算作物冠層截獲有效光合輻射來驅動，涉及水分吸收、光合過程、干物質分配和植株生長以及衰老等基本生理生態過程，模擬作物物候、葉面積變化、生長發育及生物產量形成，并在模擬過程中考慮生長過程中溫度、水分和氮素等脅迫因素影響[24-25]。模型運行需要輸入氣象數據、土壤數據、作物參數和田間管理信息4部分[3,24]。氣象數據至少包括模擬期內每日最低和最高氣溫、每日降雨量和太陽輻射量；土壤數據包括土壤容重、孔隙率、凋萎系數、田間持水率等土壤特性參數，并采用分層形式表達每一層土壤物理、水力、形態結構等特征；作物參數包括生態型、物種型和品種型3種參數，描述作物各器官光合作用、呼吸作用等生理過程，一般在模型率定調參時僅考慮品種型參數，不建議修改生態型和物種型參數[26]；田間管理信息包括種植日期、種植密度、行距、灌溉、施肥類型及施肥量等信息[27]。

本研究所需的試驗數據來源于江蘇省漣水縣水利科學研究站(119.27°E，33.78° N，海拔7.5 m)。該地區屬于暖溫帶季風性半濕潤氣候，多年平均氣溫14.8℃，日照時數2 280 h，降雨量時間分布不均，多集中于5—9月，多年平均降雨量為979 mm，多年平均蒸發量為1 056 mm。模型所需的氣象數據來源于中國氣象數據網(http:∥data.cma.cn)淮安站。氣象數據主要包括日最高氣溫、最低氣溫、降雨量及日照時數等。模型所需的太陽輻射量根據獲得的日照時數由埃斯曲郎經驗公式估算得到[28]。模型所需的土壤數據主要來源于田間試驗的實測數據及模型中相同土壤類型的默認土壤數據。供試區耕地層土壤為砂壤土，0～30 cm土壤容重為1.4 g/cm3，田間持水率為28.7%(質量含水率)。模型所需輸入的土壤參數及取值范圍如表1所示，土壤參數的上下限為試驗區不同地塊實測土壤參數最大值和最小值。

表1 選取的DSSAT模型中土壤參數Tab.1 Definition and intervals of soil parameters in DSSAT model

模型所需的作物參數來源于文獻[29]，該文獻調參并驗證了江蘇省漣水縣水利科學研究站種植的淮麥11的品種參數值，分別為春化天數(P1V)35.57 d，光周期影響因子(P1D)94.81%，灌漿期積溫(P5)707.2℃，開花期冬小麥單位冠層生物量籽粒數(G1)28.86粒/g，標準籽粒質量(G2)23.52 mg，成熟期標準條件下單蘗質量(G3)1.228 g，葉熱間隔(PHINT)60.00℃。模型所需的田間管理數據來源于試驗的觀測記錄和調查。冬小麥于2006年10月15日播種，行距15 cm，播深5 cm，播種密度為400粒/m2。田間管理措施如灌溉、施肥、除草、病蟲害防治等按當地正常田間管理進行。

本研究主要考慮模型的2類輸出，即冬小麥生產指標和氮分布指標。其中生產指標包含地上生物量(CWAM)、產量(HWAM)、收獲指數(HIAM)；氮分布指標包括表征植物各器官中氮素含量的4個輸出量，即地上生物量中氮素含量(CNAM)、籽粒中氮素含量(GNAM)、根中氮素含量(RNAD)及葉中氮素含量(LNAD)，表征土壤中氮素含量變化的4個輸出量，即作物吸收氮素量(NUCM)、土壤淋失氮素量(NLCM)、土壤中硝態氮含量(NITD)及銨態氮含量(NHTD)，表征土壤中氮素轉化的4個輸出量，即氮礦化量(NMNC)、硝化量(NITC)、反硝化量(NDNC)及氨揮發量(AMLC)。

1.2 EFAST方法

擴展傅里葉幅度檢驗(EFAST)法是一種基于方差分解的全局敏感性分析方法，結合了FAST和Sobol法的優點改進而來。其算法簡單介紹如下：

設有模型y=f(x1,x2,…,xm)，通過傅里葉轉換可將其轉換為y=f(s)，轉換函數為

(1)

(2)

(3)

(4)

式中ωi——參數xi的振蕩頻率，i=1,2,…,m

φi——每個參數xi的隨機初相位，取[0,2π]

s——標量變量，取[-π,π]

p——傅里葉變換參數

Ap、Bp——傅里葉振幅

參數xi變化引起的模型輸出方差Vi為

(5)

(6)

p∈Z={-∞,…,-1,1,…,+∞}

模型總方差為

(7)

對標量s在[-π,π]中等間隔取樣，傅里葉振幅Ap和Bp近似計算式為

(8)

(9)

式中Ns——取樣總數

sk——標量s的第k個取樣值

模型總方差可分解為

(10)

式中Vi——參數xi自身變化引起的方差

Vij——參數xi通過參數xj作用貢獻的方差

V12…m——參數x1通過其余m-1個參數相互作用貢獻的方差

通過歸一化處理后，參數xi的一階敏感性指數Si可表示其對模型輸出總方差的直接貢獻，即

(11)

總敏感性指數為

(12)

式中V-i——不包括參數xi的其他所有參數方差之和

1.3 模型參數敏感性分析方案

全局敏感性分析試驗采用SimLab(Version 2.3)進行分析，CERES-Wheat模型運行借助于RStudio進行模擬運行，具體方案如下：

(1)在SimLab中定義模型輸入參數的取值范圍及分布形式。本研究中土壤參數取值范圍如表1所示，并假設參數在范圍內服從均勻分布。

(2)在參數范圍內進行隨機取樣，生成多維參數集。EFAST法認為采樣次數大于參數個數65倍的分析結果才有效，本研究共產生1 455組參數。

(3)將生成的參數集寫入對應的CERES-Wheat模型文件中，運行模型并整理模擬結果。

(4)將模擬生成的數據整理成SimLab可識別的文本格式，通過SimLab進行分析并得到最終的敏感性分析結果。

2 結果與分析

2.1 冬小麥生產指標敏感性分析

由土壤參數對冬小麥生產指標的一階敏感性指數和全局敏感性指數分析結果(圖1)可知，對于地上生物量平均全局敏感指數最大的前3個參數為SDUL、SLHB和SLRO，分別可解釋地上生物量變化方差的70%、57%和48% (圖1a)。一階和全局敏感性分析結果表明，SLHB對CWAM的影響主要通過與其他參數的交互作用實現，而SDUL、SLRO對CWAM的影響主要是其直接效應引起的，與其他參數的交互作用對CWAM的影響較小。在參數設定范圍內，CWAM隨SDUL、SLHB、SLRO參數值增加呈現增加趨勢。除SDUL、SLRO、SLHB、SLNI 4個參數外，其余參數全局敏感性指數均小于0.10，表明對CWAM的影響較小，甚至可忽略不計。與CWAM敏感性結果類似，HWAM和HIAM對SDUL、SLRO、SLHB、SLNI較為敏感(圖1b、1c)，在參數值獲取時需要重點關注，尤其是在目標土壤存在高異質性時。

圖1 土壤參數對冬小麥生產指標的敏感性分析Fig.1 Results of global sensitivity analysis of soil parameters for production indexes of winter wheat

2.2 冬小麥氮素分布指標敏感性分析

土壤參數對冬小麥中氮素分布指標的敏感性結果(圖2)表明，SLHB是影響地上生物量中氮素含量最敏感的參數，可解釋CNAM變化方差的75%，且主要通過與其他參數的交互作用產生影響(圖2a)。在參數設定范圍內，CNAM隨SLHB參數值的增加而增加。SDUL則是對籽粒中氮素含量最敏感的參數，可解釋GNAM變化方差的70%，且其直接效應占主導(圖2b)。在參數設定范圍內，GNAM與SDUL呈正相關，GNAM隨SDUL增加而增加。對地上生物量和籽粒中氮素含量平均一階敏感性指數大于0.10的參數有SDUL、SLRO、SLHB和SLNI，表明這些參數是影響CNAM和GNAM較為敏感的土壤參數，與其余土壤參數相比，在野外獲取這4個參數時需充分考慮土壤異質性的影響。對于根中氮素含量而言，僅有SLNI是引起冬小麥根氮素含量的最關鍵土壤參數，且RNAD與SLNI呈正相關關系；其余參數全局敏感性指數均小于0.10，表明對RNAD變化影響小于10%(圖2c)。對于葉中氮素含量，平均一階和全局敏感性指數大于0.10的參數均為SLNI和SDUL，且其直接效應占主導，而與其他參數的間接效應對LNAD變化的影響較小(圖2d)。在參數設定范圍內，LNAD與SLNI和SDUL均呈正相關。

2.3 土壤氮素形態指標敏感性分析

由土壤參數對模型模擬的作物吸收氮量、氮素淋失量及土壤中硝態氮、銨態氮含量的敏感性分析結果(圖3)可知，影響作物吸收氮量較敏感的參數有SLRO、SDUL、SLNI和SLHB，且SLRO、SDUL和SLNI的直接效應占主導，與其他參數的間接效應影響較小(圖3a)。在參數設定范圍內，NUCM與SLRO、SDUL、SLNI均呈正相關關系。對氮素淋失影響最敏感的參數為SLDR，可解釋其變化方差的65%，且其直接效應占主導(圖3b)。對NLCM，平均一階敏感性指數大于0.10的參數還有SDUL、SBDM和SLNI。雖然SLCL對NLCM的直接影響較小，僅可解釋NLCM變化方差的2%，但其通過與其他參數交互作用所引起的間接影響可解釋NLCM變化方差的26%，表明在研究模擬氮素淋失時，SLCL也是需要重點測定的土壤參數之一。在CERES模型中，NLCM與SLDR、SBDM、SLNI均呈現正相關關系，而與SDUL呈負相關關系。對于土壤硝態氮含量，僅對SLNI和SBDM較為敏感，分別能解釋NITD變化方差的76%和29%(圖3c)。在參數設定范圍內，NITD分別隨SLNI和SBDM增加而增加。土壤銨態氮含量則對SLNI和SLHB變化較為敏感，分別能解釋NHTD變化方差的69%和17%，且均是直接效應占主導(圖3d)。在參數設定范圍內，NHTD與SLNI呈正相關關系，而與SLHB呈負相關關系。在研究模擬氮素分布時，需要重點考慮土壤異質性對SLNI、SLHB和SBDM的影響。

圖3 土壤參數對作物吸收氮量、氮淋失量及土壤中硝態氮、銨態氮含量的敏感性分析Fig.3 Results of global sensitivity analysis of soil parameters for nitrogen uptake and leached, content of total soil nitrate and ammonium nitrogen

2.4 土壤氮素轉化指標敏感性分析

由土壤參數對模型模擬涉及土壤氮素轉化的氮礦化、硝化、反硝化及氨揮發的一階敏感性指數和全局敏感性指數(圖4)可知，NMNC和NITC均對SLNI最敏感，分別解釋其方差變化的89%和79%，而且參數間的交互效應較小，參數的直接影響占主導。其次NMNC和NITC對SBDM也較敏感，而對其余參數的敏感性較差，甚至可忽略不計(圖4a、4b)。在參數設定范圍內，NMNC與SLNI和SBDM均呈正相關關系，與NMNC類似，NITC隨SLNI和SBDM增加而增加。對于氮反硝化量，與NITC的分析結果不同，僅有SLDR對其影響較大，可解釋其方差變化的92%，其余參數對其影響很小，全局敏感性指數均不足0.10(圖4c)。在參數設定范圍內，NDNC與SLDR呈負相關關系，NDNC隨SLDR增大而減小。對于氨揮發量直接影響較大的參數有SLNI、SLHB和SLDR，且與SLNI、SLHB呈正相關關系，而與SLDR呈負相關關系。雖然SBDM、SSKS和SRGF對AMLC的直接影響很小，但通過與其他參數的交互作用可分別解釋AMLC方差變化的12%、11%和10%(圖4d)。

圖4 土壤參數對氮素礦化量、硝化量、反硝化量及氨揮發量的敏感性分析Fig.4 Results of global sensitivity analysis of soil parameters for nitrogen mineralization, nitrification, denitrification and ammonia volatilization

3 討論

擴展傅里葉幅度檢驗法不僅考慮了單個參數對模型輸出的直接影響，而且考慮了參數間的交互作用，解決了局部敏感性分析法僅考慮參數直接效應而忽視參數間間接效應的不足[1]。從全局敏感性分析來看，冬小麥生產指標對土壤水分參數SDUL和SLRO較為敏感。田間持水率反映了土壤持水能力，在生產實踐中，常用來作為灌溉上限和計算灌水定額的指標。模型模擬土壤水分對作物生長的影響主要通過影響土壤蒸發蒸騰量來影響水分在土層中的分布，進而影響作物根系吸水及作物生長[12]。SLRO用來估算降雨徑流量，反映了冬小麥對降雨的利用效率，其取值與土地利用、坡度、土壤類型以及土壤濕度有關。模型輸出的3個生產指標敏感性分析結果存在類似趨勢，主要是因為這3個指標存在較強的相關關系。

氮素是制約作物生長和產量形成的主要因素[22]。在以模擬作物產量等生產能力指標為研究重點時，需要重點確定SDUL和SLRO取值的同時，也不能忽視SLNI和SLHB的影響。對于地上生物量和籽粒中的氮素而言，最敏感的參數是田間持水率，而對根和葉中的氮素而言，最敏感的參數則是土壤中總氮含量，這可能是因為冬小麥根、葉生長主要受氮素含量限制，而進行干物質分配時主要受土壤水分含量的影響。氮素施用對作物生長及干物質形成具有良好的作用，可增加株高、千粒質量及干物質量[21,30]。同時莊恒揚等[22]指出，開花期后作物葉、莖中的氮素開始向籽粒中轉移，而根系的吸氮能力也隨生育期進程而不斷下降。此外，在CERES模型中，模型假設葉片中的氮素優先轉移、補充營養器官[22]，這可能是造成籽粒中氮素含量對田間持水率和土壤總氮含量均敏感而根和葉中氮素含量僅對土壤總氮含量敏感的原因之一。作物吸收氮量敏感性分析結果與作物地上生物量及籽粒中含氮量的結果類似，這也反映了作物吸收的氮素大部分轉移到了籽粒中[21-22,30]。當研究側重于氮素在作物中的分布情況時，應重點考慮土壤異質性對SDUL、SLRO、SLNI和SLHB的影響，對這些參數的取值需更加慎重。

氮素淋失過程實際上是土壤氮素與降雨、徑流相互作用的過程[23,31]，氮素淋失量與土壤導流能力及土壤氮素含量有密切聯系，因此在測定土壤參數時應考慮土壤異質性對SLDR、SLRO、SDUL、SBDM和SLNI的影響。在模型模擬過程中，土壤氮素淋失量隨SLDR、SLRO、SBDM和SLNI的增加而增加，但隨SDUL的增大而減少，這可能是田間持水率增加可減少排水量，從而導致氮素淋失量降低。土壤中的含氮物質在各種微生物及一系列物理、化學過程的作用下會發生形態及狀態的轉化，氮素礦化是其中的一種生物化學過程，氮素礦化過程會受土壤質地、含氮量、溫度和濕度等諸多因素的影響[23]。在本研究中，CERES-Wheat模型土壤氮素礦化主要受土壤氮素含量及土壤容重的影響，其他土壤參數的影響效應很小。

硝態氮和銨態氮是土壤中無機氮的兩種存在形式[23,30]，土壤總氮含量是影響硝態氮和銨態氮含量最敏感因素，此外，硝態氮含量還受土壤容重的影響，而銨態氮含量還受土壤酸堿度的影響，而且隨土壤堿性增加而減少。氮素硝化反應是影響土壤硝態氮含量的重要因素，對氮素硝化量的敏感性分析結果與硝態氮含量一致，土壤總氮含量和土壤容重是影響其變化的敏感參數。而氮素反硝化量主要受排水比率的影響，這主要由于反硝化作用是反硝化細菌在缺氧條件下還原硝酸鹽的過程[32]，而模型中排水比率決定了土壤中水分排出速率及氧氣狀況。對于氨揮發量而言，除土壤總氮含量對其敏感外，與土壤含水率有關的參數和土壤酸堿度對其也比較敏感，這也與氨氣易融于水及其溶液呈堿性易與酸性溶液中和的特性相一致。在研究土壤氮素轉化及不同形態氮素分布時，需要同時考慮土壤異質性對SLDR、SLRO、SDUL、SLLL、SRGF、SSKS、SBDM、SLNI和SLHB的影響，而其他參數對土壤氮素含量影響較小，可適度降低土壤參數獲取難度。

4 結論

(1)對于冬小麥生產指標，田間持水率、徑流曲線數、土壤酸堿度及土壤總氮含量是影響冬小麥地上生物量和產量較為敏感的土壤參數，需要考慮土壤異質性的影響，而其他參數的影響較小，可將原來需測參數數量降低至4個，且冬小麥生產指標與這4個參數均呈正相關關系，降低了參數獲取難度。

(2)對于冬小麥氮素分布指標，田間持水率和徑流曲線數是影響地上生物量及冬小麥籽粒中氮素含量較為敏感的土壤參數，而土壤總氮含量則是影響冬小麥根和葉中氮素含量最敏感的參數。在研究作物氮素分布模擬時，需重點測量土壤田間持水率、徑流曲線數、土壤酸堿度及土壤總氮含量4個參數。

(3)對于氮素轉化與分布指標，田間持水率和徑流曲線數是影響作物從土壤中吸收氮素量較為敏感的參數，而排水比率是影響氮素淋失量、反硝化量最敏感的參數，土壤總氮含量則是影響土壤中硝態氮、銨態氮、氮礦化、氮硝化量最敏感的參數，土壤總氮含量、土壤酸堿度及排水比率則是影響氮素氨揮發量較為敏感的土壤參數。在模擬氮素轉化與分布時，需重點考慮土壤異質性對排水比率、徑流曲線數、凋萎含水率、飽和含水率、飽和導水率、根系影響因子、土壤容重、總氮含量及土壤酸堿度的影響，原需測量的15個土壤參數可最大簡化為9個重點測量參數。