基于GLUE的土壤溶質運移參數反演及不確定性*

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

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

（2 中國科學院大學，北京 100049）

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

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

1 材料與方法

1.1 試驗材料

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

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

1.2 靜態批量平衡試驗及阻滯系數計算

制備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+濃度。每個濃度處理設置3個重復，并做無土空白實驗。利用靜態批量平衡實驗結果，采用等溫吸附模型計算土壤對Cu2+的吸附量，并計算阻滯系數。

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

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

1.3 示蹤及溶質運移試驗

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

1.4 控制方程

忽略微生物過程，并假定液相中溶質與土壤相間是平衡交換的。考慮線性等溫吸附的一維穩態流溶質運移可以用CDE來描述：

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

1.5 NLLS方法

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

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

1.6 GLUE方法

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

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

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

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

P95CI表達式如下：

式中，NQin為落在95%置信區間內的觀測點個數，n為總觀測點數。

MNS表達式如下：

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

2 結果與討論

2.1 NLLS參數估計及數值模擬

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

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

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

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

表2 NLLS法對Br-和Cu2+穿透曲線擬合的“最優”參數及擬合效果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方法參數估計及其不確定性

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

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

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

2.4 模型輸出的不確定性

與土柱b類似，土柱a和土柱c示蹤和Cu2+運移試驗反演結果分別如表3和表4所示。可以看出，GLUE確定的v、D、Rd、μ 的置信區間均完全包含NLLS的置信區間。說明在三種土壤中，NLLS的參數置信區間均小于實際可接受的參數范圍。若僅使用NLLS算法的置信區間，將導致大量原應被接受的參數組合被舍棄。

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

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

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

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

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

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

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

3 結 論

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