基于不同地质统计方法的渗透系数场对污染物运移的影响

岳松梅，杨 蕴，吴剑锋，吴吉春

(1.表生地球化学教育部重点实验室，南京大学地球科学与工程学院，江苏 南京 210023;2.河海大学地球科学与工程学院，江苏 南京210098)

研究大尺度地下水流、污染物浓度和运动规律问题时，区域性模型构建中的一个突出问题是勘探资料分布不均或缺乏，以及含水系统本身的非均质性，造成了水文地质参数的空间变异性。在通常的地下水数值模拟中只能用简单的参数分区来描述参数的非均质性，但由于分区过大而忽略分区内部参数的差异性，往往导致模拟结果的不确定性。寻求一种尽可能利用有限的勘探资料，对未知区域内的含水层参数进行合理估值的方法，是目前大区域地下水流模拟中的关键问题［1－2］。

在影响地下水系统的所有不确定因素中，渗透系数是表征含水层特性的一个关键要素［3－5］，它在空间上的分布变化及其复杂，其不确定性主要是由自然界含水介质的非均质各向异性和取样、测试中的失真以及试验量测误差造成的。这种空间分布的不确定性使我们将其以随机变量来处理，在不可能测得含水层每一点的介质特征和水力特征的前提下，引入含水层空间随机场，将实际含水层视为随机场的一个实现。但实际情况中，渗透系数不仅具有随机性，也具有一定的结构性［6－9］。地质统计学就是研究这种具有“二重性”的区域化变量的数学工具，其以变差函数为基本工具来研究分布于空间并呈现一定结构性与随机性的自然现象［10－11］，因而用地质统计学方法来研究水文水资源水环境系统参变量的空间变异性，可以定量地揭示水文水资源水环境系统参变量在空间不同方向上的变化规律，辨识各参数的最大变异方向，尤其是它能给出区域化变量的最优估计值及估计方差［10－11］，而这种特性是其他任何一种估计方法所不曾拥有的。

本文在前人研究的基础上［12］，将地质统计学推广到大尺度的地下水污染运移的研究工作中，结合GSLIB软件库［18］，利用美国麻省军事保留区(Massachusetts Military Reservation，MMR)的场地实测数据，采用普通克里格法和指示克里格法、顺序高斯模拟法和顺序指示模拟法等四种地质统计方法，插值估测和模拟再现含水层渗透系数随机场，进而对比研究四种渗透系数场对大尺度污染物运移的影响。

1 四种地质统计方法简述

1.1 普通克里格和指示克里格方法

普通克里格法属线性平稳地质统计学范畴，是地质统计学中最基本的方法，也是最常用的方法之一。设区域化变量Z(u)满足二阶平稳假设，其数学期望为未知常数，协方差函数和变异函数存在且平稳，以任一未知点为中心点的某一邻域范围内n个样本点，赋予相应权重计算未知点的属性值，形成估计领域内n个信息值的线性组合。通过求得一个无偏的最优估计，并给出每个估值的误差方差，该方法即为普通克里格方法(ordinary kriging，OK)［10－11］。

然而，在实际研究中还会碰到采样数据中存在特异值得问题，所谓特异值是指那些比全部数值的均值或中位数高很多的值，其既非分析误差所致，也非采样方法等人为误差引起，指示克里格法(indicator kriging，IK)就是为解决上述问题而发展起来的一种非参数地统计学法。该方法在给定条件数据的前下可以估计任意位置上的条件累积分布函数(conditional cumulative distribution function，CCDF)，适用于连续型和离散型的信息预测［10－11］。

相比于普通克里格，指示克里格法是将对区域化变量Z(u)的研究转化为对其指示函数的研究，它无需假设数值来自某种特定分布的总体，也无需对原始数据进行变换，可以估计落于某些区间的每一部分的值，也可以估计比某一门坎值高(或低)的那部分值。

1.2 顺序高斯模拟和顺序指示模拟方法

顺序高斯模拟(sequential Gaussian simulation，SGSIM)也称序贯高斯模拟，是将顺序模拟思想和高斯模拟相结合的一种随机模拟方法。它假设一个高斯随机域，基于变异函数特征，对每一个待模拟的空间位置利用克里格方法求构建高斯分布函数，之后随机在其中抽取一个数值作为该点的模拟值。

顺序指示模拟(sequential indicator simulation，SISIM)将指示克里格方法和顺序模拟算法相结合，实现了非参数的连续型与离散型的分布模拟。该方法不需要对原始条件数据分布的参数形式做任何假设，而是在现有资料的基础上，通过一系列的门槛值把条件数据转化成指示数据，相应的变差函数模拟转换为指示变差函数模型，并采用指示克里格法对每个网格点处的条件概率分布进行估计［14－16］。

2 污染物运移及其结果评估

研究地下水污染问题必须用两类方程来描述:第一类方程用来描述液体的流动;第二类方程用来描述溶质即污染物质在地下水中的运移。在一般或者较复杂的水文地质条件下，污染物运移只有通过数值方法来实现，而且无论是采用有限单元或者有限差分方法，都要通过地下水流模拟地下水的流速场才能由运移模型得到污染物运移的模拟结果。在这里，分别采用常用的软件水流模拟模型MODFLOW和溶质运移模型MT3DMS来进行模拟。

2.1 渗透系数场的产生

2.1.1 原始数据的基本统计特征

包括均值、方差、标准差、变异系数、偏态系数和峰态系数的计算与分析;画出原始数据分布直方图并分析数据是否满足正态分布特征。

根据原始数据分布直方图，依据“偏度、峰度检验法”对原始数据进行正态分布检验，进而判断原始数据的正态性。

2.1.2 原始数据的稳健性处理

地质统计学必须满足数据为正态分布这一要求，如不满足正态分布，直接计算的传统地质统计学变差函数必然有偏，进而影响估计的效果，故需对数据进行稳健性处理。分别对数据先后进行去类分析、特异值处理和正态变换，使数据接近正态分布。

2.1.3 渗透系数分布场的产生

本文利用GSLIB软件库，对渗透系数进行原始数据基本参数统计和数据稳健性处理，然后采用普通克里格法和指示克里格法、高斯序列模拟法和指示序列模拟法分别对数据进行插值和条件模拟，得到相应的渗透系数(lnK)分布场。

2.2 污染羽的空间矩评估

污染物在空间上的变化特征通常可以通过刻画污染羽随时间变化的空间矩来评价。污染羽在空间上的零阶矩表示污染物的总质量;一阶矩表示污染羽的质心位置;二阶矩表示污染羽围绕质心的分布范围。

污染羽在任一时刻t关于原点的零阶矩，即污染物总质量，可以表示为

污染羽在空间上关于质心的二阶矩可以表示为协方差变量U:

本文以MMR为实例，建立水流模型和溶质运移模型，利用空间矩评估污染羽的变化范围。由于模型考虑化学反应的两种基本类型:固液表面反应(吸附作用)和一阶速率反应，污染物的质量(零阶矩)在运移过程中会发生改变，因此我们考虑污染羽的零阶矩、一阶矩和二阶矩随渗透系数场的空间变异特征。

3 实例研究

3.1 算例概况

美国麻省军事保护区位于科德角(Cape Cod)的法尔茅斯镇(Town of Falmouth)附近，始建于20世纪初，面积大约为89 km2。该区为第四纪冰川沉积物覆盖，这些沉积物组成的非承压含水层为科德角几个社区的主要供水水源。位于MMR的东南部的美国空军的给养基地，由于长期的军事活动，出现了以三氯乙烯(Trichloroethylene，TCE)为主要污染物的化学物溢漏区(Chemical－Spill 10，CS－10)，对整个区域的环境造成严重破坏，直接影响了当地居民正常的生产生活。因此，本次研究就以CS－10区为模拟计算区域(图1)。

图1 研究区位置、CS－10模拟区域及TCE污染治理系统的平面分布［19］

3.2 模型的建立

本文建立了描述研究区含水层中TCE污染物运移的数学模型，采用有限差分法求解。其中只有含水层渗透系数为一随机变量，其他模型参数均为常数。水流模型中的渗透系数K是通过地质统计学对实测的值进行插值和条件模拟得到的，变化范围较大，在泥沙含水层中可以小于3 m/d，粗砂含水层中能大于91 m/d。在假定研究区水流呈稳定态，而溶质运移呈非稳定态的基础上，整个研究区域剖分为159列、161行和21层，总模拟面积为57 km2，总模拟时间为30 a，1个应力期。水平方向步长为34m，越接近边界步长会越大，垂直方向上每一层的厚度变化较大，最小可以小于1.5 m，最大可以大于15 m。侧向边界为定水头边界，由区域水流模型计算然后插值得到;上部为定流量边界，补给量变化范围为41～86 cm/a;下部隔水为零通量边界。该区内的地下水水流主要是沿水平方向从南向西南流动，水平的水力坡度平均值为0.001。

溶质运移模型是建立在MT3DMS的基础上的。该模型采用和MODFLOW水流模型相同的网格剖分，侧向和垂向上都是定通量边界。在边界上，不考虑弥散项而只考虑对流项，如果流量为0，则溶质运移也等于0。其余参数如，有效孔隙度为0.3，横向弥散度、纵向弥散度和垂向弥散度分别为11 m、1.1 m 和 0.11 m［19］。TCE 浓度超过 5μgL－1的污染羽三维分布图见图2所示，

图2 溶质运移模型运行前CS－10 TCE污染羽浓度的三维分布图

污染物在含水层中运移前初始污染物浓度均值为16 μg/L，最高浓度值远超过 4 690 μg/L［13］。在美国，建立在保护人类健康和环境的安全饮用水法规的条件下，TCE的最大污染标准(MCL)是5 μg/L，这说明研究区域被高浓度的TCE污染。

3.3 渗透系数分布场的产生

在对研究区已知的375个观测数据进行稳健统计分析，并在此基础上构造该区域化变量的几何各向异性套合模型，最后与GSLIB数据库相结合，采用四种不同的地质统计学方法分别对数据进行插值和条件模拟，得到相应的渗透系数(lnK)分布场。对于克里格插值来说，一旦各种参数确定(变差函数模型)以后，渗透系数(lnK)分布场就确定了，对应的污染物运移结果也确定了;而对于条件模拟来说，则每一次产生的渗透系数(lnK)分布场都是随机的，相应的污染物运移结果也是随机的。在这里，分别对克里格插值和条件模拟得出的结果做出对比，条件模拟模拟次数采用300次。

3.4 模拟结果与分析讨论

普通克里格和指示克里格法所得的渗透系数(lnK)均值在3.7左右，范围在－3～7之间，标准偏差在1左右，但指示克里格值(0.607 6)比普通格里格值(1.104 9)偏小。四种地质统计方法得到的渗透系数场(lnK)估计方差的统计参数详见表1。

表1 四种地质统计方法得到渗透系数场(lnK)估计方差的统计参数

由表1可以得出，两种克里格插值估计方差均值在1左右，指示克里格法略小于普通克里格法，但标准偏差比普通克里格法略大。图3分别为污染运移5 475 d(15 a)和10 950 d(30 a)后，克里格插值下污染羽浓度的三维分布图。

图3 克里格插值生成渗透系数场所对应的15 a和30 a后污染羽分布图(μgL－1)

从图3中不难看出，普通克里格法的污染羽扩散更加均匀，指示克里格所对应的污染物浓度极值范围更大。对于污染羽的二阶空间矩(表2)评估，普通克里格插值对应的污染羽在空间上的展布范围(二阶矩)明显大于和指示克里格，这主要受渗透系数的空间变异方差的影响，渗透系数方差的偏大，污染羽在空间上的二阶矩随之升高，表明污染物围绕中心的扩散范围不断加大。这些不同之处是由于指示克里格法是用个别范围个别变异函数来刻画该区域化变量的关联性的特性，无需假设数值来自某种特定分布(如正态分布)的总体，也无需对原始数据进行变换(如对数变换)，因此指示克里格法不必去掉重要而实际存在的高值数据即可处理各种不同现象［12］，而普通克里格法表现得更圆滑效应，渗透系数估计方差大且平滑，进而污染羽扩散更加均匀。

估计法是用加权平均法对研究区域的未知量求得线性、无偏和最小方差的内插估计量及其相应的估计精度，而顺序模拟法是通过系列随机模拟对地质变量进行局部估计，克服了估计法的平滑，能够更好的再现真实曲线的波动性。对于条件模拟每一次产生的渗透系数(lnK)分布场都是随机的，相应的污染物运移结果也是随机的，这里采用300次模拟次数。

顺序高斯模拟和顺序指示模拟(模拟次数300次)所得的渗透系数(lnK)均值相同，即3.7997，范围均在－3～7之间，但标准偏差顺序指示模拟值(0.5941)比顺序高斯模拟值(0.758 2)偏小。这两种方法得到渗透系数场(lnK)估计方差均值在1左右，标准偏差顺序高斯模拟值略小，详见表1。

对比克里格插值得出的结果，四种地质统计法所得到的渗透系数(lnK)均值在3.7左右，所对应污染羽扩散范围基本一致，污染羽的质心坐标也基本相同，验证了质心位置是由渗透系数的平均值来决定的理论。图4分别为污染运移5 475 d(15 a)和10 950 d(30 a)后，条件模拟(模拟次数300次)下污染羽的浓度分布。表3为条件模拟模拟次数300次所对应的15 a和30 a后污染羽在空间上的分布特征。

综合图4和表3，不难看出，顺序指示模拟的高浓度范围明显高于顺序高斯模拟法，但污染羽的二阶空间矩略小于后者。

表2 克里格插值方法所对应的15 a和30 a后污染羽在空间上的分布特征

图4 条件模拟次数300次模拟生成渗透系数场所对应的15 a和30 a后污染羽分布图

两种模拟方法在模拟次数增加的情况下，污染羽的扩散范围与质心坐标没有明显差别，这有力的验证了质心位置是由渗透系数的平均值来决定的理论;但污染羽的空间二阶矩在某一方向上却有明显减小的现象，这与二阶空间矩主要受渗透系数的空间变异方差影响的结论相呼应。单独对比顺序高斯模拟与顺序指示模拟这两种方法，在模拟次数增加的情况下，渗透系数估计方差与污染物空间矩评估值这两方面，后者的变化程度明显高于前者。这是由于顺序指示模拟采用指示克里金来估计，指示克里金不同于其他克里金法，它并不依赖于变量的平稳性和要求变量服从某种分布。顺序指示模拟法可以模拟复杂各向异性的地质现象及连续性分布的极值，对于具有不同连续性分布的变量，可利用不同的变差函数进行表征，建立各向异性的模拟结果，进而更能体现数据的差异性。所以，采用顺序指示模拟方法模拟的污染羽分布，其高浓度范围明显高于顺序高斯模拟得出的结果，并且空间二阶矩的变化程度明显更大。

表3 条件模拟模拟次数300次所对应的15 a和30 a后污染羽在空间上的分布特征

4 结语

(1)本文利用野外场地实测数据，采用四种地质统计方法，插值估测和模拟再现随机渗透系数场，进而对比研究四种渗透系数场对大尺度污染物运移的影响。研究结果表明，污染羽一阶空间矩证明了污染羽的质心坐标基本相同，而四种地质统计法所得到的渗透系数(lnK)范围均在－3～7之间，均值在3.7左右，这有力的验证了质心位置是由渗透系数的平均值来决定的理论;通过指示克里格和顺序指示模拟法得到的渗透系数(lnK)标准偏差比其余两种更小，污染羽空间二阶矩偏小，这也与二阶矩主要受渗透系数的空间变异方差影响的结论相呼应。

(2)四种地质统计方法在处理样本数据时各有利有弊，应根据样本的特点来作出相应对策。指示克里格与普通克里格得出的结果最大不同之处在于，前者污染羽的极值特点更加明显，后者污染羽扩散更加均匀。条件模拟追求的是模拟的真实性，克服了估计法的平滑效果，能较好地再现真实曲线的波动性;随着模拟次数的增大，渗透系数(lnK)估计方差与污染羽空间二阶矩呈现变小的趋势，并且顺序指示模拟程度更加明显，这是由于顺序指示模拟在模拟复杂各向异性的地质现象及连续性分布的极值上更加突出。

(3)需要说明的是，本文进行的污染物运移模拟是建立在四种地质统计方法估计模拟渗透系数场的前提下实现的，不可避免得会出现估计误差，除此之外，在变差函数理论模型拟合中存在必然的人为误差和模型误差。用数值法计算中也存在误差，模型只能计算位于网格单元中心的结点处的浓度，而监测孔的实际位置大多不是位于有限差分网格单元的中心。

［1］Kai－ Wei Juang，Dar－ Yuan Lee，Timothy R.Ellsworth.Using rank－order geostatistics for spatial interpolation of highly skewed data in a heavy － metal contaminated site［J］.Journal of Environmental Quality，2001，30:894－903.

［2］Bastante F G，Ordóez C，Taboada J，et al.Comparison of indicator kriging，conditional indicator simulation and multiple－point statistics used to model slate deposits［J］.Engineering Geology，2008，98:50－59.

［3］Diana dell＇Arciprete，Riccardo Bersezio，Fabrizio Felletti，et al.Comparison of three geostatistical methods for hydrofacies simulation:a test on alluvial sediments［J］.Hydrogeology Journal，2012，20:299－311.

［4］Freeze R A.A stochastic conceptual analysis of one－dimensional groundwater flow in non － uniform homogeneous media［J］.Water Resources Research.1975，11(5):725－741.

［5］施小清，吴吉春，袁永生.渗透系数空间变异性研究［J］.水科学进展.2005，16(2):210－215.

［6］陈彦，吴吉春.含水层渗透系数空间变异性对地下水数值模拟的影响［J］.水科学进展.2005，16(4):482－48.

［7］刘文婷，朝伦巴根，刘艳伟，王力，宋君.地质统计学方法在渗透系数空间变异性研究中的应用［J］.水利科技与经济.2010，16(4):364－366.

［8］李少龙，张家发，杨金忠.渗透系数空间相关性对渗流场统计特征影响分析［J］.长江科学院院报.2009，26(10):17－20.

［9］闫婷婷，吴剑锋.渗透系数的空间变异性对污染物运移的影响研究［J］.水科学进展.2006，17(1):29－36.

［10］张仁铎.空间变异理论及应用［M］.北京:科学出版社.2004.

［11］刘爱利，王培法，丁园园.地统计学概论［M］.北京:科学出版社.2012.

［12］刘玲玲，吴剑锋，吴吉春.不同地质统计方法在确定渗透系数场中的对比研究［J］.水文地质工程地质.2009，5:66－71.

［13］吴剑锋，彭伟，钱家忠.基于INPGA的地下水污染治理多目标优化管理模型［J］.地质论评.2011，57(3):437－443.

［14］李少华，张昌民，王振奇.利用顺序指示模拟方法预测储集层岩性［J］.新疆石油地质.2002，23(1):59－62.

［15］宋永忠，于晶，年静波.随机顺序指示模拟技术与应用［J］.大庆石油地质与开发.2003，1:52－54.

［16］余振，何静，魏福吉.序贯指示模拟和序贯高斯模拟在某地区精细流体预测中的联合应用［J］.天然气地球科学.2012，23(6):1170－1174.

［17］施小清，吴吉春，吴剑锋.多个相关随机参数的空间变异性对溶质运移的影响［J］.水科学进展.2012，23(4):509－515.

［18］Deutsch C V，Journel A G.GSLIB－Geostatistcal Software Library and User’s Guide［M］.2nd ed.New York:Oxford University Press，1998.

［19］Zheng Chunmiao，Wang P P.2002.A field demonstration of the simulation－ optimization approach for remediation system design.Ground Water.40(3):258 ～265.