基于非均質強弱程度評價GMS求參模塊適用性

(貴州大學 資源與環境工程學院，貴陽 550025)

0 引 言

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

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

本文著重圍繞有限差分法進行，以均質模型概化為對照，模擬非均質強弱程度不同的含水層的抽水試驗，利用GMS軟件中PES Process和Pest兩個模塊交替運行，導入各條件下的抽水試驗數據進行模型識別，比較該軟件在不同條件下所得的結果，從而判別該軟件反求參數的適用精度以及適用范圍。

1 GMS求參模塊

地下水模擬系統(Groundwater Modeling System)，簡稱GMS[2]，提供3種反求參數模塊的界面，包括內部過程模塊(MODFLOW 2000 PES process)和自動參數估計的外部效應模塊(PEST、 UCODE)。

1.1 PES模塊

PES為MODFLOW的內置模塊，使用高斯牛頓法進行參數優化，采用“參數分區法”進行模型校正。

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

PES Process找到最優解，MODFLOW則會自動讀入最優解進行正演，同時輸出水頭值分布情況。

1.2 PEST模塊

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

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

2 正向模擬

2.1 滲透系數場

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

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

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

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

圖1 非均質性不同的滲透系數場

由圖1可知，各方向相關長度不變，隨著C值增大，空間變異性增強，流通性能較好的水流通道逐漸增多，滲透系數值出現異常。相鄰兩點的滲透系數值相差越大，含水層的非均質程度增強。

2.2 抽水試驗概況

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

該模擬場內有P1、P2、…、P9共9個觀測孔，圍繞抽水井規則分布。研究區抽水試驗平面布置情況見圖2。

圖2 抽水試驗布置圖

在其他水文地質條件不變的情況下，上述模型水頭計算值作為野外觀測值進行參數計算，對比計算結果與實測滲透系數差異以評估GMS軟件反求參數模塊的效果。

3 結果分析

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

3.1 均質含水層

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

圖3 均質含水層計算值與觀測值的擬合結果

3.2 非均質含水層

在0～100內按一定梯度改變Contribution(基臺值)，Range(各方向相關長度)為定值10，從而隨機生成非均質強弱程度不同的共16組數據。

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

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

表1 不同程度非均質含水層參數計算結果對比

圖4 計算值與觀測值隨非均質強弱程度改變的關系曲線

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

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

3.3 討 論

對于個別反演所得計算值異常出現最值的現象，不能說明試驗沒有價值。針對這一問題，現做以下分析：

1)隨機參數模擬本身就具有隨機性。C值較小時，計算所得等效滲透系數值是圍繞設定的平均值上下產生細微的波動，而C值的增大其波動幅度也隨之增大，由此加大了最值出現的可能。

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

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

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

4 結 論

同一地下水流模型中，變程(R)作為定量時，非均質強弱程度的整體變化趨勢是隨著基臺值(C)的增加而增強。基臺值增大，空間相關度增大，在一定程度上加大了含水層介質的連通性，導致出現數值較高的滲透系數。

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

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

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

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