基于向?qū)c(diǎn)法反演水文地質(zhì)參數(shù)

何金沙，李春光，，呂歲菊，楊佩瑤，黃傳霽

（1.寧夏大學(xué)土木與水利工程學(xué)院，銀川750021；2.北方民族大學(xué)數(shù)值計(jì)算與工程應(yīng)用研究所，銀川750021）

0 引 言

水文地質(zhì)參數(shù)在地下水資源評價(jià)、數(shù)值模擬、開發(fā)利用等方面具有十分重要的作用。地下水?dāng)?shù)值模型所能承載的某一種參數(shù)的最大維度等于該模型的網(wǎng)格剖分?jǐn)?shù)目，大多數(shù)情況下觀測數(shù)據(jù)包含的信息有限，直接對所有網(wǎng)格的參數(shù)進(jìn)行反演具有一定的難度。而參數(shù)化方法通過某種特定數(shù)學(xué)關(guān)系簡化實(shí)際參數(shù)組，無需對數(shù)值模型單個(gè)網(wǎng)格進(jìn)行參數(shù)估計(jì)，僅反演簡化后的參數(shù)組，再對已估參數(shù)點(diǎn)進(jìn)行插值，從而得到含水層參數(shù)的估計(jì)場。傳統(tǒng)的地下水參數(shù)分區(qū)人為的將模型區(qū)域劃分為若干個(gè)小分區(qū)，每個(gè)分區(qū)用一個(gè)參數(shù)表示，在參數(shù)估計(jì)過程中只反演這些分區(qū)的參數(shù)的值[1]。然而這種區(qū)域的劃分往往帶有一定的主觀性，不能很好地反映真實(shí)參數(shù)場的空間結(jié)構(gòu)。

近年來，隨著數(shù)值模擬技術(shù)的發(fā)展，向?qū)c(diǎn)法成為求解非均質(zhì)水文地質(zhì)參數(shù)的一種新途徑。向?qū)c(diǎn)法通過在模型區(qū)域中布設(shè)一定數(shù)量的點(diǎn)，估計(jì)這些點(diǎn)處的參數(shù)值并進(jìn)行插值以獲得整個(gè)參數(shù)場的估計(jì)場。該方法與地統(tǒng)計(jì)學(xué)的結(jié)合可有效克服傳統(tǒng)分區(qū)的缺點(diǎn)[2]，獲得的參數(shù)場更符合真實(shí)情況[3]。在先驗(yàn)信息（如地質(zhì)資料）豐富的情況下，可以依據(jù)當(dāng)?shù)氐牡刭|(zhì)情況布置觀測井與向?qū)c(diǎn)，獲得更加精確的水文地質(zhì)參數(shù)場。例如Young S[4]基于Edwards-Trinity 高原和Pecos 峽谷的地質(zhì)資料，通過在不同研究區(qū)域內(nèi)均勻布置向?qū)c(diǎn)，求解了水文地質(zhì)參數(shù)。在缺乏先驗(yàn)信息的情況下，以往的研究都是將向?qū)c(diǎn)與觀測井均勻布置在整個(gè)研究區(qū)域。例如姜蓓蕾[5]使用該種布置方式，討論了向?qū)c(diǎn)個(gè)數(shù)對反演結(jié)果的影響。這種布置方式雖然被廣泛采用，但是學(xué)者們并未對其展開深入研究，主要原因是向?qū)c(diǎn)的數(shù)量關(guān)系可以定量描述，但觀測井與向?qū)c(diǎn)位置描述起來不是十分方便。本文采用觀測井與向?qū)c(diǎn)分布范圍占研究區(qū)面積的比例來描述兩者的位置關(guān)系，該方法可以不受區(qū)域面積和形狀影響，不僅對于規(guī)則研究區(qū)域有效，還可以運(yùn)用到不規(guī)則的研究區(qū)域。在此基礎(chǔ)上，通過理想算例研究了觀測井與向?qū)c(diǎn)不同分布范圍以及滲透系數(shù)場初值對反演結(jié)果的影響，為向?qū)c(diǎn)法在水文地質(zhì)參數(shù)反演中的推廣提供了依據(jù)。

1 計(jì)算程序與算法

1.1 向?qū)c(diǎn)法

向?qū)c(diǎn)法起源于地質(zhì)統(tǒng)計(jì)學(xué)，通過在模型區(qū)域內(nèi)布設(shè)一定數(shù)目的點(diǎn)(向?qū)c(diǎn))，反演這些點(diǎn)處的參數(shù)值，再經(jīng)過克里金插值得到整個(gè)參數(shù)場的估計(jì)場[6]。向?qū)c(diǎn)法反演參數(shù)場的主要步驟見圖1，其中Modflow 程序調(diào)用次數(shù)表示地下水模型Modflow.exe 運(yùn)行次數(shù)，優(yōu)化迭代次數(shù)等于PEST 程序循環(huán)次數(shù)。向?qū)c(diǎn)法用于解決高維不適定問題，即向?qū)c(diǎn)個(gè)數(shù)要大于觀測井個(gè)數(shù)[7]。在先驗(yàn)信息不足的情況下，向?qū)c(diǎn)應(yīng)均勻分布在研究區(qū)內(nèi)。此外，在參數(shù)估計(jì)中不同位置的觀測井對目標(biāo)函數(shù)的貢獻(xiàn)不同，觀測井應(yīng)盡可能的包含更多的含水層信息。如果初始參數(shù)向量與最優(yōu)參數(shù)向量相差較大，參數(shù)估計(jì)問題陷入局部最優(yōu)，難以找到目標(biāo)函數(shù)的全局最優(yōu)解。

圖1 向?qū)c(diǎn)法反演參數(shù)場的主要步驟流程圖Fig.1 Flow chart of the main steps of inversion of parameter field by the pilot point method

1.2 Gauss-Marquardt-Levenberg 迭代算法

大多數(shù)情況下，模型參數(shù)與觀測值之間的關(guān)系都是非線性的，在求解過程中必須將非線性模型線性化。本文主要研究滲透系數(shù)與觀測水位之間的關(guān)系，假設(shè)二者之間的非線性關(guān)系用函數(shù)c=M(b)表示，該函數(shù)可看作是n維參數(shù)空間映射到m維觀測空間[函數(shù)c=M(b)所有待估模型參數(shù)的導(dǎo)數(shù)必須連續(xù)可導(dǎo)]。某一初始滲透系數(shù)向量b0對應(yīng)的模型觀測水位向量為c0，即：

任意向量b與其對應(yīng)向量c的關(guān)系，可用如下泰勒展開式表示：

當(dāng)b-b0足夠小時(shí)，可以忽略式(2)中二階以上的高階項(xiàng)，從而簡化為：

其中矩陣J是函數(shù)c=M(b)在b0處的雅克比矩陣(Jacobian Ma?trix)，共有m行n列，每行的元素為對應(yīng)的待估模型參數(shù)的偏導(dǎo)數(shù)。

目標(biāo)函數(shù)是模型計(jì)算水位與實(shí)際水位之差的平方和，而參數(shù)估計(jì)過程實(shí)際上是最小化目標(biāo)函數(shù)的過程。假設(shè)n階對角矩陣Q的第Z個(gè)對角元素表示的是第i個(gè)實(shí)際觀測值的權(quán)重W的平方，cr表示實(shí)測水位向量，則非線性模型參數(shù)反演問題的目標(biāo)函數(shù)表示為：

目標(biāo)函數(shù)需要反復(fù)迭代更新雅可比矩陣J與參數(shù)向量b以使b-b0到達(dá)最小，定義參數(shù)更新向量b-b0為u，則u可表示為：

假定殘差向量r=cr-c0，則式(5)可寫為：

若目標(biāo)函數(shù)Φ的梯度為g，則向量

在公式(6)中加入Marquardt參數(shù)，使估計(jì)初期參數(shù)更新向量u靠近向量g的反方向，u的表達(dá)式變?yōu)椋?/p>

式中：α為Marquardt 參數(shù)；I為n階單位矩陣；梯度向量g又可定義為g=-2JTQr。

高斯-列文伯格-馬夸爾特迭代算法是在高斯—牛頓迭代法的基礎(chǔ)上進(jìn)行改進(jìn)而來的，用于求解廣義非線性加權(quán)最小二乘最小化問題[8]。該方法在每次迭代過程中，通過對待估模型參數(shù)進(jìn)行求導(dǎo)將非線性模型線性化，同時(shí)更新參數(shù)向量，計(jì)算新的目標(biāo)函數(shù)，不斷重復(fù)迭代優(yōu)化過程，直到目標(biāo)函數(shù)滿足要求。

1.3 PEST程序

PEST( 全稱Parameter Estimate) 是由澳大利亞Water Numerical Computing 咨詢公司開發(fā)的用于參數(shù)估計(jì)與不確定性分析的程序。PEST 程序第一個(gè)版本誕生于1994年，最初只具有非線性參數(shù)估計(jì)功能，經(jīng)過20 多年的發(fā)展，目前已經(jīng)具備了強(qiáng)大的正則化功能與不確定性分析功能。PEST 程序使用高斯-列文伯格-馬夸爾特非線性參數(shù)估計(jì)方法，相比于其他非線性參數(shù)估計(jì)方法，該法可以在有效估計(jì)參數(shù)的同時(shí)節(jié)省大量運(yùn)行模型的次數(shù)。最新版本的PEST 程序支持向?qū)c(diǎn)法并且?guī)в袑ｉT針對地下水問題的程序集，因此可以很方便地與Modflow等地下水模型進(jìn)行藕合使用。

參數(shù)估計(jì)的目的是找到一個(gè)滲透系數(shù)向量，即滲透系數(shù)估計(jì)場，使得模型計(jì)算水位最大限度地?cái)M合實(shí)際觀測水位。為準(zhǔn)確評估先驗(yàn)信息精度與反演精度，采用RMSE（Root Mean Square Error，均方根誤差）作為衡量初始滲透系數(shù)場和滲透系數(shù)估計(jì)場（與實(shí)際滲透系數(shù)場）精度的評判標(biāo)準(zhǔn)，公式如下：

式中：n為參數(shù)個(gè)數(shù)；yt和ya分別表示參數(shù)的真實(shí)值和模型計(jì)算值；RMSE的值越小反演精度越高。

2 理想算例

本文采用確定性地下水?dāng)?shù)值模擬軟件Visual Modflow[9]建立了一個(gè)二維承壓含水層數(shù)值模型，如圖2所示。假定該模型中含水層水平等厚，平面是邊長為1000m 的正方形，底板高程為0 m，頂板高程為5 m，北邊界為定水頭邊界，其余邊界為隔水邊界，上部接受0.000 5 m/d 的面狀補(bǔ)給，模擬計(jì)算時(shí)間為365 d。通過克里金指數(shù)變異函數(shù)插值出地下水滲透系數(shù)對數(shù)場，作為實(shí)際滲透系數(shù)場（注：高斯變異函數(shù)在PEST 程序中運(yùn)行結(jié)果會超出向?qū)c(diǎn)值設(shè)定上下限，不建議使用[10,11]）。在研究區(qū)布置36 口井，作為地下水水位觀測井。由于向?qū)c(diǎn)法中參數(shù)反演屬于不適定問題，向?qū)c(diǎn)數(shù)要大于觀測井?dāng)?shù)，因此在研究區(qū)內(nèi)布置了64個(gè)向?qū)c(diǎn)。

圖2 理想算例中真實(shí)滲透系數(shù)場(log K)Fig.2 The real permeability coefficient field in the ideal case(log K)

在實(shí)際滲透系數(shù)場中選擇16 個(gè)點(diǎn)，這些點(diǎn)的滲透系數(shù)值是已知的，對其進(jìn)行普通的克里金插值獲得初始滲透系數(shù)場（見圖3），其與實(shí)際滲透系數(shù)場之間的RMSE記作R1，表示先驗(yàn)信息精度，此時(shí)的模型計(jì)算水位與實(shí)際水位還有較大差距。經(jīng)過PEST 程序參數(shù)反演后，模型計(jì)算水位能夠最大限度地?cái)M合實(shí)際水位[12]，得到的參數(shù)場稱為滲透系數(shù)估計(jì)場，與實(shí)際滲透系數(shù)場之間的RMSE記作R2，表示反演精度。模型計(jì)算水位與滲透系數(shù)值取自系統(tǒng)，誤差范圍為0.001。

圖3 16個(gè)點(diǎn)(三角形)位置與初始滲透系數(shù)場(log K)Fig.3 The location of 16 points(triangles)and the initial permeability coefficient field(log K)

3 分析與討論

為了研究觀測井、向?qū)c(diǎn)分布范圍以及不同滲透系數(shù)初值對反演結(jié)果的影響，結(jié)合理想案例，分別設(shè)置不同分布范圍的觀測井和向?qū)c(diǎn)模型及不同初值的初始滲透系數(shù)場模型，并對反演結(jié)果進(jìn)行分析與討論。

3.1 不同觀測井分布范圍反演結(jié)果分析

在參數(shù)估計(jì)中，觀測井所包含的含水層信息越多，模型反演精度越高。為了研究不同觀測井分布范圍對反演結(jié)果的影響，將64 個(gè)向?qū)c(diǎn)均勻分布在整個(gè)研究區(qū)，36 口觀測井均勻布置在占研究區(qū)總面積16%、36%、64%、81%和100%的正方形范圍內(nèi)（如圖4所示），正方形中心與研究區(qū)域中心重合。通過PEST 程序?qū)Τ跏紳B透系數(shù)場進(jìn)行反演求解，限于篇幅，僅給出觀測井分布范圍占研究區(qū)總面積16%、81%和100%的結(jié)果，分別見圖5～圖7?？梢钥闯?，當(dāng)觀測井布置范圍適當(dāng)時(shí)，反演結(jié)果較準(zhǔn)確地反映了真實(shí)滲透系數(shù)場的空間分布，說明向?qū)c(diǎn)法能夠有效求解地下水反演中的不適定問題。

圖4 觀測井分布范圍示意圖Fig.4 Schematic diagram of the distribution range of observation wells

圖5 分布范圍為16%滲透系數(shù)估計(jì)場（log K）Fig.5 The distribution range is 16%of the estimated field of permeability coefficient(log K)

圖6 分布范圍為81%滲透系數(shù)估計(jì)場（log K）Fig.6 The distribution range is 81%of the estimated field of permeability coefficient(log K)

圖7 分布范圍為100%滲透系數(shù)估計(jì)場（log K）Fig.7 The distribution range is 100%of the estimated field of permeability coefficient(log K)

觀測井不同分布范圍情景下反演結(jié)果見表1。盡管初始滲透系數(shù)場與真實(shí)滲透系數(shù)場之間的R1 很小，甚至小于R2，但仍然不采用初始滲透系數(shù)場直接作為滲透系數(shù)估計(jì)場。這是因?yàn)楸敬卫硐氚咐姓鎸?shí)滲透系數(shù)場數(shù)據(jù)是已知的，僅作為參照場，可以計(jì)算出R1與R2值。而實(shí)際工程中的真實(shí)滲透系數(shù)場是未知的，擬合地下水水位反演后的也只是滲透系數(shù)估計(jì)場，無法計(jì)算R1與R2的值。

表1 不同觀測井分布范圍反演結(jié)果Tab.1 Inversion results of distribution range of different observation wells

觀測井中實(shí)際水位與計(jì)算水位的RMSE變化區(qū)間為2.36×10-7～6.35×10-4，表明水位擬合十分良好。當(dāng)觀測井分布范圍為16%時(shí)，R2值為0.447，反演效果最不理想；當(dāng)觀測井分布范圍從16%增加到64%時(shí)，對應(yīng)的R2 值依次減小(0.447>0.156>0.069)，說明隨著觀測井分布范圍增加，獲取的含水層參數(shù)信息增加，反演精度逐漸提高；當(dāng)分布范圍為81%與100%時(shí)，對應(yīng)的R2 值為0.056、0.057，再次說明觀測井分布范圍增加到一定程度時(shí)，獲取的含水層信息不再有太多變化，反演結(jié)果能夠有效反映滲透系數(shù)場的空間分布，反演精度逐漸趨于穩(wěn)定。

3.2 不同向?qū)c(diǎn)分布范圍反演結(jié)果分析

為了研究不同向?qū)c(diǎn)分布范圍對參數(shù)估計(jì)的影響，將36口觀測井均勻分布在整個(gè)研究區(qū)，而64 個(gè)向?qū)c(diǎn)均勻分布在占研究區(qū)總面積的16%、36%、64%、81%和100%正方形范圍內(nèi)，正方形中心依然與研究區(qū)域中心重合。

向?qū)c(diǎn)不同分布范圍情景下反演結(jié)果見表2。地下水位觀測值與計(jì)算值的RMSE位于3.73×10-7～5.53×10-7之間，研究區(qū)水位擬合良好。當(dāng)向?qū)c(diǎn)分布范圍為16%、36%，64%時(shí)，對應(yīng)的R2值依次減小(0.194>0.065>0.045)；而當(dāng)向?qū)c(diǎn)分布范圍為81%、100%時(shí)，R2 值分別為0.052、0.057，雖比64%對應(yīng)的R2值0.045大，但增加幅度很小，基本保持穩(wěn)定，說明增加向?qū)c(diǎn)分布范圍同樣可以獲得更高精度的滲透系數(shù)場。當(dāng)向?qū)c(diǎn)集中在部分區(qū)域時(shí)，獲取的滲透系數(shù)參數(shù)信息較少，導(dǎo)致目標(biāo)函數(shù)在負(fù)方向上步長變小，下降速率變慢。在本次模擬中，雖然目標(biāo)函數(shù)初始值相同，但步長越小，達(dá)到目標(biāo)函數(shù)要求所調(diào)用Modflow程序和優(yōu)化迭代的次數(shù)就會越多[13]。

表2 不同向?qū)c(diǎn)分布范圍反演結(jié)果Tab.2 Inversion results of distribution range of different guide points

在參數(shù)化方法中，樣本點(diǎn)的選取對參數(shù)反演結(jié)果影響很大，在研究區(qū)域內(nèi)，抽樣應(yīng)該具有無偏性和一致性。當(dāng)分布范圍較小時(shí)，抽樣數(shù)據(jù)集中在較小的范圍，與真實(shí)滲透系數(shù)場和觀測數(shù)據(jù)差距較大，無法代替整個(gè)參數(shù)組與觀測組。隨著向?qū)c(diǎn)與觀測井分布范圍的增加，選擇的樣本區(qū)域擴(kuò)大，盡管選擇的樣本數(shù)量不變但有效樣本數(shù)增加，反演精度提高，R2 逐漸減小。隨著分布范圍接近整個(gè)研究區(qū)域，有效樣本數(shù)目變化不明顯，R2值基本保持穩(wěn)定。

3.3 不同滲透系數(shù)場初值反演結(jié)果分析

為研究滲透系數(shù)場初值對反演精度的影響，將36 口觀測井與64 個(gè)向?qū)c(diǎn)分別均勻布置在整個(gè)研究區(qū)域內(nèi)，在16 個(gè)已知滲透系數(shù)值點(diǎn)對數(shù)值的基礎(chǔ)上增加0%、-10%、-20%、30%、50%、100%，重新進(jìn)行克里金插值得到不同初值的初始滲透系數(shù)場，對應(yīng)R1 值分別為0.053、0.151、0.297、0.466、0.768、1.527。

在此基礎(chǔ)上進(jìn)行反演，結(jié)果如表3所示?？梢钥闯?，觀測水位與計(jì)算水位RMSE 值在3.33×10-7～4.71×10-7范圍內(nèi)，研究區(qū)水位擬合良好。R2 與R1 呈正相關(guān)，R1 值越小，R2 值也越小，擬合精度越高；反之，擬合效果越差。這是因?yàn)閷τ诜蔷€性問題，如果初始參數(shù)向量b0在參數(shù)空間中的位置與使目標(biāo)函數(shù)最小化的最優(yōu)參數(shù)向量b的位置相差較大時(shí)，參數(shù)估計(jì)問題容易陷入局部最優(yōu)，使目標(biāo)函數(shù)全局最小難以實(shí)現(xiàn)。同時(shí)由于模型要求b-b0盡可能小，目標(biāo)函數(shù)的最小化需要反復(fù)迭代更新雅可比矩陣與向量b，即模型調(diào)用Modflow 程序與優(yōu)化迭代的次數(shù)也會隨著R1 的增加而增加。因此，初始參數(shù)向量b0越接近真實(shí)情況，反演精度越高[14]。

表3 不同初值的初始滲透系數(shù)場反演結(jié)果Tab.3 Inversion results of initial permeability coefficient field with different initial values

4 結(jié) 論

本文通過PEST 程序與地下水流模擬軟件Modflow 耦合，對設(shè)定的理想算例進(jìn)行滲透系數(shù)的高維反演，研究了向?qū)c(diǎn)與觀測井分布范圍、滲透系數(shù)場初值對反演水文地質(zhì)參數(shù)的影響，得出以下結(jié)論：

（1）向?qū)c(diǎn)法是一種有效反演水文地質(zhì)參數(shù)的方法，在參數(shù)反演應(yīng)用方面具有較高的優(yōu)越性，避免了人為概化分區(qū)的局限性，反演結(jié)果能較好地反映參數(shù)場的非均質(zhì)性和空間結(jié)構(gòu)。

（2）隨著觀測井分布范圍的增加，參數(shù)反演精度逐漸提高并最終達(dá)到穩(wěn)定，在觀測井分布范圍接近整個(gè)研究區(qū)域時(shí)達(dá)到最高。

（3）隨著向?qū)c(diǎn)分布范圍增加，反演精度逐漸提高并趨于穩(wěn)定，調(diào)用Modflow 程序和優(yōu)化迭代的次數(shù)變少。向?qū)c(diǎn)分布范圍覆蓋整個(gè)區(qū)域時(shí)，模型調(diào)用次數(shù)與優(yōu)化迭代次數(shù)最少，此時(shí)該種布置方式較為簡單，推薦實(shí)際工程中使用。

（4）采用向?qū)c(diǎn)方法反演參數(shù)時(shí)，初始參數(shù)信息越接近真實(shí)情況，反演精度愈高。