基于非均质强弱程度评价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之间的某一定值，具体的临界点需要更加细化的试验数据和大量的模拟对照。