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

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

0 引 言

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

近年來(lái)，隨著數(shù)值模擬技術(shù)的發(fā)展，向?qū)c(diǎn)法成為求解非均質(zhì)水文地質(zhì)參數(shù)的一種新途徑。向?qū)c(diǎn)法通過(guò)在模型區(qū)域中布設(shè)一定數(shù)量的點(diǎn)，估計(jì)這些點(diǎn)處的參數(shù)值并進(jìn)行插值以獲得整個(gè)參數(shù)場(chǎng)的估計(jì)場(chǎng)。該方法與地統(tǒng)計(jì)學(xué)的結(jié)合可有效克服傳統(tǒng)分區(qū)的缺點(diǎn)[2]，獲得的參數(shù)場(chǎng)更符合真實(shí)情況[3]。在先驗(yàn)信息（如地質(zhì)資料）豐富的情況下，可以依據(jù)當(dāng)?shù)氐牡刭|(zhì)情況布置觀測(cè)井與向?qū)c(diǎn)，獲得更加精確的水文地質(zhì)參數(shù)場(chǎng)。例如Young S[4]基于Edwards-Trinity 高原和Pecos 峽谷的地質(zhì)資料，通過(guò)在不同研究區(qū)域內(nèi)均勻布置向?qū)c(diǎn)，求解了水文地質(zhì)參數(shù)。在缺乏先驗(yàn)信息的情況下，以往的研究都是將向?qū)c(diǎn)與觀測(cè)井均勻布置在整個(gè)研究區(qū)域。例如姜蓓蕾[5]使用該種布置方式，討論了向?qū)c(diǎn)個(gè)數(shù)對(duì)反演結(jié)果的影響。這種布置方式雖然被廣泛采用，但是學(xué)者們并未對(duì)其展開(kāi)深入研究，主要原因是向?qū)c(diǎn)的數(shù)量關(guān)系可以定量描述，但觀測(cè)井與向?qū)c(diǎn)位置描述起來(lái)不是十分方便。本文采用觀測(cè)井與向?qū)c(diǎn)分布范圍占研究區(qū)面積的比例來(lái)描述兩者的位置關(guān)系，該方法可以不受區(qū)域面積和形狀影響，不僅對(duì)于規(guī)則研究區(qū)域有效，還可以運(yùn)用到不規(guī)則的研究區(qū)域。在此基礎(chǔ)上，通過(guò)理想算例研究了觀測(cè)井與向?qū)c(diǎn)不同分布范圍以及滲透系數(shù)場(chǎng)初值對(duì)反演結(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é)，通過(guò)在模型區(qū)域內(nèi)布設(shè)一定數(shù)目的點(diǎn)(向?qū)c(diǎn))，反演這些點(diǎn)處的參數(shù)值，再經(jīng)過(guò)克里金插值得到整個(gè)參數(shù)場(chǎng)的估計(jì)場(chǎng)[6]。向?qū)c(diǎn)法反演參數(shù)場(chǎng)的主要步驟見(jiàn)圖1，其中Modflow 程序調(diào)用次數(shù)表示地下水模型Modflow.exe 運(yùn)行次數(shù)，優(yōu)化迭代次數(shù)等于PEST 程序循環(huán)次數(shù)。向?qū)c(diǎn)法用于解決高維不適定問(wèn)題，即向?qū)c(diǎn)個(gè)數(shù)要大于觀測(cè)井個(gè)數(shù)[7]。在先驗(yàn)信息不足的情況下，向?qū)c(diǎn)應(yīng)均勻分布在研究區(qū)內(nèi)。此外，在參數(shù)估計(jì)中不同位置的觀測(cè)井對(duì)目標(biāo)函數(shù)的貢獻(xiàn)不同，觀測(cè)井應(yīng)盡可能的包含更多的含水層信息。如果初始參數(shù)向量與最優(yōu)參數(shù)向量相差較大，參數(shù)估計(jì)問(wèn)題陷入局部最優(yōu)，難以找到目標(biāo)函數(shù)的全局最優(yōu)解。

圖1 向?qū)c(diǎn)法反演參數(shù)場(chǎng)的主要步驟流程圖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ù)與觀測(cè)值之間的關(guān)系都是非線(xiàn)性的，在求解過(guò)程中必須將非線(xiàn)性模型線(xiàn)性化。本文主要研究滲透系數(shù)與觀測(cè)水位之間的關(guān)系，假設(shè)二者之間的非線(xiàn)性關(guān)系用函數(shù)c=M(b)表示，該函數(shù)可看作是n維參數(shù)空間映射到m維觀測(cè)空間[函數(shù)c=M(b)所有待估模型參數(shù)的導(dǎo)數(shù)必須連續(xù)可導(dǎo)]。某一初始滲透系數(shù)向量b0對(duì)應(yīng)的模型觀測(cè)水位向量為c0，即：

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

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

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

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

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

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

若目標(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)而來(lái)的，用于求解廣義非線(xiàn)性加權(quán)最小二乘最小化問(wèn)題[8]。該方法在每次迭代過(guò)程中，通過(guò)對(duì)待估模型參數(shù)進(jìn)行求導(dǎo)將非線(xiàn)性模型線(xiàn)性化，同時(shí)更新參數(shù)向量，計(jì)算新的目標(biāo)函數(shù)，不斷重復(fù)迭代優(yōu)化過(guò)程，直到目標(biāo)函數(shù)滿(mǎn)足要求。

1.3 PEST程序

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

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

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

2 理想算例

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

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

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

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

3 分析與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

（2）隨著觀測(cè)井分布范圍的增加，參數(shù)反演精度逐漸提高并最終達(dá)到穩(wěn)定，在觀測(cè)井分布范圍接近整個(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í)該種布置方式較為簡(jiǎn)單，推薦實(shí)際工程中使用。

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