基于过程模拟及不确定性分析的污染物健康风险评价方法与流程

本发明涉及一种基于过程模拟及不确定性分析的典型有机污染物的健康风险评价方法，属于污染物健康风险评价技术领域。

背景技术

在许多水资源短缺地区，地下水是居民生活、灌溉及工业生产的重要水源。但随着社会经济的快速发展，地下水污染已成为威胁生态系统与人体健康的重要环境问题。同时，各种理化性质较为复杂的典型有机污染物出现，如四氯乙烯、多氯联苯等，这些污染物密度比水大且不溶于水，增加了人们对地下水污染监测及防治工作的难度。为了降低地下水中典型污染物对人类健康的危害，并为污染场地的治理提供决策信息，针对典型有机污染物的健康风险评价已成为地下水环境管理的重要组成部分。

人体健康风险(hhr)评价是指评估目前或者未来人体可能受到污染环境中化学物质的影响从而对健康产生危害的可能性，1983年美国科学院提出了健康风险评价的四步法：危害鉴定、剂量反应评估、暴露评估、风险表征。1989年美国环境保护署提出了《超级基金项目人体健康风险评价导则》，导则中提出的人体一生中患癌症概率增加的可能性(ilcr)被广泛用于衡量风险评价结果，最大可接受风险值为1×10^-6：

r(x)＝1-exp[-add(x)×cpf](i)

其中，r(x)代表研究区任何一个控制平面处的ilcr值；add表示日均污染物暴露量；cpf表示潜在致癌因子代谢率；add的表达式又可以写作：

其中，ir代表日饮水量；bw代表体重；ed代表暴露时长；ef代表暴露频率；at代表人的寿命，代表控制面上的平均浓度。

污染物浓度分布是人体健康风险评价的重要影响因素，常规的健康风险评价中针对污染物浓度数据多通过实地取样及实验室分析获得，这种方法简单快捷，但没有考虑到污染物在地下水中的迁移转化过程。目前一种基于过程模拟的健康风险评价方法逐渐被人们认可，这种方法耦合了地下水污染物迁移模型与健康风险评价模型，研究了污染物的时空分布规律。

地下水污染物迁移的模拟过程会受到各种不确定性因素的干扰，如模型参数的不确定性、模型结构的不确定性及观测数据的误差。这些不确定性往往会导致模拟出来的污染物浓度可靠性不高，使得最终的健康风险评价结果精度较差，极易造成地下水管理及决策的失误，浪费人力、物力资源。

技术实现要素：

针对于上述现有技术的不足，本发明的目的在于提供一种基于过程模拟及不确定性分析的污染物健康风险评价方法，以实现降低健康风险评价结果的不确定性，通过研究污染物迁移模型的不确定性来增加人体健康风险评价模型的可靠性。

为达到上述目的，本发明采用的技术方案如下：

本发明的一种基于过程模拟及不确定性分析的污染物健康风险评价方法，包括步骤如下：

1)耦合地下水污染物迁移数值模型与人体健康风险评价模型；

2)对污染物在地下水中的浓度变化进行模拟；

3)采用马尔科夫链蒙特卡洛方法处理模拟过程中的参数不确定性，将进行不确定性分析后的参数代入地下水污染物迁移数值模型与人体健康风险评价模型计算，得到最终评价结果。

优选地，所述步骤1)具体为：获得污染物在地下水中的迁移过程，将地下水污染物迁移数值模型与人体健康风险评价模型相耦合，建立地下水污染物运移模型。

优选地，所述步骤1)中的地下水污染物迁移数值模型，采用多相流达西定律建立其平衡方程：

其中，t代表时间，代表孔隙度，sβ代表β相的饱和度，ρβ代表β相的密度，代表κ质在β相中的质量分数，vn代表流域，γn代表流域面积，k代表绝对渗透率，krβ代表β相的相对渗透率，μβ代表β相的黏度；pβ代表β相的流体压力，g代表重力加速度，n代表内法向量，q^κ代表单位体积内生成或消耗热量的速率。

优选地，所述步骤2)具体为：利用tough2软件模拟非溶解相有机污染物从污染源至地下水的迁移过程，同时预测一段时间后地下水中污染物的空间分布规律。

优选地，所述步骤3)具体包括：马尔科夫链蒙特卡洛方法基于贝叶斯理论处理地下水数值模拟中的不确定性，其基础是贝叶斯公式：

式中，θ代表地下水参数中的随机变量，p(θ)代表参数的先验分布密度，p(y|θ)代表参数的似然函数，含义是当模型参数为θ时，模型输出与现有观测资料y之间的相似度；p(θ|y)代表参数的后验分布密度；

通过建立服从平稳分布π(x)的马尔科夫链，并在其分布内进行随机抽样，在整个马尔科夫链的演化过程中实现对目标函数概率分布空间的充分搜索，对关键参数进行反求。

优选地，所述步骤3)具体还包括：马尔科夫链蒙特卡洛方法通过抽样算法实验对目标系统的模拟，采用延迟拒绝自适应(dreamzs)抽样算法对模型的先验信息进行处理。

优选地，所述步骤3)具体还包括：延迟拒绝自适应抽样算法分为预热期和正式演化期两部分，其步骤如下：

预热期：

(a)生成模型参数的先验信息，通过先验分布得到一组初始样本[θⁱ，i＝1,……,n]，n表示平行马尔科夫链的数量，初始样本分别作为n条链的起点，θⁱ表示第i链的当前状态，代表第i链当前状态的第j个元素，θ＝θ[θ1,……,θd]，d代表需要进行识别的模型的参数维数；

(b)t＝1，t代表时间，crm＝m/ncr，m＝1,……,ncr，cr为子空间演化的交叉概率，ncr提前设定，pm＝1/ncr，表示crm所对应的概率，lm＝0，是更新pm的参数；

(c)选择似然函数，计算每条链起点的概率密度l[m(θⁱ|i，z)]，i＝1,…,n；

(d)为链i生成备选点xⁱ：

式中，δ为用于生成备选点的异链对数；r1(j)，r2(n)∈[1,……,n]，且r1(j)≠r2(n)≠i,j＝1,……,δ；e，ε是来自d维均匀分布和正态分布的随机数；γ(δ，d′)是跳跃尺度，取决于δ和d′，d′为子空间演化时d的替代值；

(e)基于多项式分布p(p1,……,pcn)，从1,……,ncr中抽样得到m；

(f)令交叉概率cr＝m/ncr，lm＝lm+1；

(g)在1-cr的概率下，用对的每个元素进行替代，u为[0,1]间的随机数；当cr＝1，不对进行替代：

(h)计算备选点xⁱ的概率密度和接受概率α(θⁱ，xⁱ)：

(i)判断是否接受xⁱ；若α≥u，则接受xⁱ为链i的样本，否则拒绝；

(j)计算归一化的平方跳跃距离δm：

式中，rj表示n条链的当前第j维的标准差；

(k)为链1,……,n重复步骤(d)-(j)；

(l)cr值概率分布更新：

(m)t＝t+1；

(n)重复步骤(k)–(m)直到t满足预热期时段长度；

(o)iqr统计，去除外层链；计算每条链后50％样本的后验密度对数的均值w[w1,……,wn]，iqr＝q3-q1，q1、q3分别表示n条链w的1/4与3/4分位数；当w＜q1-2iqr时，称为外层链；预热期时有链被识别为外层链，则需实施另一个预热阶段，再进行iqr检测，直到没有外层链；

正式演化期：

(p)包含上述预热期中的步骤(a)-(i)；

(q)为链1,……,n重复步骤(a)；

(r)t＝t+1；

(s)对步骤(b)-(c)重复l次，l为平行链的演化次数，提前设置；

(t)使用每条链的后50％样本进行收敛测试；

(u)若对于反演参数的每维样本均达到收敛标准，则停止演化，若未达到，回到步骤(d)。

本发明的有益效果：

本发明的方法可以从降低污染物迁移模型不确定性角度提高健康风险评价可靠性，能够更全面的刻画“源-路径-受体”这一风险评价过程，对于最终造成的风险值进行定量评价。并可以评估未来一段时间污染物对人体健康风险的影响，为污染场地管理与风险决策提供依据。

附图说明

图1为二维砂箱示意图；

图2为模拟起始时刻砂箱内污染物分布情况示意图；

图3a为渗透率参数k1后验分布密度示意图；

图3b为渗透率参数k2后验分布密度示意图；

图4为砂箱实验风险值的概率密度函数示意图。

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

实施例：借助二维砂箱实验反求介质的渗透率k，并计算砂箱右边界处的ilcr值。二维砂箱长(x)60cm，宽(y)45cm，厚(z)1.6cm。在x＝0与x＝60cm处设为定水头边界，砂箱内水流流速1m/d，污染物注入点位于x＝30cm，y＝40cm，z＝0.8cm，如图1所示。模拟开始时刻恰好污染物停止注入，分布如图2，风险控制面取最右侧靠近右边界的地方。在图中黑色虚线方框内均匀的布设了16×38(x×y)个观测点。考虑到填砂过程中的非均质性，尽管在实验设计时图中黑色实线方框并未布设透镜体，但污染物在这个位置并未沿下方透镜体右侧绕流而出，猜测黑色实线方框区域存在局部非均质性。因此针对二维砂箱反演的渗透率有两个：一个是砂箱背景质的渗透率k1，一个是黑色方框处的渗透率k2。针对二维砂箱实验的健康风险评价过程具体步骤如下：

(1)根据二维砂箱水力与介质条件，建立二维非均质承压含水层地下水动力学模型与地下水污染物迁移模型；

(2)选取渗透率k为待估参数，生成其先验信息，假设k的先验概率服从均匀分布，k1∈[1.0×10^-10，20×10^-10]，k2∈[1.0×10^-11，10×10^-11]；

(3)设置3条马尔科夫链，即n＝3。交叉概率参数ncr＝3，差别进化的异链对数δ＝3，每5次迭代，将跳跃尺度γ设为1.0；

(4)针对未知参数渗透率k，利用砂箱中的608个观测点污染物饱和度作为识别数据，采用dreamzs-mcmc方法进行模型参数反演；其中预热期所有链总迭代次数为5000，正式模拟期所有链的总迭代次数为20000；

(5)收集模拟样本，得到渗透率k的后验分布密度，如图3a及图3b，k1的分布范围集中在2.4×10^-10～2.7×10^-10，k2的分布范围与先验一致，但集中在1.0×10^-11～4.0×10^-11；我们将k1、k2的后验分布密度作为污染物迁移模型中的输入值，代入步骤(1)获得控制面处的污染物浓度后验分布密度；

(6)将浓度分布密度代入健康风险评价模型：得到风险值的概率密度函数，如图4所示，风险值的分布范围是1.356×10^-5～1.363×10^-5。

结合经验数据，可以得到砂箱中的k1、k2的经验值k1＝1.35×10^-10，k2＝3.66×10^-11，计算经验值条件下的风险值，结果为1.357×10^-6，发现其坐落在概率密度函数的左侧，而经过不确定性分析后的风险值大概率可能位于1.362×10^-6附近。显然，由参数经验值得到的结果会造成对人体健康风险的低估，从而导致决策上的失误。因此，在风险评价过程中得到确切的输入参数范围对最终的评价结果有重要意义。

本发明方法可以有效降低健康风险评价的不确定性，为水环境保护及污染场地管理决策提供支持，避免人力物力浪费。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。