基于Morris和Sobol法的茶園地上部生物量與氮磷流失模擬敏感性分析

葉志成,張莉莉,廖凱華,朱青,賴曉明,郭長強

(1.中國科學院南京地理與湖泊研究所流域地理學重點實驗室,南京 210008;2.中國科學院大學南京學院,南京 211135)

基于過程的生物地球化學模型作為模擬一系列生態現象和生態過程的有效工具,對于研究土壤-作物系統中的物質運移和元素轉化等過程有重要意義。相比于有限的點位觀測和大田試驗,生物地球化學模型具有可重復、耗時短的優點,可以有效節約人力和物力成本。但是為了模擬生態系統中的復雜過程,模型往往被構建為一個復雜的動態系統,其中包含大量的輸入參數,而這些參數均具有不同程度的不確定性,這種不確定性會嚴重影響模型輸出的準確性。因此,通過敏感性和不確定性分析識別模型輸出不確定性的來源及影響程度,被認為是應用模型的先決條件[1]。通過敏感性分析可以識別不同模型參數對于模型輸出的影響程度,對于影響程度較小的模型參數采用經驗值,著重研究分析影響程度大的模型參數,從而提高模型率定的效率與精度[2]。不確定性分析結果顯示出對模型輸出結果誤差的估計。這兩種分析方法往往結合使用,以提升模型結果的可靠性?；谶^程的DNDC模型由于其對土壤碳氮遷移過程的準確追蹤和模型的易操作性,在中國各地都得到了廣泛應用,其敏感性也得到了廣泛檢驗。但是前人的研究領域大多集中在農田[3]、果園[4]和濕地[5]等,很少有研究關注DNDC模型在東南濕潤區坡地茶園模擬領域的效果,并且前人的研究通常采用單一敏感性分析方法進行敏感性分析,很少結合多種方法。當選用的敏感性分析方法不同時,預測的參數敏感性也會發生變化。在中國,茶樹作為一種重要的經濟作物被廣泛種植,茶園為了追求經濟效益,往往施用大量化肥,導致過量的活性氮磷殘留在土壤中,造成土壤酸化和氮磷的大量淋失。同時,前人的研究也發現茶園土壤是亞熱帶地區乃至全球N2O 排放的重要來源,排放因子遠遠高于其他農田[6]。近年來,天目湖流域內大量林地被開發為茶園,茶園面積在2017年已經擴大到流域面積的10%[7],由于茶葉種植能帶來巨大的經濟效益,這種增長可能會一直持續。且流域內茶園多為坡面種植,結合流域土壤土層較薄的特點,在降雨事件發生時存在極高的面源氮磷流失風險。實地調查表明,由于茶園的迅速擴張,已經導致天目湖水質急劇惡化[8]。因此,深入研究茶園土壤的氮磷遷移、轉化過程以及土壤與作物之間的耦合關系具有重要的現實意義。本研究利用定性(Morris法)和定量(Sobol法)兩種敏感性分析方法和基于蒙特卡洛的不確定性分析方法,基于流域內徑流小區實測數據,確定DNDC模型的輸出結果(茶園的地上部生物量、氮磷淋失量、水分滲漏量以及N2O,NO 和NH3排放量)對土壤、地形和作物參數的全局敏感性,并對其進行排序確定每個模型輸出的最敏感參數,探討敏感參數出現的原因,并量化每個輸出結果的不確定性,為實現茶園土壤氮磷運移轉化過程的準確模擬提供基礎。

1 材料與方法

1.1 站點概況及參數選擇

本研究以中國科學院南京地理與湖泊研究所“天目湖流域生態觀測研究站”(北緯31°16′8″,東經119°24′32″)為例,基于站點觀測到的氣象數據、土壤數據和農業管理數據進行分析。該地區位于天目湖流域,氣候為亞熱帶季風性氣候,多年平均溫度17.3℃。多年平均降水量為1 165.3 mm,大約65%的降水集中在6—9 月份。試驗區位于低山丘陵區(高程＜200 m),坡面為直型坡,平均坡度9°,接近天目湖流域茶園平均種植坡度,坡面土壤類型為黃棕壤,滲透能力強,難以形成地表徑流,大量水肥以壤中流的形式輸移[9]。坡上等高種植茶樹,修剪、除草等管理措施均為當地常規管理方式,據調查流域內茶園一年施肥兩次,合計折氮量為400～600 kg/hm2,折磷量為100～150 kg/hm2。茶園坡地小區的坡頂、坡中和坡底分別設有土壤水樣監測點,在土層底部放入土壤滲漏液采集器,從2019年12月開始以大約每月一次的頻率采集土壤水分滲漏液,測定其溶解性總氮濃度(DTN)、溶解性總磷濃度(DTP)、硝態氮濃度(NO-3-N)和正磷酸鹽濃度(PO3-4-P)等數據。每年的3月和11月各施肥一次,肥料累計折氮量和折磷量分別為306.75,106.08 kg/hm2,并在3月追肥后對茶樹進行采摘和修剪,測定所采茶葉的水浸出物、茶多酚和游離氨基酸等指標和干茶重。

模型模擬輸入的氣象數據來自觀測站實地測量數據,選取2018—2020年3年的逐日最高氣溫、最低氣溫和日降雨量輸入模型(圖1),模型最后的輸出值取3 a的平均值。所需的土壤和地形參數主要來源于田間試驗的實測值,使用土鉆采集坡地小區不同深度土壤樣品,風干后測定其田間持水量、質地、容重等。作物參數選取DNDC模型中的默認值。各參數的取值范圍為標準值±10%,并假定所有參數的分布為均勻分布(表1)。

表1 選取的DNDC模型參數及取值范圍Table 1 Selected parameters and their value ranges in DNDC model

圖1 2018-2020年天目湖流域氣象因子時間變化Fig.1 Temporal variation of meteorological factors in Tianmu Lake Basin from 2018 to 2020

1.2 DNDC模型

DNDC(DeNitrification-DeComposition)模型是一個基于過程的生物地球化學模型,其核心是模擬農業生態系統中碳氮等元素的運移轉化等過程[13],在預測農作物產量、農田土壤碳動態、化肥對水體污染及農田溫室氣體排放等方面得到了廣泛應用。模型由兩部分組成,第一部分包括土壤氣候、作物生長和有機物分解3個子模型,通過氣候、土壤、作物和農業管理措施4個驅動因素來模擬土壤環境條件(溫度、濕度、p H 值、氧化還原電位等)、作物生長(作物對氮磷和水的吸收量和作物根、莖、葉和籽粒的生長量)和土壤碳氮庫的動態變化〔CO2的產生、NH+4-N,NO-3-N 和土壤有機碳(SOC)濃度剖面〕;第二部分由硝化、反硝化和發酵子模型組成,有機物分解產生的碳氮磷等元素被輸入這3個子模型中,進而模擬出氮磷的淋失量以及N2O,NO,N2,NH3和CH4的生產、消耗和排放。

1.3 敏感性分析方法

敏感性分析方法主要分為局部敏感性分析和全局敏感性分析兩種,局部敏感性分析通過逐一改變單個參數的值來分析單個參數變化對模型輸出變量的影響,適用于線性或接近線性的模型,而全局敏感性分析則充分考慮了多個參數變化及參數之間的交互作用對模型輸出變量的影響,適用于耦合過程較多的非線性模型。全局敏感性也可以分為定性和定量兩種,定性分析主要評價各參數對模型結果影響的相對大小,方法包括Morris參數篩選法[14]和多元回歸法等;而定量方法則可以給出每個參數對模型結果的貢獻率,包括基于方差的Sobol法[15]和擴展傅里葉幅度檢驗法。

1.3.1 Morris法 Morris法是一種是在局部敏感性分析方法一次一個變量法(one-at-a-time method)基礎上改良而成的一種全局敏感性分析方法,其原理是通過一次只改變一個參數的取值方法輪流計算參數的基本效應值(EEi),計算公式如下:

式中:EEi為第i個單位的基本效應值;f(X1,X2,…,Xi+Δ,…,Xk)為參數變化后模型的輸出;f(X1,X2,…,Xk)為參數變化前模型的輸出;Δ為參數的變化量。

Morris敏感性指數的計算公式為:

式中:平均值μi越大代表參數Xi對模型輸出的影響越大,即敏感性越強;標準差σi則代表參數之間交互作用的強弱,σi值越高,交互作用越強[14]。

1.3.2 Sobol法 相對于定性的Morris篩選法,基于方差的Sobol法可以定量求解復雜非線性模型中各個參數及參數間相互作用所產生的敏感性,其主要思想是將模型函數分解成多項遞增項之和,通過隨機采樣計算參數對模型響應的總方差(V)及各項偏方差(V1,Vij,…,V1,2,…,n),求得各參數的靈敏度。

Sobol敏感性指數的計算公式為:

一階敏感性指數

二階敏感性指數

總敏感性指數

式中:一階敏感性指數Si代表參數Xi的變化對于整個模型方差變化的貢獻,用于定量描述參數Xi對模型輸出的影響;二階敏感性指數Sij代表參數Xi和Xj之間的相互作用對于模型輸出方差產生的影響;以此類推,總敏感性參數STi為各階敏感性指數之和,不僅反映了參數單獨變化產生的影響,也考慮了參數之間的相互作用對輸出結果的影響,更詳細的計算過程請參考Sobol[15]的論文。

1.4 分析方案與數據處理

本研究基于Morris法和Sobol法對茶園小區的茶樹地上部生物量、氮磷淋失量、水分滲漏量及N2O,NO和NH3排放量進行敏感性和不確定性分析,選擇8個土壤參數、1個地形參數和11個作物參數輸入模型。敏感性和不確定性分析主要借助基于蒙特卡洛的敏感性分析軟件Simlab完成,本文的敏感性分析方案借鑒Campolongo等[16]提出的一種兼具準確性和計算成本的方法,具體方案如下:(1)對所有輸入參數采用Morris法進行分析,模型運行次數為r×(k+1)次,其中r取4,k取19,共運算80次。(2)計算敏感性均值μi和標準差σi,對計算出的μi進行分類,取的參數為Morris敏感參數。(3)對Morris法得到的敏感參數進行Sobol法分析,模型運行次數為N×(2k+2),其中N取32,k取11,共運算768次,得出一階敏感指數Si和總敏感性指數STi,將STi＞0.10作為Sobol敏感指數的取值標準。(4)對各個輸出結果的μ值與STi值進行皮爾遜相關性分析?；赟obol法產生的768次模擬,得出各輸出結果的概率統計分布并量化其不確定性。本研究所用的各項數據采用Microsoft Excel,SPSS和Origin作數據處理與作圖分析。

2 結果與分析

2.1 Morris法結果

Morris法的分析結果如圖1 所示,位置越靠右的模型參數敏感度越高,位置越靠上的模型參數與其他參數的交互作用越強。對于地上部生物量,土壤參數的敏感度極低,而作物參數P1(最大籽粒產量)和P2(成熟時籽粒占總生物量之比)的μ值遠遠大于其他參數,與其他參數交互作用較強。在模型模擬過程中,作物參數P1和P2決定了在最佳的生長條件下籽粒的最大生物量以及成熟時籽粒占全株生物量的比例,是影響地上部生物量的主要參數(圖2A)。

圖2 基于Morris法的DNDC模型敏感性分析結果Fig.2 Sensitivity analysis based on Morris method for DNDC model

DNDC模型中氮淋失量主要由土壤中的硝態氮濃度、黏粒含量和水分滲漏量3種因素共同作用,并且在模擬中黏粒含量對于氮淋失量存在指數級影響,因此試驗結果顯示S5(黏粒含量)是最主要的影響因素(圖2B)。作物參數P1,P2會影響植物的日需氮量,從而改變土層中的氮濃度。P6(葉碳氮比)和S6(孔隙度)也通過影響植物的需水量和土壤水分活動間接對氮素淋失產生影響,交互作用較強;對于磷淋失量,土壤參數S1(容重)的μ值遠高于其他參數,是影響磷淋失量的主要參數。在模型中,淋失過程只考慮溶解態磷,土壤容重通過影響土壤的吸附解吸以及有機態磷和溶解態磷之間的轉換過程影響磷淋失量(圖2C);S8(坡度)是影響水分滲漏量的最直接因素(圖2D),DNDC通過每日作物的需氮量和需水量來計算當日的蒸騰量,因此P1,P2,P11(需水系數)和S5也間接影響模擬的水分滲漏量,其中P11與其他參數的交互作用較強。

由于N2O和NO都是硝化反應和反硝化反應的中間產物,產生量的比例取決于反應的進行程度,因此二者的敏感參數存在相似性,其中S1,S7(SOC),P1,P2,S3(田間持水量)均為N2O和NO排放量的敏感參數,但是影響程度不同(圖2E,2F)。這是由于不同的環境因子對硝化反硝化過程的制約程度不同,其中影響N2O排放量的主要參數為S1,影響NO 排放量的主要參數為S2(p H)和P2。同時,S2也是影響NH3排放量的主要參數,這是由于土壤中的酸堿性對于各土層的NH+4-N濃度影響較大,進而影響NH3的產生。

2.2 Sobol法分析結果

相比于Morris分析結果,Sobol法可以更直觀地觀察出參數的敏感性大小及其相互作用強弱(圖3)。參數的總敏感性指數與一階敏感性指數之差即為該參數與其他參數的交互作用對模型輸出項的影響程度。對于地上部生物量(圖3A),敏感參數與Morris法分析結果相似,P1和P2仍是影響地上部生物量的主要參數,總敏感性指數分別為0.27,0.29,但是與Morris法不同的是,Sobol法的結果顯示S5對地上部生物量也具有很強的影響,總敏感性指數為0.24,并且與其他參數的相互作用程度較高;對于氮、磷淋失量、水分滲漏量和NH3排放量,Sobol法分析結果也與Morris法分析結果相似,影響氮淋失量的主要參數是S5,影響磷淋失量的主要參數是S1,影響水分滲漏量的主要參數是S8,影響NH3排放量的主要參數是S2,總敏感性指數分別為0.63,0.46,0.28,0.42。

圖3 基于Sobol法的DNDC模型敏感性分析結果Fig.3 Sensitivity analysis based on Sobol method for DNDC model

但對于N2O 排放量和NO 排放量,結果略有不同。Sobol法分析結果顯示參數S1,S3,P1,P2和S7對N2O 排放量的總敏感性指數較為接近,介于0.11～0.13,其中S7的總敏感性指數略大,為0.13(圖3E);影響NO 排放量的主要參數是S2,參數S1,S3,S5,P1,P2的總敏感性指數介于0.09～0.12,對于輸出結果的影響程度也較高(圖3F)。

2.3 Morris法與Sobol法敏感性分析結果比較

對Morris法和Sobol法所得出來的敏感參數,分別按照其敏感性指數(μ和STi)的大小從左到右進行排序并對兩種方法得出的參數敏感度順序做相關性分析,得出的結果見表2。除了地上部生物量和N2O 排放量,Morris法與Sobol法得出的最敏感參數均一致。對于地上部生物量、氮淋失量、磷淋失量、水分滲漏量和NH3排放量,參數敏感度順序相關性極高,平均為0.90,但N2O 排放量和NO 排放量的參數敏感度順序相關性略低,平均為0.76,總平均相關性為0.86。

表2 Morris法與Sobol法敏感參數排序和相關性比較Table 2 Morris method and Sobol method for sensitive parameter sequencing and correlation comparison

2.4 不確定性分析結果

上述結果表明,隨著模型輸出的不同,各參數的敏感性也不相同。為了進一步分析DNDC模型模擬的茶園地上生物量和氮磷流失等輸出結果,有必要進一步識別DNDC 模型輸出結果的誤差,量化各輸出結果的不確定性?；赟obol法建立的768次模擬結果,對其進行統計分析,結果見表3。其中,磷淋失量的變異系數值最大,為43.84%,說明模擬結果的離散程度較大,誤差較大。水分滲漏量的變異系數值最小,為2.36%,說明模擬的準確度較高,誤差較小。

表3 模型各輸出統計結果Table 3 DNDC model outputs statistical results

3 討論

在DNDC模型中,作物參數主要控制生物量的增長和分配以及影響作物與土壤的物質交換。P1是決定作物日潛在生長量的關鍵參數,對于作物的生物量大小起到最主要的控制作用。P2代表作物成熟時籽粒的生物量占總生物量之比,影響著生物質在籽粒、葉、莖和根之間的分配。在模型模擬的過程中,當模擬的作物還在生長期時,每日增加的生物量只在葉、莖和根中根據比例分配,而當作物的生長指數(累計積溫占總積溫的比例)超過設定的生長期比例時,新增的生物量全部分配到籽粒中[16]。因此,P2對于地上生物量大小具有很強的影響作用。本研究所采用的兩種敏感性分析方法都顯示P1和P2是影響作物地上生物量的主要參數。茶樹作為一種葉芽收獲作物,在采摘后莖葉留存率較高(本文設置為0.9),因此在收獲后會進行人工修剪(本文設置修剪比例為0.4),產生的大量凋落物繼續輸入土壤,并且作為多年生作物,其根系與土壤的物質交換活動頻繁,對土壤p H 值、有機碳含量等影響顯著。在DNDC 模型中,作物的生物量和碳氮比決定了其養分吸收量和蒸騰量,作物的根系呼吸以及產生的凋落物也會影響土壤中的有機物分解與硝化反硝化過程[17]。本研究的結果也表明作物參數P1,P2和P6對模型輸出氮淋失量、水分滲漏量、N2O 和NO 排放量均具較高的敏感性。

相比于作物參數,土壤參數對土壤中的各種過程具有更直接的影響。前人的研究表明土壤的黏粒含量與土壤的持水保肥能力息息相關[18],黏粒含量高的土壤養分保持能力較好,可以在一定程度上減少氮的淋失量,而黏粒含量較低的砂質土壤中的氮素較易流失。另外,DNDC模型主要采用土層硝態氮濃度、水分滲漏量及黏粒含量等來模擬預測氮淋失量[19],因此本文的分析結果也顯示土壤的黏粒含量是影響氮淋失量的最主要參數。DNDC 模型中的磷素在土壤中的平衡主要依賴于土壤的吸附-解吸過程,研究顯示在p H 值為4～10的土壤中,吸附過程是控制溶解態磷的主要作用[20],土壤對磷的吸附能力受到土壤膠體、土壤團聚體和黏粒含量等因素的綜合影響,而這些因素與土壤容重密切相關,因此成為影響磷淋失量的最主要參數。正是由于模擬磷淋失過程中涉及到的參數眾多,模擬的磷淋失量的不確定性也較大。坡度在模型中被用來調節降雨時的地表徑流量,同時也對地下滲漏量產生顯著影響,本文的結果也顯示S8是對水分滲漏量影響最大的參數。DNDC模型中的硝化與反硝化子模型模擬N2O,NO 和NH3三種氣體的輸出,作為土壤氮循環的重要過程,許多研究對影響硝化與反硝化過程的因素進行了深入研究。朱兆良等[21]研究發現在酸性到微堿性的p H 值范圍內,硝化作用和p H 值呈顯著正相關,而Li等[22]研究結果顯示當土壤p H 值大于8.5或小于6.5時,土壤的硝化活性受到抑制。與農田土壤相比,茶園土壤的p H 值較低,這是由于茶樹為喜鋁作物,其根際效應和富鋁凋落物的輸入,會使土壤中的鋁元素富集,導致茶園土壤偏酸。本文的結果顯示,土壤p H 值是影響NO 排放量的主要參數,而對N2O 排放量影響不明顯,這是由于模型在對反硝化速率的計算過程中,當p H 值大于5.25時,p H 值對NO 反硝化速率的影響系數會顯著加大,而對于N2O 反硝化速率,p H 值需大于6.25[23]。本文設置的p H 值范圍為5.769～7.051,在此范圍內p H 值對于NO 反硝化速率的影響程度大于對N2O反硝化速率。NH3的揮發來源于土壤中的NH+4-N 的轉化,p H 值對于土壤中NH+4-N 濃度影響顯著,因此是影響NH3排放量的主要參數。SOC含量同樣對土壤的硝化與反硝化過程存在影響,一方面土壤中的有機物可以為參與硝化或反硝化過程的微生物提供能量并且可以通過礦化分解為反應提供底物,從而促進硝化或反硝化過程,另一方面由于有機物周圍作為土壤呼吸的“熱點”,有利于形成厭氧環境,進一步促進反硝化作用[24]。研究表明,由于茶園的高施肥量,土壤中的SOC 含量也高于同一地區的稻田和林地[25-26],SOC的增加刺激了土壤中的微生物活動,促進N2O 和NO 的形成。Bollmann等[27]研究表明,由于黏粒含量高的土壤保水能力較強,氣體擴散率較低,有利于形成厭氧環境,因此黏粒含量對于土壤的反硝化勢有著顯著的正效應。這與本文的研究結果也一致,盡管順序不同,上述影響硝化反硝化過程的因素均為N2O,NO 和NH3排放量的敏感參數。同樣,這些參數自身的不確定性和模擬過程的不確定性相互疊加,導致模擬值的不確定性偏大。

以不同國家(法國、希臘、意大利、葡萄牙和西班牙)的5 個水稻區為例,Confalonieri等[28]應用Morris與Sobol兩種方法對WARM 模型模擬出的水稻產量進行敏感性分析,研究發現模型輸出對參數的敏感性在不同地點或氣候中會產生變化,并且兩種敏感性分析方法的敏感參數也不一致,因此當模型在新環境中使用時,重新檢測參數敏感性是必不可少的。但是對于同一研究區的敏感性分析結果,許多研究發現不同的敏感性方法往往表現出相似的排序方案。DeJonge等[29]利用Morris與Sobol兩種方法對CERES-Maize模型進行敏感性分析發現Morris方法的μ和Sobol方法的STi之間的r值高達0.93,表現出相當高的相關性。Confalonieri等[30]基于意大利北部稻田對WARM模型應用多種敏感性分析方法,包括Morris法、3種基于回歸的方法和兩種基于方差的方法,結果發現這幾種方法均能提供相似的敏感性排序。本文的研究結果顯示,模型輸出結果的總相關性為0.86,略低于Qin等[31]的結果,但是高于Campolongo 等[16]的結果。其中,地上部生物量、氮淋失量、磷淋失量、水分滲漏量和NH3排放量,敏感性參數的相關性較高,平均值為0.90,但是N2O 和NO 兩組的相關性略低,分別為0.79,0.73。從圖2E—F可以看出,各個參數之間的相互作用程度和非線性程度較高,而Morris法對于識別各參數間的交互作用和非線性效應效果較差,導致兩種方法的相關性較低。

兩種敏感性分析方法各有利弊,Sobol法可以量化各參數對模型輸出的貢獻率,得到的結果更加有可信度,但是運算成本較高,Morris法所需要的模型運行次數較少,但獲取的敏感性分析結果不夠精確。因此在具體應用時,可以結合兩種方法優勢,先應用Morris方法確定敏感參數的范圍,再應用Sobol法精確量化,快速準確地篩選出模型敏感參數。根據本文得到的結果,作物參數對于各個模型輸出均存在不同程度的敏感性,因此在模型率定過程中,可以先考慮調整重點作物參數,再根據敏感性分析得出的敏感參數順序,校準土壤參數。

4 結論

(1)敏感性分析結果顯示,影響地上部生物量最主要的參數為P2(成熟時籽粒占總生物量之比);影響氮淋失量最主要的參數為S5(黏粒含量);影響磷淋失量最主要的參數為S1(容重);影響水分滲漏量最主要的參數為T1(坡度);影響N2O 排放量最主要的參數為S7(初始土壤有機碳含量);影響NO 和NH3排放量最主要的參數為S2(p H 值)。

(2)Morris法和Sobol法得到的參數敏感度順序之間的平均相關性為0.86,相關性較高,兩者具有一定可替代性。Morris法優在方便省時,Sobol法優在準確全面,將兩者結合可以優化敏感性分析過程。另外,不確定性分析結果顯示磷淋失量的不確定性最大,變異系數為43.84%;水分滲漏量的不確定性最小,變異系數為2.36%。

(3)在天目湖流域對DNDC 模型進行參數率定時,建議在識別出敏感參數的基礎上,先重點校正作物參數,然后校正土壤參數,可以有效提高模型的校正效率。針對不確定性較大的模型輸出(如磷淋失量),在未來校正模型時,建議在本文的基礎上細化參數的范圍,進而進一步提高敏感性分析結果的準確性,降低模型輸出的不確定性。未來可以考慮設置不同的氣候條件、土壤類型和施肥情景來分析模型的敏感性和不確定性,以提高分析結果的適用范圍。