閆一凡 李曉鵬 張佳寶 劉建立?

（1 中國科學(xué)院南京土壤研究所，南京 210008）

（2 中國科學(xué)院大學(xué)，北京 100049）

對流-彌散方程（Convection-dispersion equation，CDE）可用于表征溶質(zhì)的時空動態(tài)特征并揭示溶質(zhì)在多孔介質(zhì)中的遷移轉(zhuǎn)化機理，是最常用的描述土壤中溶質(zhì)運移過程的模型之一［1］。由于溶質(zhì)運移模型的關(guān)鍵參數(shù)（如彌散系數(shù)等）存在較強的非線性或尺度依賴性，通常難以通過實驗直接準確測定［2］。利用監(jiān)測數(shù)據(jù)通過數(shù)值模型反演來識別或優(yōu)化這些參數(shù)，是應(yīng)用CDE解決實際問題的常規(guī)做法。傳統(tǒng)的參數(shù)估計方法主要包括非線性最小二乘法（Nonlinear least squares，NLLS）、高斯-牛頓法（Gauss-Newton）［3］、模擬退火法［4］等。但這些方法容易陷入局部最優(yōu)，無法應(yīng)對“異參同效”并進行參數(shù)和模擬結(jié)果的不確定性分析，在一定程度上制約了模型這一預(yù)測和決策工具的潛力發(fā)揮［5］。

1992年，Beven和Binley［6］提出基于貝葉斯方法的廣義似然不確定性估計（Generalized likelihood uncertainty estimation，GLUE）方法。它是一種基于蒙特卡羅模擬，且可識別多種潛在“異參同效”現(xiàn)象的參數(shù)優(yōu)化方法，同時也是不確定性分析的一種重要手段［5］。不同于NLLS這類傳統(tǒng)尋優(yōu)算法，GLUE方法認為存在一個似然度函數(shù)，符合一定似然度的參數(shù)組合均可接受，且離實測值越近該參數(shù)組合的似然度越大，可信度越高。GLUE方法自提出后，已被廣泛應(yīng)用于水文模型、降雨—徑流模型以及作物生長等各個領(lǐng)域的模型參數(shù)及響應(yīng)界面的不確定性和敏感性分析中［7-11］，在實際應(yīng)用中表現(xiàn)優(yōu)良。其似然函數(shù)、采樣方法及校正標準的選擇及其對結(jié)果的影響也得到了眾多學(xué)者的探索［12-15］。與傳統(tǒng)尋優(yōu)算法相比，GLUE方法的優(yōu)勢在于：1）適用于非線性系統(tǒng)；2）無需假設(shè)誤差的概率分布［12］；3）可綜合反映所有來源的誤差［6］；4）可定量分析不確定性。因此，本文選擇GLUE方法來探討數(shù)值反演CDE模型參數(shù)及模型預(yù)測結(jié)果的不確定性，期望為數(shù)值模擬結(jié)果的風(fēng)險評估提供一定參考。

1 材料與方法

1.1 試驗材料

試驗所用的土壤采自中國科學(xué)院封丘農(nóng)業(yè)生態(tài)國家試驗站（114.51°E～114.60°E，34.98°N～35.06°N），為黃河沖積成因的典型潮土。自表層向下，土壤質(zhì)地類型（美國制）分別為砂壤土（0～20 cm）、壤砂土（20～40 cm）和砂黏壤土（40～60 cm），每種質(zhì)地采集三個重復(fù)，共9個土樣。實驗所用土柱均為高15 cm、直徑9 cm的透明有機玻璃柱。采樣時，在土柱內(nèi)壁涂抹少許凡士林潤滑，以減少采樣阻力、緩減對原狀結(jié)構(gòu)的破壞。在采集原狀土柱同一地點，取適量擾動土樣帶回實驗室，自然風(fēng)干后用激光粒度儀（Mastersizer 3000，英國）測定機械組成，烘干法測定容重，電位法測定pH。土壤樣本理化性質(zhì)如表1。

表1 土壤基本理化性質(zhì)Table 1 Physical and chemical properties of soil sample

1.2 靜態(tài)批量平衡試驗及阻滯系數(shù)計算

制備5、10、15、20、25 mg·L-1的硝酸銅（Cu(NO3)2）溶液備用。在50 mL的離心管中分別加入2.00 g土壤，加入20 mL制備好的不同濃度梯度的硝酸銅溶液，加塞密封。恒溫箱振蕩器保持溫度23℃±1℃，震蕩24 h。然后以4 000 r·min-1的速度離心30 min，過濾分離出上清液，用電感耦合等離子光譜儀（PerkinElmer Optima 8000，美國）測定溶液中的Cu2+濃度。每個濃度處理設(shè)置3個重復(fù)，并做無土空白實驗。利用靜態(tài)批量平衡實驗結(jié)果，采用等溫吸附模型計算土壤對Cu2+的吸附量，并計算阻滯系數(shù)。

采用線性吸附方程，分別對三種質(zhì)地土壤中Cu2+靜態(tài)批量平衡試驗的數(shù)據(jù)進行擬合并計算阻滯系數(shù)，其公式如下：

式中，Rd為阻滯系數(shù)；ρ為土壤干容重，g·cm-3；Kd為分配系數(shù)，L·mg-1；θ為土壤體積含水量，cm3·cm-3。

1.3 示蹤及溶質(zhì)運移試驗

示蹤試驗在三種不同質(zhì)地（a砂壤土，b壤砂土，c砂黏壤土）的土柱中進行。試驗過程保持室溫20℃±1℃。用蠕動泵自下向上緩慢飽和土柱，待形成穩(wěn)定流場時，輸入一個孔隙體積（Pore volume，PV）pH為5.5、濃度為0.05 mol·L-1的溴化鉀（KBr）惰性示蹤溶液，然后連續(xù)洗脫30 h。示蹤試驗完畢后，以同樣的泵入速度向土柱泵入pH為3.5、濃度為0.5 mol·L-1的Cu(NO3)2溶液1 PV，然后連續(xù)洗脫30 h。出流溶液用自動部分收集器采集。出流液中Br-濃度由Br-選擇電極測定，Cu2+濃度由高效液相色譜-等離子體質(zhì)譜儀（Agilent 7700x，澳大利亞）測定。

1.4 控制方程

忽略微生物過程，并假定液相中溶質(zhì)與土壤相間是平衡交換的?？紤]線性等溫吸附的一維穩(wěn)態(tài)流溶質(zhì)運移可以用CDE來描述：

式中，Rd為阻滯系數(shù)；C為液相中溶質(zhì)的濃度，mg·L-1；x為距離，cm；t 為時間，min；μ 為匯項（μ C）的速率系數(shù), min-1，用于描述土壤溶液中基于一階衰減的不可逆化學(xué)持留；v為平均孔隙流速，cm·min-1；D為水動力彌散系數(shù)，cm2·min-1。

1.5 NLLS方法

NLLS算法是以實測值和模型計算值間的誤差平方和最小為準則來識別非線性靜態(tài)模型參數(shù)的一種參數(shù)估計方法。目前應(yīng)用最廣的基于NLLS的土壤溶質(zhì)運移參數(shù)反演程序是由美國鹽土實驗室開發(fā)的CXTFIT。CXTFIT采用Levenberg- Marquardt算法，通過匹配實測溶質(zhì)穿透曲線（Breakthrough curve, BTC）來反演彌散系數(shù)、孔隙水流速等溶質(zhì)運移模型參數(shù)，并預(yù)測溶質(zhì)濃度隨時間和空間的分布規(guī)律。其優(yōu)化目標是尋求一組參數(shù)使實測和計算值之間的決定系數(shù)R2最大化、殘差平方和（SSQ）最小。目標函數(shù)為決定系數(shù)R2，如下式：

式中，Ci和 fi分別為觀測值和模型計算值；C 為全部觀測值的平均值；n為總觀測點數(shù)。R2值越接近1，表明擬合效果越好。

1.6 GLUE方法

GLUE方法既考慮到“最優(yōu)”這一直觀事實，又避免了采用單一“最優(yōu)”參數(shù)組合帶來的預(yù)測風(fēng)險?？紤]到每個參數(shù)可能的空間變異和測量誤差，將平均孔隙水流速度（v）的初始取值范圍定為其測量值ve± 0.5vecm ·min-1，彌散系數(shù)（D）定為0.000 1～0.05 cm2·min-1，阻滯系數(shù)（Rd）定為1.0～3.0，匯項速率系數(shù)（μ）定為 0.001～0.1 min-1。蒙特卡羅采樣策略中采用拉丁超立方采樣方法（Latin hypercube sampling, LHS）生成隨機參數(shù)組合。

在溶質(zhì)運移參數(shù)反演時，為避免數(shù)值求解困難，一般會利用示蹤試驗結(jié)果先確定v和D兩個參數(shù)，然后將其固定再擬合其他參數(shù)。為此，本文中對GLUE方法的應(yīng)用也分為兩個階段，階段一：通過對示蹤試驗數(shù)據(jù)的擬合確定可接受的參數(shù)v和D；階段二：為參數(shù)組（v，D）隨機搭配100組經(jīng)過LHS采樣得到的（Rd，μ）重新組合，生成全新的參數(shù)組合（v，D，Rd，μ）。

為方便與NLLS方法比較，選用Nash–Sutcliffe函數(shù)（形同決定系數(shù)R2）作為似然函數(shù)，定量描述模擬結(jié)果與觀測值間擬合的優(yōu)劣程度。

式中，Lu,t和Ll,t分別為95%置信區(qū)間的出流濃度上下限，mg·L-1；n為總觀測點數(shù)，Ro,t為出流濃度的觀測值，mg·L-1。

P95CI表達式如下：

式中，NQin為落在95%置信區(qū)間內(nèi)的觀測點個數(shù)，n為總觀測點數(shù)。

MNS表達式如下：

式中，i為可接受的采樣編碼，N為可接受的采樣數(shù)量。模型擬合結(jié)果越好，ARIL值越接近于0，P95CI越接近于100% ，MNS越接近于1。

2 結(jié)果與討論

2.1 NLLS參數(shù)估計及數(shù)值模擬

用平衡模型擬合Br-穿透曲線，同時擬合參數(shù)v 和D。由于示蹤劑Br-不發(fā)生吸附解析、降解沉淀等物理化學(xué)反應(yīng)，其阻滯系數(shù)Rd=1，沉淀項速率系數(shù) μ=0 min-1?；贜LLS平衡模型對Br-穿透曲線的擬合結(jié)果對應(yīng)的最優(yōu)參數(shù)及決定系數(shù)見表2。由表2可知，NLLS反演的“最優(yōu)”參數(shù)組合對示蹤離子穿透曲線（BTC）擬合效果極佳，R2均大于0.98，均方根誤差RMSE＜0.046。

由于Cu2+運移試驗條件均保持與示蹤試驗時一致，在擬合Cu2+穿透曲線時，仍選用第一階段識別的v和D值，僅對Rd和μ進行優(yōu)化，結(jié)果如圖1和表2。觀察可知，擬合曲線的峰值略低于觀測值，這可能是因為土柱與管壁之間所形成的微弱優(yōu)勢流所致，也可能是實驗中存在著一定的測量誤差導(dǎo)致的。林青［17］在使用兩點非平衡吸附研究重金屬運移的過程中也出現(xiàn)類似情況。但總體而言，模擬的結(jié)果符合觀測值的趨勢，R2均大于0.93且RMSE也保持在10-2數(shù)量級（表2）。

圖1 Cu2+穿透曲線（BTCs）觀測值與非線性最小二乘法（NLLS）擬合值的對比Fig. 1 Comparison between measured breakthrough curves ( BTCs) of Cu2+and non-linear least square (NLLS) fitted BTCs

2.2 第一階段GLUE方法參數(shù)估計

以土柱b（壤砂土）的Br-穿透曲線擬合為例，來說明GLUE方法的反演結(jié)果。圖2中每一個散點均代表模型運行一次的結(jié)果，所以可視作模型響應(yīng)界面的投影。從指定區(qū)域10 000次蒙特卡羅采樣中，有818個參數(shù)組合的似然值＞ 0.9，被判定為“行為集”。由圖2和表2可知，最大似然值MNS=0.987 2時的參數(shù)（v，D）=（0.031 7 cm·min-1，0.003 9 cm2·min-1），與NLLS得到的唯一“最優(yōu)”解（v，D）=（0.032 1 cm·min-1，0.004 2 cm2·min-1），非常接近。且NLLS方法擬合的R2為0.987 4，與最大似然值（MNS）幾乎相同?？梢灶A(yù)知，蒙特卡羅采樣次數(shù)越多，模型運行的次數(shù)越多，GLUE方法得到的MNS與NLLS尋優(yōu)方法得到的“最優(yōu)解”越接近。

與參數(shù)的初始取值范圍相比，“行為集”的參數(shù)取值范圍明顯縮小。平均孔隙流速v的初始和更新后取值范圍分別為［0.015 6，0.046 8］ cm·min-1、［0.028 6，0.035 4］ cm·min-1；彌散系數(shù)D的初始和更新后取值范圍分別為［0.000 1，0.050 0］cm2·min-1、［0.000 1，0.022 1］ cm2·min-1，區(qū)間寬度分別縮小了55.0%和55.8%。但NLLS方法確定的95%置信限（表2），如圖2中散點圖頂端的線段所示，仍較更新后的GLUE的參數(shù)取值范圍窄很多，說明NLLS方法在Br-運移模型參數(shù)反演時存在“異參同效”現(xiàn)象，且其參數(shù)v和D的置信限遠小于實際可被接受的參數(shù)范圍，導(dǎo)致大量可接受參數(shù)被忽略。

表2 NLLS法對Br-和Cu2+穿透曲線擬合的“最優(yōu)”參數(shù)及擬合效果Table 2 Optimum parameters of NLLS fitting Br- and Cu2+ BTCs and fitting effect

圖2 GLUE方法獲得的似然值 ＞0.9的v和D散點圖Fig. 2 Dotty plots of v and D of 0.90 acquired by GLUE for bromide BTCs

2.3 第二階段GLUE方法參數(shù)估計及其不確定性

第二階段是通過GLUE方法對Cu2+穿透曲線的擬合進行不確定性分析（此處仍以土柱b為例）。為818組0.9的“行為集”分別搭配100組通過拉丁超立方采樣（LHS）得到的隨機參數(shù)Rd和 μ，組成81 800個新的參數(shù)組合。逐次運行模型，得到403個參數(shù)組合滿足＞0.9，其對應(yīng)參數(shù)的散點圖如圖3所示。其他兩種質(zhì)地的土柱也分別得到了448和648組“行為集”。

圖3 GLUE獲得的似然值＞0.9的參數(shù)的散點圖Fig. 3 Dotty plots of the parameters of＞ 0.90 acquired by GLUE fitting Cu breakthrough

GLUE方法受到爭議的原因之一在于選擇似然函數(shù)和閾值的主觀性［19-20］。區(qū)分“行為集”和“非行為集”的閾值選擇的確存在主觀性，但這種選擇是建立在對未來模型預(yù)測的有效性基礎(chǔ)之上的，并非純粹的主觀判斷［21］。似然函數(shù)的選擇亦是如此，本研究中選用的似然函數(shù)與NLLS方法的決定系數(shù)公式形式相同，目的是便于NLLS和GLUE兩種方法的結(jié)果對比。近年來，關(guān)于如何選定合理的似然函數(shù)也得到了更多的關(guān)注。如Zhang等［22］為了強調(diào)模型在實際應(yīng)用中的不確定性引入了多種似然函數(shù)，并基于概率計算來選定最合理的似然函數(shù)。Freni等［15］用多種形式的似然函數(shù)估計了城市洪水模型結(jié)果的不確定性并對模型進行了校正。這些研究表明，Nash–Sutcliffe函數(shù)適用于數(shù)據(jù)量不大和需要掌握單一參數(shù)對模擬結(jié)果影響的情況。當觀測值不斷增多且旨在調(diào)整模型適用性時，指數(shù)類型的似然函數(shù)則更合適。

2.4 模型輸出的不確定性

與土柱b類似，土柱a和土柱c示蹤和Cu2+運移試驗反演結(jié)果分別如表3和表4所示?？梢钥闯?，GLUE確定的v、D、Rd、μ 的置信區(qū)間均完全包含NLLS的置信區(qū)間。說明在三種土壤中，NLLS的參數(shù)置信區(qū)間均小于實際可接受的參數(shù)范圍。若僅使用NLLS算法的置信區(qū)間，將導(dǎo)致大量原應(yīng)被接受的參數(shù)組合被舍棄。

表3 NLLS和GLUE反演得到的Br-運移模型參數(shù)及其95%置信區(qū)間Table 3 Parameter of the Br- transport model and 95% confidence intervals obtained using NLLS and GLUE

表4 NLLS和GLUE反演得到的Cu2+運移模型參數(shù)及其95%置信區(qū)間Table 4 Parameter of the Cu2+ transport model and 95% confidence intervals obtained using NLLS and GLUE

GLUE方法對應(yīng)MNS的參數(shù)組合與NLLS“最優(yōu)”參數(shù)組合的模型預(yù)測結(jié)果如圖4所示。由圖可知，在土柱b中，MNS參數(shù)組合（v，D，Rd，μ）=（0.034 9 cm·min-1，0.005 3 cm2·min-1，1.112，0.003 2 cm-1）與NLLS“最優(yōu)”參數(shù)組合（v，D，Rd，μ ）=（0.032 1 cm·min-1，0.004 2 cm2·min-1，1.034，0.003 1 cm-1），對觀測值的擬合極佳，其R2均為0.937。但MNS與NLLS“最優(yōu)”的參數(shù)組合并不一致，表明不同的參數(shù)組合可以達到類似的模擬結(jié)果，即“異參同效”現(xiàn)象。在土柱a和土柱c中亦存在相同現(xiàn)象。GLUE方法獲取的95%置信區(qū)間較NLLS方法獲取的置信限寬很多，其覆蓋的觀測點比例P95CI的均值為84.30%，而NLLS方法則僅為46.05%，這表明NLLS方法獲取的置信限無法良好預(yù)測Cu2+穿透曲線，特別是峰值區(qū)。GLUE方法的95%置信區(qū)間對試驗觀測點的高覆蓋率還表明了模型定義及模型結(jié)構(gòu)滿足了溶質(zhì)運移模擬的需求［23］。

圖4 NLLS和GLUE兩種方法預(yù)測Cu2+ 穿透曲線的不確定性結(jié)果對比Fig. 4 Comparison between NLLS and GLUE used in predicting Cu2+ BTC in uncertainty

三個定量化指標平均相對區(qū)間長度（ARIL）、MNS和置信區(qū)間內(nèi)的觀測點比例（P95CI）的統(tǒng)計結(jié)果如表5所示。砂壤土、壤砂土和砂黏壤土中，通過NLLS方法獲取的“最優(yōu)解”的決定系數(shù)R2（分別為0.960，0.937，0.954）與GLUE方法獲取的MNS的似然值（分別為0.960，0.937，0.953）高度一致，這在圖4也得到直觀體現(xiàn)。GLUE與NLLS獲得的ARIL均值分別為79.66和135.5，其中壤砂土顯著高于其他兩種質(zhì)地土壤。這主要是由于壤砂土內(nèi)Cu2+穿透曲線的相對濃度在0.000 1以下時，其對應(yīng)的置信區(qū)間和置信限的寬度較其他質(zhì)地更寬所致?？傮w而言，GLUE方法對應(yīng)的ARIL值顯著低于NLLS，這說明MNS對應(yīng)的參數(shù)組合對模型的擬合效果可媲美NLLS的“最優(yōu)”參數(shù)組，且GLUE計算的置信區(qū)間對觀測值的覆蓋度更高，平均相對區(qū)間寬度也更窄，因此，在參數(shù)及模型輸出不確定性分析兩個方面GLUE方法均優(yōu)于NLLS方法。

表5 NLLS和GLUE方法擬合Cu2+ 穿透曲線的不確定性定量化結(jié)果Table 5 Quantification of the uncertainties in fitting Cu2+ BTCs using NLLS and GLUE

NLLS方法計算得到的“最優(yōu)”參數(shù)組合（表2）及其對實測值的擬合結(jié)果（圖1）均表現(xiàn)優(yōu)秀。但NLLS方法的殘差并非白噪聲，這表示該方法的參數(shù)估計可能存在偏差［13］。本文的研究結(jié)果也證實了該結(jié)論。盡管“最優(yōu)”參數(shù)及其置信限的確包含在通過GLUE方法確認的95%置信區(qū)間內(nèi)，但圖2和圖3均顯示，在NLLS法確定的置信限以外仍有大量參數(shù)組合可極好地擬合觀測值。這也證實了NLLS尋優(yōu)的結(jié)果——“最優(yōu)”參數(shù)組合并無看似的高“魯棒性”。再者，NLLS方法對BTC的95%置信限僅覆蓋了低于65%的觀測點（圖3和表5），該結(jié)果再次驗證了以上結(jié)論。由此可知，試驗觀測值所能提供的信息可能并不足以識別真實的參數(shù)值。在多個可接受的參數(shù)向量中，僅通過部分觀測值，通常無法判定哪一個參數(shù)向量才是參數(shù)的真值，對于數(shù)據(jù)稀缺的野外試驗情況更是如此［24］。但依然應(yīng)該關(guān)注和改進傳統(tǒng)的尋找單一最優(yōu)解的算法。其一，對于設(shè)計性實驗，其理論完備且邊界條件等可被嚴格控制，此時傳統(tǒng)尋優(yōu)算法的單一“最優(yōu)解”仍有重要意義。其二，GLUE等不確定性分析方法對計算能力要求較高，若計算所有可行性參數(shù)組合并運行模型驗證超出計算能力，此時就需要單一“最優(yōu)解”來進行預(yù)測。因此，傳統(tǒng)尋優(yōu)算法和GLUE這類承認多解性的算法亦可相輔相成、適當結(jié)合。

3 結(jié) 論

溶質(zhì)運移模型的參數(shù)估計過程和模型預(yù)測本身就蘊含不確定性，無論觀測數(shù)據(jù)量多少、數(shù)據(jù)質(zhì)量高低均無法完全避免此問題。本文利用傳統(tǒng)尋優(yōu)方法（NLLS）和不確定性分析方法（GLUE）對土壤溶質(zhì)運移模型的不確定性進行了研究，結(jié)果表明，NLLS作為一種參數(shù)反演方法簡單明了，易于操作，對Br-及Cu2+穿透曲線的擬合效果較好。但“異參同效”現(xiàn)象的存在說明NLLS法的“最優(yōu)”參數(shù)應(yīng)用于預(yù)測時的“魯棒性”不及預(yù)想。GLUE方法清晰直接地表明有多個參數(shù)組合能滿足擬合要求，其獲取的對應(yīng)于最大似然值MNS的參數(shù)向量亦可較好擬合Br-及Cu2+的穿透曲線。由似然值大于0.9的“行為集”計算出各參數(shù)的后驗取值范圍顯著大于且包括NLLS方法獲取的參數(shù)取值范圍。NLLS方法較窄的參數(shù)置信區(qū)間導(dǎo)致許多可被接受的參數(shù)組合被摒棄，僅使用單一“最優(yōu)”解的尋優(yōu)方法預(yù)測溶質(zhì)運移存在極大的不確定性。