ILUES算法融合ERT数据反演污染源参数与渗透系数场

张瑞城，周念清，江思珉，夏学敏

（同济大学土木工程学院，上海 200092）

随着现代工农业生产的快速发展和污染物排放量的不断增加，各种污染物通过不同的途径入渗进入到地下水体中，导致地下水污染问题日趋严重。在污染场地研究中，除了采用传统的监测、取样、分析和评价外，通过建立数值模型对地下水污染物的迁移转化过程进行模拟预测已成为研究地下水污染问题的有效途径和重要手段。准确获取模型参数是进行数值模拟预测的关键，然而，由于参数存在空间变异性的特点，难以利用稀疏的观测信息来准确描述非均质场地的模型参数［1］。近些年来，将地球物理方法（如高密度电阻率成像法、地质雷达法等）和间接观测手段引入污染水文地质研究中，其应用越来越广泛［2-4］。

求解地下水逆问题通常需要进行数据同化，利用浓度、水头值以及采用地球物理方法获取的间接观测值来对模型参数进行反演，从而对系统状态提供更科学可靠的估计［5］。其中，基于集合的数据同化算法有集合卡尔曼滤波器（EnKF）和集合平滑器（ES）等，在地下水逆问题中已被广泛使用［6-8］，二者的区别是EnKF需要同时更新模型参数和状态，而ES仅进行参数更新，就能有效避免EnKF中参数和状态的不一致问题［9］。为了提高计算的运行效率和准确性，有学者提出了一种迭代局部更新集合平滑（iterative local updating ensemble smoother，ILUES）算法，通过采用局部更新策略，实现了对高维强非线性问题的求解［9-10］。

国内外有关地下水污染参数反演研究，大多集中于污染源参数（源强和位置）或场地水文地质参数（渗透系数场、孔隙度等）等单一类型参数的识别，而对污染源源强和渗透系数场的联合反演少有研究［11-12］。与此同时，现有数据同化方法主要通过融合水头与污染物浓度观测值进行参数反演，融合地球物理数据和水头浓度数据的研究相对较少［1，13］。有鉴于此，本文综合考虑污染源参数与场地水文地质参数，将污染源源强和渗透系数场作为未知模型参数，构建基于迭代局部更新平滑器算法（ILUES）的数据同化框架，并融合ERT（高密度电阻率法）数据进行联合参数反演，通过数值算例验证了该方法的有效性。

1 研究方法

1.1 地下水流与溶质运移模型

对于饱和多孔介质中的三维地下水流运动，有如下基本微分方程：

式中：Kx，Ky，Kz分别为渗透系数在X，Y，Z方向上的分量，m·d-1；H为水头，m；W为单位时间从单位体积含水层中流入或流出的水量，d-1；Ss为的贮水率，m-1；t为时间，d。

针对饱和多孔介质中的溶质运移问题，其基本控制方程如下：

式中：C为浓度，g·m-3；t为时间，d；Dij为水动力弥散系数，m2·d-1；vi为含水介质中的实际水流速度，m·d-1；qs为单位体积含水层的源汇项流量，d-1；Cs为源汇项中溶质浓度，g·m-3；∑Rn为化学反应项总和，g·m3·d。

在建立地下水流三维运动方程和溶质运移模型的基础上，分别采用MODFLOW程序［14］求解地下水流模型，用MT3DMS程序［15］求解溶质运移方程。

1.2 地层电阻率与溶质浓度关系的构建

溶质具有导电性，不同的溶质其导电性与电阻率差异比较大。Archie公式用于刻画溶液电阻率与地层电阻率以及含水饱和度之间关系［16］，其有效性得到了许多学者的证实和采用，基本形式如下：

式中：φ为有效孔隙度；Sw为含水饱和度；Rw为溶液电阻率，Ω·m；Rt为地层总电阻率，Ω·m；m为黏结指数；n为饱和度指数。

溶液电阻率与温度和溶液浓度等因素的有关，Sen［12］综合考虑了这两种因素对溶液电阻率的影响，将物理化学理论与室内实验相结合，提出了如下模型公式。

式中：T为溶液的温度，根据场地调查数据取为15℃，其余参数同前。

1.3 地球物理反演模型

高密度电阻率法以岩土体导电性差异为基础，通过供电电极（A、B）向地下空间输入稳定的电流I，然后根据测量电极（M、N）的电位差分析电场的分布，进而推算出视电阻率的空间分布。视电阻率正演模型可表示如下：

地球物理反演模型是利用反演方法通过对ERT数据（观测值）的拟合实现正演模型参数（地层电阻率）的反演识别，可以利用开源程序Pygimli构建地球物理模型并求解电阻率的反演问题［17］。

1.4 基于ILUES算法的数据同化框架

ILUES算法是一种迭代形式的ES算法，其协方差矩阵扰动策略来源于Emerick和Reynolds［18］所提出的ES-MDA算法，通过添加扰动的观测误差协方差矩阵对观测数据进行多次同化，从而实现对参数样本的迭代更新。

为了更好地解决高维非线性参数反演问题，ILUES算法没有采用传统集合平滑方法的全局更新策略，而是使用了一种优化的局部更新策略，即只对每个参数样本的局部样本集合进行更新。在ILUES算法中，对于集合中的每一个样本，其局部样本集合通过以下距离评价因子进行筛选：

本文构建了融合ILUES算法、Archie公式、Sen模型和溶质运移模型的数据同化框架，其浓度和水头及地球物理数据（ERT数据）作为数据同化的观测值。算法流程如图1所示。

图1 基于ILUES算法的数据同化框架Fig.1 The framework of data assimilation based on ILUES algorithm

2 算例研究

研究区为矩形区域（二维XZ剖面，100m×40m），概化为非均质各向异性的承压含水层，用边长为2m的正方形网格将研究区剖分为20行50列有限差分网格（如图2）。假设含水层水流运动为稳定流，研究区下边界为隔水边界，左右边界为定水头边界（左边界水头为50.0m，右边界水头为49.0m），上边界无源汇项。含水介质孔隙度为0.30，纵向弥散度为10m，垂直横向弥散度为3m。初始时刻，研究区无污染物。含水层渗透系数K满足对数正态分布，其均值（ˉˉˉˉˉlnK）与方差（σ2lnK）分别为2.0和1.0，x，z方向的相关长度分别为60m和24m，变差函数为指数型，真实渗透系数场见图2。

图2 污染场地概念模型Fig.2 The conceptual model of contaminated site

假定研究区存在3个潜在点源污染，分别标记为S1、S2和S3（S2处实际上无污染物释放），污染物进入含水层的释放强度如表1，场地内布置有12个监测点（图2）。观测点的浓度观测值是利用地下水流和溶质运移模型正演计算得到。模拟总时长为200d，均匀离散为10个应力期，且只在前6个应力期内发生污染物释放，其污染源强度见表1。由正演模拟值和正态分布的测量噪声计算得到观测值。

表1 各应力期的污染源源强真实值Tab.1 The reference true values of pollution sources strength in different stress periods

式中：Csimu表示模拟值；Cobs为观测值；L为噪声水平（取10%）；δ为满足标准高斯分布的随机偏差。

高密度电阻率成像法采用温纳装置进行探测，沿地表布置一排电极，参见图2，共获取了10组ERT观测数据，ERT测量值相对误差取10%。为了研究利用场地观测数据进行污染溯源和渗透系数场估计的可行性，设计了算例1和算例2，分别通过同化传统观测数据（水头与浓度）、地球物理数据（ERT数据）来反演未知模型参数。算例1中选取12个观测点（图2中三角形）的水头和t＝［20，30，40，50，60，70，80］day的浓度作为观测值；算例2选取t＝［20，40，60，80］d的ERT监测数据作为观测值，算例的具体设置见表2。

表2 算例设置Tab.2 The setting of case studies

考虑对空间上连续变化的渗透系数场进行反演时，渗透系数场的维度与模型的网格剖分程度有关（本算例为20×50＝1 000维）。高维地下水逆问题的求解效率较低，一般需要对渗透系数场进行降维，采用Karhunen-Loève展开［19］的方法对其进行降维。

式中：ξi是独立标准高斯分布的随机数；τi和si(x，y)是特征值和特征函数；NKL是KL展开的项数。

为了定量评估ILUES算法同化2类不同观测数据的参数反演效果，引入均方根误差（root-meansquare error，RMSE）作为评价标准，该值反映了由ILUES算法反演得到的模型参数估计值与其真实值之间的差距，RMSE的值越小（趋近于0），说明参数反演的精确度越高。RMSE的计算如下：

式中：ŷn为参数反演得到的模型参数估计值；yn为模型参数的真实值；Ne为ILUES算法中集合的大小。

3 结果与讨论

在2个算例中，渗透系数场的KL展开项数NKL＝60，可以保留lnK场约95%的变异性，从而对渗透系数场（20×50＝1 000维）的识别降维为对KL展开的60个高斯随机数的反演。对于潜在污染点源S1、S2和S3，污染物释放发生在前6个应力期，共有18个污染源强度参数需要识别，其先验分布为［0，20］的均匀分布。因此，对于2个算例总共需要反演识别的未知模型参数数量均为78个。为了获得准确的参数反演结果，根据以往数值算例的经验，将ILUES算法的参数设置为：集合数Ne＝1 000，迭代次数NIter＝10，α＝0.1。

3.1 污染源源强反演结果

图3为算例1和算例2的污染源源强的反演结果。图3a～3c是利用水头和［20，30，40，50，60，70，80］d的浓度观测值，得到的源强反演结果，经10次迭代后基本收敛，但与各应力期的源强真值仍有较大的偏差；相较而言，图3d～3f融合［20，40，60，80］d的ERT时变数据反演得到的污染源源强，经7次迭代后已经收敛于真实值，其中S2各应力期的源强值均接近于0（即无污染发生）。

图3 污染源源强的反演结果Fig.3 The inversion results of contaminant source strength

因此，融合ERT数据的ILUES算法可以更好地估计污染源源强。值得注意的是，算例2中ERT数据测量时间为［20，40，60，80］d，而算例1中污染物浓度观测时间为［20，30，40，50，60，70，80］d，在较少时间点采集的ERT数据的反演结果反而更好，其原因在于ERT等地球物理方法虽然测量精度不如水头、浓度数据，但可以便捷地获取大量具有空间连续性的观测数据，其采集的观测数据量（620个）远多于算例1中的水头与浓度数据（96个）。

为了更好地量化和对比使用ERT数据和水头浓度数据的污染源源强反演效果，将算例1和算例2中源强反演结果的RMSE值列于表3。由表3可知，使用ERT数据作为观测值的算例2，任一污染源的任一应力期对应的均方根误差均小于使用水头与浓度数据的算例1，这进一步说明了利用ILUES算法融合ERT数据进行参数反演可以获得更加精确的结果。

表3 污染源参数反演结果的均方根误差Tab.3 Root-mean-square error of Case 1 and Case 2

3.2 渗透系数场的反演结果

图4和图5是经过10次迭代后lnK的反演结果。对于算例1和算例2，作为参照场的lnK真实场如图4a和图5a；初始时刻随机生成的lnK场如图4b和图5b，为了有效对比算例1与算例2的反演精度，两算例采用相同的初始随机lnK场。分别随机选取2个后验估计，见图4c～4d和图5c～5d，lnK的估计均值场与方差场（反映集合的离散型和估计的不确定性）分别见图4e、4f和图5e、5f。对比污染物分布形态及高、低值区域，发现算例2拟合程度更好，且对渗透系数场的刻画更为精细。

在使用水头与浓度数据的算例1中（见图4），其2组后验估计场在轮廓上基本反映了lnK场的形态及高、低值区域，但是在细节刻画上远不如使用ERT数据的算例2，其反演得到的lnK场与参照场仍有不小的差距，同时如图4f所示算例1的估计方差已然趋近于0，说明在现有基础上继续进行迭代对反演结果的提升非常有限。在使用ERT数据的算例2中（见图5），两组后验估计场与参照场更为接近，方差场（图5f）在研究区左右两侧较高，其原因是利用温纳装置获得的呈倒梯形分布的ERT数据无法覆盖到相关区域。由此可知，融合ERT数据的ILUES方法可以更加精确地估计污染场地的渗透系数场。

图4 渗透系数场的反演结果（算例1）Fig.4 The inversion results of hydraulic conductivity in Case 1

图5 渗透系数场的反演结果（算例2）Fig.5 The inversion results of hydraulic conductivity in Case 2

在此基础上，将算例1与算例2反演所得ln K场代入地下水流与溶质运移模型，获取其末时刻的污染物浓度分布，并与作为参照的真实浓度分布进行对比，如图6。图6a反映了研究区末时刻真实的污染物浓度分布；图6b是使用初始随机lnK场正演得到的末时刻浓度分布；图6c、6d分别为算例1和算例2反演所得lnK场代入地下水流与溶质运移模型中所获取的污染物浓度分布情况。经过对比可以得到以下结论：虽然初始随机生成的lnK场与真实场有较大的差距，但是经过算例1与算例2的参数反演后，其分布基本接近于真实场，证明了ILUES算法在求解参数反演问题上的有效性；算例2末时刻的浓度分布情况与真实浓度分布更为接近，在细节的刻画上更加精确，进一步证明了融合ERT数据的ILUES方法的优越性。

图6 ln K正演所得末时刻浓度分布Fig.6 The final concentration distribution from forward modeling of ln K

4 结论

基于ILUES方法构建融合溶质运移模型和地球物理模型的地下水污染溯源的数据同化框架，并通过数值算例分析，得到以下结论：

（1）使用ILUES算法同化水头与浓度数据和ERT数据，可以实现对各应力期污染源源强和渗透系数场的联合反演，并且随着迭代次数的增加，反演结果越来越趋近于真实值，证明基于ILUES算法的数据同化框架对于求解地下水参数反演问题的有效性。

（2）在ILUES数据同化框架下，设计了两个数值算例对比同化浓度与水头数据（算例1）和ERT数据（算例2）在相同条件下对于反演结果精度的影响。在对污染源源强的反演上，算例2结果的RMSE值在3个潜在污染源的全部应力期均小于算例1。并且在对渗透系数场的反演上，算例2对于lnK场细节的刻画上明显优于算例1。充分说明了融合ERT数据的ILUES方法在求解地下水参数反演问题上的优越性。

（3）算例2通过在较少的观测时间点采集信息，得到的反演结果精度反而优于观测时间点较多的算例1，其原因在于使用高密度电法一次可以获取大量的具有空间连续性的观测信息，而传统的浓度与水头数据，只能在不连续的少量观测井处获取。在获取观测信息的数量与规模上，结合地球物理方法的算例2要明显优于使用传统观测方法的算例1。并且由于利用高密度电法等地球物理方法进行观测的便利性和观测数据的连续性，其在实际场地应用中将更有优势。

（4）本文仅研究了对于非均质渗透系数场的反演，而在实际应用中，孔隙度、弥散度等水文地质参数也存在着不同程度的非均质性，且各参数之间存在着一定的因果关系和相关性。后续研究计划结合机器学习和深度学习等方法同时反演多种水文地质参数。

作者贡献声明：

张瑞城：论文的主要撰写，数据处理与分析。

周念清：指导论文的总体框架，协助修改与定稿。

江思珉：协助完善数值方法的理论部分。

夏学敏：协助完善图表的绘制与处理。