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

统的参数分区方法很难准确地刻画低渗裂隙介质场的非均质性，因此需要利用有效的方法优化数值模型参数，使 其计算结果更加准确、可靠$建立了烟台某地下水封洞库的低渗介质地下水渗流场数值模型，并对模型中的关键

参数一&岩体渗透系数进行了参数反演，得到了模型中两个关键层位的岩体渗透系数场$结果表明：研究区岩体

渗透系数表现出强烈的非均质性，丁烷洞库在一94〜一100 m深度上形成了强渗透带，造成丁烷洞库地下水涌水量

偏大$参数优化后的数值模型的计算精度进一步提高，能够准确地刻画地下水封洞库的运行状态$关键词：地下水封洞库；地下水渗流场;参数反演；岩体渗透系数；非均质性；数值模型中图分类号：X143；P641. 2

文献标识码：A

DOI： 10. 13578/j. cnki. issn. 1671-1556. 2019. 04. 001Parameter Inversion Study of Low Permeability Media SeepageNumerical Model for Groundwater-sealed CavernTANG Liansong,CHEN Gang, HU Cheng(School of Environmental Studies , China University of Geosciences (Wuhan ) ,Wu#a( 430074 , China}Abstract:In the groundwater numerical simulation work,the model parameters are the key factors affecting

the accuracy of numerical simulation. The traditional parameter partitioning method can hardly describe the

anisotropism of low permeability fracture medium field,so an effective parameter optimizing method is ex- pectedtoincreasetheveracityoNthenumericalmodel.ThispaperestablishesaseepagemodeloNaground-

water-sealed cavern in Yantai City,and optimizes the permeability coefficients of two crucial layers by the

parame;erinversion me;hod.Theresul;sindica;e;ha;;hepermeabili;ycoeficien;ofrock massiss;rongly anisotropic in the study area,and a high permeability zone exists in —94〜一100 m depth of Butane cavern,

which makes the water inflow of Butane cavern more than that of others. After the parameter optimization

analysis, heseepagemodelismoreaccura;ewhichcandescribe;heopera;ionalcondiionof;hegroundwa-

;er-sealedcavernprecisely.Key words： groundwater-sealed cavern； groundwater seepage field； parameter inversion； permeability coef-

6entofrock mass(an6sotrop6sm(numer6calmodel在地下水封洞库的建设和运营中，地下水压力的 分布直接影响洞库的水圭寸效果，因此准确地刻画地下

模拟的过程中，模型参数的选取是决定模型计算准确 度的重要影响因素之一。以往的数值模拟工作中，往

水渗流场对于保证地下水封洞库的安全和稳定运行 往是依据当地地质条件将研究区简单地划分为若干

至关重要1+$由于地下水封洞库工程的复杂性，在其

*均质区，分区内部为均质，分区间为非均质，以整个分 区为单个参数进行参数估计，这种方法被称为传统参

设计、施工和运营阶段,一般都需要借助数值模拟方 法对地下洞库所在水文地质单元的地下水渗流场特 征、水位动态变化及涌水量进行模拟和预测。在数值收稿日期=2019-03-15

数分区法或常规分区法②，该方法虽然可以极大地减

少模型参数的数量，但无法精确地刻画含水层的渗透

修回日期：2019-07-05基金项目：江苏省地质勘查专项资金项目(016046322)作者简介：唐连松(1993—)，男，硕士，主要研究方向为地下水数值模拟。E-mail： tlsl11406@gmail. com2安全与环境工程第26卷性分布，尤其是对于强烈非均质各向异性的裂隙岩体 类型为松散岩类孔隙水、浅层基岩网状裂隙水、深层

含水层，由于裂隙岩体含水层的渗透系数表现出强烈 的空间变异性，仅利用几个钻孔试验数据进行常规统 计，根本无法全面地反映出整个研究区含水层渗透性

基岩脉状裂隙水。研究区地下水主要接受大气降水 的入渗补给，沿地下裂隙网络向中西部低洼处径流,

排泄于第四系残坡积或冲洪积层中，或被人工开采。

的分布特征。而采用参数反演分析法对影响地下水 研究区地下水动态监测数据表明，浅层地下水随季

渗流场模拟的重要参数一&渗透系数进行优化可以 节动态变化明显，深层裂隙水长期保持相对稳定。 研究区的基本地质概况见图1$0

提高模拟的准确度。反演分析法实际上是一种参数

优化方法,通过拟合地下水水位与实际地下水水位的 接近程度来确定裂隙岩体含水层的渗透系数「3+，接近

1 000 2 000 m程度越高，说明拟合得到的渗透系数越精确。参数反

演分析法分为直接反演法和间接反演法两种「4+。直

接反演法就是将实际的地下水水位分布作为已知数 （已有观测值），将含水层渗透系数作为未知数进行求

解，当模拟得出的地下水水位分布与实际地下水水位 分布很接近时，求解完成。间接反演法是在前期水文

地质勘探及调查的基础上，综合研究区岩体结构和水 文地质特性，从总体上确定裂隙岩体含水层的相对渗 透性,并圈定渗透系数的大致范围，再以此作为已知 条件去拟合地下水水位分布，并通过比较地下水水位

的计算值和实际值，逐步修正含水层渗透系数以达到 较好的模拟效果「5+。本文利用FEFLOW软件建立了烟台某地下水封

洞库运营阶段的研究区地下水渗流场数值模型，由于 花岗岩体含水层的强烈非均质性，其模拟结果与实测 数据存在一定的误差，尤其是洞库地下水涌水量模拟

值与实测数据有较大偏差，因此，进一步以现场水力

试验数据作为先验信息,利用地下压力计的压力监测 数据及工程结构内部的流量监测数据，对建立的数值 模型参数进行了反演，在提高模型模拟精度的同时刻

画了强烈非均质含水介质的渗透系数场。1研究区概况研究区位于山东半岛北部、烟台市西侧，为某公 司地下水封洞库所在场地。该洞库工程包括1个丙 烷洞库、1个丁烷洞库％个LPG洞库，3个洞库在

平面上呈“品-字形分布。其中，位于北部的丙烷洞

库总库容为50万m3，主洞室所在高程为一120〜— 146 m（黄海高程），水幕巷道所在高程为一94〜— 100 m（黄海高程）；位于西侧的丁烷洞库和位于

东侧的LPG洞库总库容均为25万m3，主洞室所在

高程为一90〜一116 m（黄海高程），水幕巷道所在 高程为一64〜一70 m（黄海高程\"研究区含水介质主要为第四系覆盖层和中生代 燕山早期中粗粒黑云母二长花岗岩，地下水的主要

洞库

水幕巷道

冲积神負砾混變鑫炉石、含

风积相中细砂、粉砂图1研究区地质概况Fig. 1 Geology introduction of the study area2数值模型的建立2・1模型的边界和数学模型根据研究区的地质条件，将研究区概化为非均

质各向异性三维潜水流，模型上边界为入渗补给、蒸 发排泄边界；下边界（标高为一176 m）为隔水边界,东部和西部边界是沿流域分水岭划分的边界，在模

型中概化为隔水边界;南部边界部分为隔水边界，部

分为第二类边界（定流量边界），流量根据区域地下

水径流模数计算给出；北部为黄海，是研究区的最低 排泄基准面，定为第一类边界$洞室和水幕系统等 工程结构按第一类边界处理$依照上述条件建立数

学模型如下：I （\"「b2

H y\" [W = 3H一8 !d \"（x，y，;）

，t$09（8幕）+牛\"H\\ 亠\"y ）+。\"\"H 、）+ #yy H（X，y，;） # So，$0H （x，y，;,0） = H0H （x，y，;，t ） = H18 --\"n

S2 =q（，y，;，t） 第）期唐连松等：封洞库低渗介质地下水渗流场数值模型的参数反演3式中：Q为研究区域;8竝、心、8；分别为9%;主 方向的渗透系数（m/J）;：o为初始地下水水位

（m）；Hi为洞库或水幕水位（m）;S°为地下水自由

表面;Si为第一类边界，表示黄海、河流、洞库和水

幕的位置;宀为单位储水系数（1/m）;出为给水度;# 为源汇项;W为潜水面上的蒸发排泄，降雨入渗补 给,井的抽水量和泉的排泄量等;S&为第二类边界;

＜为第二类边界的水流通量（m/d） $2.2模型的网格划分和参数赋值2.2.1 模型的网格划分采用三角网格对模型进行剖分，对洞室及水幕系

统进行局部加密，在平面上将4 993.7 mX5 165 m的 研究区剖分为49 308个有限单元格，在垂向上按现

场水力试验分段（现场压水试验按10 m 一段进行）及

洞库工程结构将研究区共划分为18层，得到研究区 地下水渗流场数值模型的三维结构见图2。图2研究区地下水渗流场数值模型三维结构Fig・ 2 Three-dimensional structure of the model

in the study area2.2.2 模型的参数赋值根据研究区地质图、水文地质图，并结合现场水

文地质试验结果，对模型参数进行赋值，其中第四系 覆盖层和岩体强风化层即模型的第1层和第2层按

参数分区法赋值,其参数分区情况见图3。模型第1 层为第四系覆盖层，在第1层分区中，区域1、为低

山区域且第四系覆盖层较薄，区域3、4为残坡积物 堆积形成的第四系覆盖层，区域5、6为九曲河冲积 物形成，区域7为九曲河河道，区域8、9为水田；模 型第2层为岩体强风化层，其中区域3的风化程度

稍强于区域1、2。模型中其他各层的参数利用现场 钻孔压水试验结果的平均值进行赋值。由此得到研

究区地下水渗流场数值模型各层各分区具体参数的

初始值见表1$图3研究区地下水渗流场数值模型第1、2层参数分区图

Fig. 3 Distribution of hydrogeological parameters in layer1andlayer2ofthenumericalmodelofthe groundwater seapage field in the study area表1研究区地下水渗流场数值模型各层各分区

参数的初始值Table 1 Initial value of each layer parameter of the

numerical model of the groundwater seapagefieldinthestudyarea分层分区号水平渗透系 垂直渗透系

单位储水 数/（m- d—1）数/（m- d-1）给水度系数/m—111.555200.15552021.296000.1296000.101X10—40.10312.960001.2960001X10—4第四系覆48.640000.8640000. 201X10—4盖层（模50.728000.0728000. 201X10—40.型第1 层）62.592000.259200 201X10—4786.400008.6400000. 251X10—4817.280001.7280000. 301X10—40. 2599.750000.9750001X10—40. 251X10—4强风化层10.355600.0355605X10—（模型第20.112320.0.011232 0252 层）30.864000.0864000. 025X10—50. 055X10—5丁烷、LPG 洞库水幕 层（模型 0. 003 00 0.000 300 7X10+7第6层） 丙烷洞库 水幕层（模 0.002 05 0.000 205 7X10+7型第10层）2.3地下洞库及其周围区域地下水压力及水位的

监测

在地下洞库运营过程中，为了保证洞库的水封 条件长期有效，需要对地下洞库及其周围区域的地 下水压力进行全方位监测,因此在部分水幕孔中安

装有地下压力计。地下压力计安装位置为3个洞库 的水幕层，其中丁烷、LPG洞库水幕层（标高为

—64〜一70 m）位于模型第6层，丙烷洞库水幕层

（标高为一94〜一100 m）位于模型第10层，模型在 这两个层位的初始水文地质参数见表1$同时，将

部分勘察钻孔改造为地下水水位监测井，对地下外

4安全与环境工程第26卷围的地下水水位进行长期监测$2.4模型的模拟精度分析工程场地的地层岩性较为单一，而且工程建设会避开 导水性较强或规模较大的地质构造，因此传统的参数

利用读数较为稳定、数据较为可靠的地下压力 分区法很难应用于这种情况。为此，本文采用PEST

计及地下水水位监测井(共31个监测点位)获得的 数据对地下水渗流场数值模型计算结果进行校验#

程序对影响数值模型计算结果的主要参数一&渗透 系数进行参数反演，以提高模型模拟的精度。PEST参数优化算法是一种较为成熟、应用较

其结果见图4$模拟地下水水位/m50「 40 为广泛的参数优化方法$最初PEST是由Doher­

ty™ 发出的一种通用参数优化算法，之后Tokm

等⑺提出了一种结合Tikhonov正则化及T-SVD

30 20 10■-------------------------------‘ 丿#-50-40 -30 -20> -10,•苗0

10 20 30 40 50实测地下水水位/m护 -20-30 -40 -50图4地下水水位模拟值与实测值的拟合情况Fig. 4 Fitting of the simulated and measured

groundwater level由图4可见,受研究区花岗岩含水层强烈非均质

性的影响，数值模型的模拟精度一般，地下水水位31 个监测点位模拟值与实测值的残差平方和为264. 86, 地下水水位模拟值与实测值的最大误差为7. 4 m$地下洞库地下水涌水量模拟值与实测值的对

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

(单位:mD/d)Tab:e2 Comparisonoftheamountofsimu:atedwaterinflow and measured water inflow(unit： m3 /d)洞库涌水量模拟值涌水量实测值丙烷洞库227336丁烷洞库441696LPG洞库360432由表2可知，地下洞库地下水涌水量的模拟值

与实测值也存在一定的误差，尤其是丁烷洞库地下 水涌水量的模拟值与实测值的误差较大，说明丁烷

洞库处可能存在渗透系数较大的强渗透带，这种情

况下需要对模型参数进行进一步的优化以提高模型 模拟的精度$3模型参数反演3.1 PEST参数反演方法由于该地下洞库修建于较为完整的花岗岩体中，

正则化方法优点的混合正则化方法(SVD-Assist),

进一步优化了 PEST算法，使其更加完善$ PEST

算法在国内外地下水数值模型参数优化工作中取得 了较好的效果，董艳辉等*〕利用并行化的PEST算 法对大尺度的地下水数值模型进行了参数优化，得 到了较为理想的结果；王礼恒等凹利用假想模型进

一步验证了 PEST算法在优化地下水流动模型参数 上的可靠性；姜蓓蕾等*10+对PEST算法所采用的向 导点-正则化方法进行了深入的研究，证明了 PEST

是一种有效的参数反演方法，同时证明该方法中并 不是向导点越多反演精度越高，向导点过多反而会 影响计算精度且增加计算量$PEST 程序采用的是 Gauss-Marquardt-Leven-

berg非线性参数估计算法，该算法是一种改进的 Gauss-Newton法*11+$该算法在每一次迭代过程

中，计算模型结果对各待估参数的导数，利用这些导 数将非线性模型线性化，并更新一次参数向量，计算

新的目标函数后，再进入下一次迭代过程，通过反复 迭代直至目标函数达到最小值$实际工作中，大部分模型都是非线性的，假设模 型的n个初始参数值储存在向量X。中，模型的m个

实际观测值储存在向量J0中，则X0与J0的关系为\"0==(X。)

(1)根据泰勒公式，任意一组参数向量X与其对应

的观测向量y之间的关系可以表示为y= M (x) = M (x0 ) + =' (x0 ) (x — x0 ) +

='(x°) M\" (x°) ()2! (x —X0)十…十n! (x —X0)”当x)X0时，忽略高阶无穷小项，则式(2)可简 化为

\"=\"。+丿(X —X0)

(3)式中J为m行n列的偏导数雅克比矩阵：Ji，+=^MxT (x

)

⑷(、由此可得出非线性模型参数反演过程中的目标 函数为第4期唐连松等：封洞库低渗介质地下水渗流场数值模型的参数反演5$= （y —y° — J（X —X0））TQ（y —y° — J（X —X°））（5）

式中：Q为m行n列的实际观测值的权重矩阵，参 数优化过程就是不断迭代更新参数向量X使得目标

函数①最小化的过程「12+。3.2参数反演范围由于前期勘察工作及洞库运营阶段的监测工作 都是围绕着洞库所在区域开展的，缺少洞库工程外 围的实测资料，因此参数反演的范围为以洞库区域

为中心的1 000 mX1 000 m正方形区域（见图5） $ 地下水水位监测点和岩体渗透性原位试验点主要集图5模型参数反演范围Fig. 5 Range of parameter inversion中在丁烷、LPG洞库水幕层（一64〜—70 m）和丙烷

洞库水幕层（一94〜— 100 m）这两个关键层位，这 是由于在洞库运营过程中，水幕巷道及水幕孔的持 续供水对保证洞库水封效果起着至关重要的作用,

而且在深部修建这种特殊的工程结构也会使这两个 层位的地下水动态变得更加复杂，对这两个层位的 模型参数进行优化是客观评价洞库水封效果、提高

水幕系统工作效率的必要工作，因此本次以岩体渗

透系数原位试验点和地下水水位监测点的监测数据 对丁烷、LPG洞库水幕层和丙烷洞库水幕层这两个

关键层位的岩体渗透系数（8）进行参数反演，丁烷、

LPG洞库水幕层和丙烷洞库水幕层的原位试验点、

监测点及水位分布情况见图6和图7。3.3参数反演结果和模型优化结果运行PEST程序对丁烷、LPG洞库和丙烷洞库 水幕层的岩体渗透系数进行参数反演后，得到洞库

岩体渗透系数的反演结果（以lgK表示）以及地下 水水位模拟值与实测值的拟合结果见图8和图9。

丁烷、LPG洞库地下水水位模型模拟值与实测值的0 100 200 300 400 500 600 700 800 900 1 000X/m上监测点 •原位试验点一水幕巷道 等水位线图6 丁烷、LPG洞库水幕层的原位试验点、监测点

及水位分布Fig.6 In-situtestpoints，monitoringpointsandground-

water:eve:distributionin:ayerofwatercurtains ofButanecavernandLPGcavern0 100 200 300 400 500 600 700 800 900 1 000X/m上监测点•蠶试-LPG、丁烷洞库—霹黯 談图7丙烷洞库水幕层的原位试验点、监测点及水位分布 Fig.7 In-situtestpoints，monitoringpointsandground-

water:eve:distributionin:ayerofwatercurtain

ofPropanecavern残差平方和为39. 82,丙烷洞库地下水水位模拟值 与实测值的残差平方和为19. 15$由图8和图9可见，经参数反演后的丁烷、LPG 洞库水幕层的岩体渗透系数最大值为0.006 1 m/d,6安全与环境工程第26卷模拟地下水水位/m(a)岩体渗透系数反演结果(b)地下水水位模拟值与实测值拟合结果图8 丁烷丄PG洞库水幕层岩体渗透系数反演结果（以lgK表示）以及地下水水位模拟值与实测值的拟合结果Fig. 8 Parameter inversion result of layer of water curtains of Butane cavern and LPG cavern and the fittingresult of the groundwater level simulated values and the meosured values模拟地下水水位/m1.01.5 2.02.53.03.54.0X/m(a)岩体渗透系数反演结果(b)地下水水位模拟值与实测值拟合结果图9丙烷洞库水幕层岩体渗透系数反演结果（以lgK表示）以及地下水水位模拟值与实测值的拟合结果Fig.9 Parame;erinversionresul oflayerof;hewa;ercur;ainofPropanecavernand;hefiingresul;of;he

groundwa;erlevelsimula;edvaluesand;hemeasuredvalues最小值为0. 000 03 m/d,丙烷洞库水幕层的岩体渗透 系数最大值0. 059 m/d,最小值为0. 000 105 m/d,其

层两个层位岩体渗透系数的反演结果与现场实测数

据具有相似的统计特征，说明其反演结果能够较为 客观地反映研究区的实际情况$根据概率密度函数 可知，丁烷、LPG洞库水幕层岩体的渗透系数以

0.000 5 m/d为主，丙烷洞库水幕层的渗透系数以

渗透系数在空间上表现出强烈的非均质性。假设岩体渗透系数服从对数正态分布，将其反

演结果绘制成岩体渗透系数概率密度曲线，同时将

现场16个钻孔中对应层位的压水试验数据与岩体 渗透系数概率密度曲线进行对比，其结果见图10。0. 003 9 m/d 为主。利用参数优化后的数值模型计算得到的各洞库 地下水涌水量见表3。由图10可见，丁烷、LPG洞库和丙烷洞库水幕

第4期唐连松等：封洞库低渗介质地下水渗流场数值模型的参数反演7方法来刻画地下水封洞库特定层位的岩体渗透系数，

经过参数优化后,数值模型地下洞库拟合精度得到提 高,模拟结果表明：丁烷、LPG洞库水幕层的岩体渗透

系数以0.000 5 m/d为主，渗透性最强处岩体渗透系 数为0. 006 1 m/d,但这一层位上没有形成明显的强 渗透带;丙烷洞库水幕层的岩体渗透系数以0.003 9

m/d为主，渗透性最强渗透系数为0. 059 m/d,这一

层位上出现了明显的强渗透带,这也解释了为什么丁

-5.5 -4.5 -3.5 -2.5

-1.5 -0.5

0.5 lgK烷洞库从开始运行起就出现了地下水涌水量偏大的 (b)丙烷洞库水幕层图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参数优化后洞库地下水涌水量模拟值与实测值的对比(单位:m3/d)Table3 CompaHison of the amount of simulated wateH

inflowand measuHed wateHinflow afteH

parameter optimizing洞库涌水量模拟值涌水量实测值丙烷洞库339336丁烷洞库LPG 洞库670696401432由表3可知，与初始数值模型的计算结果(见表2)相比较，参数优化后的数值模型计算结果与实测

数据更加相符，说明经过参数优化后的数值模型具 有更高的拟合精度。4结论与展望结合研究区水文地质条件分析和现场水文地质

试验数据,建立了烟台某地下水封洞库运营阶段研究 区地下水渗流场数值模型，通过对比监测数据发现该

均质数值模型的计算结果与实测数据有较大偏差，本 研究进一步采用参数反演与现场试验数据相结合的

现象。由此说明，经参数反演后的地下水渗流场数值 模型能够更加客观、准确地评价水封洞库的运行状

况。由于工程实际情况较为复杂，本文只对岩体渗透

系数这一个参数进行了反演分析，实际上数值模型的

计算结果往往是由模型计算所用到的所有参数共同 决定的，因此数值模型的参数反演工作还有待进一步

完善,对数值模型所用到的所有参数进行综合分析与 协同反演将是未来研究的发展趋势。参考文献：*1+李术才，张立，马秀媛，等.大型地下水封石油洞库渗流场时空演

化特征研究[J+.岩土力学,2013,34(7) ：1979-1986.[2+刘玉珍，王本德，姜英震.基于可变水文地质参数的地下水系统

数学模型[J+.水科学进展，2009,20(3) ：398-402.[3+刘迎曦，李守巨，李正国，等.岩体渗透系数反演的数值方法及其

适定性*+•辽宁工程技术大学学报,2000,19(4) ：375-378,*4+纪成亮，李晓昭，王驹，等.裂隙岩体渗透系数确定方法研究[J + .

工程地质学报，2010,18(2) ：235-240.*5+孙蓉琳.玄武岩渗透系数尺度效应及顺序指示模拟[D+.武汉：中

国地质大学(武汉),2006.[6+ Doherty J. Ground water model calibration using pilot points and

regularization]J+. Ground Water ,2003,41(2) : 170-177.*7+ Tokin M,Doherty J. A hybrid regularized inversion methodology

for highly parameterized environmental models [J]. Water Re­

sources Research , 2005 , 41, W10412. doi: 10. 1029/2005WR00

3995[8+董艳辉，李国敏，郭永海，等.应用并行PEST算法优化地下水模

型参数*+•工程地质学报,010,18(1):140-144.*+王礼恒，董艳辉，李国敏，等.基于PEST的地下水数值模拟参数

优化的应用[J+.工程勘察,2014,42(3):38-42.*10]姜蓓蕾，冯燕生，施小清，等.向导点-正则化方法在水文地质参

数反演中的应用[J + .南京大学学报(自然科学)，2016,52 (3): 438-447*11]Maneta M P, Wallender W W. Pilot-point based multi-objective

calibrationinasurface-subsurfacedistributedhydrologicalmodel

[J]. Hydrological Sciences Journal,2013,58(2) ：390-406.*12]Skahill B E, Baggett J S, Frankenstein S, et al. More efficient

PEST compatible model independent model calibrationCJ]. Eni-

ronmental Modelling and Software ,2008,24(4) ：436-451.通讯作者：陈 刚(1967—)男，博士，副教授，主要从事工程水文地

质等方面 的研究 $ E-mail: chengang@ cug. edu. cn