考慮相關(guān)長(zhǎng)度誤差的含水層滲透系數(shù)反演研究

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

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

0 引 言

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

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

1 模型描述

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

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

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

1.2 反演結(jié)果的評(píng)價(jià)準(zhǔn)則

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

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

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

2 數(shù)值實(shí)驗(yàn)

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

2.1 正演計(jì)算

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

圖1 K的參考場(chǎng)Fig.1 The material field of K

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

2.2 SLE參數(shù)反演

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

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

3 結(jié)果與討論

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

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

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

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

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

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

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

3.2 誤差分析

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

隨著反演相關(guān)長(zhǎng)度數(shù)值在正演相關(guān)長(zhǎng)度基礎(chǔ)上依次的增加時(shí)，參數(shù)K的反演結(jié)果誤差與K實(shí)測(cè)值數(shù)據(jù)量之間變化的關(guān)系如圖6所示。由此分析：隨著反演相關(guān)長(zhǎng)度依次增加10%、20%、30%、40%直到100%時(shí)，以豎直黑實(shí)線為基準(zhǔn)，當(dāng)橫坐標(biāo)K實(shí)測(cè)值數(shù)在30 個(gè)以內(nèi)時(shí)，縱坐標(biāo)L1、RMSE值越來(lái)越小，Cor、R2值越趨近于1；當(dāng)橫坐標(biāo)K實(shí)測(cè)值數(shù)大于30 個(gè)時(shí)，縱坐標(biāo)L1、RMSE、Cor、R2均變化幅度較小，越來(lái)越平緩。由表2可知，當(dāng)K實(shí)測(cè)值數(shù)據(jù)量由0、16、32、48、64、依次增加至80 個(gè)，反演相關(guān)長(zhǎng)度增加10%時(shí)，決定系數(shù)R2的值分別為0.192、0.383、0.474、0.522、0.565、0.583，參數(shù)K反演精度分別提高了101%、145%、173%、193%、203%。當(dāng)反演相關(guān)長(zhǎng)度增加100%時(shí)，決定系數(shù)R2的值分別為0.311、0.392、0.461、0.542、0.523、0.541，參數(shù)K反演精度分別提高了25%、48%、74%、67%、77%。當(dāng)反演相關(guān)長(zhǎng)度增加200%時(shí)，決定系數(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可知，無(wú)論增加還是減小反演相關(guān)長(zhǎng)度，當(dāng)K實(shí)測(cè)值數(shù)據(jù)僅有16 個(gè)，即K實(shí)測(cè)值數(shù)據(jù)量較少，相關(guān)長(zhǎng)度誤差存在的情況下，反演結(jié)果誤差都是比較大的。當(dāng)橫坐標(biāo)K實(shí)測(cè)值數(shù)據(jù)量開始依次增加時(shí)，縱坐標(biāo)R2、Cor的值開始向1 靠近；L1、RMSE值越來(lái)越小。分析可知，以豎直黑實(shí)線為基準(zhǔn)，隨著橫坐標(biāo)K實(shí)測(cè)值數(shù)據(jù)量的增加，當(dāng)K實(shí)測(cè)值數(shù)據(jù)大于30 個(gè)時(shí)，各個(gè)相關(guān)長(zhǎng)度誤差下的反演結(jié)果誤差值變化幅度都比較小，雖然呈遞增趨勢(shì)，但是增加的程度越來(lái)越平緩。由此可知，橫坐標(biāo)K實(shí)測(cè)值數(shù)據(jù)量的增加能夠有效地提高含水層滲透系數(shù)的反演精度，但是持續(xù)的增加K實(shí)測(cè)值數(shù)，直到一定程度時(shí)，縱坐標(biāo)反演結(jié)果誤差的變化異常緩慢，并沒有很大的改變反演精度。比如，減小反演相關(guān)長(zhǎng)度下，以取值2 700 m 為例，相關(guān)長(zhǎng)度誤差減小10%，當(dāng)橫坐標(biāo)K實(shí)測(cè)值從0增加到32 個(gè)時(shí)，參數(shù)K的反演精度提高了185%；由此對(duì)比，當(dāng)橫坐標(biāo)K實(shí)測(cè)值從32 增加到80 個(gè)時(shí)，參數(shù)K的反演精度只提高了26%，并沒有前期精度提高的那么明顯。同理，增大反演相關(guān)長(zhǎng)度，以取值6 000 m 為例，相關(guān)長(zhǎng)度誤差增加100%，前期反演精度提高了48%，后期反演精度只提高了17%。因此，通過(guò)數(shù)值實(shí)驗(yàn)表明，在相關(guān)長(zhǎng)度誤差存在的基礎(chǔ)上，K實(shí)測(cè)值的數(shù)據(jù)量會(huì)對(duì)反演結(jié)果有影響，K實(shí)測(cè)值的數(shù)據(jù)量越多，反演結(jié)果越好。但是過(guò)多的K實(shí)測(cè)值數(shù)據(jù)量會(huì)造成數(shù)據(jù)過(guò)剩，當(dāng)K實(shí)測(cè)值數(shù)達(dá)到一定程度之后，觀測(cè)信息的冗余性使得反演精度提升的越來(lái)越慢。因此，對(duì)于本文設(shè)置的數(shù)值實(shí)驗(yàn)而言，K實(shí)測(cè)值數(shù)據(jù)量應(yīng)設(shè)置在30個(gè)以內(nèi)。

圖5 減小反演相關(guān)長(zhǎng)度下K實(shí)測(cè)值數(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)長(zhǎng)度下K實(shí)測(cè)值數(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)在進(jìn)行反演計(jì)算時(shí)，若K實(shí)測(cè)值較難獲得大量準(zhǔn)確數(shù)據(jù)，改變反演相關(guān)長(zhǎng)度的數(shù)值會(huì)增加含水層滲透系數(shù)的反演結(jié)果誤差。

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

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