基于贝叶斯模型的地下水风险源污染概率估计方法研究

李 璐, 殷乐宜, 牛浩博, 刘伟江, 陈 坚*

1.生态环境部环境规划院, 长江经济带生态环境联合研究中心，北京 100012 2.生态环境部土壤与农业农村生态环境监管技术中心，北京 100012

我国地下水环境风险源点多面广，且地下水污染具有一定的隐蔽性，一旦发生污染，修复治理难度大且成本一般较高，但风险源周边地下水监测水平较低，因此，如何准确识别地下水污染来源对于地下水污染修复治理责任认定具有重要意义. 地下水污染源识别一般是基于一定的调查或监测数据，推断污染源的位置或污染路径的过程[1]. 如何通过较少的监测数据来识别污染来源，已经成为近年来国内外学者重点研究的领域.

目前，地下水污染源识别的主要方法包括地统计学法[2]、地球物理探测法[3-5]、同位素溯源法[6-7]和模型反演法等[8]. 其中，地统计学法要求多个位置甚至是长时间序列的监测采样结果，识别成本较高，一般用于已知污染源的地下水污染程度分析；地球物理探测法是基于地下水污染前后场地物性差异导致的物理场观测值的变化进行污染分布情况分析的方法，但一般仅适用于与岩层存在较大物性差异以及30 m深度以内的污染探测，且随着精度要求越高，可探测范围越浅；同位素溯源法目前已被广泛应用于环境污染事件的仲裁、环境污染物的来源分析与示踪中，但该方法只能对存在一定程度的特定无机物或有机物污染信息的污染源进行识别，且需要大量同位素和水化学监测数据；模型反演方法中应用比较广泛的求解方法主要包括模拟-优化法[9-11]及不确定性分析方法等，其中，不确定性分析方法因对监测数据需求量较小且可获得相对可靠的反演结果，被越来越多的专家学者用于环境污染问题的反演识别，主要包括随机游走粒子方法[12]、最小相对熵方法[13]、伴随状态方法[14]、贝叶斯统计方法[15-18]等.

贝叶斯统计方法以从监测值中获取参数信息为目标，将参数先验概率密度函数与样本似然函数结合在一起，形成一套非常灵活、直观、易于理解的统计方法；此外，通过与其他分析方法联用，还可以解决污染源贡献率计算、反演精度提升等问题[19-20]. 张双圣等[21-22]将贝叶斯模型与地下水二维水质对流-扩散方程相耦合，建立了地下水污染源参数反演模型，实现了基于8个监测点对1个污染源排放过程的反演. 然而，现实条件下很大概率上存在监测数据的严重缺失，甚至只有单个监测点一期数据的情况，为解决在监测数据较少但存在多个地下水风险源的情景下，如何判定监测数据异常的来源及风险源污染概率如何计算的问题，该研究提出了基于贝叶斯公式的地下水风险源污染概率估计研究方法，并开展案例分析，以期提供一种地下水污染来源判别方法的有效途径，对地下水污染风险防控具有重要科学意义.

1 研究方法

1.1 贝叶斯概率计算公式

该研究利用贝叶斯公式，计算造成监测点指标异常是由某一个风险源造成的概率，并将计算得到的概率最大值作为最优解，从而形成地下水风险源污染概率估计方法.

贝叶斯模型可以表达为

P(A|B)=P(A)×P(B|A)/P(B)

(1)

式中：P(A|B)为后验概率，代表在B监测点m指标异常是由A风险源造成的概率；P(B|A)为似然度，代表A风险源造成B监测点m指标异常的概率；P(A)为A风险源的先验概率，代表A风险源发生m指标异常的概率；P(B)为标准化常量，代表B监测点观测到m指标出现异常的概率.

1.2 先验概率修正

对于某个风险源来说，其是否发生污染取决于风险源是否排放目标污染物，而与B监测点指标是否异常无关，在相关参数未知的情况下，风险源排放目标污染物和不排放目标污染物的概率相等，即可将初始先验概率〔p0(Si)〕设为0.5. 然而，实际上各污染源性质并不相同，若确定风险源排放目标污染物时，可将p0(Si)设置为1.0；若确定风险源不排放目标污染物时，可将p0(Si)设置为0；对于不同目标污染物来说，p0(Si)也各不相同. 除此之外，还需要利用风险源的原材料、生产工艺、产品、存在时间、废水排放特征、环境保护设施运行情况等相关软数据信息进行先验概率修正. 对于Si风险源来说，修正后的先验概率〔p(Si)〕与初始先验概率〔p0(Si)〕、污染释放可能性系数(Li)和污染物可能释放量系数(Qi)有关，即

p(Si)=p0(Si)×Li×Qi

(2)

式(2)中先验概率修正时所需数据通常来源于场地调查、环境统计、专家建议和文献资料等. 该研究基于《地下水污染防治分区划分工作指南》[23]中所述污染源荷载评价办法，提出了利用风险源存在时间和废水排放量等进行Li和Qi估算的方法，其中，Li与风险源存在时间有关，存在时间为(0,5]、(5,10]、(10,20]、(20,30]、>30 a时，Li可分别按照0.1、0.2、0.5、0.8、1.0取值；Qi与废水排放量有关，废水排放量为(0,1×104]、(1×104,1×105]、(1×105,5×105]、(5×105,1×106]、>1×106m3/a时，Qi可分别按照0.2、0.4、0.6、0.8、1.0取值.

若风险源为废弃地块，其废水排放量可能为0，无法按照上述计算方法进行先验概率修正，因而提出了基于《地下水污染源强评价、分类与防控技术研究》[24]中点源源强计算方法以及风险源特征的Li和Qi的计算方法，其中，Li与风险源的防渗措施运行情况有关. 若风险源已开展防渗，防渗措施运行(0,1]、(1,5]和>5 a时，Li可分别按照0.2、0.6、0.8取值，若未开展防渗的可将Li设置为1.0；Qi与渗水面积有关，渗水面积为(0,1×103]、(1×103,1×104]、(1×104,1×105]、(1×105,1×106]、>1×106m2时，Qi可分别按照0.2、0.4、0.6、0.8、1.0取值.

1.3 似然度计算

似然度的计算一般涉及风险源与监测点的向量距离、水头、流向、孔隙度、渗透系数、弥散度、渗漏浓度、监测点观测值等参数. 似然度也是一种对计算值和观测值之间接近程度的度量，计算值和观测值越接近，则似然度的值就越大，反之亦然. 若已知多点或多次观测值，可通过极大似然思想和溶质运移方程求解似然函数[25-26]，但受监测条件限制而无法准确求解溶质运移方程时，也可通过对流弥散方程对似然度进行简化求解.

假设研究区地下水流可概化为一个等厚的各向同性均质含水层稳定流，根据流场特征可知，似然度p(Si,m)可以表示为

p(Si,m)∝cosα×Δhi/Li2

(3)

也可以表示为

p(Si,m)/W=cosα×Δhi/Li2

(4)

式中：α为Si风险源和监测点连线与地下水流向的夹角；Δhi为水头差，m；Li为Si风险源与监测点的距离，m；W为常数.

1.4 标准化常量及后验概率计算

标准化常量为归一化的积分常数，可以表示为

(5)

也可以表示为

(6)

则后验概率公式可以表示为

(7)

根据式(7)计算得到的后验概率最优解则为可能造成目标污染物超标的污染源[1].

2 案例分析

2.1 研究区概况

该研究选取石家庄市某工业集聚区作为研究对象，该区域位于滹沱河冲洪积扇平原地带，地势平坦，属华北暖温带半湿润地区，受大陆性季风气候影响，年均气温12.2 ℃，年均降水量474.0 mm. 该工业集聚区于2000年底批复建设的装备制造业及重化工产业基地，入驻了包括无机盐制造、化学农药制造、其他合成材料制造、其他基础化学原料制造、有机化学原料制造等行业类别在内的多家企业，目前部分企业已停产.

通过2019年4月现场调查采样发现，该区域东南部某一农灌井中，Cr6+含量超过GB 5749—2006《生活饮用水卫生标准》指标限值及GB/T 14848—2017《地下水质量标准》Ⅲ类标准限值3倍多，超过GB 5084—2005《农田灌溉水质标准》指标限值1倍多；CHCl3有检出，但其含量未超过GB/T 14848—2017《地下水质量标准》Ⅲ类标准限值. 为查找该监测点Cr6+和CHCl3含量异常的原因，拟根据笔者所建立的风险源污染概率估计方法进行分析.

基于收集到的该工业集聚区内企业的建成时间、废水排放量、原材料、产品等相关软数据，得到研究区域内8个可能造成污染事件的风险源(记为S1～S8)，风险源及监测点(记为O1)分布位置见图1. 研究区含水层岩性主要为中砂、黏土互层，渗透系数为15～20 m/d，含水层平均厚度为18 m，地下水由西北流向东南，与区域多年水位方向一致.

2.2 风险源先验概率计算

通过对风险源相关软数据进行分析，与Cr6+含量和CHCl3含量具有极强相关性的风险源为S6和S3，其中，S6风险源的原材料包括铬矿，铬矿及其矿渣中所含Cr6+极易通过雨水淋滤作用进入地下水，而S3风险源的原材料中包括CHCl4，在进入土壤或地下水中极易降解为CHCl3等. 8个风险源中，由于S6风险源废水排放量为0，需按废弃地块进行先验概率计算. 8个风险源的行业类别、存在时间、废水排放量等软数据信息及先验概率计算结果见表1.

2.3 污染事件似然度计算

根据研究区水文地质条件，笔者将地下水流概化为均质各向同性的潜水平面二维稳定流. 根据风险源与监测点的相对位置、水头差等数据，利用式(4)计算污染事件发生的似然度与标准化常量的比值〔p(Si,m)/W〕，结果见表2.

图1 地下水风险源及监测点位置及地下水位等值线示意Fig.1 Sketch of the groundwater risk sources location, the monitoring point location and the groundwater table

表1 研究区域内风险源相关信息及先验概率值

2.4 后验概率计算结果分析

基于2.2节和2.3节计算得到的先验概率和似然度，按照式(7)计算得到不同指标异常时不同风险源的后验概率，计算结果如表3所示.

由表3可知，O1监测点Cr6+含量异常由S6风险源造成的后验概率最大，即O1监测点Cr6+含量异常由S6风险源造成的概率为76.2%，由此表明，监测点Cr6+含量异常最有可能由某无机盐制造业污染源造成；O1监测点CHCl3含量异常由S1和S3风险源造药制造业污染源造成.

表2 似然度计算结果

成的后验概率分别为32.7%和23.6%，由此表明，监测点CHCl3含量异常最有可能由一个或两个化学农

表3 后验概率计算结果

进一步调研发现，S6风险源于2014年5月停产，铬渣贮存库房顶棚、围墙破损严重，库房内大量铬渣堆存尚未处置，2016年5月已对其实施挂牌督办. 但此次研究结果显示，Cr6+污染羽已超过1 500 m，存在巨大安全隐患，亟需治理，进一步佐证了监测点Cr6+污染最有可能由S6风险源造成. 但S1风险源和S3风险源从生产工艺过程中涉及的原辅料和产物来说，S3风险源主要以吡啶、四氯化碳、氰化钠、百草谷二氯盐原药等生产有机农药，S1风险源主要以阿维菌素、吡虫啉、粉剂、乳剂等生产农药制品，但是粉剂和乳剂的成分不清，需要进一步补充分析S1风险源是否释放或产生CHCl3，因此，S1和S3风险源是否造成CHCl3含量异常仍需多个监测点或多期监测数据进行反演验证.

2.5 讨论

该研究建立的地下水污染来源识别方法通过将风险源与监测点之间污染物迁移的相关关系进行了简化和求解，在此过程中筛选了关键参数进行分析，主要目的是在有限数据中挖掘最可能解，在地下水污染发现初期锁定最可能的污染来源，以期减少常规地下水污染调查的工作量，也可以使后期的地下水污染调查更具有针对性，为地下水风险管控提供了有效方法和手段.

但是，在研究过程中还发现一些可能影响判定结果的情景，需要下一步重点研究：①当存在行业相似、后验概率结果无显著差异的两个或多个风险源时，需要视现场条件开展少量补充监测，并利用监测结果进行后验概率反演，进一步利用不确定分析理论甄别污染来源，在监测数据充足时还可实现对污染释放过程的反演[27-28]. ②当存在水位波动较大等情况时，利用稳定流进行污染来源判定可能造成准确性下降，需收集区域内多年水位变化情况，分时段进行概率分析再进行加权概率叠加. ③当区域空间异质性显著或存在优势通道时，则无法将研究区水流简单概化为均质稳定流进行风险源污染概率计算，需要通过大量水文地质资料及水位监测数据进行研究区概念模型更新，并利用溶质运移方程得到浓度计算值与实测值的拟合结果，通过不断迭代优化似然函数[21,26]，理论上可以实现非均质条件下的污染溯源，这也是目前和未来一段时间相关研究的重要方向，同时进一步提升不确定性分析方法在实际案例中的应用水平.

3 结论

a) 利用已知的多个风险源的建成时间、废水排放量、原材料、产品等相关软数据信息修正先验概率，并基于对流弥散方程提出似然度的简化求解方法，提出了基于贝叶斯模型的地下水风险源污染概率估计方法，用于计算得到监测点指标异常来源于每个风险源的后验概率，通过结果判定监测值异常的最可能来源，为我国目前监测水平较低现状下污染来源识别问题的解决提供了有效途径.

b) 对石家庄市某工业集聚区内8个风险源下游一个农灌井中Cr6+含量和CHCl3含量异常事件的分析发现，Cr6+含量异常来源于S6风险源的后验概率为76.2%，即Cr6+含量异常最有可能由某无机盐制造业污染源造成；CHCl3含量异常来源于S1和S3风险源的后验概率分别为32.7%和23.6%，监测点CHCl3含量异常最有可能由一个或两个化学农药制造业污染源造成.

c) 该研究仅针对有限信息下的污染来源进行概率判断，风险源是否造成地下水污染仍需多个监测点数据进行验证分析，并结合溶质运移方程继续开展特殊情景下的方法优化研究，使其更加满足地下水污染溯源分析和责任认定应用需求.