數值反演模型模擬層狀土壤水力性質*

周曉冰，欒永霞，林 青，徐紹輝?

（1. 青島大學環境科學與工程學院，山東青島 266071；2. 青島市城陽區水資源服務中心，山東青島 266109）

土壤水力性質是水分運動方程的重要參數，主要包括水力傳導率K（θ）、土壤水分特征曲線θ（h）等，是建立土壤水分運動數學模型的基礎。參數的可靠性直接影響水分運動模型預測的準確性。盡管有許多實驗室和田間方法可以用來確定包氣帶中以土壤水分特征曲線和水力傳導率為代表的土壤水力性質，但大多數方法需要靜態或穩態流條件來滿足相應分析方法的假設條件或需要較精密的儀器，這使得測量耗時耗力，可行性不高，特別是當土壤的空間變異性很強且區域范圍較大時，更增加了測量的難度。因此，人們開始發展間接方法來估計土壤水力性質。間接方法有土壤轉換函數法、物理-經驗法、分形幾何法、土壤形態學網絡模型和數值反演方法等[1]。其中數值反演方法由于具有以下優點受到研究者的關注[2-4]：（1）無需精密的測量儀器；（2）不需要達到穩態水流，在初始和邊界條件上較直接方法更加靈活；（3）可以從單個瞬態流實驗同時估算土壤水分特征曲線和水力傳導率；（4）為優化參數提供置信區間。但數值反演方法在確定土壤水力性質時通常存在解的非唯一性問題。

關于如何減少數值反演方法解的非唯一性和提高參數估計的準確性，國內外學者做了許多相關研究[5-7]。Toorman等[8]借助響應曲面方法對目標函數進行分析，使用出流量、壓力水頭、含水量的不同組合作為目標函數的測量變量，結果表明，如果在目標函數中包含壓力水頭值，對參數估計的敏感性將一定程度有所提高。此外，壓力水頭的測量位置對參數估計也有影響，可以據此確定需要測定的變量及位置。van Dam等[9]通過出流實驗模擬結果，說明僅采用累計出流量數據容易導致解的非唯一性，而加入θ（h）信息可以得到可靠的參數解。Eching等[10-11]在多步出流實驗的累積出流量數據基礎上添加測定的土壤壓力水頭數據，以三次不同初始參數反演結果的變異系數表征解的唯一性，在反演土壤水分特征曲線和非飽和水力傳導率時得到了較好的擬合效果，驗證了增加壓力水頭信息可以提高解的唯一性。?imunek 和 van Genuchten[12]應用了張力圓盤入滲數據進行數值反演，求解水力特征參數，通過參數與目標函數的響應曲面分析使用哪些入滲數據可以提高解的唯一性，認為可以適當增加剖面的壓力水頭信息和含水量信息來提高計算準確性，而加入累積入滲量信息不能改善擬合結果。Zijlstra和Dane[13]將反演方法應用于層狀土壤，土壤剖面由一個、兩個或三個不同的層位組成，設定水分含量測量值為時間和深度的函數，下邊界設為壓力水頭關于時間的函數，以最小化目標函數法優化參數，結果認為隨著優化參數數量的增加，反演問題可能變得不適定，反演的不適定性取決于輸入數據集的大小。張俊等[14-15]利用數值反演方法耦合瞬時出流實驗推求土壤水力學參數，結果顯示，單步出流實驗在目標函數中增加壓力水頭信息、多步出流實驗聯合累積出流量和壓力水頭信息，均可以減小反演問題的非唯一性。

綜上所述，利用數值反演方法推求土壤水力性質方法雖然有多種，前人也做了關于降低數值反演方法解的非唯一性和提高參數估計準確性的相關研究，但通常只考慮單層土壤，而對層狀土壤反演其土壤水力性質的研究還較少，需要深入開展這方面的工作。本研究結合室內層狀土柱入滲試驗，運用 Hydrus-1D[16]中的反演（inverse solution）模塊求解層狀土壤水力性質，并分析目標函數中包含累積入滲量、土壤壓力水頭及兩者組合對解的唯一性及準確性的影響。以此探究獲得可靠層狀土壤水力性質需要測量的試驗變量。

1 材料與方法

1.1 供試土壤

土樣取自青島市即墨區移風店鎮農田，土壤類型為棕壤。經測定土壤的基本物理性質如表1所示。

表1 供試土壤的理化性質 Table 1 Physico-chemical properties of tested soils

由表1結果可知，采樣點的土體構型相對屬于“上粗下細型”。土壤質地分為兩種：0～50 cm為質地較粗的砂壤，其砂粒含量明顯高于下層，滲透性較好；50 cm以下為粉砂壤，下部土壤黏粒含量較高，滲透性較差。

1.2 試驗方法

試驗采用室內層狀土柱定水頭積水入滲方法，有機玻璃柱填裝土柱，馬氏瓶供水，使馬氏瓶內管底端與液面相平以維持2 cm恒定水頭（圖1）。有機玻璃柱規格：內徑19 cm，高80 cm，一側開六孔安裝張力計，上端設進水口，下端設出水口。馬氏瓶規格：內徑9 cm。土柱分三層按田間土層排列順序填裝，從上至下，土層厚度分別為20 cm、20 cm和25 cm。由于田間原狀土容重較大，室內土柱無法還原田間狀態，故試驗按照設計容重填裝，從上至下容重分別為1.57 g·cm–3、1.45 g·cm–3、1.40 g·cm–3。土柱表層和底層分別鋪設2 cm、5 cm石英砂，由于石英砂質地與表層土相差較小，與底層相差較大，為了減小模擬誤差，上部土層厚度概化為22 cm，下部土層厚度概化為27 cm，土柱總厚度為69 cm。6個張力計從上至下分別放置在9、19、29、39、49和59 cm深度處。土柱底端張力計讀數開始變化時停止通水。試驗全程監測各張力計讀數和馬氏瓶讀數。

1.3 土壤水分運動模型的建立

Hydrus-1D模型是由國際地下水模型中心公布的數值模擬軟件，用來模擬一維飽和-非飽和多孔介質中水、熱、溶質的運移過程。

土壤水分以垂向運動為主，因此模型采用一維垂直入滲Richards方程描述土壤水分運動：

土壤水力函數用van Genuchten-Mualem 模型表示：

式中，C（h）為比水容量（cm–1），Ks為飽和水力傳導率（cm·min–1），K（h）為非飽和水力傳導率（cm·min–1），Se為飽和度，θs為土壤飽和含水量（cm3·cm–3），θr為殘余含水量（cm3·cm–3），α為進氣吸力相關的參數（cm–1），m和n為形狀系數，h為壓力水頭（cm）；z為土壤深度（cm）；t為時間（min）。其中，m=1–1/n，l=0.5。

以上方程包含了5個獨立參數，θr，θs，α，n和Ks，未知參數通過最小化目標函數獲取[15]。參數的優化計算采用Hydrus-1D，version 4.16（US Salinity Laboratory）。Hydrus-1D中目標函數最小化的方法采用LM法，該方法結合了牛頓法和最速下降法，為優化參數提供置信區間，并作為一種標準方法在水力學中求非線性擬合最小平方和而被普遍采用[17-19]。

初始條件：

上邊界：

下邊界：

式中，h0（z）為模擬開始時的土壤壓力水頭值，監測點處為實測值，其余節點采用插值法計算得到，H0為土壤表面積水深度（2 cm）。

時間和空間離散：對試驗條件下的一維水分運動進行數值模擬，模擬時長為T=1 220 min，對應于59 cm深度處水勢傳感器讀數開始變化的時間。計算過程采用可變時間步長?？臻g步長為0.5 cm，依據土柱總厚度69 cm將土壤剖分為138個單元，共139個節點。

初始參數的選?。篐ydrus-1D模型各層土壤水力學參數的選取至關重要，雖然最小化目標函數方法計算復雜度低，但對參數初始值敏感，且得到的優化參數可能是局部最優，而非全局最優，因此需要嘗試不同的初始參數值。對于van Genuchten- Mualem中的六個參數θr，θs，α，n，Ks，l，其中l通常設為0.5。θr采用風干土的質量含水量與體積含水量的關系計算得到，θs通過試驗時土柱內各層土壤容重計算出的孔隙度獲得。為降低反演的復雜程度，反演計算時固定參數θr和θs（表2），主要對3層土壤的α，n，Ks共9個參數進行優化。Carsel 和 Parrish[20]根據概率分布的方法分析大量數據得到了不同土壤類型的van Genuchten非飽和水力參數，選取相應土質的α，n，Ks參數值作為反演的初值，并根據反演效果調整反演時α，n，Ks的初值，以達到更好反演效果，調整后的參數見表3。

表2 未參與反演的部分土壤水力學參數 Table 2 Soil hydraulic parameters not involved in inversion

1.4 評價指標

為了探究層狀土壤一維水分運動中測量變量的類型和數量對數值反演參數唯一性和準確性的影響，分別將累積入滲量（Q）、不同位置壓力水頭（H）及兩者組合（Q+H）作為目標函數中的測量變量，對土壤水力學參數進行優化求解。每次反演過程使用3組不同初始值（表3）進行運算。對目標函數為Q、H和Q+H的9次水分運動模擬過程分別用情景（Ⅰ，Ⅱ，Ⅲ）、（Ⅳ，Ⅴ，Ⅵ）和（Ⅶ，Ⅷ，Ⅸ）表示。采用R2，RMSE指標對擬合效果進行定量化評價。R2也稱為決定系數，代表了回歸方程對實測值的擬合度，值越接近1，擬合效果越好。RMSE為均方根誤差，值越小代表擬合效果越好。對于解的非唯一性用變異系數（CVs）進行評價，并結合優化參數的相關性矩陣以及優化參數的標準誤差對反演效果進行分析。

表3 各情景下的參數初始值 Table 3 Initial values of parameters relative to scenario

2 結果與討論

2.1 不同目標函數的擬合結果

圖2表示目標函數中測量變量為累積入滲量時擬合結果，圖3表示目標函數中測量變量為壓力水頭時擬合結果，圖4和圖5為目標函數中測量變量為累積入滲量與壓力水頭組合時擬合結果。由圖2—圖5可以看出模擬值與實測值符合程度較好，但19 cm處壓力水頭的模擬較差。導致該結果的可能原因如下：（1）19 cm處張力計讀數受氣泡影響讀數上升較慢，未能及時反映出壓力水頭的變化，即由于儀器引起的測量誤差；（2）19 cm處張力計位置靠近中層土壤上界面，土體構型為上粗下細，進 水受空氣阻力影響，水分在19 cm處運移較慢，土壤達到飽和需要較長時間，而模型未考慮進氣阻力致使模型默認水分到19 cm處后很快達到飽和，因此模擬值與測量值吻合度較差。

2.2 不同目標函數的模擬精度

表4為9種情景的反演效果定量化評價，總體擬合效果良好，R2均在0.855～0.999之間，RMSE值在0.005 50 cm～0.038 4 cm之間。三種目標函數下R2大小順序為Q >Q+H>H，RMSE 值大小順序為Q

表4 各情景下的模擬效果評價 Table 4 Evaluation of the simulation in effect relative to scenario

2.3 不同目標函數反演結果的變異系數

由表5可以看出，目標函數僅有累積入滲量時，變異系數總體最小，目標函數僅有壓力水頭時變異系數總體最大。這可能由于累積入滲量實驗誤差較小，而19 cm處壓力水頭模擬值與實測值相差較大，導致目標函數為壓力水頭時變異系數較大。當目標函數為壓力水頭和累積入滲量的組合時，變異系數總體較目標函數只有壓力水頭時小。變異系數越小，解的非唯一性越小。說明在壓力水頭數據基礎上加入累積入滲量數據可以減小解的非唯一性。

表5 不同目標函數反演結果的變異系數 Table 5 Coefficient of variation of the inversion relative to objective function

2.4 不同目標函數反演參數的相關性

參數相關矩陣反映了兩個參數值之間的相關性。值為–1或+1表示完全線性相關，而0表示不相關。參數之間的高度相關性會導致低估參數的不確定性，降低收斂速度并增加非唯一性[21]。根據Hydrus-1D軟件模擬結果反演信息（Inverse solution information）中的相關性矩陣（correlation matrix），分析不同目標函數在9種情景下數值反演的優化參數之間的相關系數，相關系數絕對值大于0.9的數值表示相關性顯著。通過對比不同目標函數下參數間相關性顯著的數量，分析參數的不確定性和非唯一性。

三種目標函數情況下反演參數相關性顯著的總數量，目標函數為Q、H和Q+H時分別為24、6和7個。由此可見，當目標函數中只有Q時，反演參數相關性顯著的總數明顯多于其他兩種情況。說明目標函數只包含累積入滲量時，參數之間相關性較強。此時模型建議固定相關性較強的兩個參數中的一個，以便于對另一個參數更準確求解，否則在此情況下得出的參數優化結果不確定性和非唯一性很強，易出現異參同效現象。這主要是由于層狀土壤水力性質參數較多和水分運動模型的高度非線性造成的。而目標函數為壓力水頭及累計入滲量和壓力水頭組合時，參數相關性顯著的數量基本相等，且大幅減少，這是由于目標函數中含壓力水頭時，提供了不同深度處土壤壓力水頭隨時間變化的信息，限定了更多約束條件，減少了異參同效現象。因此，目標函數中包含壓力水頭時對降低反演過程參數非唯一性有明顯效果。

2.5 不同目標函數反演參數的標準誤差

Hydrus-1D提供了關于擬合參數的統計信息，其中參數的標準誤差是度量擬合結果精確度的指標，由最優值附近的95%置信區間確定，標準誤差越小，參數預測的95%置信區間越窄，擬合結果確定性越高。對三種目標函數下優化參數的標準誤差進行分析，9種情景下優化參數的標準誤差列在表6中，9個參數的標準誤差中，除α1之外，其余參數的標準誤差大小順序均為Q+H

表6 各情景優化參數的標準誤差 Table 6 Standard error of optimized parameters relative to scenario

2.6 驗證反演參數的準確性

為進一步比較目標函數中不同測量變量對數值反演求得的參數準確度的影響，對累積入滲量反演優化的參數結果（表7）進行正演，預測各觀測點壓力水頭值隨時間變化，并與實測值進行對比（圖6），驗證參數的準確性[15]。同樣方法對壓力水頭反演優化的參數結果用累計入滲量預測值與實測值進行對比，得出結果如圖7所示。

表7 各情景下參數的優化結果 Table 7 Optimization of the parameters relative to scenario

由圖6可以看出目標函數只包含累積入滲量時所求參數對壓力水頭的預測值與實測值相差較大，三次模擬R2分別為0.334、0.338、0.332，RMSE分別為0.096 1 cm、0.095 6 cm、0.096 5 cm。由圖7可知目標函數值包含壓力水頭時所求參數對累積入滲量的預測值與實測值符合程度較高，三次模擬R2分別為 0.996、0.994、0.993，RMSE 分別為0.028 76 cm、0.034 02 cm、0.038 56 cm。由此可見運用累積入滲量數據對土壤水力學參數進行數值反演的準確性較差，而目標函數中包含壓力水頭時反演結果準確性較好，基本能反映累計入滲量和壓力水頭隨時間的變化關系。

3 結 論

目標函數僅含有累積入滲量時反演模型模擬值與實測值符合度雖然較高，決定系數R2達到0.999，均方根誤差RMSE小于0.005 62 cm；但用壓力水頭實測值驗證參數準確性的效果較差，R2最高僅為0.338，RMSE小于0.096 5 cm。使用壓力水頭作為目標函數進行參數優化求解，擬合度R2達到0.855以上，RMSE小于0.038 4 cm；對反演得到的參數用累計入滲量實測值進行驗證，R2可達到0.99以上，RMSE小于0.038 56 cm，效果較好。累積入滲量和壓力水頭聯合反演求土壤水力性質擬合度較高，R2達到0.905以上，RMSE小于0.035 4 cm。雖然累積入滲量和壓力水頭聯合反演的R2和RMSE指標次于累積入滲量，但累積入滲量和壓力水頭聯合反演可 以一定程度降低優化參數的相關性和標準誤差，從而降低解的非唯一性，提高解的準確性。綜上所述，由于層狀土壤水力性質的參數較多，數值反演方法求參數的不確定性和非唯一性較強，若單獨使用一種測量數據進行反演計算，很難確定參數的準確性和唯一性。此外，若這一測量變量誤差較大，將嚴重降低結果準確性。因此數值反演過程應盡量使用不同類型和較多數量的測量變量以增加約束條件，提高數值反演解的唯一性、確定性和準確性。