考慮相關長度誤差的含水層滲透系數反演研究

李株丹，劉登川，申同慶，朱 磊

（1.寧夏大學土木與水利工程學院，銀川750021；2.旱區現代農業水資源高效利用教育部工程研究中心，銀川750021）

0 引 言

詳細準確的描述非均質含水層水文地質參數的空間分布特征，是準確預測地下水流動和污染物運移過程的基礎[1]。Yeh提出了一種基于地質統計的參數估計反演方法，稱為連續線性估計算法（Successive linear estimator），即SLE 算法[2,3]。該算法能夠獲取含水層水力參數的空間分布信息，是一種可靠而有效地刻畫含水層非均質性的反演方法。在實際工程中，滲透系數作為水文地質數值模擬的主要參數，準確描述其空間分布特征具有非常重要的意義[4]。

許多學者在提高滲透系數反演精度的研究上做了大量的工作，Carrera 和Van 等[5,6]分析了觀測數據的質量和數量在反演計算時對反演結果的影響，利用有效的數據去做反演計算時，容易得到比較好的反演結果。王文娟等[7]通過研究不同觀測井、抽水井間距、抽水速率大小和采樣時間等對反演精度的影響，提出了抽水試驗的優化設計方案，反演精度明顯提高。Mu 等[8]通過數值實驗分析得到了改變觀測井數和抽水試驗數會對滲透系數反演結果造成的影響，當增加觀測井數和抽水試驗數時，能夠減小滲透系數模擬值與真實值之間的誤差。蔣立群等[9]通過改變抽水試驗次數分析了滲透系數反演精度的提升程度，隨著抽水試驗數的增加，反演精度也在逐步提升，當抽水實驗數增加到一定程度時，即觀測信息過多時，反演精度的提升越來越小。Yu 等[10]研究了在非均質性很強的土壤中收集觀測資料時，觀測數據的獲取較困難。通常，在假定含水層均質各向同性、等厚和無限延伸等條件下，野外現場抽水試驗是能夠有效地獲取觀測信息，但是這些假定與實際情況往往不符。現場抽水試驗也受試驗區條件、經費等各種因素限制，有限的觀測數據很難準確地描述含水層空間分布。因此，滲透系數實測值很難準確而有效地獲得，目前也尚缺乏滲透系數實測值與反演結果誤差的關系。本文將基于情景模擬，通過連續線性估計算法研究在相關長度誤差影響下的含水層滲透系數的反演精度，并結合數值實驗得到滲透系數實測值數據量與反演結果誤差的關系。因此，該研究對實際工程中選取較適宜的滲透系數實測值數據量具有現實和指導意義。

1 模型描述

1.1 連續線性估計方法(SLE)

SLE 算法是將水力參數(即滲透系數)與水頭在估計過程中的非線性關系結合起來，通過使用線性估值來處理，最大化提高數據的有用性，將均值和協方差矩陣輸入模型初始值作為先驗信息，以滿足控制方程的求解條件[11,12]。SLE 算法中的迭代過程可以概化為貝葉斯迭代，同時更新均值和協方差矩陣。在穩定流條件下，用水頭值估計滲透系數可以用以下公式迭代式求解：

式中：上標r表示迭代次數；(xm)表示在xm點第r次迭代的模擬值；hj*(xj)表示在xj點處水頭觀測值；hj(r)(xj)表示在xj點處水頭模擬值；ωmj(r)為權重項，表示第r次迭代時，在xj點處觀測到的條件水頭與模擬水頭之差[即hj*(xj)和hj(r)(xj)]對點xm處估計值的貢獻。

1.2 反演結果的評價準則

在本文中使用一階矩L1，相關系數Cor，決定系數R2，均方根誤差RMSE來評價SLE 反演結果優劣的標準[13]。L1，RMSE，Cor的數學表達式為：

式中：n表示數據的個數；和Xi分別表示第i個觀測井的水頭觀測值和估計值；和分別指觀測值和估計值的均值；和σXi分別表示觀測值與估計值的標準差。

其中，一階矩L1 指絕對誤差；均方根誤差RMSE值大于0，指估計值和觀測值偏差的平方與n比值的平方根；相關系數Cor取值在?1 到1 內，表示觀測值與估計值之間的相關程度；決定系數R2在0到1之間取值，是相關系數的平方。當L1和RMSE數值越小，Cor和R2數值越大時，則觀測值與模擬值之間的誤差越小，此時模擬結果越好。

2 數值實驗

本研究將采用基于SLE 算法的VSAFT2 軟件進行一個綜合的數值算例，該軟件能夠進行2D 地下水水分運動、溶質運移的正演及含水層地質參數的反演計算，計算結果可用Tecplot軟件處理和顯示。

2.1 正演計算

設定水力參考隨機場為一個（10 000×10 000）m2的水平二維飽和非均質穩定流含水層，包含10 000 個（X、Y軸分別為100 個）單元，設定滲透系數的均值為6.5 m/d，方差為4，相關長度X、Y軸均為3 000 m，具有指數協方差函數。含水層的左右邊界均設定為隔水邊界，上、下邊界設定為定水頭邊界，水頭值分別取5 m、4 m。模型共設置一組抽水試驗，每次抽水試驗共設置9口井，其中包括1口抽水井P1，既抽水又觀測，8口觀測井，抽水流量為3 000 m3/d。含水層滲透系數K的參考場見圖1，正演計算水頭值見圖2。

圖1 K的參考場Fig.1 The material field of K

圖2 水頭等值線圖Fig.2 Head contour map

2.2 SLE參數反演

反演計算設定參數除相關長度外均與正演計算時設定數值相同，計算時改變反演相關長度，通過依次增加K實測值數據量計算得到滲透系數估計值的空間分布。相關長度、K實測值數具體取值可見表1。

表1 相關長度、K實測值數取值Tab.1 Value of correlation scale and K field measured volume

3 結果與討論

3.1 滲透系數反演結果

當滲透系數K反演的相關長度取值和正演相關長度相同，均取值為3000m，且K實測值數據量沒有增加時，滲透系數K反演結果如圖3所示。

當減小滲透系數K反演的相關長度，以1 500 m 為例，K實測值數據量分別設置為0、16、32、48、64、80個時，滲透系數K反演結果如圖4所示。

圖3中反演相關長度為3 000 m，圖4反演相關長度為1 500 m，且K實測值均為0 個時，比較可知，當反演相關長度為1 500 m 時，參數K的模擬值越偏離真實值，誤差較大。因此表明，當減小反演相關長度數值時，會對參數K的反演結果有影響，反演結果誤差會隨著相關長度數值的減小而增大。

圖3 相關長度取值3 000 m時K反演結果(R2=0.171)Fig.3 K inversion result at correlation scale of 3 000 m(R2=0.171)

圖4 相關長度取值1 500 m時K反演結果Fig.4 K inversion result at correlation scale of 1 500 m

分析圖4散點圖可知：當K實測值數為0 個時，散點圖中的大部分估計值較分散，偏離45°線較多；當K實測值數為16個時，左右兩側的估計值開始向45°線靠近；當K實測值數為32 個時，邊緣散開的估計值開始趨近于45°線；隨著K實測值數依次增加為48、64、80 個時，由散點圖可知，大部分點越來越趨近于45°線，說明估計值越來越接近于真實值，當K實測值數為80 個時，散點圖圍繞45°線均勻分布，反演結果最好。由此可知：隨著橫坐標K實測值數據量的增加，含水層滲透系數的估計值越來越接近于真實值，即估計值與真實值之間的誤差越來越小，滲透系數反演結果越好，反演精度越來越高。

3.2 誤差分析

隨著反演相關長度數值在正演相關長度基礎上依次減小時，參數K的反演結果誤差與K實測值數據量之間變化的關系如圖5所示。由此分析：隨著反演相關長度依次減小10%、20%、30%、40%直到100%時，以豎直黑實線為基準，當橫坐標K實測值數在30 個以內時，縱坐標L1、RMSE值越來越小，Cor、R2值越來越大；當橫坐標K實測值數大于30 個時，縱坐標L1、RMSE值依然在減小，只是減小的幅度越來越平緩，同理Cor、R2值增加的程度也越來越小。由表2可知，當K實測值數據量由0、16、32、48、64、依次增加至80 個，反演相關長度減小10%時，決定系數R2的值分別為0.161、0.382、0.463、0.522、0.564、0.580，參數K反演精度分別提高了136%、180%、223%、251%、260%；當反演相關長度減小50%時，決定系數R2的值分別為0.071、0.321、0.413、0.472、0.522、0.561，參數K反演精度分別提高了355%、480%、570%、642%、700%；當反演相關長度減小90%時，決定系數R2的值分別為0.052、0.062、0.113、0.174、0.203、0.262，參數K反演精度分別提高了21%、142%、241%、302%、420%。

隨著反演相關長度數值在正演相關長度基礎上依次的增加時，參數K的反演結果誤差與K實測值數據量之間變化的關系如圖6所示。由此分析：隨著反演相關長度依次增加10%、20%、30%、40%直到100%時，以豎直黑實線為基準，當橫坐標K實測值數在30 個以內時，縱坐標L1、RMSE值越來越小，Cor、R2值越趨近于1；當橫坐標K實測值數大于30 個時，縱坐標L1、RMSE、Cor、R2均變化幅度較小，越來越平緩。由表2可知，當K實測值數據量由0、16、32、48、64、依次增加至80 個，反演相關長度增加10%時，決定系數R2的值分別為0.192、0.383、0.474、0.522、0.565、0.583，參數K反演精度分別提高了101%、145%、173%、193%、203%。當反演相關長度增加100%時，決定系數R2的值分別為0.311、0.392、0.461、0.542、0.523、0.541，參數K反演精度分別提高了25%、48%、74%、67%、77%。當反演相關長度增加200%時，決定系數R2的值分別為0.353、0.384、0.494、0.513、0.495、0.523，參數K反演精度分別提高了9%、40%、46%、40%、49%。

表2 決定系數R2增加值Tab.2 Determination coefficient R2 increment

3.3 討 論

分析圖5和圖6可知，無論增加還是減小反演相關長度，當K實測值數據僅有16 個，即K實測值數據量較少，相關長度誤差存在的情況下，反演結果誤差都是比較大的。當橫坐標K實測值數據量開始依次增加時，縱坐標R2、Cor的值開始向1 靠近；L1、RMSE值越來越小。分析可知，以豎直黑實線為基準，隨著橫坐標K實測值數據量的增加，當K實測值數據大于30 個時，各個相關長度誤差下的反演結果誤差值變化幅度都比較小，雖然呈遞增趨勢，但是增加的程度越來越平緩。由此可知，橫坐標K實測值數據量的增加能夠有效地提高含水層滲透系數的反演精度，但是持續的增加K實測值數，直到一定程度時，縱坐標反演結果誤差的變化異常緩慢，并沒有很大的改變反演精度。比如，減小反演相關長度下，以取值2 700 m 為例，相關長度誤差減小10%，當橫坐標K實測值從0增加到32 個時，參數K的反演精度提高了185%；由此對比，當橫坐標K實測值從32 增加到80 個時，參數K的反演精度只提高了26%，并沒有前期精度提高的那么明顯。同理，增大反演相關長度，以取值6 000 m 為例，相關長度誤差增加100%，前期反演精度提高了48%，后期反演精度只提高了17%。因此，通過數值實驗表明，在相關長度誤差存在的基礎上，K實測值的數據量會對反演結果有影響，K實測值的數據量越多，反演結果越好。但是過多的K實測值數據量會造成數據過剩，當K實測值數達到一定程度之后，觀測信息的冗余性使得反演精度提升的越來越慢。因此，對于本文設置的數值實驗而言，K實測值數據量應設置在30個以內。

圖5 減小反演相關長度下K實測值數據量與反演結果誤差的關系Fig.5 The relationship between the measured value of K and the error of the inversion result with reducing the inversion correlation scale

圖6 增加反演相關長度下K實測值數據量與反演結果誤差的關系Fig.6 The relationship between the measured value of K and the error of the inversion result with increasing the inversion correlation scale

4 結 論

(1)在進行反演計算時，若K實測值較難獲得大量準確數據，改變反演相關長度的數值會增加含水層滲透系數的反演結果誤差。

(2)相關長度誤差較大時，改變K實測值數據量會對反演結果造成影響。當依次增加K實測值數據量時，滲透系數反演結果越來越好，從而提高了反演精度，但是增加過多的K實測值到某一程度時，反而會使反演精度提升幅度變得越來越緩慢。

(3)由于過多的K實測值數據會造成觀測信息的冗余性，隨著K實測值數據量的增加，反演精度提升非常緩慢，對于本文設置的數值實驗，K實測值數應小于30 個，所以在實際工程中可以根據水文地質條件選取較適宜的K實測值數據量。