梁　浩，胡克林，李保國（中國農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院，農(nóng)業(yè)部華北耕地保育重點(diǎn)實(shí)驗(yàn)室，北京 100193）

?

基于PEST的土壤-作物系統(tǒng)模型參數(shù)優(yōu)化及靈敏度分析

梁浩，胡克林※，李保國
（中國農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院，農(nóng)業(yè)部華北耕地保育重點(diǎn)實(shí)驗(yàn)室，北京 100193）

摘要：農(nóng)業(yè)生產(chǎn)管理系統(tǒng)模型輸入?yún)?shù)多，參數(shù)率定過程十分耗時(shí)費(fèi)力，大大限制了其推廣應(yīng)用。該研究以華北平原2 a的冬小麥-夏玉米田間試驗(yàn)觀測數(shù)據(jù)為基礎(chǔ)，使用PEST（parameter estimation）參數(shù)自動(dòng)優(yōu)化工具對土壤-作物-大氣系統(tǒng)水熱碳氮過程藕合模型（soil water heat carbon and nitrogen simulator，WHCNS）的土壤水力學(xué)參數(shù)、氮素轉(zhuǎn)化參數(shù)和作物遺傳參數(shù)進(jìn)行自動(dòng)尋優(yōu)，同時(shí)計(jì)算分析模型參數(shù)的相對綜合敏感度，并將優(yōu)化結(jié)果與土壤實(shí)測水力學(xué)參數(shù)和試錯(cuò)法的模擬結(jié)果進(jìn)行比較。參數(shù)敏感度分析結(jié)果表明，18個(gè)模型參數(shù)的相對綜合敏感度較高，其中土壤水力學(xué)參數(shù)普遍具有較高的敏感度，以飽和含水率敏感度最高；作物參數(shù)中，作物生長發(fā)育總積溫和最大比葉面積具有較高的綜合敏感度；而氮素轉(zhuǎn)化參數(shù)的敏感度遠(yuǎn)低于土壤水力學(xué)參數(shù)和作物參數(shù)。評價(jià)模型模擬效果的統(tǒng)計(jì)性指標(biāo)（均方根誤差、模型效率系數(shù)和一致性指數(shù)）表明，PEST法比實(shí)測水力學(xué)參數(shù)的模擬精度有所提高，其中土壤含水率、土壤硝態(tài)氮含量、作物產(chǎn)量和葉面積指數(shù)的均方根誤差分別降低了61.8%、23.5%、73.6%和23.3%。同時(shí)PEST法比試錯(cuò)法對土壤水分和作物產(chǎn)量的模擬精度也有較大提高，但對土壤氮素和葉面積指數(shù)的模擬精度提高不明顯。由于該方法大大節(jié)約了模型校準(zhǔn)時(shí)間，在較短的時(shí)間內(nèi)獲得了明顯高于試錯(cuò)法的模擬精度，因此PEST軟件在WHCNS模型參數(shù)自動(dòng)優(yōu)化中是一個(gè)值得推廣的工具。

關(guān)鍵詞：靈敏度分析；作物；模型；PEST；參數(shù)優(yōu)化；氮循環(huán)；WHCNS

梁浩，胡克林，李保國. 基于PEST的土壤-作物系統(tǒng)模型參數(shù)優(yōu)化及靈敏度分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2016，32（3）：78－85.doi：10.11975/j.issn.1002-6819.2016.03.012http://www.tcsae.org

Liang Hao, Hu Kelin, Li Baoguo. Parameter optimization and sensitivity analysis of soil-crop system model using PEST[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(3): 78－85. (in Chinese with English abstract)doi：10.11975/j.issn.1002-6819.2016.03.012http://www.tcsae.org

Email：checkmyself@qq.com

Email：hukel@cau.edu.cn

0　引　言

土壤-作物系統(tǒng)模型在作物估產(chǎn)、水肥優(yōu)化管理和區(qū)域農(nóng)業(yè)生產(chǎn)決策等方面應(yīng)用非常廣泛[1]，而建立在動(dòng)力學(xué)基礎(chǔ)之上的農(nóng)業(yè)生產(chǎn)系統(tǒng)模型具有輸入?yún)?shù)多、土壤和作物參數(shù)的時(shí)空變異性強(qiáng)、參數(shù)存在較大的不確定性等特點(diǎn)。在模型應(yīng)用之前需要對輸入?yún)?shù)進(jìn)行率定，而參數(shù)之間往往存在著非線性關(guān)系。因此，土壤-作物系統(tǒng)模型參數(shù)的校驗(yàn)一直被認(rèn)為是模型研究和應(yīng)用的難點(diǎn)[2-3]。目前大多數(shù)模型工作者采用傳統(tǒng)的試錯(cuò)法來調(diào)試模型參數(shù)，該方法需要根據(jù)研究者自身的知識和對模型的了解程度對輸入?yún)?shù)進(jìn)行主觀的調(diào)試，不僅耗時(shí)費(fèi)力，而且模型的模擬精度也不是令人十分滿意[4-5]。另外，由于尺寸效應(yīng)的存在，實(shí)驗(yàn)室測得的土壤水力學(xué)等參數(shù)很難直接運(yùn)用到模型中。這些因素大大限制了模型的推廣和應(yīng)用。

近年來，發(fā)展了一些參數(shù)自動(dòng)優(yōu)化工具，如HYDRUS-1D[6]、universal inverse code (UCODE)[7]、parameter estimation（PEST）[8]，shuffle complex evolution method（SCE-UA）算法[9]、個(gè)體種算法（improved adaptive genetic algorithm， IAGA）[10]等，其中PEST采用基于改進(jìn)的非線性最小二乘優(yōu)化算法，已經(jīng)在水文學(xué)模型中得到了廣泛運(yùn)用[11-13]。如Elci等[11]使用PEST工具自動(dòng)率定了地下水模型（modular three-dimensional finite-difference ground-water flow model，MODFLOW）的水力學(xué)參數(shù)，并分析了土耳其伊茲密爾地區(qū)地下水污染物的擴(kuò)散情況。董艷輝等[12]應(yīng)用并行PEST算法優(yōu)化了MODELOW模型的參數(shù)，結(jié)果表明并行化的PEST可以將地下水模型參數(shù)優(yōu)化效率提高3.7倍。Kim等[13]將PEST工具用于率定降雨-徑流模型（solute transport with infiltration and runoff，STWIR），研究發(fā)現(xiàn)PEST對降雨產(chǎn)流參數(shù)有較高的率定效率。

參數(shù)自動(dòng)優(yōu)化方法在土壤-作物系統(tǒng)模型方面的運(yùn)用并不多見，現(xiàn)有的研究主要是對土壤參數(shù)或作物參數(shù)的反演[10,14-17]，如Fang等[14]利用PEST程序估算了root zone water quality model（RZWQM）模型土壤水力學(xué)參數(shù)，發(fā)現(xiàn)PEST軟件在土壤水力學(xué)參數(shù)尋優(yōu)方面表現(xiàn)出了明顯的優(yōu)勢和較高的效率。莊嘉祥等[10]使用優(yōu)勢遺傳算法（individual advantages genetic algorithm，IAGA）成功估算了ORYZA2000水稻模型中的作物遺傳參數(shù)，表明IAGA方法對模型作物參數(shù)的估算具有較高的準(zhǔn)確性。Nolan等[15]將PEST程序用于RZWQM模型中，估計(jì)了土壤水力學(xué)參數(shù)和氮素轉(zhuǎn)化參數(shù)。Rashad等[3]利用PEST軟件優(yōu)化了生物地球化學(xué)模型DayCent的氮素轉(zhuǎn)化參數(shù)，并利用PEST計(jì)算得到的敏感度雅克比矩陣分析了溫室氣體排放的不確定性，認(rèn)為PEST能較好的反演DayCent模型的氮素轉(zhuǎn)化參數(shù)。上述研究分別對土壤水力學(xué)、作物遺傳和氮素轉(zhuǎn)化參數(shù)進(jìn)行了優(yōu)化，目前將土壤水力學(xué)參數(shù)、氮素轉(zhuǎn)化參數(shù)和作物參數(shù)同時(shí)進(jìn)行優(yōu)化的研究還不多見。

課題組最近開發(fā)了土壤-作物-大氣系統(tǒng)耦合模型（soil water heat carbon nitrogen simulator，WHCNS），該模型已經(jīng)在農(nóng)業(yè)生產(chǎn)系統(tǒng)水氮管理[18-19]和有機(jī)無機(jī)配施方面[20]得到了一定的應(yīng)用，但是由于該模型參數(shù)較多，模型校驗(yàn)過程十分復(fù)雜，大大限制了該模型的推廣應(yīng)用。另外，目前還未對該模型參數(shù)進(jìn)行全面的靈敏度分析。因此，本研究以華北平原2 a的冬小麥-夏玉米田間試驗(yàn)觀測數(shù)據(jù)為基礎(chǔ)，利用PEST參數(shù)優(yōu)化工具對土壤水力學(xué)參數(shù)、氮轉(zhuǎn)化參數(shù)和作物遺傳參數(shù)同時(shí)進(jìn)行優(yōu)化，并對模型參數(shù)進(jìn)行綜合敏感度分析，旨在為土壤-作物-大氣系統(tǒng)模型參數(shù)的自動(dòng)優(yōu)化提供一種便捷方法，為WHCNS模型的快速推廣應(yīng)用提供技術(shù)支持。

1　材料與方法

1.1田間試驗(yàn)和數(shù)據(jù)測定

于2009年10月10日－2011年9月28日在山東省泰安市山東農(nóng)業(yè)大學(xué)農(nóng)學(xué)實(shí)驗(yàn)站（36°09′37″N、117°09′18″E，海拔130 m）進(jìn)行田間試驗(yàn)。試驗(yàn)地冬小麥-夏玉米一年兩作，本研究選取了3塊代表不同土地生產(chǎn)力水平的農(nóng)田，各田塊面積大小為30 m×10 m，土壤基本理化性質(zhì)如表1?；谛←湺嗄昶骄a(chǎn)量將這3個(gè)田塊分為高T1、中T2、低產(chǎn)田T3，分別對應(yīng)平均產(chǎn)量＞9 000 kg/hm2、約7 500 kg/hm2和約6 000 kg/hm2）。

冬小麥（泰農(nóng)18）在10月上旬播種，6月中旬收獲，基本苗為240萬株/hm2，耕作方式為深耕，秸稈還田。夏玉米（鄭單958）在6月下旬播種，9 月下旬收獲；密度為6.6萬株/hm2，耕作方式為免耕。3塊農(nóng)田水肥及耕作管理完全一致。具體的田間水肥管理見文獻(xiàn)[18]。

土壤0～160 cm水分?jǐn)?shù)據(jù)每10 d利用中子儀（CNC503B，北京核子儀器有限公司）觀測1次，測定時(shí)每20 cm 1層。土壤無機(jī)氮每月測定1次，用土鉆取0～160 cm鮮土樣，每20 cm 1層，用流動(dòng)分析儀（AA3，Bran and Luebbe）測定。土壤飽和導(dǎo)水率采用定力水頭方法測定[21]，土壤水分特征曲線采用砂箱和壓力薄膜法測定[21]，采用van Genuchten模型[22]擬合得到水力學(xué)參數(shù)，具體參數(shù)見表1。作物葉面積每月測定1次，采用葉片長×寬×0.75計(jì)算得到；地上部干物質(zhì)也是每月測定1次，將作物樣在75℃下烘干稱質(zhì)量獲得；作物成熟后劃定3 m×3 m區(qū)域，測定該區(qū)域內(nèi)產(chǎn)量,具體方法見文獻(xiàn)[23]。氣象數(shù)據(jù)（日最高、最低和平均氣溫、降雨量、日照時(shí)數(shù)、風(fēng)速和相對濕度）由當(dāng)?shù)貧庀笳咎峁?/p>

表1　農(nóng)田土壤剖面基本理化性質(zhì)Table 1　Basic soil physical and chemical properties of soil profile for three fields

1.2土壤水熱碳氮模型

研究采用本研究小組研發(fā)的WHCNS模型，該模型可模擬土壤水分運(yùn)動(dòng)及碳氮循環(huán)的關(guān)鍵過程，包括地表徑流、土壤溫度、剖面水分動(dòng)態(tài)、土面蒸發(fā)和作物蒸騰、有機(jī)質(zhì)周轉(zhuǎn)、氮素遷移與轉(zhuǎn)化各項(xiàng)、作物生長等過程，可以模擬農(nóng)田水肥管理、作物輪作、耕作等措施對土壤水氮運(yùn)移和作物產(chǎn)量的影響，該模型原理及詳細(xì)過程見參考文獻(xiàn)[18]。

1.3PEST軟件優(yōu)化算法及參數(shù)敏感度

1.3.1參數(shù)優(yōu)化算法流程

對于非線性系統(tǒng)，PEST軟件基于Gauss-Marquardt-Levevberg （GML）算法對目標(biāo)函數(shù)求解最小值，該算法是介于牛頓法與梯度下降法之間的一種非線性優(yōu)化方法，同時(shí)具備高斯牛頓法的快速收斂性和梯度下降法的全局搜索性，能夠在多維的參數(shù)空間內(nèi)優(yōu)化模型輸入?yún)?shù)，迭代逐步逼近目標(biāo)函數(shù)最小值，具體算法參見文獻(xiàn)[8]。

PEST算法的核心是對目標(biāo)函數(shù)求解最小值，目標(biāo)函數(shù)值一般是由實(shí)驗(yàn)觀察值與模擬值的差值組成，根據(jù)觀測值的種類和模型模擬的需求，可以將觀測值分組并附以相應(yīng)權(quán)重，以期望達(dá)到實(shí)測值與模擬值最好的匹配，研究包括土壤含水率、土壤硝態(tài)氮含量、作物產(chǎn)量、作物葉面積指數(shù)4組目標(biāo)函數(shù)，其中冬小麥和夏玉米同時(shí)優(yōu)化等同對待，目標(biāo)函數(shù)不區(qū)分C3或C4作物，其形式如下：

式中Φ是目標(biāo)函數(shù)；下標(biāo)sw、sn、y和l分別對應(yīng)土壤含水率、土壤硝態(tài)氮含量、作物產(chǎn)量、作物葉面積指數(shù)；p分別是各對應(yīng)目標(biāo)的觀測值個(gè)數(shù)；Oi和Si分別代表第i個(gè)實(shí)測值和第i個(gè)模擬值；ω分別是各對應(yīng)目標(biāo)的權(quán)重系數(shù)，其取值方法見2.1小節(jié)。

1.3.2參數(shù)靈敏度分析

在PEST迭代優(yōu)化過程中，計(jì)算雅克比矩陣J占據(jù)了其大多數(shù)的運(yùn)行時(shí)間，每一次迭代至少調(diào)用模型N次（優(yōu)化參數(shù)個(gè)數(shù)），研究設(shè)置最大迭代次數(shù)為10次。PEST在優(yōu)化時(shí)會估計(jì)各個(gè)參數(shù)的相對綜合敏感度（RCS，relative composite sensitivity），RCS反映的是輸入?yún)?shù)對目標(biāo)函數(shù)值的敏感程度，其值越高代表參數(shù)對目標(biāo)函數(shù)的響應(yīng)程度越高，公式如下：

式中si為第i個(gè)參數(shù)的綜合敏感度；m是觀測值的加權(quán)平均值。

1.3.3PEST軟件的參數(shù)優(yōu)化步驟

首先，選擇WHCNS模型需要優(yōu)化的參數(shù)，以及各個(gè)參數(shù)取值范圍、初始值；然后，建立PEST軟件能夠識別的模型參數(shù)輸入文件，觀測數(shù)據(jù)文件和控制文件，以實(shí)現(xiàn)PEST軟件與WHCNS模型之間的數(shù)據(jù)交換；其次，在DOS窗口下，執(zhí)行控制文件，計(jì)算機(jī)根據(jù)控制文件中使用者給定的路徑地址自動(dòng)調(diào)用WHCNS模型，PEST在優(yōu)化運(yùn)行前會計(jì)算出各個(gè)組分的殘差，通過反復(fù)調(diào)試權(quán)重系數(shù)使得各組分殘差相等；最后，設(shè)置最大迭代次數(shù)，開始優(yōu)化。

1.4模擬效果評價(jià)

選擇3種統(tǒng)計(jì)參數(shù)：均方根誤差（root mean square error，RMSE）、模型模擬效率（Nash-Sutcliffe modeling efficiency，E）和一致性指數(shù)（agreement index，d）來評價(jià)模型的模擬值與實(shí)測值的吻合程度[18]。

2　結(jié)果與分析

2.1參數(shù)靈敏度分析

RCS反映的是輸入?yún)?shù)對目標(biāo)函數(shù)值的敏感程度，其值越高代表參數(shù)對目標(biāo)函數(shù)的響應(yīng)程度越高，而目標(biāo)函數(shù)由土壤含水率、土壤硝態(tài)氮含量、作物產(chǎn)量和作物葉面積指數(shù)四部分構(gòu)成，這4組觀測值在數(shù)量級上差異很大，PEST在優(yōu)化運(yùn)行前會計(jì)算出各個(gè)組分的殘差，通過反復(fù)調(diào)試權(quán)重系數(shù)可以使得4個(gè)組分殘差相等。經(jīng)過調(diào)試，土壤含水率、土壤硝態(tài)氮含量、作物產(chǎn)量和作物葉面積指數(shù)的權(quán)重系數(shù)分別設(shè)置為40、0.3、5和0.01時(shí)系統(tǒng)給出的各組分殘差相同。圖1給出了18個(gè)參數(shù)的RCS值，整體上看，相比氮素轉(zhuǎn)化參數(shù)，土壤水力學(xué)參數(shù)和作物遺傳參數(shù)普遍具有較高的敏感度。土壤水力學(xué)參數(shù)中飽和含水率有最高的敏感度程度，其敏感程度順序?yàn)棣萻＞n＞α＞θr＞Ks；作物遺傳參數(shù)方面，冬小麥和夏玉米的生長發(fā)育積溫Tsum和最大比葉面積SLAmax都具有較高的敏感度，作物各個(gè)階段的作物系數(shù)的敏感度其次，對于模型的率定可以優(yōu)先考慮敏感度較高的參數(shù)。PEST計(jì)算的氮素轉(zhuǎn)化參數(shù)普遍敏感度較低，其中與硝化作用有關(guān)的參數(shù)較敏感，最大硝化速率常數(shù)Vn*敏感度最高，其次是硝化半飽和常數(shù)Kn，一些研究結(jié)果也表明土壤氮素轉(zhuǎn)化參數(shù)的敏感度較低[15,17]。

圖1　WHCNS模型各輸入?yún)?shù)相對綜合敏感度Fig.1　Relative composite sensitivity of each input parameter of WHCNS model

2.2基于PEST的參數(shù)優(yōu)化

參數(shù)自動(dòng)優(yōu)化的初始土壤水力學(xué)參數(shù)設(shè)置為實(shí)測值，而作物遺傳參數(shù)和碳氮循環(huán)參數(shù)來源于模型默認(rèn)值，土壤基本物理參數(shù)采用實(shí)測數(shù)據(jù)（土壤飽和導(dǎo)水率、飽和含水率、殘余含水率、水分特征曲線擬合參數(shù)a和n），共選取18個(gè)水力學(xué)參數(shù)和作物遺傳參數(shù)作為自動(dòng)優(yōu)化對象，其中包括5個(gè)土壤水力學(xué)參數(shù)、5個(gè)氮素轉(zhuǎn)化參數(shù)和8個(gè)作物遺傳參數(shù)，具體的參數(shù)選擇、各個(gè)參數(shù)取值范圍及優(yōu)化結(jié)果見表2。

表2　用于PEST優(yōu)化的水力學(xué)參數(shù)和作物遺傳參數(shù)初值和取值范圍Table 2　Initial value and range of soil hydraulic parameters and crop genetic parameters use for PEST

2.3模型校準(zhǔn)與驗(yàn)證

2.3.1模型校準(zhǔn)

本研究采用T2處理進(jìn)行模型參數(shù)率定，采用其他2個(gè)處理校驗(yàn)?zāi)Ｐ?。模型?jì)算時(shí)段為2009年10月10日－2011年9月28日，期間連續(xù)運(yùn)行，以實(shí)測土壤水力學(xué)參數(shù)作為模型初始值，包括土壤飽和導(dǎo)水率、飽和含水率、殘余含水率、水分特征曲線擬合參數(shù)a和n，然后使用PEST優(yōu)化軟件自動(dòng)率定表2中的參數(shù)，在取值范圍內(nèi)自動(dòng)尋優(yōu)，使其最大程度地吻合T2處理實(shí)測的土壤含水率、土壤硝態(tài)氮含量、葉面積指數(shù)和作物產(chǎn)量。參數(shù)率定后輸入T1和T3處理進(jìn)行模型校驗(yàn)。最后與實(shí)測土壤水力學(xué)參數(shù)的模擬結(jié)果進(jìn)行對比。從表3可以看出，PEST軟件對T2處理的率定結(jié)果較好，特別是對土壤含水率和作物產(chǎn)量的率定效果有較大改進(jìn)，含水率RMSE從0.085 cm3/cm3下降到0.032 cm3/cm3，E值從?0.33提高到0.80，d值增加了0.16；作物產(chǎn)量RMSE降低73%，E值提高到0.93，d值從0.81增加到0.97；土壤硝態(tài)氮和作物葉面積指數(shù)模擬效果也有相應(yīng)提高。同時(shí)將PEST方法與試錯(cuò)法的結(jié)果也進(jìn)行了比較（表3），在土壤水分和作物產(chǎn)量模擬方面，PEST方法比試錯(cuò)法有更高的精度，這與一些學(xué)者的研究結(jié)果是一致的[15,17]；而土壤氮素和葉面積指數(shù)的模擬效果并沒有顯著的提升?？傮w上來看，PEST方法在參數(shù)率定中不僅大大節(jié)約了調(diào)試時(shí)間，而且取得了高于試錯(cuò)法的模擬精度，提高了模型工作效率。

注：試錯(cuò)法結(jié)果來源于文獻(xiàn)[18]；RMSE是均方根誤差，土壤含水率單位為cm·cm，硝態(tài)氮單位為mg·kg，產(chǎn)量單位為kg·hm，E是模型模擬效率，d是一致性指數(shù)。Note: Results of trial-and-error method come from reference[18]; RMSE is root mean square error, cm3·cm-3for soil moisture, mg·kg-1for nitrate concentration, kg·hm-2for yield, E is Nash-Sutcliffe modeling efficiency, d is agreement index.

2.3.2土壤含水率驗(yàn)證

表3中給出了T1和T3處理土壤含水率驗(yàn)證的統(tǒng)計(jì)指標(biāo)（RMSE、E和d值），其RMSE介于0.032～0.037 cm3/cm3之間，E值均大于0.8，d值均大于0.9，說明PEST方法能夠得到較高的模擬精度。從3個(gè)統(tǒng)計(jì)指標(biāo)可以看出，與實(shí)測值參數(shù)法模擬結(jié)果相比，PEST優(yōu)化算法模擬的土壤含水率有更高的精度。以T1處理土壤各層含水率（10、30、70和110 cm）PEST方法和實(shí)測參數(shù)法模擬值為例（圖2），土壤含水率在70 cm以上土層變化比較劇烈，隨著深度的增加而變化比較平緩，主要是由于上層土壤受到降雨、灌溉和大氣條件的影響比較大，與實(shí)際情況相符。深層土壤含水率變幅很小，說明已接近田間持水量。2種方法對表層土壤含水率的模擬吻合程度都較好，而PEST法對底層含水率（70和110 cm）的模擬精度有明顯的提升。

圖2　T1處理土壤剖面體積含水率模擬值與實(shí)測值比較Fig.2　Comparion of simulated and measured soil water content in different soil depths for T1

2.3.3土壤硝態(tài)氮含量驗(yàn)證

表3中給出了T1和T3處理土壤硝態(tài)氮含量驗(yàn)證的RMSE、E和d值，其RMSE介于8.6～8.9 mg/kg之間，E值范圍為?0.15～0.04，d值范圍為0.68～0.79，由于氮素轉(zhuǎn)化過程十分復(fù)雜，加上硝態(tài)氮具有較大的空間變異性，許多研究認(rèn)為負(fù)的E值是可以接受的[24]。從3個(gè)統(tǒng)計(jì)值的變化上可以看出，與實(shí)測參數(shù)法相比，PEST法模擬的土壤硝態(tài)氮含量精度明顯較高。圖3顯示了PEST法與實(shí)測參數(shù)法模擬的T1處理土壤各層硝態(tài)氮含量（10、30、70、和110 cm）與實(shí)測值的對比情況，可以看出表層土壤硝態(tài)氮變化劇烈，而底層相對平緩，模擬值與實(shí)測值的變化趨勢基本一致，反映了硝態(tài)氮在土體中的實(shí)際變化情況，其中PEST優(yōu)化算法總體上提高了土壤硝態(tài)氮的模擬精度（表3），特別是對于土壤表層硝態(tài)氮的模擬效果明顯好于實(shí)測參數(shù)法。

圖3　T1處理土壤剖面硝態(tài)氮模擬值與實(shí)測值比較Fig.3　Comparion of simulated and measured soil nitrate concentration at observed depths for T1

2.3.4作物生長驗(yàn)證

表3中給出了T1和T3處理作物產(chǎn)量驗(yàn)證的RMSE、 E和d值，其RMSE介于140～227 kg/hm2之間，E值均大于0.7，d值均大于0.9，2種方法都對冬小麥和夏玉米的產(chǎn)量有很高的模擬精度。從3個(gè)統(tǒng)計(jì)指標(biāo)可以看出，與實(shí)測值參數(shù)法相比，PEST優(yōu)化算法模擬的作物產(chǎn)量有更高的精度，其中T1和T3處理的RMSE分別降低64% 和72%。圖4給出3個(gè)處理冬小麥和夏玉米葉面積指數(shù)的模擬與實(shí)測值的對比結(jié)果。

圖4　3種處理葉面積指數(shù)擬值與實(shí)測值比較Fig.4　Comparison of simulated and measured crop leaf area index for 3 treatments

由于PEST對每個(gè)葉面積指數(shù)觀測值具有相同的權(quán)重，故PEST優(yōu)化算法對于實(shí)測值越多的時(shí)間段其優(yōu)化的精度也越高[25]。由于第1季田間管理較差，加上取樣點(diǎn)較少，故第1季小麥的葉面積指數(shù)模擬效果相對較差，而其余作物生長季的模擬效果均較好?？傮w來看，兩種方法對玉米葉面積指數(shù)的模擬吻合程度都較高，特別是對于T2和T3第2年玉米的模擬。與實(shí)測值參數(shù)法相比，PEST法明顯提高了小麥葉面積指數(shù)的模擬精度，特別是對第2年3個(gè)處理葉面積指數(shù)的模擬有較大的改善，T1、T2和T3處理的RMSE分別降低了6%、20%和15%，E 和d值都有所提高，其中T3處理E值從0.65提高到0.74。

3　討　論

從模型參數(shù)敏感度分析來看，土壤水力學(xué)參數(shù)和作物遺傳參數(shù)普遍具有較高的敏感度，這與孫美等[17]對RZWQM模型參數(shù)的靈敏度分析結(jié)果相似。在冬小麥和夏玉米生長過程中，上層土壤（通常0～50 cm）的根系密度比較大與土壤交互作用也最強(qiáng)，導(dǎo)致其水力學(xué)參數(shù)也有較高的敏感度，這與Nolan等[15]的研究結(jié)果相似。對于模型參數(shù)的率定可以優(yōu)先考慮敏感度較高的參數(shù)，作物遺傳參數(shù)方面，冬小麥和夏玉米的生長發(fā)育積溫Tsum和最大比葉面積SLAmax都具有較高的敏感度，作物各個(gè)階段的作物系數(shù)的敏感度其次；2種作物最大同化率AMAX參數(shù)的敏感程度差異最大，該參數(shù)代表了作物同化CO2的能力，而C3作物（如小麥）和C4作物（如玉米）有不同的光合反應(yīng)途徑，在模型中使用不同的光合作用公式[18]，這可能是導(dǎo)致其敏感度不同的原因。此外，冬小麥和夏玉米的最大根深Rmax敏感程度也有較大差異，中國華北平原冬小麥生育期內(nèi)土壤較為干旱[21]，而根系深度參數(shù)控制著作物吸水的最大深度，這是導(dǎo)致冬小麥Rmax參數(shù)較為敏感的原因。PEST計(jì)算的氮素轉(zhuǎn)化參數(shù)普遍敏感度較低，其中與硝化作用有關(guān)的參數(shù)較敏感，最大硝化速率常數(shù)Vn*敏感度最高，其次是硝化半飽和常數(shù)Kn，一些研究結(jié)果也表明土壤氮素轉(zhuǎn)化參數(shù)的敏感度較低[15,17]。通常對模型參數(shù)進(jìn)行敏感性分析可以減小需要優(yōu)化的參數(shù)，本研究對所有的模型參數(shù)均進(jìn)行了相對綜合敏感度分析，發(fā)現(xiàn)只有18個(gè)模型參數(shù)對模型輸出有影響，故對該18個(gè)模型參數(shù)均進(jìn)行了優(yōu)化。從參數(shù)優(yōu)化的結(jié)果來看（表2），仍然存在一些不太敏感的參數(shù)（如飽和導(dǎo)水率、殘余含水率、硝化半飽和常數(shù)、根系最大深度等），它們的優(yōu)化結(jié)果基本沒有變化。因此，如何設(shè)定相對綜合敏感度的閾值，剔除一些不敏感性的參數(shù)，利用敏感參數(shù)進(jìn)行分層調(diào)參的結(jié)果對比，從而提高優(yōu)化效率，這需要進(jìn)一步的研究。

在模型的驗(yàn)證過程中發(fā)現(xiàn)氮素模擬統(tǒng)計(jì)值E值出現(xiàn)了負(fù)值，由于氮素轉(zhuǎn)化過程十分復(fù)雜，加上硝態(tài)氮具有較大的空間變異性，許多研究認(rèn)為負(fù)的E值是可以接受的[24]。從統(tǒng)計(jì)值上看PEST的優(yōu)化可以達(dá)到較高的模擬精度，但同時(shí)也存在一定的局限，特別對于作物葉面積的模擬。王禮恒等[25]利用PEST軟件優(yōu)化地下水?dāng)?shù)值模型時(shí)發(fā)現(xiàn)觀測資料越豐富，優(yōu)化結(jié)果與真實(shí)值越接近，反之優(yōu)化結(jié)果誤差較大，本研究中葉面積指數(shù)觀測值較少，特別是第1年的冬小麥生育期內(nèi)觀測值只有3個(gè)，使得PEST優(yōu)化對第1季小麥的葉面積模擬產(chǎn)生了較大的誤差。

最大迭代次數(shù)的設(shè)置對模型參數(shù)的優(yōu)化精度和效率有比較大的影響。圖5顯示PEST優(yōu)化次數(shù)與殘差平方和（sum of squared weighted residuals，SSWR）之間的關(guān)系，可以看出SSWR值隨迭代次數(shù)的增加而減少，當(dāng)?shù)螖?shù)達(dá)到7次以后，SSWR值不再減小，因此，本研究認(rèn)為設(shè)置最大迭代次數(shù)為10次可以取得較為合理的結(jié)果，這也與前人的研究結(jié)果相似[17]。一些研究表明，PEST軟件改進(jìn)了GML算法的空間尋值方法，優(yōu)化初值的選取會影響優(yōu)化時(shí)的迭代次數(shù)和模型調(diào)用次數(shù)，但對優(yōu)化結(jié)果并無影響[25]。由于模型參數(shù)種類的復(fù)雜性，本文并沒有對不同初值情況下PEST的優(yōu)化效果進(jìn)行分析，這需要進(jìn)一步研究。

圖5　PEST參數(shù)優(yōu)化收斂曲線Fig.5　Convergence curve of PEST parameter optimization

4　結(jié)　論

本文使用PEST參數(shù)自動(dòng)優(yōu)化工具對土壤-作物系統(tǒng)水熱碳氮模型（WHCNS）的土壤水力學(xué)參數(shù)、氮素轉(zhuǎn)化參數(shù)和作物遺傳參數(shù)進(jìn)行了自動(dòng)尋優(yōu)和靈敏度分析。

參數(shù)敏感度分析結(jié)果表明，作物參數(shù)的敏感度最高（其中生長發(fā)育積溫和最大比葉面積），其次是土壤水力學(xué)參數(shù)（其中飽和含水率最大），而氮素轉(zhuǎn)化參數(shù)最低。作物參數(shù)中，作物生長發(fā)育總積溫和最大比葉面積有較高的敏感度，土壤水力學(xué)參數(shù)中飽和含水率的敏感度最高。

模型統(tǒng)計(jì)性檢驗(yàn)指標(biāo)表明，PEST優(yōu)化算法較實(shí)測參數(shù)法對土壤含水率、葉面積指數(shù)和作物產(chǎn)量的模擬精度有較大的提高，但對土壤硝態(tài)氮的模擬精度提高不明顯，其中土壤含水率、土壤硝態(tài)氮含量、作物產(chǎn)量和葉面積指數(shù)的RMSE分別降低了61.8%、23.5%、73.6%和23.3%，同時(shí)PEST法較試錯(cuò)法對土壤水分和作物產(chǎn)量的模擬精度也有較大程度提高，然而對土壤硝態(tài)氮含量和葉面積指數(shù)的模擬精度提高不明顯，但是該方法大大節(jié)約了模型校驗(yàn)時(shí)間，在較短的時(shí)間內(nèi)就取得了高于“試錯(cuò)法”的模擬精度。因此，PEST軟件在WHCNS模型參數(shù)自動(dòng)優(yōu)化中是一種值得推廣的工具。

[參考文獻(xiàn)]

[1] Ahuja L R, Rojas K W, Hanson J D, et al. The Root Zone Water Quality Model[M]. Highlands Ranch, Ohios: Water Resources Publications LLC, 2000.

[2] Stella T, Frasso N, Negrini G, et al. Model simplification and development via reuse, sensitivity analysis and composition: A case study in crop modeling[J]. Environmental Modelling & Software, 2014(59): 44－58.

[3] Rashad R, Sandeep K, Luo Y Q, et al. An algorithmic calibration approach to identify globally optimal parameters for constraining the DayCent model[J]. Ecological Modelling, 2015(297): 196－200.

[4] 李穎，王康，周祖昊. 基于SWAT模型的東北水稻灌區(qū)水文及面源污染過程模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2014，30(7)：42－53. Li Ying, Wang Kang, Zhou Zuhao. Simulation of drainage and agricultural non-point source pollutions transport processes in paddy irrigation district in North-East China using SWAT[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(7): 42－53. (in Chinese with English abstract)

[5] 郝遠(yuǎn)遠(yuǎn)，徐旭，任東陽，等. 河套灌區(qū)土壤水鹽和作物生長的HYDRUS-EPIC 模型分布式模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2015，31(11)：110－116. Hao Yuanyuan, Xu Xu, Ren Dongyang, et al. Distributed modeling of soil water-salt dynamics and crop yields based on HYDRUS-EPIC model in Hetao Irrigation District[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(11): 110－116. (in Chinese with English abstract)

[6] ?im?nek J, ?ejna M, Saito H, et al. The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media[M]. V4.0. California: University of California Riverside, 2008.

[7] Poeter E P, Hill M C. UCODE, a computer code for universal inverse modeling[J]. Computer Geosciences, 1999, 25: 457－462.

[8] Doherty J. PEST Model-Independent Parameter Estimation User Manual [M]. 5th Edition. Brisbane, Australia: Watermark Numerical Computing, 2004.

[9] Li Q, Chen X, Luo Y, et al. A new parallel framework of distributed SWAT calibration[J]. Journal of Arid Land, 2015, 7(1): 122－131.

[10] 莊嘉祥，姜海燕，劉蕾蕾，等. 基于個(gè)體優(yōu)勢遺傳算法的水稻生育期模型參數(shù)優(yōu)化[J]. 中國農(nóng)業(yè)科學(xué)，2013，46(11)：2220－2231. Zhuang Jiaxiang, Jiang Haiyan, Liu Leilei, et al. Parameters optimization of rice development stages model based on individual advantages genetic algorithm[J]. Scientia Agricultura Sinica, 2013, 46(11): 2220－2231. (in Chinese with English abstract)

[11] Elci A, Karadas D, Fistikoglu. The combined use of MODFLOW and precipitation-runoff modeling to simulate groundwater flow in a diffuse-pollution prone watershed[J]. Water Science & Technology, 2010, 62(1): 180－188.

[12] 董艷輝，李國敏，郭永海，等. 應(yīng)用并行PEST算法優(yōu)化地下水模型參數(shù)[J]. 工程地質(zhì)學(xué)報(bào)，2010，18(1)：140－145. Dong Yanhui, Li Guoming, Guo Yonghai, et al. Optimization of model parameters for groundwater flow using parallelized PEST method[J]. Journal of Engineering Geology, 2010, 18(1): 140－145. (in Chinese with English abstract)

[13] Kim K, Whelan G, Purucker S T, et al. Rainfall–runoff model parameter estimation and uncertainty evaluation on small plots[J]. Hydrological Processes, 2014, 28(20): 5220－5235.

[14] Fang Q X, Green T R, Ma L, et al. Optimizing soil hydraulic parameters in RZWQM2 under fallow conditions[J]. Soil Science Society of America Journal, 2010, 74(6): 1897－1913.

[15] Nolan B T, Puckett L J, Ma L W, et al. Predicting unsaturated zone nitrogen mass balances in agricultural settings of the United States[J]. Journal of Environmental Quality, 2010(39): 1051－1065.

[16] 房全孝. 根系水質(zhì)模型中土壤與作物參數(shù)優(yōu)化及其不確定性評價(jià)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2012，28(10)：118－123. Fang Quanxiao. Optimizing and uncertainty evaluation of soil and crop parameters in root zone water quality model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(10): 118－123. (in Chinese with English abstract)

[17] 孫美，張曉琳，馮紹元，等. 基于PEST軟件RZWQM2模型參數(shù)優(yōu)化與驗(yàn)證[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2014，45(11)：146－153. Sun Mei, Zhang Xiaolin, Feng Shaoyuan et al. Parameter optimization and validation for RZWQM2 model using PEST method[J]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(11): 146－153. (in Chinese with English abstract)

[18] 梁浩，胡克林，李保國，等. 土壤-作物-大氣系統(tǒng)水熱碳氮過程耦合模型構(gòu)建[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2014，30(24)：54－66. Liang Hao, Hu Kelin, Li Baoguo, et al. Coupled simulation of soil water-heat-carbon-nitrogen process and crop growth at soil-plant-atmosphere continuum system and its validation[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(24): 54－66. (in Chinese with English abstract)

[19] Li ZhouJing, Hu Kelin, Li Baoguo, et al. Evaluation of water and nitrogen use efficiencies in a double cropping system under different integrated management practices based on a model approach[J]. Agricultural Water Management, 2015(159): 19－34.

[20] 梁浩，胡克林，李保國. 不同施肥管理模式下農(nóng)田氮素淋失及水氮利用效率模擬分析[J]. 農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào)，2015，34(7)：1317－1325. Liang Hao, Hu Kelin, Li Baoguo. Simulation and analysis of nitrate leaching and water and nitrogen use efficiencies in farmland under organic and mineral fertilizer management practices[J]. Journal of Agro-Environment Science, 2015, 34(7): 1317－1325. (in Chinese with English abstract)

[21] 劉光崧. 土壤理化分析與剖面分析描述[M]. 北京：中國標(biāo)準(zhǔn)出版社，1996.

[22] Van Genuchten M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5): 892－898.

[23] 劉海濤，胡克林，李保國，等. 土壤剖面基礎(chǔ)性質(zhì)差異對農(nóng)田水氮過程和產(chǎn)量的影響[J]. 中國農(nóng)業(yè)科學(xué)，2015，48(7)：1348－1360. Liu Haitao, Hu Kelin, Li Baoguo, et al. Effects of soil profile basic properties on water and nitrogen movement and crop yield[J]. Scientia Agricultura Sinica, 2015, 48(7): 1348－1360. (in Chinese with English abstract)

[24] Kersebaum K C, Hecker J M, Mirschel W, et al. Modelling water and nutrient dynamics in soil-crop systems[M]. Netherlands: Springer, 2007.

[25] 王禮恒，董艷輝，李國敏，等. 基于PEST的地下水?dāng)?shù)值模擬參數(shù)優(yōu)化的應(yīng)用[J]. 工程勘察，2014(3)：38－42. Wang Liheng, Dong Yanhui, Li Guomin, et al, Application of groundwater numerical simulation for parameter optimization based on PEST[J]. Geotechnical Investigation & Surveying, 2014(3): 38－42. (in Chinese with English abstract)

Parameter optimization and sensitivity analysis of soil-crop system model using PEST

Liang Hao, Hu Kelin※, Li Baoguo
(College of Resources and Environmental Sciences, China Agricultural University, Key Laboratory of Arable Land Conservation in North China, Ministry of Agriculture, Beijing 100193, China)

Abstract:Agricultural production management system models usually require numerous input parameters and the calibration and validation of the parameters are time-consuming, which significantly limit the use of models. This study aimed at improving the efficiency and accuracy of a soil-crop system model (soil water heat carbon and nitrogen simulator, WHCNS) using a model-independent optimization tools (parameter estimation, PEST) and data from field experiments. A two-year field experiment was conducted from October 2009 to October 2011 in Tai’an City, Shandong Province in North China Plain. The crop rotation was winter wheat-summer maize, and three fields with high, middle and low productivity levels based on the wheat yields (named T1, T2 and T3 treatments, respectively) were selected to test the WHCNS model. The dynamics of soil water content and soil nitrate concentration in different soil depths were monitored, crop dry matter and leaf area index at the key crop growth stages and yield data were measured. PEST was used to optimize model parameters and to calculate the relative composite sensitivity (RCS) of each input parameter for WHCNS model. The optimization parameters involved the majority modules of the model, such as soil water dynamic, nitrogen transformation and crop growth. The objective function of the optimization model were consist of four different groups of field data, including soil water content, soil nitrate concentration, crop yield and leaf area index. And the inverse solution was obtained through minimizing the object function using PEST program base on Gauss-Marquardt-Levevberg algorithm. The results of PEST were then compared with the simulations based on measured soil hydraulic parameters and the trial-and-error method. The statistical analysis (root mean square error, model efficiency, and agreement index) indicated that the PEST optimization method provided better accuracy and efficiency than the other two methods. For example, PEST method significantly decreased RMSE of soil water content, nitrate concentration, crop yield and leaf area index by 61.8%, 23.5%, 73.6% and 23.3%, respectively. Furthermore, the accuracy of simulated water content, nitrate concentration and crop yields were significantly improved by using PEST method. However, there were no significant improvements for the soil nitrogen concentrations and leaf area index, compared to the trial-and-error method. With sensitivity analysis, we identified 18 key parameters that had relatively higher sensitivity. Among these 18 parameters, soil water hydraulic parameters and crop genetic parameters had higher sensitivity than soil nitrogen transformation parameters. Among soil water hydraulic parameters, the soil saturated water content had the highest sensitivity; among crop parameters, the total cumulative available temperature and maximum specific leaf area showed the highest sensitivity; and among soil nitrogen transformation parameters, the maximum soil nitrification rate showed the highest sensitivity. Overall, the sensitivity of nitrogen transformation parameters was generally lower compared with those of soil hydraulic parameters and crop parameters. The sensitivity of crop parameters was significantly different between wheat (C3 crop) and maize (C4 crop), e.g., the maximum root depth and the maximum assimilation rate for maize showed a higher sensitivity than those of wheat, suggesting that model calibration and validation should be crop specific. The PEST method not only greatly saved time for model calibration, but also achieved significant higher simulation accuracy than that by trial-and-error method. In conclusion, the PEST parameter optimization program is a useful tool and should be adopted in calibration and application of soil-crop models.

Keywords:sensitivity analysis; crops; models; PEST; parameter optimization; nitrogen cycle; WHCNS

通信作者：※胡克林，男，湖北鐘祥人，博士，教授，主要從事土壤空間變異及溶質(zhì)運(yùn)移的研究。北京中國農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院，100193。

作者簡介：梁浩，男，湖北仙桃人，博士生，主要從事資源環(huán)境系統(tǒng)模型的研究。北京中國農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院，100193。

基金項(xiàng)目：國家自然科學(xué)基金項(xiàng)目（41171184、51139006）；長江學(xué)者和創(chuàng)新團(tuán)隊(duì)發(fā)展計(jì)劃項(xiàng)目（IRT0412)

收稿日期：2015-09-08

修訂日期：2015-12-10

中圖分類號：S152

文獻(xiàn)標(biāo)志碼：A

文章編號：1002-6819(2016)-03-0078-08

doi：10.11975/j.issn.1002-6819.2016.03.012