(貴州大學(xué) 資源與環(huán)境工程學(xué)院，貴陽 550025)

0 引 言

地下水?dāng)?shù)值模擬是評(píng)價(jià)人類工程地質(zhì)活動(dòng)對(duì)地下水的影響、模擬自然界水文地質(zhì)過程發(fā)生和發(fā)展的主要方法和手段之一[1]。為建立可以較真實(shí)地反映研究區(qū)水文地質(zhì)條件和地下水流運(yùn)動(dòng)基本特征的地下水流數(shù)值模型，確定模型中的水文地質(zhì)參數(shù)是必不可缺的研究工作。而伴隨著地下水研究尺度的擴(kuò)大、精度的提高以及復(fù)雜性的增加，使得地下水模型的設(shè)定條件趨于復(fù)雜，導(dǎo)致自動(dòng)優(yōu)化的求解效率和計(jì)算精度下降，因此選用合理科學(xué)的參數(shù)估計(jì)程序顯得十分重要。

GMS軟件中一系列的反演求參數(shù)工具，采用自動(dòng)交替迭代的方式，不論是PES Process模塊的內(nèi)部調(diào)用還是PEST模塊的外源調(diào)用都能夠在更短的時(shí)間內(nèi)達(dá)到比人工-試估法更高的精度。

本文著重圍繞有限差分法進(jìn)行，以均質(zhì)模型概化為對(duì)照，模擬非均質(zhì)強(qiáng)弱程度不同的含水層的抽水試驗(yàn)，利用GMS軟件中PES Process和Pest兩個(gè)模塊交替運(yùn)行，導(dǎo)入各條件下的抽水試驗(yàn)數(shù)據(jù)進(jìn)行模型識(shí)別，比較該軟件在不同條件下所得的結(jié)果，從而判別該軟件反求參數(shù)的適用精度以及適用范圍。

1 GMS求參模塊

地下水模擬系統(tǒng)(Groundwater Modeling System)，簡稱GMS[2]，提供3種反求參數(shù)模塊的界面，包括內(nèi)部過程模塊(MODFLOW 2000 PES process)和自動(dòng)參數(shù)估計(jì)的外部效應(yīng)模塊(PEST、 UCODE)。

1.1 PES模塊

PES為MODFLOW的內(nèi)置模塊，使用高斯牛頓法進(jìn)行參數(shù)優(yōu)化，采用“參數(shù)分區(qū)法”進(jìn)行模型校正。

對(duì)于建立完全包含觀測層數(shù)據(jù)的地下水流模型，具體操作[3]如下：①選擇PES為其反向模塊；②對(duì)參數(shù)區(qū)進(jìn)行劃分，給模型的每個(gè)區(qū)域設(shè)置參數(shù)分區(qū)的關(guān)鍵值，包括滲透系數(shù)、大氣降水入滲率等，并將其映射到MODFLOW的網(wǎng)格單元；③為PES賦初始水頭，以便于模型迅速收斂；④編輯各分區(qū)的參數(shù)，設(shè)置最大迭代次數(shù)、參數(shù)估計(jì)截止數(shù)等；⑤經(jīng)過計(jì)算機(jī)檢驗(yàn)和運(yùn)行MODFLOW，自動(dòng)迭代求得觀測值和計(jì)算值的差值函數(shù)的最小值。

PES Process找到最優(yōu)解，MODFLOW則會(huì)自動(dòng)讀入最優(yōu)解進(jìn)行正演，同時(shí)輸出水頭值分布情況。

1.2 PEST模塊

PEST是利用現(xiàn)有模型對(duì)不便于直接測量的參數(shù)進(jìn)行推測估值的計(jì)算方法，通常將其用途分為三大類：數(shù)據(jù)解譯、模型校準(zhǔn)、預(yù)測分析。PEST 通過不斷地調(diào)整模型參數(shù)的值，直至模型的計(jì)算值與實(shí)驗(yàn)室或野外觀測獲取的相應(yīng)值盡可能地靠近，以此獲得模型參數(shù)的估計(jì)值[4-5]。

PEST的核心算法是用Gauss-Marquardt-Leven-berg算法對(duì)目標(biāo)函數(shù)求解最小值，即模型參數(shù)的計(jì)算值與實(shí)際觀測值的差異函數(shù)，模型的計(jì)算值主要依賴于模型參數(shù)[6]。對(duì)線性模型來說，只需一步就可以完成參數(shù)優(yōu)化；對(duì)于非線性模型，則需要不斷地迭代以求解參數(shù)的最優(yōu)值。而每次迭代開始前，都采用當(dāng)前狀態(tài)下的最優(yōu)參數(shù)集，通過泰勒級(jí)數(shù)展開，使得模型參數(shù)與模型觀測值之間線性化。如此反復(fù)運(yùn)行，并通過比較之前迭代的參數(shù)值與當(dāng)前迭代的參數(shù)值和目標(biāo)函數(shù)估計(jì)值的變化，由PEST選擇是否繼續(xù)進(jìn)行迭代優(yōu)化。

2 正向模擬

2.1 滲透系數(shù)場

給定一個(gè)滲透系數(shù)(K)為36 m/d的均質(zhì)承壓含水層，正向模擬得到抽水試驗(yàn)的觀測孔水頭值的數(shù)據(jù)。

使用GMS軟件內(nèi)置的高斯模擬法生成等概率的滲透系數(shù)的隨機(jī)空間分布。具體應(yīng)用生成滲透系數(shù)場時(shí)，只需輸入滲透系數(shù)的地質(zhì)統(tǒng)計(jì)參數(shù)，包括其取值范圍、均值、方差、變異函數(shù)類型、各方向相關(guān)長度等，快速生成需要的滲透系數(shù)隨機(jī)場[7-8]?？臻g相關(guān)度與基臺(tái)值的取值有關(guān)，各方向相關(guān)長度(R)一定時(shí)，基臺(tái)值(C)較小，則參數(shù)分布的集中趨勢較明顯，出現(xiàn)異常值的概率較低；反之亦然。

非均質(zhì)承壓含水層，以Range(變程)為定量10，改變Contribution(基臺(tái)值)，從而改變含水層非均質(zhì)強(qiáng)弱的程度。本文運(yùn)用高斯模擬改變C、R值，依次生成C5R10、C10R10、…、C100R10總共16組非均質(zhì)性強(qiáng)弱不同的試驗(yàn)數(shù)據(jù)，并通過利用該數(shù)據(jù)進(jìn)行正向模擬得到16組非均質(zhì)程度不同的含水層抽水試驗(yàn)觀測孔的數(shù)據(jù)。

改變C值生成的滲透系數(shù)場見圖1。圖1中僅給出C5R10、C10R10、C20R10、C40R10、C60R10和C100R10這6組具有代表性的滲透系數(shù)分布圖。其中，若C=5，R=10，則以C5R10的形式表示這一含水層的非均質(zhì)程度，且以C表示基臺(tái)值，R表示各方向相關(guān)長度，該表達(dá)方式只應(yīng)用于本文。

圖1 非均質(zhì)性不同的滲透系數(shù)場

由圖1可知，各方向相關(guān)長度不變，隨著C值增大，空間變異性增強(qiáng)，流通性能較好的水流通道逐漸增多，滲透系數(shù)值出現(xiàn)異常。相鄰兩點(diǎn)的滲透系數(shù)值相差越大，含水層的非均質(zhì)程度增強(qiáng)。

2.2 抽水試驗(yàn)概況

利用GMS軟件的MODFLOW模塊建立理想的地下水流模型，流域面積為791×478 m2，采用規(guī)則矩形網(wǎng)格剖分：水平向300×200，共計(jì)60 000個(gè)單元格，垂向上為單層。設(shè)定無窮邊界以降低邊界對(duì)流場的影響，使抽水井以定流量300 m3/d抽水能夠達(dá)到相對(duì)穩(wěn)定流狀態(tài)。頂板標(biāo)高20 m，底板標(biāo)高0 m，含水層厚度20 m，初始水頭20 m。該含水層為承壓含水層，左右兩邊為定水頭邊界，左邊界水頭值為30 m，右邊界水頭值為25 m。

該模擬場內(nèi)有P1、P2、…、P9共9個(gè)觀測孔，圍繞抽水井規(guī)則分布。研究區(qū)抽水試驗(yàn)平面布置情況見圖2。

圖2 抽水試驗(yàn)布置圖

在其他水文地質(zhì)條件不變的情況下，上述模型水頭計(jì)算值作為野外觀測值進(jìn)行參數(shù)計(jì)算，對(duì)比計(jì)算結(jié)果與實(shí)測滲透系數(shù)差異以評(píng)估GMS軟件反求參數(shù)模塊的效果。

3 結(jié)果分析

均質(zhì)含水層反演求參數(shù)值較為簡單，模型收斂精度高，用程序反向模擬得出的計(jì)算值誤差小，與觀測值擬合情況好。因此，本文側(cè)重于研究非均質(zhì)強(qiáng)弱程度不同的含水層的數(shù)據(jù)反演所得結(jié)果的誤差變化，而只設(shè)置一組均質(zhì)承壓含水層的數(shù)據(jù)反演作對(duì)照試驗(yàn)。

3.1 均質(zhì)含水層

地下水流模型中給定含水層的滲透系數(shù)為K=36，通過程序反向模擬求得的等效滲透系數(shù)參數(shù)值為HK=36.003 4，誤差僅為0.003 4。該含水層中抽水試驗(yàn)的觀測孔的數(shù)據(jù)與GMS軟件反演所得計(jì)算值相吻合，波動(dòng)微弱到幾近可以忽略不計(jì)，擬合結(jié)果見圖3。

圖3 均質(zhì)含水層計(jì)算值與觀測值的擬合結(jié)果

3.2 非均質(zhì)含水層

在0～100內(nèi)按一定梯度改變Contribution(基臺(tái)值)，Range(各方向相關(guān)長度)為定值10，從而隨機(jī)生成非均質(zhì)強(qiáng)弱程度不同的共16組數(shù)據(jù)。

非均質(zhì)含水層的模擬結(jié)果相較于均質(zhì)含水層而言要復(fù)雜很多，不僅出現(xiàn)了模擬所得計(jì)算值徒增突降的現(xiàn)象，而且個(gè)別模擬所得的結(jié)果大幅波動(dòng)，誤差變動(dòng)較大，一定程度上影響該試驗(yàn)整體的精度判斷。特別是C25R10、C35R10和C50R10，出現(xiàn)了滲透系數(shù)設(shè)定范圍內(nèi)的最大值100。

C5R10、…、C20R10這4組數(shù)據(jù)反向模擬所求得的等效滲透系數(shù)值與模型正演給定的滲透系數(shù)值大致吻合，誤差值較小。在這一范圍內(nèi)，GMS軟件反求參數(shù)的精度是適用的，且具有一定的參考價(jià)值。見表1和圖4。

表1 不同程度非均質(zhì)含水層參數(shù)計(jì)算結(jié)果對(duì)比

圖4 計(jì)算值與觀測值隨非均質(zhì)強(qiáng)弱程度改變的關(guān)系曲線

非均質(zhì)含水層中，R=10、C值取值小于20時(shí)，利用該組數(shù)據(jù)進(jìn)行隨機(jī)參數(shù)模擬所得的值與反向模擬求出的等效滲透系數(shù)值之間誤差較小，具有參考價(jià)值；R=10、C值取值大于20時(shí)，給定的滲透系數(shù)值與計(jì)算的等效滲透系數(shù)值之間誤差有逐漸增大的趨勢。

由表1和圖4提供的數(shù)據(jù)可以看出，隨著C值的增大，隨機(jī)參數(shù)模擬所得數(shù)據(jù)的隨機(jī)性增大，含水層的非均質(zhì)性不斷增強(qiáng)，反演求出的等效滲透系數(shù)值與給定的滲透系數(shù)值之間的誤差也隨之增大。

3.3 討 論

對(duì)于個(gè)別反演所得計(jì)算值異常出現(xiàn)最值的現(xiàn)象，不能說明試驗(yàn)沒有價(jià)值。針對(duì)這一問題，現(xiàn)做以下分析：

1)隨機(jī)參數(shù)模擬本身就具有隨機(jī)性。C值較小時(shí)，計(jì)算所得等效滲透系數(shù)值是圍繞設(shè)定的平均值上下產(chǎn)生細(xì)微的波動(dòng)，而C值的增大其波動(dòng)幅度也隨之增大，由此加大了最值出現(xiàn)的可能。

2)滲透系數(shù)可以用對(duì)數(shù)正態(tài)分布來描述其空間變異性，需考慮隨機(jī)模擬引起的不確定性和隨機(jī)模型中參數(shù)的不確定性的影響。

3)一對(duì)C值與R值進(jìn)行隨機(jī)參數(shù)模擬時(shí)僅生成一個(gè)滲透系數(shù)場，用于反演求參數(shù)，這也增大了誤差加大的可能性。實(shí)際工作中，應(yīng)利用由高斯模擬法生成多組等概率的滲透系數(shù)的隨機(jī)空間分布，進(jìn)行反向模擬，取其平均值作對(duì)比，結(jié)果會(huì)更加精確。

值得注意的是，兩邊界都設(shè)置成定流量邊界，使得抽水井無法對(duì)流場產(chǎn)生影響，導(dǎo)致模型反演不收斂，無論滲透系數(shù)為何值都可能與該觀測孔的數(shù)據(jù)擬合良好，形成解的不唯一性問題。同時(shí)，數(shù)值模擬求取的水文地質(zhì)參數(shù)受多個(gè)因素影響(如模擬驗(yàn)證期的長短、觀測資料的豐富程度與準(zhǔn)確度等)，只能算作曲線擬合參數(shù)或等效參數(shù)，與實(shí)際水文地質(zhì)參數(shù)具有一定的出入。

4 結(jié) 論

同一地下水流模型中，變程(R)作為定量時(shí)，非均質(zhì)強(qiáng)弱程度的整體變化趨勢是隨著基臺(tái)值(C)的增加而增強(qiáng)?；_(tái)值增大，空間相關(guān)度增大，在一定程度上加大了含水層介質(zhì)的連通性，導(dǎo)致出現(xiàn)數(shù)值較高的滲透系數(shù)。

根據(jù)C25R10、C35R10、C50R10等3組數(shù)據(jù)異常的討論分析，由于參數(shù)隨機(jī)模擬的隨機(jī)性，可暫時(shí)性地忽略該數(shù)據(jù)對(duì)整體結(jié)果的影響。由此得出以下結(jié)論：

1)均質(zhì)含水層中，模型收斂精度高，GMS軟件的反求參數(shù)模塊是適用的，且精度較高。

2)對(duì)于非均質(zhì)含水層，基臺(tái)值(C)取小于20的值，變程(R)不變，該軟件的精度有所保證，反演結(jié)果與給定值之間誤差較小，具有一定的參考價(jià)值。

3)R=10時(shí)，程序適用范圍的界限初步估算在C值取20～25之間的某一定值，具體的臨界點(diǎn)需要更加細(xì)化的試驗(yàn)數(shù)據(jù)和大量的模擬對(duì)照。