封洞库低渗介质地下水渗流场数值模型的参数反演

唐连松，陈 刚，胡 成

(中国地质大学(武汉)环境学院，湖北 武汉 430074)

在地下水封洞库的建设和运营中，地下水压力的分布直接影响洞库的水封效果，因此准确地刻画地下水渗流场对于保证地下水封洞库的安全和稳定运行至关重要[1]。由于地下水封洞库工程的复杂性，在其设计、施工和运营阶段，一般都需要借助数值模拟方法对地下洞库所在水文地质单元的地下水渗流场特征、水位动态变化及涌水量进行模拟和预测。在数值模拟的过程中，模型参数的选取是决定模型计算准确度的重要影响因素之一。以往的数值模拟工作中，往往是依据当地地质条件将研究区简单地划分为若干均质区，分区内部为均质，分区间为非均质，以整个分区为单个参数进行参数估计，这种方法被称为传统参数分区法或常规分区法[2]，该方法虽然可以极大地减少模型参数的数量，但无法精确地刻画含水层的渗透性分布，尤其是对于强烈非均质各向异性的裂隙岩体含水层，由于裂隙岩体含水层的渗透系数表现出强烈的空间变异性，仅利用几个钻孔试验数据进行常规统计，根本无法全面地反映出整个研究区含水层渗透性的分布特征。而采用参数反演分析法对影响地下水渗流场模拟的重要参数——渗透系数进行优化可以提高模拟的准确度。反演分析法实际上是一种参数优化方法，通过拟合地下水水位与实际地下水水位的接近程度来确定裂隙岩体含水层的渗透系数[3]，接近程度越高，说明拟合得到的渗透系数越精确。参数反演分析法分为直接反演法和间接反演法两种[4]。直接反演法就是将实际的地下水水位分布作为已知数(已有观测值)，将含水层渗透系数作为未知数进行求解，当模拟得出的地下水水位分布与实际地下水水位分布很接近时，求解完成。间接反演法是在前期水文地质勘探及调查的基础上，综合研究区岩体结构和水文地质特性，从总体上确定裂隙岩体含水层的相对渗透性，并圈定渗透系数的大致范围，再以此作为已知条件去拟合地下水水位分布，并通过比较地下水水位的计算值和实际值，逐步修正含水层渗透系数以达到较好的模拟效果[5]。

本文利用FEFLOW软件建立了烟台某地下水封洞库运营阶段的研究区地下水渗流场数值模型，由于花岗岩体含水层的强烈非均质性，其模拟结果与实测数据存在一定的误差，尤其是洞库地下水涌水量模拟值与实测数据有较大偏差，因此，进一步以现场水力试验数据作为先验信息，利用地下压力计的压力监测数据及工程结构内部的流量监测数据，对建立的数值模型参数进行了反演，在提高模型模拟精度的同时刻画了强烈非均质含水介质的渗透系数场。

1 研究区概况

研究区位于山东半岛北部、烟台市西侧，为某公司地下水封洞库所在场地。该洞库工程包括1个丙烷洞库、1个丁烷洞库、1个LPG洞库，3个洞库在平面上呈“品”字形分布。其中，位于北部的丙烷洞库总库容为50万m3，主洞室所在高程为-120～-146 m(黄海高程)，水幕巷道所在高程为-94～-100 m(黄海高程)；位于西侧的丁烷洞库和位于东侧的LPG洞库总库容均为25万m3，主洞室所在高程为-90～-116 m(黄海高程)，水幕巷道所在高程为-64～-70 m(黄海高程)。

研究区含水介质主要为第四系覆盖层和中生代燕山早期中粗粒黑云母二长花岗岩，地下水的主要类型为松散岩类孔隙水、浅层基岩网状裂隙水、深层基岩脉状裂隙水。研究区地下水主要接受大气降水的入渗补给，沿地下裂隙网络向中西部低洼处径流，排泄于第四系残坡积或冲洪积层中，或被人工开采。研究区地下水动态监测数据表明，浅层地下水随季节动态变化明显，深层裂隙水长期保持相对稳定。研究区的基本地质概况见图1。

图1 研究区地质概况Fig.1 Geology introduction of the study area

2 数值模型的建立

2. 1 模型的边界和数学模型

根据研究区的地质条件，将研究区概化为非均质各向异性三维潜水流，模型上边界为入渗补给、蒸发排泄边界；下边界(标高为-176 m)为隔水边界，东部和西部边界是沿流域分水岭划分的边界，在模型中概化为隔水边界；南部边界部分为隔水边界，部分为第二类边界(定流量边界)，流量根据区域地下水径流模数计算给出；北部为黄海，是研究区的最低排泄基准面，定为第一类边界。洞室和水幕系统等工程结构按第一类边界处理。依照上述条件建立数学模型如下：

式中：Ω为研究区域；Kxx、Kyy、Kzz分别为x、y、z主方向的渗透系数(m/d)；H0为初始地下水水位(m)；H1为洞库或水幕水位(m)；S0为地下水自由表面；S1为第一类边界，表示黄海、河流、洞库和水幕的位置；μs为单位储水系数(1/m)；μd为给水度；ε为源汇项;W为潜水面上的蒸发排泄，降雨入渗补给，井的抽水量和泉的排泄量等；S2为第二类边界；q为第二类边界的水流通量(m/d)。

2. 2 模型的网格划分和参数赋值

2.2.1 模型的网格划分

采用三角网格对模型进行剖分，对洞室及水幕系统进行局部加密，在平面上将4 993.7 m×5 165 m的研究区剖分为49 308个有限单元格，在垂向上按现场水力试验分段(现场压水试验按10 m一段进行)及洞库工程结构将研究区共划分为18层，得到研究区地下水渗流场数值模型的三维结构见图2。

图2 研究区地下水渗流场数值模型三维结构Fig.2 Three-dimensional structure of the model in the study area

2.2.2 模型的参数赋值

根据研究区地质图、水文地质图，并结合现场水文地质试验结果，对模型参数进行赋值，其中第四系覆盖层和岩体强风化层即模型的第1层和第2层按参数分区法赋值，其参数分区情况见图3。模型第1层为第四系覆盖层，在第1层分区中,区域1、2为低山区域且第四系覆盖层较薄，区域3、4为残坡积物堆积形成的第四系覆盖层，区域5、6为九曲河冲积物形成，区域7为九曲河河道，区域8、9为水田；模型第2层为岩体强风化层，其中区域3的风化程度稍强于区域1、2。模型中其他各层的参数利用现场钻孔压水试验结果的平均值进行赋值。由此得到研究区地下水渗流场数值模型各层各分区具体参数的初始值见表1。

图3 研究区地下水渗流场数值模型第1、2层参数分区图Fig.3 Distribution of hydrogeological parameters in layer 1 and layer 2 of the numerical model of the groundwater seapage field in the study area

分层分区号水平渗透系数/(m·d-1)垂直渗透系数/(m·d-1)给水度单位储水系数/m-111.555 200.155 5200.101×10-421.296 000.129 6000.101×10-4312.960 001.296 0000.201×10-4第四系覆48.640 000.864 0000.201×10-4盖层(模50.728 000.072 8000.201×10-4型第1层)62.592 000.259 2000.251×10-4786.400 008.640 0000.301×10-4817.280 001.728 0000.251×10-499.750 000.975 0000.251×10-4强风化层10.355 600.035 5600.025×10-5(模型第20.112 320.011 2320.025×10-52层)30.864 000.086 4000.055×10-5丁烷、LPG洞库水幕层(模型第6层)0.003 000.000 3007×10-7丙烷洞库水幕层(模型第10层)0.002 050.000 2057×10-7

2.3 地下洞库及其周围区域地下水压力及水位的监测

在地下洞库运营过程中，为了保证洞库的水封条件长期有效，需要对地下洞库及其周围区域的地下水压力进行全方位监测，因此在部分水幕孔中安装有地下压力计。地下压力计安装位置为3个洞库的水幕层，其中丁烷、LPG洞库水幕层(标高为-64～-70 m)位于模型第6层，丙烷洞库水幕层(标高为-94～-100 m)位于模型第10层，模型在这两个层位的初始水文地质参数见表1。同时，将部分勘察钻孔改造为地下水水位监测井，对地下外围的地下水水位进行长期监测。

2. 4 模型的模拟精度分析

利用读数较为稳定、数据较为可靠的地下压力计及地下水水位监测井(共31个监测点位)获得的数据对地下水渗流场数值模型计算结果进行校验，其结果见图4。

图4 地下水水位模拟值与实测值的拟合情况Fig.4 Fitting of the simulated and measured groundwater level

由图4可见，受研究区花岗岩含水层强烈非均质性的影响，数值模型的模拟精度一般，地下水水位31个监测点位模拟值与实测值的残差平方和为264.86，地下水水位模拟值与实测值的最大误差为7.4 m。

地下洞库地下水涌水量模拟值与实测值的对比，见表2。

表2 洞库地下水涌水量模拟值与实测值的对比(单位：m3/d)

由表2可知，地下洞库地下水涌水量的模拟值与实测值也存在一定的误差，尤其是丁烷洞库地下水涌水量的模拟值与实测值的误差较大，说明丁烷洞库处可能存在渗透系数较大的强渗透带，这种情况下需要对模型参数进行进一步的优化以提高模型模拟的精度。

3 模型参数反演

3.1 PEST参数反演方法

由于该地下洞库修建于较为完整的花岗岩体中，工程场地的地层岩性较为单一，而且工程建设会避开导水性较强或规模较大的地质构造，因此传统的参数分区法很难应用于这种情况。为此，本文采用PEST程序对影响数值模型计算结果的主要参数——渗透系数进行参数反演，以提高模型模拟的精度。

PEST参数优化算法是一种较为成熟、应用较为广泛的参数优化方法。最初PEST是由Doherty[6]开发出的一种通用参数优化算法，之后Tokin等[7]提出了一种结合Tikhonov正则化及T-SVD正则化方法优点的混合正则化方法(SVD-Assist),进一步优化了PEST算法,使其更加完善。PEST算法在国内外地下水数值模型参数优化工作中取得了较好的效果，董艳辉等[8]利用并行化的PEST算法对大尺度的地下水数值模型进行了参数优化，得到了较为理想的结果;王礼恒等[9]利用假想模型进一步验证了PEST算法在优化地下水流动模型参数上的可靠性;姜蓓蕾等[10]对PEST算法所采用的向导点-正则化方法进行了深入的研究，证明了PEST是一种有效的参数反演方法，同时证明该方法中并不是向导点越多反演精度越高，向导点过多反而会影响计算精度且增加计算量。

PEST程序采用的是Gauss-Marquardt-Levenberg非线性参数估计算法，该算法是一种改进的Gauss-Newton法[11]。该算法在每一次迭代过程中，计算模型结果对各待估参数的导数，利用这些导数将非线性模型线性化，并更新一次参数向量，计算新的目标函数后，再进入下一次迭代过程，通过反复迭代直至目标函数达到最小值。

实际工作中，大部分模型都是非线性的，假设模型的n个初始参数值储存在向量x0中，模型的m个实际观测值储存在向量y0中，则x0与y0的关系为

y0=M(x0)

(1)

根据泰勒公式，任意一组参数向量x与其对应的观测向量y之间的关系可以表示为

(2)

当x→x0时，忽略高阶无穷小项，则式(2)可简化为

y=y0+J(x-x0)

(3)

式中:J为m行n列的偏导数雅克比矩阵:

(4)

由此可得出非线性模型参数反演过程中的目标函数为

Φ=(y-y0-J(x-x0))TQ(y-y0-J(x-x0))

(5)

式中:Q为m行n列的实际观测值的权重矩阵，参数优化过程就是不断迭代更新参数向量x使得目标函数Φ最小化的过程[12]。

3. 2 参数反演范围

图5 模型参数反演范围Fig.5 Range of parameter inversion

由于前期勘察工作及洞库运营阶段的监测工作都是围绕着洞库所在区域开展的，缺少洞库工程外围的实测资料，因此参数反演的范围为以洞库区域为中心的1 000 m×1 000 m正方形区域(见图5)。地下水水位监测点和岩体渗透性原位试验点主要集中在丁烷、LPG洞库水幕层(-64～-70 m)和丙烷洞库水幕层(-94～-100 m)这两个关键层位，这是由于在洞库运营过程中，水幕巷道及水幕孔的持续供水对保证洞库水封效果起着至关重要的作用，而且在深部修建这种特殊的工程结构也会使这两个层位的地下水动态变得更加复杂，对这两个层位的模型参数进行优化是客观评价洞库水封效果、提高水幕系统工作效率的必要工作，因此本次以岩体渗透系数原位试验点和地下水水位监测点的监测数据对丁烷、LPG洞库水幕层和丙烷洞库水幕层这两个关键层位的岩体渗透系数(K)进行参数反演，丁烷、LPG洞库水幕层和丙烷洞库水幕层的原位试验点、监测点及水位分布情况见图6和图7。

图6 丁烷、LPG洞库水幕层的原位试验点、监测点 及水位分布Fig.6 In-situ test points,monitoring points and ground- water level distribution in layer of water curtains of Butane cavern and LPG cavern

图7 丙烷洞库水幕层的原位试验点、监测点及水位分布Fig.7 In-situ test points,monitoring points and ground- water level distribution in layer of water curtain of Propane cavern

3. 3 参数反演结果和模型优化结果

运行PEST程序对丁烷、LPG洞库和丙烷洞库水幕层的岩体渗透系数进行参数反演后，得到洞库岩体渗透系数的反演结果(以lgK表示)以及地下水水位模拟值与实测值的拟合结果见图8和图9。丁烷、LPG洞库地下水水位模型模拟值与实测值的残差平方和为39.82，丙烷洞库地下水水位模拟值与实测值的残差平方和为19.15。

由图8和图9可见，经参数反演后的丁烷、LPG洞库水幕层的岩体渗透系数最大值为0.006 1 m/d，最小值为0.000 03 m/d，丙烷洞库水幕层的岩体渗透系数最大值0.059 m/d，最小值为0.000 105 m/d，其渗透系数在空间上表现出强烈的非均质性。

图8 丁烷、LPG洞库水幕层岩体渗透系数反演结果(以lgK表示)以及地下水水位模拟值与实测值的拟合结果Fig.8 Parameter inversion result of layer of water curtains of Butane cavern and LPG cavern and the fitting result of the groundwater level simulated values and the meosured values

图9 丙烷洞库水幕层岩体渗透系数反演结果(以lgK表示)以及地下水水位模拟值与实测值的拟合结果Fig.9 Parameter inversion result of layer of the water curtain of Propane cavern and the fitting result of the groundwater level simulated values and the measured values

假设岩体渗透系数服从对数正态分布，将其反演结果绘制成岩体渗透系数概率密度曲线，同时将现场16个钻孔中对应层位的压水试验数据与岩体渗透系数概率密度曲线进行对比，其结果见图10。

由图10可见，丁烷、LPG洞库和丙烷洞库水幕层两个层位岩体渗透系数的反演结果与现场实测数据具有相似的统计特征，说明其反演结果能够较为客观地反映研究区的实际情况。根据概率密度函数可知，丁烷、LPG洞库水幕层岩体的渗透系数以0.000 5 m/d为主，丙烷洞库水幕层的渗透系数以0.003 9 m/d为主。

图10 丁烷、LGP洞库和丙烷洞库水幕层岩体渗透系数 反演结果与压水试验数据的对比Fig.10 Comparison of inversion results of rock mass permeability coefficients of Butane,LGP and Propane cavern with pressure water test data

利用参数优化后的数值模型计算得到的各洞库地下水涌水量见表3。

表3 参数优化后洞库地下水涌水量模拟值与实测值的对比(单位：m3/d)

由表3可知，与初始数值模型的计算结果(见表2)相比较，参数优化后的数值模型计算结果与实测数据更加相符，说明经过参数优化后的数值模型具有更高的拟合精度。

4 结论与展望

结合研究区水文地质条件分析和现场水文地质试验数据，建立了烟台某地下水封洞库运营阶段研究区地下水渗流场数值模型，通过对比监测数据发现该均质数值模型的计算结果与实测数据有较大偏差，本研究进一步采用参数反演与现场试验数据相结合的方法来刻画地下水封洞库特定层位的岩体渗透系数，经过参数优化后，数值模型地下洞库拟合精度得到提高，模拟结果表明：丁烷、LPG洞库水幕层的岩体渗透系数以0.000 5 m/d为主，渗透性最强处岩体渗透系数为0.006 1 m/d，但这一层位上没有形成明显的强渗透带；丙烷洞库水幕层的岩体渗透系数以0.003 9 m/d为主，渗透性最强渗透系数为0.059 m/d，这一层位上出现了明显的强渗透带，这也解释了为什么丁烷洞库从开始运行起就出现了地下水涌水量偏大的现象。由此说明，经参数反演后的地下水渗流场数值模型能够更加客观、准确地评价水封洞库的运行状况。由于工程实际情况较为复杂，本文只对岩体渗透系数这一个参数进行了反演分析，实际上数值模型的计算结果往往是由模型计算所用到的所有参数共同决定的，因此数值模型的参数反演工作还有待进一步完善，对数值模型所用到的所有参数进行综合分析与协同反演将是未来研究的发展趋势。