[发明专利]基于过程模拟及不确定性分析的污染物健康风险评价方法

权利要求书

1.一种基于过程模拟及不确定性分析的污染物健康风险评价方法，其特征在于，包括步骤如下：

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

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

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

将浓度分布密度代入人体健康风险评价模型，得到风险值的概率密度函数；

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

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

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

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

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

2.根据权利要求1所述的基于过程模拟及不确定性分析的污染物健康风险评价方法，其特征在于，所述步骤1)具体为：获得污染物在地下水中的迁移过程，将地下水污染物迁移数值模型与人体健康风险评价模型相耦合，建立地下水污染物运移模型。

3.根据权利要求1所述的基于过程模拟及不确定性分析的污染物健康风险评价方法，其特征在于，所述步骤2)具体为：利用TOUGH2软件模拟非溶解相有机污染物从污染源至地下水的迁移过程，同时预测一段时间后地下水中污染物的空间分布规律。

4.根据权利要求1所述的基于过程模拟及不确定性分析的污染物健康风险评价方法，其特征在于，所述步骤3)具体还包括：马尔科夫链蒙特卡洛方法通过抽样算法实验对目标系统的模拟，采用延迟拒绝自适应抽样算法对模型的先验信息进行处理。

5.根据权利要求4所述的基于过程模拟及不确定性分析的污染物健康风险评价方法，其特征在于，所述步骤3)具体还包括：延迟拒绝自适应抽样算法分为预热期和正式演化期两部分，其步骤如下：

预热期：

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

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

(c)选择似然函数，计算每条链起点的概率密度L[M(θⁱ|I，Z)]，i＝1,……,N；

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

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

(e)基于多项式分布p(p₁,……,p_cn)，从1,……,n_CR中抽样得到m；

(f)令交叉概率CR＝m/n_CR，L_m＝L_m+1；

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

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

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

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

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

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

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

(m)t＝t+1；

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

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

正式演化期：

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

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

(r)t＝t+1；

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

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

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