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

（1.寧夏大學(xué)土木與水利工程學(xué)院，銀川750021；2.旱區(qū)現(xiàn)代農(nóng)業(yè)水資源高效利用教育部工程研究中心，銀川750021）

0 引 言

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

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

1 模型描述

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

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

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

1.2 反演結(jié)果的評價準則

在本文中使用一階矩L1，相關(guān)系數(shù)Cor，決定系數(shù)R2，均方根誤差RMSE來評價SLE 反演結(jié)果優(yōu)劣的標準[13]。L1，RMSE，Cor的數(shù)學(xué)表達式為：

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

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

2 數(shù)值實驗

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

2.1 正演計算

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

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

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

2.2 SLE參數(shù)反演

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

表1 相關(guān)長度、K實測值數(shù)取值Tab.1 Value of correlation scale and K field measured volume

3 結(jié)果與討論

3.1 滲透系數(shù)反演結(jié)果

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

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

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

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

圖4 相關(guān)長度取值1 500 m時K反演結(jié)果Fig.4 K inversion result at correlation scale of 1 500 m

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

3.2 誤差分析

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

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

表2 決定系數(shù)R2增加值Tab.2 Determination coefficient R2 increment

3.3 討 論

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

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

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

4 結(jié) 論

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

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

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