地下水污染强度及渗透系数场的反演识别研究.pdf

在制定地下水污染修复方案时，污染源参数和渗透系数场是最重要的地下水数值模型参数，但前人研究多集中于单一类型参数的识别。文章中采用地下水污染物运移模型（MT3DMS）和数据同化方法（迭代局部更新集合平滑器，ILUES）构成地下水污染源识别的求解框架，并利用 Karhunen-Love 展开技术实现渗透系数场的参数降维，最后通过同化水头与浓度数据实现地下水污染源强和渗透系数场的联合反演。结果表明：（1）ILUES 算法能精确识别污染源参数和渗透系数场，并且具有很高的普适性；（2）精确表征渗透系数在空间上呈现出的非均质性，是预测污染物迁移路径、反演污染强度的关键；（3）ILUES 算法参数影响着反演效果，综合考虑计算效率和计算精度等，可以得到算例的最佳样本集合大小（Ne=4 000）和ILUES 算法最佳参数组合（局部临近样本集合占比=0.4，相对权重 b=4）。但在实际工程案例中，如果对精度的要求不是过高，经验组合（=0.1，b=1）更值得推荐。研究结果对于区域地下水资源调查、评价和管理等工作具有较强的实践意义，并可为后期地下水污染预测及地下水监测井网优化提供技术支撑。关键词：地下水污染；参数反演；数据同化；集合平滑器中图分类号：X523 文献标志码：A 文章编号：1000-3665（2023）04-0193-11Identification of groundwater pollution intensity and hydraulicconductivity fieldWU Yanhao1，JIANG Simin1，WU Zijun2（1.Department of Hydraulic Engineering,School of Civil Engineering,Tongji University,Shanghai200092,China；2.State Key Laboratory of Marine Geology,(Tongji University),Shanghai200092,China）Abstract：The pollution source parameters and hydraulic conductivity field are the most important parameters ofgroundwater numerical models when making groundwater pollution remediation plans.However,previous studiesfocused mainly on the identification of single type parameters.The groundwater pollutant transport model(MT3DMS)and data assimilation method(iterative local updating ensemble smoother,ILUES)are used to form asolution framework for groundwater pollution source identification,and Karhunen-Love expansion technique isused to realize parameter dimension reduction of the hydraulic conductivity field.The joint inversion ofgroundwater pollution source intensity and hydraulic conductivity field are also realized by assimilating hydraulicheads and concentration data.The results show that(1)the ILUES algorithm can accurately identify pollution 收稿日期：2022-08-17；修订日期：2022-09-24投稿网址：基金项目：国家自然科学基金项目（41976057；41676061）第一作者：吴延浩（1999-），男，硕士研究生，主要研究方向为污染水文地质学。E-mail：通讯作者：江思珉（1980-），男，博士，副教授，主要从事地下水数值模拟等方面研究。E-mail： 第 50 卷 第 4 期水文地质工程地质Vol.50 No.42023 年 7 月HYDROGEOLOGY&ENGINEERING GEOLOGYJul.，2023source parameters and permeability coefficient field,and it is of high universality.(2)Accurate characterization ofspatial heterogeneity of the coefficient of permeability is the key to predict pollutant migration path and inversionof pollution intensity.(3)The ILUES algorithm parameters affect the inversion results.By considering thecomputational efficiency and accuracy,the optimal sample set size (Ne=4 000)and the optimal parametercombination of ILUES algorithm(=0.4,b=4)can be obtained.However,in practical engineering cases,theempirical combination(=0.1,b=1)is more recommendable if the requirement for accuracy is not too high.Theresults of this study have strong practical significance for regional groundwater resources investigation,evaluationand management,and can provide technical support for later groundwater pollution prediction and optimization ofgroundwater monitoring well networks.Keywords：groundwater pollution；parameter inversion；data assimilation；ensemble smoother 地下水污染的隐蔽性对污染修复带来了巨大挑战1 2。为了应对这一挑战，地下水污染源识别（gro-undwater contaminant source identification，GCSI）受到了越来越多的关注3 4。GCSI 能够通过有限的地下水实际观测数据（包括水头和污染物浓度）和其他先验信息来获取地下水污染源的有关信息5。此外，地下介质以及污染物自身化学反应的空间异质性也给 GCSI模型参数带来了新的挑战6。GCSI 模型对含水层性质（如孔隙度、渗透系数等）及生化参数十分敏感，所以在地下水污染的实际应用中需要对渗透系数等空间异质性参数进行准确估计7 8。因此地下水污染源源强和渗透系数的同步识别在 GCSI 中是亟需解决的问题。目前，常用的地下水污染溯源及参数识别方法分为以下 2 种。第 1 种是模拟-优化（S-O）方法9，包括遗传算法10 11、模拟退火算法12 13和粒子群优化算法14 15。然而，对于观测数据点较少、未知参数较多的情况，问题的不确定性大大增加，S-O 方法效率可能较低16。第 2 种方法则是数据同化（DA）方法，包括集合卡尔曼滤波器（EnKF）和集合平滑器（ES）算法，它不仅可以考虑未知参数的不确定性，还考虑了模型结构和模型输入的不确定性17。如 Xue 等18将 EnKF 集合到贝叶斯模型平均框架，实现了地下水流模型边界和渗透系数场不确定条件下的污染源识别；Li 等19在基于正态得分集合卡尔曼滤波（NS-EnKF）的基础上提出基于正态得分集合平滑器（NS-ES），高效反演出同精度的多峰非高斯渗透系数场。EnKF 虽被广泛应用，但其对模型参数和状态是同时进行更新的，在诸如地下水模型这样非线性的模型中，就会引起参数和状态的不一致，并且 EnKF 需要频繁的同化数据。相比较而言，ES 的批处理方法只更新参数，不更新状态，并且一次性使用所有观测数据进行同化，所以ES 可能比 EnKF 更适用于地下水模型的参数反演。为了进一步提高 ES 对求解非线性问题的能力，许多学者对该方法进行了改进。Zhang 等20提出了一种迭代局部更新集合平滑器（ILUES），通过对每个样本使用 ES 进行局部更新并进行多次同化，提高了 ES 求解性能的同时，也验证了其能量化复杂水文模型的参数不确定性；贺楠楠21在 ILUES 的基础上提出多残差深度卷积神经网络耦合迭代局部更新集合平滑器反演框架（MRDCN-ILUES）反演框架，以更低的计算代价，实现了地下水中铬污染源的泄漏速率反演。基于以上分析，本研究中采用 ILUES 算法作为反演框架来估计污染源特征和渗透系数场。在以往的研究中，大多集中于污染源源强或渗透系数等单一类型参数的识别，而对二者的同步反演则少有研究。本文首先将渗透系数视为已知参数，说明了渗透系数场的准确刻画在污染物信息识别的重要性，在此基础上利用 ILUES 算法实现了地下水污染源源强和渗透系数的同步反演，并确认了最佳参数组合。1 研究方法 1.1 地下水流与溶质运移模型（1）地下水流动方程根据 Darcy 定律以及质量守恒定律，在不考虑水的密度变化条件下，三维空间孔隙介质中地下水流动的偏微分方程：x(Kxhx)+y(Kyhy)+z(Kzhz)+W=SSht（1）式中：Kx、Ky、Kz分别为渗透系数在 x、y、z 方向上 的分量/（md1）；W源汇项，用以表示流进汇或来自源的水量/d1；194 水文地质工程地质第 4 期h含水层水位/m；SS含水介质的贮水率/m-1；t时间/d。（2）污染物迁移的数学表达式xi(DijCpxj)xi(viCp)+qsCps+Rn=(Cp)t（2）式中：有效孔隙度；Dij水动力弥散系数张量/（m2d1）；Cp污染质 p 的质量浓度/（kgm3）；vii 分量的渗流速度/（md1）；Cps源汇污染质 p 的质量浓度/（kgm3）；qS源汇单位面积上通量/d1；Rn化学反应项/（kgm3d1）；xi、xj距离在相应坐标方向上的分量/m。在建立地下水流三维运动方程和溶质运移模型的基础上，分别采用 MODFLOW 程序22求解地下水流模型，用 MT3DMS 程序23求解溶质运移方程。1.2 ILUES 算法相比于 ES，由 Zhang 等20提出的 ILUES，可通过利用扰动后的观测误差的协方差矩阵对观测数据进行多次同化，解决了当参数的后验分布不服从高斯分布时，ES 方法不适用的问题，而且提高了反演高维参数的计算效率。ILUES 中采用的协方差矩阵扰动策略来源于 Emerick 等24所提出的 ES-MDA 算法。扰动协方差矩阵 CD通常由下式得到：CD=NiterCD（3）式中：CD观测误差的原始协方差矩阵；NiterILUES 算法的迭代次数。mfj(j=1,Ne)传统的集合平滑器方法采用集合内所有参数样本来更新样本，为了更好的解决高维、强非线性反演问题，ILUES 方法与传统集合平滑方法不同，采用了一种局部更新策略，仅使用局部样本来更新整体参数样本25。具体来说，对于地下水模型中的参数 m，不同于 ES 直接更新参数，ILUES 是先更新参数的 Ne 个样本的局部集合：J(m)=J1(m)Jmax1+bJ2(m)Jmax2（4）J1(m)=f(m)dTC1Df(m)d（5）J2(m)=mmfmTC1MMmmfm（6）式中：J1(m)用来度量模型预测值 f(m)与观测值 d 之间的距离；mfmJ2(m)用来度量模型参数 m 与对应样本之 间的距离；b表示两个距离度量 J1(m)和 J2(m)在综合度量标准中的相对权重，b(0,)；Jmax1Jmax2、分别表示 J1(m)和 J2(m)的最大取值；CMM维度为 NmNm的模型参数集合 M 的自 协方差矩阵。mloc,fj具体在更新局部样本集合时，ILUES 首先选择J(m)值最低的样本 Nloc=Ne，（0,1，作为样本 m 的局部临近样本集合，再用集合平滑器更新每个样本的局部临近样本集合：maj,k=mfj,k+Cloc,fMD(Cloc,fDD+CD)1dk f(mfj,k)（7）maj,k式中：更新后的局部集合；mfj,k更新前的局部集合；k选取的 J(m)最低的样本个数;Nloc局部临近集合中样本的个数；mloc,fjm局部临近样本集合占全局样本集合 的比例；Cloc,fMDmloc,fjDloc,fj=f(mfj，1),f(mfj，2),f(mfj,Nloc)维 度 为 NmNd的 模 型 参 数 局 部 临 近 样本集合与模型预测值样本集合 的 互 协方差矩阵；Cloc,fDDDloc,fj维度为 NdNd的的自协方差矩阵；dk表示在观测值 d 基础上引入随机观测误差扰动的样本；f(mfj,k)预测值。Ma=mloc,a1,mloc,a2,mloc,aNe最终得到经过更新的全部样本集合 2 算例研究本研究中，在稳态地下水流条件下，对流和弥散是二维污染物迁移系统中的主导过程。假设含水层为 200 m150 m 的非均匀各向异性承压含水层，具有二维稳态地下水流，初始模型网格尺寸为 30 行 40 列，单元尺寸为 5 m。南北边界为二类隔水边界，西部和东部边界是 20 m 和 15 m 的恒定水头分边界，含水层厚度为 1 m（图 1）。含水层参数值见表 1。lnK2lnK初始时刻，研究区无污染物。含水层渗透系数（K）满足对数正态分布，其均值（）与方差（）分别为 3.0 和 0.5，沿 x 和 y 方向的相关长度均为 50 m，变差函数为指数型，真实渗透系数场见图 2（a）。模拟总时长为 800 d，共 8 个应力期（每个应力期 100 d）。为了验证 ILUES 算法反演的有效性，本算例中考虑 3 个2023 年吴延浩，等：地下水污染强度及渗透系数场的反演识别研究 195 SjSi污染源，位置分别为 S1、S2和 S3。其中 S1在应力期内无规律性释放；S2经过人工修复后，释放强度不断减小；而 S3为新污染源，释放强度持续增大。污染源参数由 24 个参数进行量化（3 个污染源在 8 个应力期的污染释放强度），表示为，其中 i=1,2,3，代表第 i 个污染源，j=1,2,8，代表第 j 个应力期，8 个应力期所指示的时间为 100（j-1）,100 j。假设表征污染源释放强度的 24 个参数的先验分布是均匀分布，其先验区间和真实值见表 2。表 2 污染源释放强度参数的真实值及先验区间Table 2 True values and a priori interval of the pollution sourcerelease intensity parameter 应力期先验区间污染源源强/（kgd1）S1S2S3SP120,10098.054.030.0SP220,10087.078.039.0SP320,10074.064.046.0SP420,10061.085.057.0SP520,10047.071.066.0SP620,10039.061.071.0SP720,10031.045.079.0SP820,10024.089.083.0 为了估计上述未知参数，在研究区域布设 16 个观测井（图 1），并由 MODFLOW 和 MT3DMS 正演得到水头和浓度的模拟值，由正态分布的测量噪声计算得到观测值。C=CS+NCS（8）式中：C模拟值；CS观测值；N观测噪声（均值为 0，标准差分别为 0.005 m和 0.005 kg/m3)的高斯分布）；满足标准高斯分布的随机偏差。lnK/（md1）lnK/（md1）5.04.54.03.53.02.52.01.51.05.04.54.03.53.02.52.01.51.00200200长度/m长度/m宽度/m宽度/m（a）真实渗透系数场（b）克里格插值场1501500图 2 含水层渗透系数场Fig.2 Hydraulic conductivity field of the aquifer 表 1 含水层的水文地质特征Table 1 Hydrogeological features of aquifers 参数名称单位参数值含水层厚度m1.0孔隙度无量纲0.3纵向弥散度m2.0横向弥散度m0.6 污染源位置观测井位置定水头边界隔水边界15002001.01.52.02.53.03.54.04.55.0隔水边界长度/m宽度/m定水头边界lnK/（md1）S1S2S3图 1 含水层结构平面图（S1、S2和 S3为 3 个污染源位置）Fig.1 Planar structure of the aquifer(S1,S2 and S3 are the threecontamination source locations)196 水文地质工程地质第 4 期3 结果 3.1 算例场景设置在影响地下水系统的所有不确定因素中，渗透系数是表征含水层特性的一个关键要素，其变异性主要是由自然界含水介质的非均质各向异性、取样和测试的失真以及试验量测误差造成的，而渗透系数是影响污染物在地下水中运移产生不确定性的重要原因。因此，渗透系数的准确获取对于解决污染物反演问题至关重要。在实际场地中，对于含水层的渗透系数场分布的认知往往是零散的，即仅有从井中测量得到的该点处的渗透系数值。为了研究 ILUES 算法对地下水污染强度识别的准确性，对渗透系数场分别采用 3种方式进行处理（对应于 3 种场景），观测数据都是真实渗透系数场条件下获取的。lnK=31）场景 1。渗透系数场设置为恒定值是最常用的简化方法。实际上，渗透系数是通过拟合试验结果得到的，但是样本数量有限，无法覆盖整个场地。一般情况下，利用有限的样本得到的渗透系数值的平均值便成为场地渗透系数的代表性定值。本场景案例的渗透系数场取真实渗透系数场的均值（）。2）场景 2。采用 Kriging 方法对渗透系数进行插值，得到具有空间变异性的渗透系数场。受技术和经济条件制约，现实中获取的渗透系数数据往往较少，不能充分反映整个含水层渗透系数空间分布特征。因此，大多数情况下需要根据有限的数据对渗透系数进行空间插值。普通克里格法属线性平稳地质统计学范畴，是地质统计学中最基本的方法，也是最常用的方法之一。设区域化变量满足二阶平稳假设，其数学期望为未知常数，协方差函数和变异函数存在且平稳，以任一未知点为中心点的某一邻域范围内的 n 个样本点，赋予相应权重计算未知点的属性值，形成估计领域内 n 个信息值的线性组合，该方法即为普通克里格方法（ordinary kriging，OK）。本场景中，已知研究区观测位置的 K 值，利用普通克里格插值方法得到整个研究区的渗透系数场，如图 2（b）所示。3）场景 3。引入数据降维的方法表征具有空间变异性的渗透系数场。渗透系数场的维度与模型的网格剖分程度有关（本算例为 3040=1 200 维），通过降维，可以寻找到高维数据的低维子空间，并保留原始数据的大部分变异性。Karhunen-Love 展开（以下简称 KL 展开）具有离散效率高、得到的随机变量数较少等优点,是一种较为理想的降维方法，因此在这里采用 KL 展开对地下水模型的高维渗透系数场进行降维：lnK(x,y)=lnK(x,y)+NKLi=1iiSi(x,y)（9）式中：i独立标准高斯随机变量；i相关函数的特征值；Si(x,y)相关函数的特征函数；NKLKL 展开项的数目。本算例中选择 NKL=200，可以保留渗透系数场约 95.20%的变 异性。3.2 不同场景下源强的反演结果S1S1、S4S1、S7S1S1S2、S4S2、S7S2S1S3、S4S3、S7S3在上述 3 种场景中，分别获得整个研究区的渗透系数场，然后写入 MODFLOW 模型中的 BCF（包括其中 LPF、HUF 等）模块并利用 ILUES 算法反演污染强度，因此本节中，渗透系数场不会被视为一个不确定的模型参数。为了对参数不确定性进行充分估计，ILUES 算法参数集合大小 Ne=500；为了充分考虑解的收敛性，迭代次数 Niter=10。同时根据张江江26的建议，局部集合参数=0.1，距离权重参数 b=1。图 3 分别是 3 种情境下，随着迭代次数的增加，3 个污染源释放强度反演结果的变化趋势图（从表征污染源信息的24 个未知参数分别选取了 S1的，S2的和 S3的）。S1S1、S4S1、S7S1S7S3S1S1S4S3在场景 1 中，通过 ILUES 算法得到准确的污染源特征识别结果，得到的污染源特征参数反演结果偏离于参考值较多，较为明显的参数有和，如图 3（a）所示，主要原因在于在该情境下，仅利用恒定的渗透系数值来描述整个场地渗透系数的分布情况，无法充分表征污染物运移的不确定性。在场景2 中，9 个污染源参数的识别结果，如图 3（b）所示，明显优于场景 1 的结果；仅在个别参数（如和）与真实值之间出现明显偏差，其余 7 个参数识别结果接近于真实值，但是整体识别的准确性和集合样本的稳定性略逊于场景 3 的识别结果，如图 3（c）所示。可能的原因是在场景 2 中，利用若干已知渗透系数（30 个观测点）通过 Kriging 插值方法得到场地渗透系数场的分布情况，在一定程度上考虑了污染物运移的不确定性，但作用是有限的。综上所述，利用 ILUES 作为反演算法且选用相同的算法参数，从准确性和稳定性角度看，场景 3 的污染源强度识别结果都优于其他 2 种场景。结果表明：ILUES 是能够解决此类反演问题的有效算法；精确表2023 年吴延浩，等：地下水污染强度及渗透系数场的反演识别研究 197 605550S1s24540908580S4s27570504540S7s235301009896S1s1929490706560S4s1555040383634S7s13230403530S1s3252060585654S4s35250908580S7s37570605550S1s24540908580S4s27570504540S7s235301009896S1s1929490706560S4s1555040383634S7s13230403530S1s3252060585654S4s35250908580S7s37570605550S1s24540908580S4s27570504540S7s235301009896S1s1929490706560S4s1555040383634S7s13230403530S1s3252060585654S4s35250908580S7s37570先验样本1次迭代2次迭代4次迭代6次迭代8次迭代10次迭代3次迭代5次迭代7次迭代9次迭代真值（a）场景 1（b）场景 2（c）场景 3图 3 污染源释放强度反演图Fig.3 Inverted intensity of pollution source release 198 水文地质工程地质第 4 期征渗透系数在空间上呈现出的非均质性，是预测污染物迁移路径、反演污染强度的关键，在此基础上才能制定有效的场地修复措施。4 讨论 4.1 集合数量对渗透系数场同步反演的影响在对空间上高变异性的渗透系数场进行同步反演时，对于不同于污染强度的单一反演，将渗透系数同样视为未知参数，同时更新污染源源强和渗透系数两个参数，因此面临更高维地下水逆问题的求解。场景 3 中的降维方法可以有效减少参数量，在后续反演中只需反演降维后的参数，即可估计整个渗透系数场，显著降低了计算压力。因此本节在场景 3 的基础上，反演 24 个污染源参数和 200 个含水层参数。由于未知参数较多，往往会选择更大集合样本数（Ne）进行准确估计，但更大的 Ne 值计算时会更多的调用原始模型，显著增加运算时间和结果存储空间。为了平衡计算效率与计算精度，选用不同的 Ne 值（Ne=1 000，2000，3 000，4 000，5 000）并比较其在 ILUES 算法下反演的结果。由图 4 可知迭代次数 Niter=8 时基本收敛，因此为了减短计算时间，本节中取 Niter=8，其他参数不变（=0.1 和 b=1）。由于 Ne=500 时 24 个污染源参数的反演已经比较精确，因此下面仅讨论 200 个含水层参数（渗透系数场）的反演精度。（a）真实渗透系数场（b）反演场（Ne=1 000）（c）反演场（Ne=2 000）（d）反演场（Ne=3 000）（e）反演场（Ne=4 000）（f）反演场（Ne=5 000）200150000000150150150150150200长度/m长度/m长度/m长度/m长度/m长度/m宽度/m宽度/m宽度/m宽度/m宽度/m宽度/m2002002002005.04.54.03.52.53.02.01.51.05.04.54.03.52.53.02.01.51.05.04.54.03.52.53.02.01.51.05.04.54.03.52.53.02.01.51.05.04.54.03.52.53.02.01.51.05.04.54.03.52.53.02.01.51.0lnK/（md1）lnK/（md1）lnK/（md1）lnK/（md1）lnK/（md1）lnK/（md1）图 4 不同 Ne 值下渗透系数场的估计结果Fig.4 Estimation of the hydraulic conductivity field at different Ne values 图 4 为不同 Ne 值下的渗透系数场的估计结果，由图可以看出，Ne 值从 1 000 增加至 4 000 的过程中，ILUES 算法估计的渗透系数场与参考场的相似度逐渐上升，见图 4（a）（d）。但当 Ne 值从 4 000 继续增2023 年吴延浩，等：地下水污染强度及渗透系数场的反演识别研究 199 加到 5 000 时，渗透系数场的估计水平没有明显的提升，见图 4（e）（f），可能原因是 Ne=4 000 和 5 000 经8 次迭代，渗透系数场估计结果的不确定性都达到很低的水平。进一步定量评价 ILUES 方法的参数估计结果，引入均方根误差（RMSE）和平均集合离散度（AES）作为定量评价指标。其中，RMSE 用来度量待估计参数的估计值和真实值之间的匹配程度，AES 用来度量参数估计值的不确定性或置信度，2 者值越低代表参数反演效果越好。RMSE=vt1NNi=1(SiMi)2（10）AES=vt1NNi=11NeNej=1(SijbSi)2（11）式中：Si 第 i 个参数的估计值;Mi第 i 个参数的真实值;Sij第 i 个参数在第 j 个集合中的估计值;bSi第 i 个参数集合的均值，N 为参数的总个数。图 5 为不同 Ne 值下渗透系数场反演结果的评价指标 RMSE 和 AES 趋势图。可以看出，Ne 由 1 000 到4 000 时，AES 值在这个过程中由 0.012 5 下降至 0.009，降幅为 28%，RMSE 值由 0.549 6 下降至 0.211 5，降幅为 61.5%，说明用 ILUES 方法估算出的渗透系数场与真实场的相似性逐步增大。但当 Ne 由 3 000 到 4 000的过程中，AES 值和 RMSE 值的变化幅度很小，降幅分别为 9.1%和 4.4%，证实了 Ne 值到达 4 000 及以上时候，反演结果的不确定性已到达一个较低的水平。为了使运算的效率和准确性达到均衡，认为本场景中 Ne=4 000 为最佳集合样本大小。0.600.550.500.450.012 50.549 60.426 50.338 60.0130.0120.0110.0100.0090.0080.211 50.009 00.008 61 0002 0003 000集合样本数/个4 0005 0000.009 80.010 80.192 20.40RMSERMSEAESAES0.350.300.250.200.15图 5 不同 Ne 下的 RMSE 和 AES 的变化趋势图Fig.5 Trends in RMSE and AES under different Ne values 4.2 最优参数值的选取在上述模拟中，根据经验选择 ILUES 算法的系数mfj=0.1，b=1（经验组合），因为在前人的研究中它具有最好的综合效果。但在实际的情况中，由于场地的不同，导致最佳模拟结果的参数组合往往并不是经验组合。为了讨论本算例中最优参数组合，本节基于上节的场景，固定 Ne=4 000 并设置=（0.01，0.1，0.2，0.3，0.4，0.5）和 b=（0.01，0.1，1，2，3，4，5），组合后进行 67=42 次同步反演，并计算渗透系数场估计结果的RMSE 值，结果如图 6 所示。由图可知，当 0.1 时，若 b 取值较小（例如 0.01），ILUES 算法可以得到相对较好的拟合结果，但当选用 0.1 时，往往需要相对较大的 b 值（例如 b1），才能得到更好的估计结果，这与张江江26得出的结论一致，因为当 大于等于0.1 时，局部集合的样本数量也会比较大，因而更可能保留接近观测值的样本，这时保留和距离更近的样本，即选择较大的 b 值，会更有利于参数反演27。特别地，最佳参数组合（=0.4，b=4）也体现了上述规律。b120.10.010.373 40.546 60.592 70.573 80.699 10.598 30.667 50.598 90.554 20.211 60.714 10.493 20.513 10.643 90.489 10.414 20.437 70.209 50.413 00.376 70.365 80.578 50.679 20.598 30.313 70.404 50.298 80.358 70.678 90.587 70.697 10.483 30.497 70.168 30.374 40.714 3 1.021 0.774 30.693 50.548 20.617 70.664 13451.00.90.80.70.60.50.40.30.20.010.100.200.300.400.50图 6 不同（,b）参数组合下的估计渗透系数场的 RMSE 值Fig.6 RMSE values of an estimated hydraulic conductivity fieldunder different(,b)parameter combinations 进一步类比最佳组合（=0.4，b=4）与经验组合（=0.1，b=1）的反演结果，分别进行场景 3 中渗透系数场估计。考虑到 2 种参数组合下反演精度较高，因此重点比较高低渗透系数区的估计精度。图 7 为 2 种参数组合下的渗透系数场的估计结果，从形态上看两者均得到与真实场相近的估计结果。对比估计场和真实场的高、低值区，最佳组合对低渗透性区域刻画更好，即图 7 中黑色虚线圈定的区域，并且最佳组合的估计方差场也相对较小，说明估计的渗透系数场离散型较小，因此对于本算例而言，最佳组合（=0.4，b=4）即为最优参数组合。虽然本算例中最优组合为（=0.4，b=4），但在实际工程案例中，由于无法获取研究区域真实的渗透系数 200 水文地质工程地质第 4 期场分布，也就无法比较不同的参数组合。鉴于本文中经验组合（=0.1，b=1）同样取得了不错的精度，并在过往研究中均具有最好的综合效果26 27，因此在实际应用过程中，如果对精度的要求不是过高，经验组合的普适性或许更值得推荐。5 结论（1）使用 ILUES 算法同化水头与浓度数据，可以在 8 个迭代步骤内实现多点污染源源强的精确反演，这表明 ILUES 算法在识别污染源方面具有很高的普适性。但同时精确表征渗透系数在空间上呈现出的非均质性，是预测污染物迁移路径、反演污染强度的关键。（2）本研究中的场地参数大小为 3040=1 200 维，如果直接反演渗透系数场将造成极大的运算负荷。通过场景 3 的测试，依托于 KL 展开方法，可以将渗透系数场的参数维度降低至 200 维的同时保留大部分的含水层参数变异性（95.20%），在后续反演中只需反演降维后的参数，即可估计整个渗透系数场，显著降低了计算压力。（3）ILUES 算法能同步精确识别污染源参数和渗透系数场，其中集合样本数量 Ne 和 ILUES 算法参数 和 b 影响着反演效果。综合考虑计算效率和计算精度，本研究中选择 Ne=4 000 为最佳样本集合大小，5.04.54.03.52.53.02.01.51.05.04.54.03.52.53.02.01.51.00000.100.080.060.040.0200200长度/m长度/m长度/m长度/m长度/m宽度/m宽度/m宽度/m宽度/m宽度/m2002002002001501501505.04.54.03.52.53.02.01.51.0150lnK/（md1）lnK/（md1）lnK/（md1）150样本方差样本方差0.100.080.060.040.0200（a）真实渗透系数场（b）反演场（=0.1，b=1.0）（d）反演场（=0.4，b=4.0）（c）反演方差场（=0.1，b=1.0）（e）反演方差场（=0.4，b=4.0）图 7 2 种（，b）参数组合下渗透系数场的估计结果Fig.7 Estimation of the hydraulic conductivity field under the combination of two(,b)parameters2023 年吴延浩，等：地下水污染强度及渗透系数场的反演识别研究 201 （=0.4，b=4）为算法最佳参数组合。但在实际工程案例中，如果对精度的要求不是过高，经验组合（=0.1，b=1）更值得推荐。（4）实际上，水文地质问题不仅包括渗透系数的不确定性，还包括着弥散度、孔隙度等参数的不确定性。另外预测污染物运移不能忽略污染物自身的变化过程，因此化学反应参数（如速率常数、分配系数等）的准确获取对于地下水溶质运移模型同样十分重要。未来识别地下水污染源的研究可以集中在这些特征的影响上，进一步验证数据