地下水污染物溯源的数学模拟方法研究进展

龙玉桥，崔婷婷，李 伟，李砚阁，吴春勇

(1. 南京水利科学研究院，江苏 南京 210029；2. 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098)

地下水污染物溯源的数学模拟方法研究进展

龙玉桥1, 2，崔婷婷1, 2，李 伟1，李砚阁1，吴春勇1

(1. 南京水利科学研究院，江苏 南京 210029；2. 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098)

在当前我国地下水资源供需矛盾突出，地下水环境不断恶化的情况下，加强地下水污染物溯源的有关研究对于防治地下水污染、保护地下水资源、保障城乡居民饮水安全、保证经济和社会又好又快发展有着重要的实际意义。探索地下水污染物溯源的新理论和方法，在推动不适定问题的求解、不确定性研究，完善我国地下水污染理治技术与方法方面具有明显的理论意义。根据国外近三十年的研究成果，分析地下水污染物溯源问题的特点，将常用数学模拟方法归纳为模拟-优化法、解析和回归方法、直接法、随机理论法，讨论现有数学溯源方法的优缺点，认为应在复杂溯源问题、物理化学及生物作用的表达、模型计算效率和不确定性方面加强溯源数学方法的研究。

地下水；污染物；溯源；数学模拟方法

2008年，我国总供水量5 910亿 m3，其中，地下水源供水量占18.3%[1]。我国有400多个城市开采地下水，华北、西北地区城市利用地下水比例分别高达72%和66%，地下水往往是部分城市和农村唯一的供水水源[2]。地下水是我国经济和社会发展以及人民生活所必需的、不可替代的重要资源[3]。

近20年来，在我国人口较为密集、人类活动干扰大、工农业生产发达的平原地区，由于工业废水和生活污水的排放，大面积、超量化肥和农药的使用，垃圾场的淋滤和地下油罐的渗漏等原因，地下水正遭受着越来越严重的污染。据全国118座大城市浅层地下水的调查，97.5%的城市受到不同程度的污染，其中40%的城市受到严重污染[2]。在全国水资源调查评价的197万 km2平原区浅层地下水中，Ⅰ类和Ⅱ类水质区的面积仅为总流域面积的4.98%，Ⅲ类面积为35.53%，Ⅳ、Ⅴ类面积高达59.49%[4]。太湖、辽河、海河、淮河等流域地下水污染最为严重，劣于Ⅲ类水质的水质区面积占各相应流域面积的91.49%，84.55%，76.40%和67.78%[4]。

为保护地下水资源，我国在地下水污染防治方面做了大量的实践工作，如水利部开展的全国地下水保护行动、国土资源部开展的全国地下水污染调查和环保部制定的地下水污染防治规划。此外，针对地下水中硝酸盐和石油类污染物的治理技术与方法也涌现出大量的研究成果[5-10]。然而，国内的实践和研究的成果主要集中于污染防治的政策与手段、污染治理的技术与方法方面，在地下水污染物溯源方面的研究却鲜有报道。

地下水污染物溯源(地下水污染源解析)是通过有限的观测数据，查明污染源的位置及污染物迁移转化的历史，是地下水污染治理的首要步骤[11-13]。它主要包括三个方面的研究[11, 14, 15]：追溯污染物释放历史、确定污染源位置及排放时间、判断污染物在释放初期的分布。地下水污染物溯源的研究有助于制定和选择高效、经济的污染治理策略与方法，有助于判别不同污染源或污染肇事者的责任大小，也有助于确定污染治理成本在不同污染肇事者间的分配比例[16]。

地下水污染物溯源的方法主要有两类[14]：地球化学足迹法(geochemical fingerprint techniques)和数学模拟法(mathematical and simulation approaches)。地球化学足迹法依据于实地测量的污染物数据，通过同位素、地质统计分析等方法确定污染源的位置[17][18]。单独使用地球化学足迹法不能完全解决追寻污染源位置及污染物随时间演变过程等问题[19, 20]。溯源是一种典型的逆问题(不适定问题)，该类问题是难以求解的。目前，在合适的假设条件下，数学模拟法已从解一维均质的理想算例[16, 21, 22]，逐步发展到求解二维非均质实际溯源问题[23, 24]，但受到含水层参数、污染物浓度及模型自身的不确定性等因素的限制，尚不能完全满足实际应用的要求。因此，地下水污染物溯源的理论与方法还有较大的发展空间。

在当前我国地下水资源供需矛盾突出，地下水环境不断恶化的情况下，加强地下水污染物溯源的有关研究对于防治地下水污染、保护地下水资源、保障城乡居民饮水安全、保证经济和社会又好又快发展有着重要的实际意义。探索地下水污染物溯源的新理论和方法，在推动不适定问题的求解、不确定性研究，完善我国地下水污染理治技术与方法方面具有明显的理论意义。为此，本文讨论地下水污染物溯源问题的特点、常用数学模拟方法的研究进展和存在问题。

2 污染运移理论及溯源问题的特点

考虑对流、弥散、流体的汇(源)、平衡吸附作用和一级不可逆速率化学反应的溶质迁移方程为[25]：

(1)

解对流-弥散方程需要求解同时含有双曲项(与对流有关)和抛物项(与弥散有关)的偏微分方程，至今这类问题仍未完全解决，因此求解对流-弥散方程一直是地下水模拟领域的难题。溯源问题要求通过有限的观测数据，查明污染源的位置或追溯污染物迁移转化的历史，要逆时间求解对流-弥散方程[19]，难以满足解唯一性或稳定性，是一种典型的逆问题(不适定问题)。因此溯源问题的数学模拟难度比常见的溶质迁移的模拟难度更大。

弥散是污染物在地下水中迁移的主要形式，但弥散过程具有不可逆性[19]，这使溯源问题更难求解。求解溯源问题的难易程度还取决于观测数据的详实程度和要确定的系统输入的多少[26]。水动力弥散、模型不确定性、测量误差和观测数据缺失，都会增加溯源问题的求解难度[14]。溯源问题的解对输入数据的误差有极强的敏感性，即便是很小的污染羽的测量误差也能使计算得出的污染羽的历史发生巨大的变化[16]。

Liu和Ball(1999)[27]根据求解溯源问题的数学方法，大致将溯源问题分成全估计问题和函数拟合问题。全估计问题是指判定某位置上污染源随时间的变化历史或是判定某时间点上污染物的空间分布，然而却不知道这些分布的函数形式。全估计问题是在无限维空间上求最小值的问题。函数拟合问题是指已知需要估计的浓度或通量函数的一些信息，是在有限维空间上求最小值的问题。函数估计问题的待估计参数越多(如求不同时间步长对应的浓度时)，函数估计问题与全估计问题间的差别就越小。

3 数学模拟方法

地下水污染物溯源的数学方法研究已有近30 a的历史。不同的学者对数学方法的分类有不同的见解。Atmadja和Bagtzoglou(2001)将地下水污染溯源的数学模拟方法归纳为四类[19]：优化方法、解析方法、直接法、随机理论和地质统计学方法。Sun等(2006)将地下水污染溯源的数学模拟方法归纳为三类[11]：优化方法、随机理论方法、逆时求解对应弥散方程。文中采用Atmadja和Bagtzoglou(2001)的分类方法，归纳地下水污染物溯源的数学方法。

3.1 模拟-优化方法

模拟-优化方法将溯源问题转化成优化问题，即寻求污染源的位置、浓度等变量或它们的组合，使模拟模型的输出结果最接近观测数据。此方法常将正向的地下水流运动和污染物迁移的模拟模型与优化模型耦合在一起，形成模拟-优化模型。在模拟-优化模型中，正向的模拟模型被反复调用，各次模拟的结果将与观测数据进行比较，优化模型则用于寻找使模拟结果最接近观测数据的解(污染源的位置、浓度等变量)。按照优化方法的差异，可大致将模拟-优化方法进一步分为非启发式优化算法和启发式优化算法。

3.1.1 模拟优化模型的耦合方法

地下水流运动和污染物迁移的模拟模型与优化模型耦合方法主要有响应矩阵法和嵌入法。响应矩阵法[11, 28]利用响应系数矩阵描述地下水系统内脉冲(抽注水量或污染物质量)与响应(地下水位的升降或污染物浓度变化)的关系，并借助此矩阵实现模拟模型与优化模型的耦合。该方法可以显著地减少优化模型中决策变量和约束方程的个数，但它的理论基础要求地下水运动方程和边界均要具备线性齐次条件，这在实际应用中却较难满足[29]。嵌入法[12, 24, 30]是将地下水流运动和污染物迁移的模拟模型作为优化模型的一部分约束条件，实现模拟模型与优化模型的耦合。嵌入式的模拟-优化模型要将各时刻所有空间离散点对应的水头和浓度作为优化模型的决策变量，因此模型的计算量很大[31]。Datta等(2009)[31]将水流及溶质模拟模型作为外部模块，利用Jacobian矩阵将优化模型和模拟模型连接在一起并控制优化模型的搜索方向，同时确定污染源位置和含水层参数，这一方法显著减少了决策变量数，提高了计算非线性优化问题的可行性。

3.1.2 优化算法

1) 非启发式优化方法

Gorelick等(1983)[28]较早使用优化方法研究了理想的二维含水层中的污染物溯源问题，他们采用线性规划法，多元线性回归，并假设算例含水层参数不存在不确定性。Mahar和Datta(1997)[12]首先初步识别污染源的位置，并根据结果优化观测井网，最后利用优化后的观测井网得到更准确的源位置，利用非线性规划法确定污染源位置。随后，Mahar和Datta(2001)[30]将文献[12]的模拟优化方法拓展到同时估计含水层参数和识别污染源位置(约束是线性的非线性规化问题的解法是减化梯度算法和拟牛顿算法的混合方法，非线性约束问题的解法是投影增广拉格朗日算法)。Sun等(2006)[11]利用二次规划软件包SeDuMi求解了二维理想算例的溯源问题。他们考虑了模型本身的不确定性对溯源结果的影响，使用了带约束的鲁棒最小二乘方法恢复了二维理想算例中的污染物排放历史，他们的方法可减小最优解对模型和误差引起的扰动的敏感性，且易于与MT3DMS等溶质迁移模型结合。Ghafouri和Darabi(2007)[24]使用增广拉格朗日乘子法求解模拟-优化模型，从而确定污染源的位置和范围，他们将其方法应用于确定伊朗西南部的Ramhormoz含水层中的污染源位置。增广拉格朗日乘子法减少了优化问题对罚函数的系数的依赖性，也使原约束优化问题转变为无约束优化问题。他们的方法中所有的空间离散点都可能作为污染源，其模型将低浓度的点除去，从而大幅度削减了计算量。Ghafouri和Darabi的方法中还考虑的吸附作用，这是之前研究中未有考虑的。

2) 启发式算法

地下水污染物溯源问题的优化问题常具有非线性、解空间非凸的特点，这是传统优化方法难以应对的[32]。因此，启发式算法越来越多地被应用于这一领域。

遗传算法是最为常见的启发式算法。Aral等(2001)[33]将进步式遗传算法应用于解决二维含水层中的污染物溯源问题，提高了模型的计算效率。他们的方法在测量误差较大的情况下也可取得令人较为满意的结果，污染源位置的起始搜索位置对该方法的搜索效果影响不大。如果观测数据满足：(1)在一个时间步长内至少有一个观测点有观测数据；(2)独立的观测数据的数量多于待确定变量的个数，那么观测数据的缺少不会对溯源结果造成过坏的影响。遗传算法的全局搜索能力强，但局部搜索能力有所欠缺，为此Mahinthakumar和Sayeed(2005)[34]将遗传算法和局部搜索算法相结合，形成GA-LS方法，并将其应用于三维溯源问题。GA用于确定单一非点状污染源的位置，局部搜索算法则用于细化GA算法的识别解果。Singh和Datta(2006)[35]采用外部连接的形式将遗传算法和地下水流及溶质迁移模型耦合在一起，形成模拟优化模型，并将其应用于二维溯源问题。他们的方法可用于确定多个污染源的位置及排放历史，即使存在观测数据存在测量误差也可得到较好的溯源结果。当污染源个数增加时，溯源问题的复杂程度也随之增加。污染源如果不局限于几个可能位置，而是散布在一个区域内时，测量误差对溯源结果的误差影响更大。这种情况下，GA的方法要优于ANN方法。

Yeh等(2007)[36]利用禁忌搜索-模拟退火算法，结合三维污染物迁移模型，研究了理想算例中的污染源位置及污染物排放浓度与时间。禁忌搜索算法负责从疑似污染区中确定污染源位置，而模拟退火算法负责计算污染物排放浓度与时间的试探值。他们的方法具有对污染源的初始估计值的依赖程度低的优点。其成果表明至少需要6个监测点(每个监测点有四个浓度间隔的监测值)，才能正确地估计出污染源的排放信息。Yeh等的研究中的污染源是恒定源，其污染物排放的浓度和时间是不随时间变化的，然而对于间断排放的污染源，Yeh的方法仍有待于验证。

Mirghani等(2009)[26]采用进化算法研究了三维理想算例中的溯源问题。他们的研究表明含水层渗透系数的非均质性对算法的搜索效率没有决定性的影响。他们对算法的可靠性和一致性进行了较为深入的分析，并发现当采用更多的决策变量刻画污染源特征时，问题的复杂度加大，使预测误差和求解误差的非一致性变得更为显著。

Bharat等(2009)[37]以粒子群算法作为寻优工具，追溯了一维地下水流中污染羽的演变历史。Bharat等(2009)使用的算例与Skaggs和Kabala(1994)[16]相同。方法计算过程表明不同的解能得到相同的污染物排放曲线。且这些曲线呈现出陡峭的峰，需要施加全局约束平滑这些曲线。Vesselinov和Harp(2010)[38]将自适应粒子群算法和麦夸特搜索方法相结合，得到了SQUADS方法，这一方法兼顾了全局最优和局部最优解，能够高效地在解空间中搜索最优解，模拟结果与观测浓度的一致性较高。

上述各研究在解决溯源问题时，都假设已知污染源的个数，然而实际情况下污染源的数量也可能是未知的。Ayaz(2010)[39]引入和声算法求解溯源问题，并在和声算法过程中加入了一个隐含求解过程，实现在确定污染源位置和排放历史的同时，确定污染源的个数。虽然该项研究中假定边界条件、渗透系数、弥散系数不存在不确定性，但Ayaz认为考虑这些不确定性是重要的。测量误差能影响污染物的排放量的确定，不影响污染源的定位。污染源与观测井间的距离影响污染源特征的解析。当污染源的数量较多时，模型的计算量增大，这需要采用并行计算手段，从而提高计算效率。

3.1.3 计算效率的提高方法

在模拟-优化模型中，正向的模拟模型被反复调用，各次模拟的结果将与观测数据进行比较，优化模型则用于寻找使模拟结果最接近观测数据的解(污染源的位置、浓度等变量)。由于需要多次运行模拟模型，所以这一方法的计算量是非常巨大的。目前，削减计算量的方法通常有并行计算方法(Code parallelization)[26]和替代模型法(Surrogate modeling)[40]。

并行计算指同时使用多种计算资源解决计算问题的过程，是快速解决大型且复杂的计算问题的有效方法之一。并行计算方法有效地提高了计算效率，显著地缩短了计算时间，促进了溯源问题向三维问题的发展。Mahinthakumar和Sayeed(2005)[34]、Mirghani等(2009)[26]先后将并行计算引入了三维溯源问题中。Mirghani等(2009)[26]系统地分析了并行方法的效果，他们的研究表明这一方法能够胜任算例中的溯源工作，并行计算的方法使计算速度提高了约100倍。然而，并行计算是一个专业的工作，对使用者的专业素质(命令与代码操作)要求高，这限制了它在普通研究人员或用户中的推广。

相对于并行计算方法，替代模型方法则更易于被使用者所掌握。溯源问题涉及的替代模型方法主要以人工神经网络(Artificial Neural Network，ANN)为主。利用模拟模型的输出结果训练ANN，并将训练好的ANN作为替代模型用于污染物溯源问题。Singh和Datta(2004)[41]认为ANN是最优统计范型识别方法[42]的拓展与延续，他们将ANN方法用于二维理想溯源问题，在确定污染源排放量的同时，还进行了含水层参数的识别，并提出了评价模型计算效果及ANN结构的指标。Singh和Datta(2007)[43]进一步研究了ANN方法在无观测误差、有观测误差和观测数据缺失三种情况下的效果，并指出虽然这一方法能够解决理想的溯源问题，但是ANN在用于溯源问题时，存在以下局限性：(1)对ANN的训练决定溯源结果的优劣；(2)研究中使用的是均质含水层。随后Bashi-Azghadi等(2010)(23)将替代模型方法运用到真实的污染物溯源研究中。他们利用非支配排序遗传算法优化观测井位，后利用优化后的观测井位结合概率支持向量机和概率神经网络估计未知污染源的主要特征，并确定了德黑兰精炼厂的潜水含水层中的石油类碳氢化合物的污染位置和排放量。

3.2 解析方法和回归方法

解析方法的计算效率要高于模拟优化方法，适宜在水文地质条件及污染物迁移过程简单的情况下使用。Ala和Domenico(1992)[44]采用解析法研究了Otis空军基地的地下水污染强度和大小、污染物锋面的对流位置。污染物有氯化物、硼、三氯乙烯、四氯乙烯、可生物降解和不可生物降解的污染组份。Butcher和Gauthier(1994)[45]利用简化的解析解估计了Massachusetts东北部一工厂旧址的DNAPL污染物残留量，并通过假设条件，利用最小二乘算子求得简化解析解中参数。这一方法能判别残留的DNAPL源是否存在，但不能准确地估计残留污染源的体积和溶解时间。Sidauruk等(1998)[46]用解析法求解了污染物迁移的反问题，只需知道浓度分布，就能估计弥散系数、流速、污染物量及初始位置、时间起点，最小化线性回归的相关系数。由于使用了解析法，该方法只能用于几何形状简单的2维均质含水层和水流运动。

Alapati和Kabala(2000)[47]利用非线性最小二乘方法(NLS)追溯了污染物排放历史。他们的方法中没有利用正则化技术，他们的研究表明当污染物逐渐释放的情况下，NLS方法的结果对误差干扰十分敏感；当污染物大量释放的情况下，即使存在较大的误差干扰，NLS方法也能得到较好的结果，此时NLS方法可以替代Tikhonov正则化方法。

3.3 直接法

提克洛夫正则化(Tikhonov Regularization，TR)方法的主要思想是利用对解和数据误差的先验估计可以将问题的求解限定在某个较小范围内，对问题的提法进行适当的改造后，原本不适定的问题就可以转化为适定的最优化问题求解，而且先验估计表明在一定精度下用正则化方法求得的解是合理的。正则化技术需要污染源函数特征的先验知识，正则化项的权重较难确定[27]，也较难确保得到的解是最优解[34]。

Skaggs和Kabala(1994)[16]利用Tikhonov正则化追溯污染羽发展历史，并研究了污染羽的测量误差、迁移参数误差和数值稳定性对追溯污染羽的影响，他们的方法对舍入误差不敏感，但是污染羽的测量误差对计算精度影响较大。方法中的正则化项的权重对计算结果的精度影响较大，权重过小时，正则化的稳定作用失效，而权重过大时，反问题的解本质上是被人工平滑了。随后，Liu和Ball(1999)[27]在Skaggs和Kabala(1994)的研究基础上[16]，假定污染物的迁移方式以弥散为主，利用观测数据，估计了Dover空军基地低渗透性隔水层中三氯乙烯和四氯乙烯的污染羽演变历史。方法中将污染羽的演变历史假设为一个未知形式的函数，并利用Tikhonov正则化方法将溯源问题转化为一个最小化问题，该问题以估计边界上浓度函数为目的。

Neupauer等(2000)[48]比较了TR方法和最小相对熵(Minimum Relative Entropy，MRE)方法，他们的研究表明两种方法都能建立光滑的污染源函数，污染源演化历史是无误差的分段函数时，MRE要优于TR方法，而当已知的数据存在误差时，TR算法要比MRE更稳定。但是TR和MER方法的解的质量取决于具体的溯源问题和模拟者主观的输入[34]。

拟可逆方法(Quasi-Reversibility，QR)原本用于求解逆时间弥散问题[49]，Skaggs和Kabala(1995)[50]通过在移动坐标系内求解QR算子，将QR方法拓展到对流-弥散条件下的逆问题中，反演了一维含水介质中污染羽的变化过程。QR方法计算精度要劣于TR方法，但其计算效率要高于TR方法，此外，QR方法在应用于非均质介质时较TR方法容易[50]。目前，QR和TR方法尚无应用于非均质含水层的研究实例[19]。

伴随法是变分数据同化方法中的一种方法，伴随法具有数值稳定的优点，且能有效的减少模拟次数[14]，曾被用于含水层参数的估计问题[51, 52]和敏感性分析[52]。Neupauer和Wilson(1999)[53]指出逆时的污染物的位置概率和迁移时间概率是顺时的污染物的质量-体积浓度的伴随状态，并采用伴随方程代替了溶质迁移的控制方程，伴随状态替换了迁移方程中的因变量(溶质浓度)，研究了一维逆时对流-弥散方程的求解方法。随后，Neupauer和Wilson(2001)[54]又将此方法拓展到三维逆时对流-弥散方程的求解，并反演了一个二维理想算例中污染源位置。

常用的数值方法通常是进步式的，当时间步长是负值时，这些方法是不稳定的，因此难以在逆时间维的情况下追溯污染物的迁移历史[55]。Atmadja和Bagtzoglou(2001)[55]先将污染物迁移的控制方程转化为Jury Backward Beam Equation(JBBE)，然后使Jury方法和Marching方法相结合，得到Marching-Jury Backward Beam Equation(MJBBE)法。该方法不仅在求解一维非均质含水层中污染羽的空间分布与时间演化中取得了良好的效果，其计算效率也约是JBBE方法的计算效率的3倍[55]，但这种方法对噪声是敏感的[19]。

Milnes和Perrochet(2007)[14]在转换函数理论[56]的基础上，将观测值作为初始条件，反向模拟地下水流场，通过污染质浓度等值线的变化，确定污染源的位置和污染排放时间。他们的方法被用于二维非均质的理想算例中，但是方法基于的假设是已知含水层及迁移参数，已知地下水水流的流场，局限于单一不反应污染质的溯源问题。

3.4 随机理论法

Datta等(1989)[42]基于统计范型识别算法建立了追溯污染源的专业系统。他们研究了参数不确定性、测量误差、数据缺失对溯源结果的影响。他们的创立的方法，在数据缺失的情况下是较为有效的[23]。在获得相近的结果时，统计范型识别方法需要的数据要少于优化方法[16, 57]。

Wagner(1992)[15]采用非线性最大似然估计法同时估计了含水层参数、污染物的排放量，并利用一阶不确定性分析方法评价了最大似然估计法的可靠性及模拟模型的准确性。Wagner的研究中，假设含水层的参数分区和污染源的位置是已知的，因此更适用于对含水层参数和污染源位置有一定了解的溯源问题。

Bagtzoglou等(1991，1992)[13, 58]用随机走方法模拟理想二维非均质含水层中的逆时间溶质迁移方程，采用逆向的速度场达到逆时间的效果，污染物迁移模型中的对流项是逆时的，但弥散项是不变的。还采用地质统计学的方法评价了各污染源的相对重要程度。他们的方法可用于已知准确的渗透系数值及空间分布和渗透系数值及空间分布存在不确定性的情况。

Wilson和Liu(1994)[59]采用逆时随机微分方程求解溶质迁移方程。方程中的对流项是逆时的，但弥散项是不变的。他们指出在给定合适的边界条件后，这一方法能拓展到含有线性动态吸附和一阶衰变过程的溯源问题。Liu和Wilson(1995)[60]又将这一方法拓展到2维非均质含水层中。

Woodbury和Ulrych(1996)[21]利用最小相对熵研究了一维稳定流算例中污染羽的变化历史。污染源函数被离散成多个状态，每个状态被当作是一个统计变量，原问题被转换成估计这些变量的后验联合概率密度函数[27]。Woodbury和Ulrych(1996)指出这一方法不仅可以追溯污染历史，也可用于预测污染的未来的发展状态。随后Woodbury等(1998)[61]在三维的污染羽观测数据的基础上，利用该方法反演了Gloucester填埋场中1，4-二恶烷的排放历史。这一研究中，Gloucester填埋场的地下水流被概化成一维稳定渗流。Woodbury等(1998)的研究表明，早期的污染源排放历史的恢复效果较差，如果想反演出整个排放历史，污染羽的监测数据的时间系列则要尽可能的长。

Snodgrass和Kitanidis(1997)[22]将贝叶斯理论和地质统计技术相结合，将待估计的源函数离散成多个单元，每个单元都有给定的随机结构和未知的随机参数。这一方法的解更具有普适性，且不需要对未知污染源的特性和结构做出盲目的假设，但是此方法的限制是潜污染源的位置必须是先验的[19]。

地质统计反向分析法(geostatistical inversion approach)[22, 62, 63]也是一种常用溯源方法。方法假设待求的污染源的排放历史是一个未知函数，通过随机理论估计出这一函数，从而得到反问题的解。求解过程中，未知函数将被离散成多个点，每个点对应于一个观测值，对于每个离散点都要计算一次敏感性矩阵[62]。敏感性矩阵的计算量很大，多维情况下就更大[62]。随机理论方法的研究报道主要集中于一维介质中的污染质分布、点状污染源或层间污染源的排放历史，且较少涉及非均质介质中的污染物迁移问题。Michalak和Kitanidis(2004)[62]将地质统计分析与伴随法相结合，建立了追溯三维非均质污染物迁移的方法，并用其反演了二维非均质介质的理想算例中污染物的演变历史。研究中，地质统计分析法负责计算某给定点污染质的分布，伴随法则承担提高模型计算效率的重任。以往的地质统计分析法忽略了模型的不确定性，为此Sun(2007)[64]提出了鲁棒地质统计法(Robust Geostatistical approach)，在考虑参数不确定性和测量误差的前提下，将问题转化为求最小值的优化问题，最终采用半定规划方法求解。

4 结语

地下水污染物溯源是地下水污染治理的首要步骤。地下水污染物溯源工作是制定和选择污染治理策略与方法的基础，可为判定不同污染源或污染肇事者的责任大小提供依据，也可为确定污染治理成本在不同污染肇事者间的分配比例提供帮助。数学理论上，地下水污染物溯源问题本质上是一种典型的逆问题(不适定问题)，它具有解的不唯一性和解的不稳性，是国内外研究的重要难点。物理过程上，地下水污染物溯源问题涉及水流运动、溶质迁移、物理化学反应等方面，侧重于多学科方法的综合应用和交叉研究。地下水污染物溯源的数学方法的研究已有近30年的历史，目前主要的数学方法可分为模拟-优化方法、解析方法、直接法和随机理论方法，然而每种方法都有各自的优点与局限性，今后应进一步加强以下几方面的研究：

(1) 在复杂溯源问题中研究与应用。现有研究多集中于简单的理想算例，多数研究均针对一维或二维均质各向同性含水层中的污染物溯源问题，污染源多为点状源[62]。虽然这些研究中的方法在应用于简单算例时，取得了较好的效果，但却较难应对实际的工程问题。

(2) 溯源问题中物理化学及生物作用的表达。目前，溯源研究主要以单一的不反应污染物为主，少数研究中考虑了线性吸附作用，少数实例研究中考虑了有机污染物的变化。然而，地下水中污染物可能存在多种，也可能同时发生不同的物理化学及生物作用，这增加了溯源的困难程度，也使这一研究更有挑战性。

(3) 溯源模型计算效率的提高。溯源模型计算量巨大，计算耗时，尤其是采用优化方法时，计算时间长的问题更为突出。虽然并行技术和替代模型方法可以有效地削减计算时间，但是并行技术对使用者的专业素质要求高，这限制了它在普通研究人员或用户中的推广，而替代模型法不能体现溯源问题的物理基础，在使用时也受到限制。

(4) 多数溯源方法对观测数据的不确定性考虑较多，而对模型本身的不确定性考虑不足[11]。在模型参数确定的前提下，观测数据的误差是溯源结果好坏的主要决定因素。然而，实际问题中模型参数，如渗透系数和弥散系数，也存在不确定性，此时溯源问题将受到观测数据和模型本身不确定性的双重影响，异参同效的现象将更为显著。如果排除不确定性的困扰，获得正确的溯源结果，是研究人员面临的难题之一。

[1]中华人民共和国水利部. 2008年中国水资源公报[EB/OL]. http://www.mwr.gov.cn/zwzc/hygb/szygb/qgszygb/201001/t20100119_171051.html.

[2]姜建军, 文冬光. 合理开发利用地下水缓解水资源紧缺状况[J]. 中国水利.中国水利杂志专家委员会会议暨节水型社会建设高层论坛专辑. 2005，(13):36-39.

[3]薛禹群, 张幼宽. 地下水污染防治在我国水体污染控制与治理中的双重意义[J]. 环境科学学报.2009,29(3):474-481.

[4]唐克旺, 吴玉成, 侯杰. 中国地下水资源质量评价(II)-地下水水质现状和污染分析[J]. 水资源保护.2006,22(3):1-8.

[5]秦传玉, 赵勇胜, 郑苇, 等. 空气扰动技术对地下水中氯苯污染晕的控制及去除效果[J]. 吉林大学学报(地球科学版).2010,40(1):164-168.

[6]徐绍辉, 朱学愚. 地下水石油污染治理的水力截获技术及数值模拟[J]. 水利学报.1999(1):71-76.

[7]王志强, 武强, 邹祖光, 等. 地下水石油污染曝气治理技术研究[J]. 环境科学.2007,28(4):754-760.

[8]姜建军. 中国地下水污染现状与防治对策[J]. 环境保护.2007,38(10):16-17.

[9]陈秀成, 曹瑞钰. 地下水污染治理技术的进展[J]. 中国给水排水.2001,17(4):23-26.

[10]毕晶晶, 彭昌盛, 胥慧真. 地下水硝酸盐污染与治理研究进展综述[J]. 地下水.2010,32(1):97-102.

[11]Sun A Y, Painter S L, Wittmeyer G W. A constrained robust least squares approach for contaminant source release history identification[J]. Water Rescources Research, 2006,42, W04414, doi: 10.1029/2005WR004312.

[12]Mahar P S, Datta B. Optimal monitoring network and ground-water pollution source identification[J]. Journal of Water Resources Planning and Management, 1997,123(4):199-207.

[13]Bagtzoglou A C, Dougherty D E, Tompson A F B. Application of particle methods to reliable identification of groundwater pollution sources[J]. Water Resources Management, 1992,6(1):15-23.

[14]Milnes E, Perrochet P. Simultaneous identification of a single pollution point-source location and contamination time under known flow field conditions[J]. Advances in Water Resources, 2007,30(12):2439-2446.

[15]Wagner B J. Simultaneous parameter estimation and contaminant source characterization for coupled groundwater flow and contaminant transport modeling[J]. Journal of Hydrology, 1992,135(1-4):275-303.

[16]Skaggs T H, Kabala Z J. Recovering the release history of a groundwater contaminant[J]. Water Resources Research, 1994,30(1):71-79.

[17]Mansuy L, Philip R P, Allen J. Source identification of oil spills based on the isotopic composition of individual components in weathered oil samples[J]. Environmental Science and Technology, 1997,31(12):3417-3425.

[18]Rachdawong P, Christensen E. Determination of PCB sources by a principal component method with non-negative constraints[J]. Environmental Science and Technology, 1997,31(9):2686-2691.

[19]Atmadja J, Bagtzoglou A C. State of the Art Report on Mathematical Methods for Groundwater Pollution Source Identification[J]. Environmental Forensics, 2001,2(3):205-214.

[20]National-Research-Council. Groundwater Models, Scientific and Regulatory Applications[M]. Washington, D.C: National Academy Press, 1990.

[21]Woodbury A D, Ulrych T J. Minimum relative entropy inversion: theory and application to recovering the release history of a groundwater contaminant[J]. Water Resources Research, 1996,32(9):2671-2681.

[22]Snodgrass M F, Kitanidis P K. A geostatistical approach to contaminant source identification[J]. Water Resources Research, 1997,33(4):537-546.

[23]Bashi-Azghadi S N, Kerachian R, Bazargan-Lari M R, et al. Characterizing and unknown pollution source in groundwater resources systems using PSVM and PNN[J]. Expert Systems with Applications, 2010,37(10):7154-7161.

[24]Ghafouri H R, Darabi B S. Optimal Identification of Ground-Water Pollution Sources[J]. International Journal of Civil Engineering, 2007,5(2):144-154.

[25]Zheng Chunmiao, Bennett Gordon D. 地下水污染物迁移模拟[M]. 孙晋玉, 卢国平, 译. 第二版. 北京: 高等教育出版社.2009.

[26]Mirghani B Y, Mahinthakumar K G, Tryby M E, et al. A parallel evolutionary strategy based simulation-optimization approach for solving groundwater source identification problems[J]. Advances in Water Resources, 2009,32(9):1373-1385.

[27]Liu C, Ball W P. Application of inverse methods to contaminant source identification from aquitard diffusion profiles at Dover AFB, Delaware[J]. Water Resources Research, 1999,35(7):1975-1985.

[28]Gorelick S M, Evans B, Ramson I. Identifying sources of groundwater pollution: an optimization approach[J]. Water Resources Research, 1983,19(3):779-790.

[29]陈葆仁, 吴吉春, 刘淑芸. 地下水管理模型在我国实践中存在问题的讨论[J]. 水文地质工程地质.1994(6):36-39.

[30]Mahar P S, Datta B. Optimal identification of ground-water pollution sources and parameter estimation[J]. Journal of Water Resources Planning and Management, 2001,127(1):20-29.

[31]Datta B, Chakrabarty D, Dhar A. Simultaneous identification of unknown groundwater pollution sources and estimation of aquifer parameters[J]. Journal of Hydrology, 2009,376(1-2):48-57.

[32]McKinney D C, Lin M D. Genetic algorithm solution of groundwater management models[J]. Water Resources Research, 1994,6(30):1897-1906.

[33]Aral M M, Guan J, Maslia M L. Identification of contaminant source location and release history in aquifers[J]. Journal of Hydrological Engineering, 2001,6(3):225-234.

[34]Mahinthakumar G, Sayeed M. Hybrid genetic algorithm - local search methods for solving groundwater source identification inverse problems[J]. Journal of Water Resources Planning Management, 2005,131(1):45-57.

[35]Singh R M, Datta B. Identification of groundwater pollution sources using GA-based linked simulation optimization model[J]. Journal of Hydrological Engineering, 2006,11(2):101-109.

[36]Yeh H D, Chang T H, Liang Y C. Groundwater contaminant source identification by a hybrid heuristic approach[J]. Water Resources Research, 2007,43, doi:10.1029/2005WR004731.

[37]Bharat T V, Sivapullaiah P V, Allam M M. Swarm Intelligence Based Inverse Model for Characterization of Groundwater Contaminant Source[J]. Electronic Journal of Geotechnical Engineering, 2009,14(B):1-14.

[38]Vesselinov V V, Harp D R. Contaminant source identification using adaptive hybrid optimization of inverse groundwater transport model[J]. Water Resources Research, 2010.

[39]Ayvaz M T. A linked simulation-optimization model for solving the unknown groundwater pollution source identification problems[J]. Journal of Contaminant Hydrology, 2010,117(1-4):46-59.

[40]Singh R M, Datta B, Jain A. Identification of unknown groundwater pollution sources using artificial neural networks[J]. Journal of Water Resources Planning and Management, 2004,130(6):506-514.

[41]Singh R M, Datta B. Groundwater pollution source identification and simultaneous parameter estimation using pattern matching by artificial neural network[J]. Environmental Forensics, 2004,5(3):143-159.

[42]Datta B, Beegle J E, Kavvas M L, et al. Development of an expert system embedding pattern-recognition techniques for pollution-source identification, Technical Report: PB-90-185927/XAB, OSTI ID: 6855981[R]. Department of Civil Engineering, California University, Davis, CA (USA), 1989.

[43]Singh R M, Datta B. Artificial neural network modeling for identification of unknown pollution sources in groundwater with partially missing concentration observation data[J]. Water Resources Management, 2007,21(3):557-572.

[44]Ala N K, Domenico P A. Inverse analytical techniques applied to coincident contaminant distributions at Otis Air Force Base, Massachusetts[J]. Ground Water, 1992,30(2):212-218.

[45]Butcher J B, Gauthier T D. Estimation of residual dense NAPL mass by inverse modeling[J]. Ground Water, 1994,32(1):71-78.

[46]Sidauruk P, Cheng A H D, Ouazar D. Ground water contaminant source and transport parameter identification by correlation coefficient optimization[J]. Ground Water, 1998,36(2):208-214.

[47]Alabati S, Kabala Z J. Recovering the release history of a groundwater contaminant via the non-linear least-squares estimation[J]. Hydrological Processes, 2000,14(6):1003-1016.

[48]Neupauer R M, Borchers B, Wilson J L. Comparison of inverse methods for reconstructing the release history of a groundwater contamination source[J]. Water Resources Research, 2000,36(9):2469-2475.

[49]Lattes R, Lions J L. The method of Quasi-Reveribility, Applications to Partial Differential Equations[M]. New York, USA: Elsevier, 1969.

[50]Skaggs T H, Kabala Z J. Recovering the release history of a groundwater contaminant plume: method of quasi-reversibility[J]. Water Resources Research, 1995,31(11):2669-2673.

[51]Chavent G, Dupuy M, Lemommier P. History matching by use of optimal theory[J]. Society of Petroleum Engineers Journal, 1975,15(1):74-86.

[52]Yeh W W G, Sun N Z. Variational sensitivity analysis, data requirements, and parameter identification in a leaky aquifer system[J]. Water Resources Research, 1990,26(9):1927-1938.

[53]Neupauer R M, Wilson J L. Adjoint method for obtaining backward-in-time location and travel probabilities of a conservative groundwater contaminant[J]. Water Rescources Research, 1999,35(11):3389-3398.

[54]Neupauer R M, Wilson J L. Adjoint-derived location and travel time probabilities for a multi-dimensional groundwater system[J]. Water Resources Research, 2001,37(6):1657-1668.

[55]Atmadja J, Bagtzoglou A C. Pollution source identification in heterogeneous porous media[J]. Water Resources Research, 2001,37(8):2113-2125.

[56]Jury W A, Roth K. Transfer functions and solute movement through soil: theory and applications[M]. Basel, Boston, Berlin: Birk user Verlag, 1990.

[57]Datta B, Beegle J E, Kavvas M L, et al. Development of an Expert System Embedding Pattern Recognition Technique for Groundwater Pollution Source Identification[R]. Springfield, Virginia: National Technical Information Service, 1989.

[58]Bagtzoglou A C, Tompson A F B, Dougherty D E. Probabilistic simulation for reliable soulute source identification in heterogeneous porous media[J]. Water Resources Engineering Risk Assessment, NATO ASI Series, 1991,G29:189-201.

[59]Wilson J L, Liu J. Backward tracking to find the source of pollution[J]. Waste Management from Risk to Remediation, 1994,1:181-199.

[60]Liu J, Wilson J L. Modeling travel time and source location probabilities in two-dimensional heterogeneours aquifer, Las Cruces, New Mexico, 1995[C].

[61]Woodbury A D, Sudicky E, Ulrych T J, et al. Three-dimensional plume source reconstruction using minimum relative entrop inversion[J]. Journal of Contaminat Hydrology, 1998,32(1-2):131-158.

[62]Michalak A M, Kitanidis P K. Estimation of historical groundwater contaminant distribution using the adjoint state method applied to geostatistical inverse modeling[J]. Water Rescources Research, 2004,40, W08302, doi: 10.1029/2004WR003214.

[63]Butera I, Tanda M G. A geostatistical approach to recover the release history of groundwater pollutants[J]. Water Rescources Research, 2003,39, (12), 1372, doi: 10.1029/2003WR0023.

[64]Sun A Y. A robust geostatistical approach to contaminant source identification[J]. Water Resources Research, 2007,43, W02418, doi: 10.1029/2006WR005106.

A Review of Mathematical simulation methods for groundwater pollution source identification

LONG Yu-qiao1，2，CUI Ting-ting1，2，LI Wei1，LI Yan-ge1，WU Chun-yong1

(1.Department of Hydrology and Water Resources, Nanjing Hydraulic Research Institute, Nanjing 210029, China；2.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Nanjing 210098, China)

China has to face up to the groundwater resource crisis and the deteriorating groundwater environment. Reinforcing the studies on groundwater pollution source identification (GPSI) could be an important support of contaminate removing, groundwater protecting, potable water security, and society and economy development. Exploring the new theory and method of GPSI could push the studies on ill-posed problems, and improve the techniques of contaminate removing. GPSI has been studied for thirty years, and a brief review is given to conclude the characteristics of GPSI problems. The mathematical simulation method could be classified into four types: simulation-optimization method, analytical and regression method, direct method, and stochastic method. Each method has its advantages and disadvantages, further researches may focus on more complex GPSI problem, expressing physical chemistry and biological process, improving GPSI modeling efficiency, and model uncertainty.

groundwater；pollution；source identification；mathematical simulation methods

2016-09-28

国家自然科学基金项目(51409161；51509157)；江苏省自然科学基金项目(BK20140080)

龙玉桥(1984 -)，男，广东英德人，博士，主攻方向：地下水数值模拟研究。

P641.12

A

1004-1184(2017)01-0001-07