确定某一污染羽分布的佳地下水监测网络.doc

2007年第3期 确定某一污染羽分布的最佳地下水监测网络 25确定某一污染羽分布的最佳地下水监测网络*Kyung-Ho Kim Kang-Kun Lee（School of Earth and Environmental Sciences,Seoul National University, Seoul 151-742, South Korea）翻译：葛秀珍；校对：苑惠明【摘 要】 本文提出了一种新的优化地下水监测网络，确定污染羽范围的方法。将监测井安装后量化的污染物存在的不确定性期望值降低到最大的位置定为监测井的最佳位置。本项研究中，将水力传导系数作为诱发不定性的因素。使用连续的随机添加(SRA)法生成水力传导系数的随机场。随着监测网络范围的增大，污染羽分布的不确定性降低，根据这个降低的量，评价某一污染羽存在的信息判定的期望值。选择采集监测井最大信息量的最小系统为优选的监测网络。为了量化污染羽分布的不确定性，在定义域范围内，针对所有的产生的污染羽的实现值，编制污染物存在的概率图。不确定性定义为，污染物存在的概率或者不存在的概率单元的总和。在非均质水力传导系数场，本文给出了确定最佳监测网络的数值试验的结果。【关键词】 优选 监测网 不确定性 污染羽 1 引言制定污染的地下水的修复计划需要有关污染羽分布的资料。如果确定的污染羽不合适，则很难制定出节省成本的修复方案，但是，如果观察资料不足或者过多，也很可能造成设计设备的安全系数过大或者经费超支。*Optimization of groundwater-monitoring networks for identification of the distribution of a contaminant plume以前就用各种方法监测过污染羽，并且进行过一些监测井的优选安装工作。Meyer和Brill(1988)在不确定的条件下，使用地下水监测网络中定井方法监测污染物。Hudak和Loaiciga(1992)提出了定位-模拟的方法，该方法利用覆盖区增大了某一地下水监测网。Hudak和Loaiciga(1993)针对多层地下水流系统，提出了一种优化的方法。James和Gorelick(1994)评估了安装监测井，修复污染的含水层的成本期望值，为减少总的修复成本，提出了最佳位置和需要监测井的数量。Storck等(1997)提出了确定优选三维非均质的含水层监测井的位置的方法。Warrick等(1998)对监测地下污染羽的地下水监测系统的效果进行了概率分析。Bogaert 和Russo (1999)根据最小二乘法，提出了优化变量评估的空间采样设计的选择方法。Reed等(2000)研制了降低长期监测地下水污染场地经费的方法。Wu等（2005）开发了大规模污染羽长期监测的有成本效益的采样网络。Bierkens(2005)利用随机模拟地下水流和对流传输，介绍了监测地下水污染的监测网络的优选方法。以前已经讨论了监测网的优选问题，主要是最有效的监测潜在的污染源附近的污染物。然而，当监测到污染物时，列举的使用方法在编制修复污染场地的计划中是受限制的，因为优化了的监测网络可以有效地监测污染物但不能描述污染羽的位置和边界。本研究的目的是提出和论证优化监测网络的方法，最有效地降低地下水系统污染羽分布的不确定性。因此，利用地下水流和污染物传输模拟，在随机传导场产生可能的污染羽实现值。编制污染物存在的概率图，查明最佳的位置，最大地减少不确定性。2 方法2.1 序列安装模型根据最近安装井得到的信息，利用以前安装井修改的资料，序列安装模型可以确定监测井的最佳位置。这个序列是可以重复的。当根据新监测井的安装获得的量化信息量小于以前的标准设置，这个网络的扩充就要停止。这个模型的流程图见图1。26 水文地质工程地质技术方法动态 2007年第3期图 1 序列安装模型的流程图 第一步是采集预安装井的数据和评估野外参数。野外参数可以根据采集的数据和其他目标区域有关的文献资料进行评估。以前序列中获得的数据将用于参数场数字生成的野外参数。第二步是生成污染羽实现值。在这个过程中，用产生的参数场模拟地下水流和污染传输。第三步编制污染物存在的概率图。使用这份概率图，可以确定某一新的监测井的最佳位置。在这个位置，如果采集信息的期望值(EVSI)高于标准值，那么就可以安装新的监测井。否则，整个过程，包括井的安装就要被停止。在以前的实例中，数据是采自某一新近安装的井，使用这样的数据，要排除与新井污染数据相矛盾的污染实现值。接下来，要匹配井的数据，修改污染位置的概率图。这个顺序要重复进行。更进一步的工作在下一节中叙述。2.1.1 生成水力传导系数场在地下水流模拟中，存在许多不确定的参数。本项研究中，水力传导系数是就一种不确定是参数。根据均方差和协方差结构生成非均质水力传导系数场。这些非均质水力传导系数场可以从某一参数评估过程获得。连续的随机添加(SRA)法用于生成随机场(Kentwell等1999)。2.1.2 生成污染羽实现值通过蒙特卡洛法模拟，使用以前过程中形成的几组水力传导系数场得到污染羽实现值。首先，在生成的水力传导系数场上模拟地下水流场。其次，用水流模拟计算的速度场模拟溶解物污染浓度的分布。本项研究中，由于引起不确定性的单一因素是水力传导系数，所以假定在水流和溶解物传输模拟的主要区域，边界条件和野外参数值是确定的。然后，经过蒙特卡洛模拟，在任一一个水力传导系数场的任一一个位置获得浓度分布数据。第三，将浓度数据转换为二进制数据。因此，修改的设计者必须建立转换的标准浓度。如果在某一位置的浓度高于建立的标准。那么这个位置将被假定为存在污染物；否则就认为这个位置没有污染物。使用这些方法，在任一一个场的任一一个位置，浓度数据都能够转换为二进制数据。2007年第3期 确定某一污染羽分布的最佳地下水监测网络 272.1.3 绘制污染物存在的概率图在这个过程中，转换的污染物存在的二进制数据将用于生成污染物存在的概率图。这个概率图不仅能用于污染物分布分析而且还能计算某一监测井的最佳位置（见下一节的描述）。污染物在某一位置（x，y）存在的概率定义为：式中，p(x, y)是位置（x，y）污染物存在的概率；n(x, y)是在某一位置（x，y）污染物存在场地的数量；N是生成污染羽实现值的总数。用方程（1）计算污染物存在概率时，如果修复设计需要单一的污染羽分布，那么p(x, y)=0.5的等位线就可以充当这条表示的线。2.1.4 计算某一监测井最佳位置利用以前的方法绘制的污染物存在的概率图能通过下面方程获得某一监测井最佳的位置。用这个方程计算最佳位置，应该量化不确定性。某一监测井的最佳位置假定为某一监测井安装后，将量化的不确定性降低到最小的位置。这个量化的不确定因素f定义为：式中A表示区域的面积，f(x,y)表示定量污染物分布不确定性的函数。这个函数可以用因子的各种形式表示，即进行不确定性计算。本研究中，用两个模型确定不确定性。第一个模型是面积模型，设想污染羽分布的面积是不确定的。换言之，用模型1中的f确定污染物存在的概率既不是0，也不是1的场地面积。在这个模型中，确定污染物分布f(x, y)的不确定性的函数用下式表示：将方程（2）中的f（x, y)带入方程(3)，能得出f的值。图2说明了直角坐标系中，区域的面积（A）和概率分布函数p(x, y)。如果修复设计者欲忽略污染羽实现值的极端的情况，那么某一置信区间就可以用极限概率定义。这个排除区域超出了置信区间。如果把这个概率应用到方程（3），那么就能表示如下： （4）图 2 直角坐标系区域的面积（A）和概率分布函数p(x, y)，暗区是高概率区，亮区是低概率区式中C1 和C2 是置信区间边界值的概率。第二个模型（量化污染物分布不确定性）是加权区域模型。这个模型产生一个确定污染物分布不确定性的函数f(x, y)和加权因子，该加权因子使得大的不确定因素区域量化的不确定性增大。由于在概率p(x, y)为0.5的位置，污染存在的不确定性为最大，在概率为0或者1的位置为最小，所以用加权区模型确定污染分布不确定性的函数定义如下：28 水文地质工程地质技术方法动态 2007年第3期因为污染物存在的不确定性是根据污染物存在的概率p(x, y)而定，所以，方程（5）与p(x, y)有关联的。当p(x, y)0.5时，函数值最大，而且这个函数是两个线性函数的组合，这两个线性函数当概率为0.5时，都有一个峰值。因此，如果使用这个模型，当概率大约等于0.5时，能得出着重说明的不确定性的结果。用同样的方式使用面积模型，还能忽略加权模型中污染羽实现值的极端情况。方程（6）提出了应用于方程（5）的可靠性。综上所述，作者提出了确定污染物分布不确定性的两个模型，修复设计者可以根据情况选择适当的模型。利用选择量化不确定性的计算模型，就可以用下面方程得到采集信息的期望值(EVSI)：式中，x为位置矢量，定义域为xC，C是候选井点的组数。如果C定义为C=x|xA，那么目标函数可能就没有唯一的解。因此，C应该是x|xA和n(C)的真子集；f为方程（2）计算的不确定的因子，这一项表示使用安装井信息的不确定性的量；f为新的监测井在X处安装后由方程（2）计算的不确定性因子；P(x)为x处污染存在的概率；i(x)为说明x处污染物存在或者不存在的函数。假定f的采样信息与f的相同。换言之，对新近安装的监测井，在以前采样和下一次采样之间的这段时间内，假定污染羽是停滞的。所以，在f和f的采集时间间隔与污染羽迁移的速度相比相对较短的情况下，这个方法是适用的。恒等式i(x)=表示位置x处某一监测井安装之后污染物存在，i(x)=0表示上述情况污染物不存在；因此，f(i(x)=1)表示仅在位置x处，某一监测井安装之后，污染物存在的不确定因子。所以这个f(i(x)=1)可以根据污染物分布的概率图，用污染物存在的位置x对应的污染羽实现值获得。同样，f(i(x)=0)也可以根据概率图，利用污染物不存在的位置x对应的污染羽实现值获得。E（f）是根据x处监测井安装之后的数据信息得到的不确定因子f的期望值。根据贝叶斯定理，用概率的总数乘以随机变量f得到E（f）。因为EVSI在不同获得途径（根据现在给定的信息获得，或者根据监测井安装后期望的不确定因子获得）上是有差别的，EVSI值是x处（安装监测井）不确定性期望值降低的量。因此，EVSI最高的位置就是某一新监测井安装的最佳位置，因为这个位置要求不确定性降低的量最大。如果定义域是不连续的，最佳的位置可以通过简单搜索定义域内所有的节点找到。按照这样结果进行，监测井的最佳位置就可以确定。然而，如果EVSI低于修复设计者以前设置的标准值，那么新的监测井就不能安装，而且整个过程都要停止。另一方面，如果不是这样，EVSI高于标准值，那么新的监测井就可以在EVSI最大的位置安装。接下来，从新的位置进行数据采集，删除与数据采集结果不符合的污染羽实现值。利用剩余的污染羽实现值，修改污染物存在的概率图。这个过程在以前的部分提到过，后面的过程通过图1中流程图重复。2.2 同步安装模型同步安装的模型利用从预安装井获得的信息和修改最近安装的监测网得到的信息确定优选的监测网络。这个模型不同于那些安装方法中重复的序列模型。模型使用的流程见图3。这个模型的安装顺序比序列安装模型简单，但是EVSI的计算函数不同。2.2.1 优选计算网络的计算在同时进行的多级井安装的工作中，根据贝叶斯定理，用复合的方法推广了EVSI的计算。本文有两个单井的实例：污染物存在或者不存在，而n眼井的实例数量为2n。多级井的EVSI用综合公式计算。EVSI的函数定义如下：2007年第3期 确定某一污染羽分布的最佳地下水监测网络 29式中的k为表示定义域内位置节点数；xk为xk处污染物存在的事件中，k个节点数所在位置的位置矢量。利用方程（8），EVSI用下式表示：ij，XiXj(i,j = 1,n) 图 3 同步安装模型的流程图式中x1, x2, x3,. , xk 为定义域内或者x1, x2,x3,., xk; xiC 处污染物存在的位置矢量；C为候选井位置的集合。如果C定义为C=xi|xi 2A，则目标函数可能没有唯一解。因此，C应该是xi|xiA和(C)的真子集。f为方程（2）计算的不确定因子；这项是根据安装井的信息计算的不确定性的量；f是监测井安装之后，x1, x2, x3,. ,xk处，由方程（2）计算出的不确定因子。P(x)和i(x)与序列安装模型中的一样。在方程（9）中，EVSI最大时对应的位置是监测网的最佳位置。计算各种监测网的EVSI并且利用遗传运算法则(GAs) (Holl和1975)获得最佳位置。分别用MODFLOW-96(Harbaugh和McDonald 1996)和MT3DMS (Zheng和Wang1999)进行地下水流和溶解物模拟。3 数字示例3.1 数字模拟的定义域利用数值试验，本研究中提出的方法应用于两个简单的实例问题。假定实例问题的定义域为2D非承压含水层系统（图4）。假定2D定义域范围由35列和25行组成；大小为420m300m，30 水文地质工程地质技术方法动态 2007年第3期每一个网格的规格为12m12m。左边(hL=26m)和右边(hR=25m)设为定水头。上边和下边设为不渗透边界。含水层的厚度和有效空隙度分别为30和0.3m。固定含量污染源位于中央区(x,y)= 8,13)。释放周期是10年，浓度为100 mg/l。图4中的X表示污染源的位置。假定这种污染是保守的。3.2 实例: 序列安装模型为提供蒙特卡洛模拟的污染羽实现值，产生了一组水力传导系数随机场。使用SRA方法随机场生成了2500个水力传导系数场。水力传导系数的对数、标准偏差和H分别是3.4log(cm/s), 0.5log(cm/s)和0.75（图5）。图 5 水力传导系数随机场集合暗区和亮区分别表示低和高的水力传导系数在水力传导系数随机场生成了地下水流矢量场。本项研究中，利用MODFLOW-96(Harbaugh 等1996)进行地下水流数字模拟。利用传输模拟方法模拟各个水流场污染浓度分布。纵向弥散率设定为6.344m。Xu和Ecksteins(1995)关系式如下：式中aL为分散率差，L是野外规模。本研究中，横向弥散性与纵向弥散性的比为0.1。忽略了分子弥散。用MT3DMS(Zheng和Wang 1999)模拟溶解物传输（图6b）。根据污染物存在的概率图，通过溶解物传输，用污染羽实现值将每一个污染浓度分布图都转换为二进制数据。本项研究的传输的标准浓度为1mg/l。换言之，假定在污染物存在、污染浓度大于1mg/l的地方和污染物不存在且污染物浓度小于1mg/l的地方进行传输数据的二进制转换(图6c)。图 6 污染羽实现值生成和转换过程 a 地下水流模拟.暗区表示地下水位高区，亮区表示地下水位低区.实线表示地下水位等值线.b 溶解物模拟.暗区表示高浓.亮区表示低浓度.c 浓度数据转换位二进制.暗区表示污染物存在区域，亮区表示污染物不存在区域.2007年第3期 确定某一污染羽分布的最佳地下水监测网络 31图 4 数字模拟的定义域。X表示水源的位置 (x,y)=(8,13)利用浓度分布图转换二进制数据，绘制污染物存在的概率分布图；图7b就说明了这个结果。 本研究中，将从2500个随机场中选择的随机场假定为有效的污染分布，根据这个假定，模拟某一污染羽(图7a)并且依据这个有效的分布，每一个优选的监测网络安装之后，在每一眼井的位置确定污染物的存在。利用某一绘制的概率分布图，计算定义域内每一个网格单元的EVSI，然后确定EVSI最大的最佳的监测井的位置。利用SRA方法生成随机场，图7中表示了污染物存在的概率分布，根据评估的有效的污染分布，用50很大的概率发散度等位线估算污染羽。然而，安装5个井之后，由于不确定性降低量很大，50的概率等位线与假定的有效污染分布非常相似(图7c)。在图7d的图表中，随着安装井的数量的增多，EVSI最大值便快速降低。图 7 a假定的有效污染分布.浓度数据转换为二进制形式.暗区表示污染物存在区域，亮区表示污染物不存在区域；b 监测井安装前，污染物存在概率图，彩色定标线条表示概率值，深黑色实线是50概率的等值线；c序列安装5眼监测井之后，污染物存在的概率；.彩色定标线条表示概率值，深黑色实线是50概率的等值线；5个黑点是监测井的位置；d根据监测井是数量得到的最大的EVSI和污染羽实现值的数量，点线表示标准的EVSI线；当最大的EVSI低于这个线时，另外的监测井的安装就要停止32 水文地质工程地质技术方法动态 2007年第3期图 8 根据SRA方法获得的，为适合某一非均质水力传导系数集合而设计的优选监测网和监测井 安装以前污染物存在的概率图。彩色定标线条表示概率的值，深黑色实线是50概率等值 线，4个黑点表示监测井选择的优选的位置，X表示污染源的位置；a 方程（3）确定的 不确定性的实例，b 方程（4）确定的不确定性的实例， c 方程（5）确定的不确定性的实例， d 方程（6）确定的不确定性实例3.3 实例: 同步安装模型在同步安装的模型中，采用序列安装模型的定义域、条件和方法进行数字试验(图5,6)。用同步安装模型计算EVSI，根据生成的污染羽实现值，利用污染存在的概率图查明最佳的监测井的位置。对实验情境中的每一个实例，监测井的数量设为4。图8说明了最佳监测井的位置。利用4个不同的方法计算污染物存在的不确定性。将图8b中的d、c1=0.02 和c2 =0.98分别代入方程(4)和(6)。在图8a、b是用面积模型确定不确定性的实例，因为污染物存在的概率低于10的区域范围非常大，4个优选监测井中有3个和2个区域分别位于概率低于10处。然而，图8c、d中使用的是加权面积模型。在这些例子中，4个监测井中的两个和一个分别位于污染物存在概率高于10处。监测井位于概率高于面积模型的区域。在图8a、c中，最佳的监测井所在位置的污染物存在概率只比图8b、d高一点，可靠性也是一样。监测井最佳的排列分布在能获得各个方向污染物分布信息的各个方向上。在下节中，通过与各种监测网对比，评估监测网排列的效率（表1）。表 1 每一个监测网确定的不确定性对比2007年第3期 确定某一污染羽分布的最佳地下水监测网络 33序列模型同步模型网络 #1网络 #2网络#3网络#4面积模型34.22863.75486.602110.463104.73174.553面积模型 （可靠性使用）26.40035.79843.93542.58743.78840.341加权面积模型4.9416.1647.7787.3777.2