CN104820745A - 地表水环境介质中有机化学品暴露水平预测方法 - Google Patents

Abstract

本发明涉及环境生态风险评价领域，具体来说，涉及一种地表水环境下基于质量平衡和动态消减的地表水环境介质中有机化学品暴露水平预测模型和数据分析系统。随着现代化工业生产废水经过处理后最终排放进入周围水体环境，并随着水体挥发进入大气，通过吸附作用或沉降作用进入水体底部沉积物中。本发明建立了一个模型结构简单、预测规则透明、易于理解、有效的暴露预测模型系统，并且能够根据研究的目的和要求选择合适级别的模型，对有机化学品的环境暴露和安全管理以及环境生态风险评价具有重要意义。

Description

地表水环境介质中有机化学品暴露水平预测方法

技术领域

本发明涉及环境生态风险评价领域，具体来说，涉及一种地表水环境下基于质量平衡和动态消减的多介质暴露预测模型和数据分析方法。

背景技术

近年来，我国化学工业迅猛发展，有机化学品的种类、产量、使用量以及排放量急剧增加，而产生的废水排放进入污水处理厂，并最终进入周围水体环境，进入水环境系统的有机化学品在各种跨介质迁移过程的作用下，能够从水体挥发进入大气，可以通过吸附作用富集在悬浮颗粒物中，并通过颗粒物沉降作用进入水体底部沉积物中，因此，在进行有机化学品的生态环境风险评价时，应充分考虑有机化学品在环境各介质中的暴露水平，这就使得有机化学品在多介质环境中的跨介质迁移与转化过程显得非常重要，但是对数量众多的有机化学品在各种环境介质中的时空分布进行研究是不可行的。由于计算机软硬件技术的快速发展，基于建立的多介质暴露水平预测模型，通过计算机语言研发出预测评估软件，可以很大程度上提高效率和准确率，地表水环境介质中有机化学品暴露水平预测方法能够综合、细致和真实的描述污染物在实际水环境中的分配、转移与转化过程，将对已经存在的或新合成的有机化学品的环境影响提供早期评价，这对有机化学品的环境暴露和安全管理以及环境生态风险评价具有重要意义。

目前，已有研究者将多介质暴露预测模型广泛的应用于污染物归趋行为的研究。如文献“Chemosphere,1983,12:981-997”构建以水相为主体，包括沉积物相、生物相和气相的多介质暴露预测模型，描述了污染物在湖泊中的稳态行为；紧接着，文献“Chemosphere,1985,12:1193-1208”将河流分成一系列连接湖段，并假设每一湖段混合均匀并作为独立的水体，用多介质暴露预测模型模拟不挥发有机物在河流中的归趋；为了更加深入的了解污染物在环境中的动态行为，文献“environmental science&technology,1995,29:1200-1209”建立了连续时间变化的多介质预测模型，模拟了持久性有机污染物在英国南部地区空气-土壤两相中的长期变化；随后，文献“environmental science&technology,2000,34:2373-2379”中引入矩阵来描述污染物在不同环境相中环境行为，从而构建了以逸度为基础的水环境下多介质暴露水平预测方法，将质量平衡方程以矩阵的形式联立起来。但是目前多介质暴露预测方法形式单一，无法提取预测规则和预测能力可能会出现过拟合的矛盾问题，可理解性差，不利于模型应用和行为解释，数据分析难等问题。因此，亟需建立一个结构简单、预测规则透明、易于理解、有效的暴露水平预测方法，并且能够根据研究的目的和要求选择合适级别的模型，然后通过Crystal Ball数据系统分析输入参数对输出结果的影响程度，得到影响模型输出结果的关键参数，并通过进一步优化，使预测值与实测值有较好的吻合度；采用Monte Carol仿真功能协助分析，给出模型的不确定性信息，最后提供一个能够较准确的模拟出有机化学品在地表水环境介质中暴露水平预测方法，为环境管理与决策提供有价值的信息。

发明内容

本发明的目的是针对上述不足之处，建立一个结构简单、预测规则透明、易于理解、有效的暴露预测模型系统，为环境管理与决策提供有价值的信息。

地表水环境介质中有机化学品暴露水平预测方法，其特征在于：

步骤一、确定待测地表水环境系统的主环境相及其子相，测量并获得该地表水系统的主要环境参数；确定目标污染物的理化参数及背景排放参数。

步骤二、将所获得的数据输入到该系统中，首先以Mackay多介质逸度模型方法为基础获得各个环境相的逸度容量、降解反应消减参数、相间质量交换参数，并根据步骤一中污染物的背景排放参数，建立以环境相化合物的质量变化＝进入环境相的化合物量-离开环境相的化合物量-降解反应消减的化合物量的方程，该系统将根据已编入的矩阵、逆矩阵和矩阵乘积运算MMULT函数的语法，得到目标污染物在地表水多介质环境稳态系统下各个环境相的逸度和浓度；同时该系统可以将降解反应消减参数、相间质量交换参数和步骤一中污染物的排放参数整合到不同的标志单元格中，并根据研究目的和实际情况确定研究区域在时间上的大小，在Microsoft Windows环境下实现Crystal Ball软件与Matlab R2009b进行连接，实现数据的交换和同步更新，采用龙格-库塔函数对微分方程组进行计算，得到以时间为基础的暴露水平预测模型的逸度，再根据逸度计算得到该地表水环境体系中各个环境相的逸度、浓度和质量分布，可以较为准确地描述残留污染物在地表水环境体系中长期的浓度变化规律，并预测了受污染河流的恢复时间；可计算区域范围内各个境相中目标污染物的动态质量交换和动态富集情况，可为环境管理与决策提供有价值的信息。

步骤三、将地表水环境多介质暴露水平预测方法嵌入到Crystal Ball软件中，并可以根据生成的变量灵敏度结果，得到影响模型输出结果的关键参数，并通过进一步优化，使预测值与实测值有较好的吻合度；采用Monte Carol仿真功能协助分析，给出模型的不确定性信息。

具体地，步骤一所述待测地表水环境系统的主环境相及其子相可以划分为空气相、水相和沉积物相，其中空气相包括气相和气溶胶相，水相包括纯水相和悬浮颗粒物，沉积物相包括沉积物固相和水相；

具体地，步骤二所述地表水环境介质各主环境相及其子相的逸度容量分别为：

气相：Z_A＝1/(RT)

气溶胶相：Z_Q＝Z_A*6*10⁶/P_S

空气相：Z_TA＝V_Q*Z_Q+V_A*Z_A(V_Q+V_A＝1)

纯水相：Z_W＝S/Ps

悬浮颗粒物：Z_P＝Z_W*ф_P*ρ_P*K_OC

水相：Z_TW＝V_W*Z_W+V_P*Z_P(Vw+Vp＝1)

沉积物固相：Z_S＝Z_W*ф_sed*ρ_sed*K_OC

沉积物相：Z_TS＝Vw*Z_W+V_S*Z_S(Vw+VS＝1)

所述Z为逸度容量(mol·m^-3·pa)，R为气体常数(pa·m³·mol^-1·k^-1)，T为环境温度(k)，Ps为过冷蒸汽压(pa)；V_Q为空气中气溶胶微粒的含量(-)，V_A为空气中气体的含量(-)，S为水中的溶解度(mol/m³)，P_s为饱和蒸汽压(pa)，ф_P为悬浮颗粒物的有机碳含量(-)，ρ_P为悬浮颗粒物的密度(kg/m³)，Koc为有机碳分配系数(-)，Vw为水中纯水的含量(-)，Vp为水中悬浮颗粒物的含量(-)，ф_sed为沉积物的有机碳含量(-)，ρ_sed为沉积物的密度(kg/m³)，Vs为沉积物中固体相的含量(-)；

具体地，步骤二所述各个环境相的降解反应消减参数分别为：

空气相：D_RA＝k_RA*V_A*Z_A

$k_{\mathrm{RA}}=\frac{\ln 2}{t_{1/2A}}$

水相：D_RW＝k_RW*V_W*Z_W

$k_{\mathrm{RW}}=\frac{\ln 2}{t_{1/2W}}$

沉积物相：D_RS＝k_RS*V_S*Z_S

$k_{\mathrm{RS}}=\frac{\ln 2}{t_{1/2S}}$

其中，k_RA、k_RW和k_RS分别为目标污染物在空气相、水相和沉积物相中的降解速率(s^-1)，t_1/2A、t_1/2W和t_1/2S为化合物在空气相、水相和沉积物相中的半衰期(s)；

具体地，步骤二所述各个环境相间的质量交换参数分别为：

水向空气的挥发：D_W-A＝1/(1/(K_VA*Aw*Z_A)+1/(K_VW*Aw*Z_W))

空气与水的交换：D_A-W＝D_W-A+D_RW+D_QD+D_QW

D_RW＝A_A-W*U_R*Z_W

D_QW＝A_A-W*U_R*Q*V_Q*Z_Q

D_Qv＝A_A-W*U_Q*V_Q*Z_Q

空气的平流流入：D_A＝G_A*Z_A

水的平流流入：D_W＝G_w*Z_w

沉积物与水间的交换：D_S-W＝D_Y+D_RS

D_Y＝1/﹛1/(k_SW*A_S-W*Z_W)+Y₄/B_MW*A_S-W*Z_W)﹞﹜

D_RS＝U_RS*A_S-W*Z_S

水与沉积物间的交换：D_W-S＝D_Y+D_DS

D_DS＝U_DP*A_W-S*Z_P

所述K_VA水-气界面气侧传质系数(m/h)，K_VW为水-气界面中水侧传质系数(m/h)，A_A-W和A_S-W分别被为空气-水界面和沉积物-水界面的面积(m²)，U_R为年平均降雨量(m/h)、Q为清除速率、G_A为空气的平流流入速率(m³/h)，G_w为水径流量(m³/h)，k_sw为底泥-水界面水侧传质系数，Y₄为底泥中扩散路径长度(m)，B_MW为水中分子扩散系数(m²/h)，U_RS为底泥再悬浮速率(m/h)，U_DP为底泥沉降速率(m/h)；

具体地，步骤二所述各个环境相中目标污染物的质量变化微分方程组分别为：

空气相：

V_AZ_Adf_A/dt＝E_A+G_A*C_A+f_W*D_W-A-f_A*(D_A-W+D_RA+D_A)，

水相：

V_WZ_Wdf_W/dt＝E_W+G_W*C_W+f_A*D_A-W+f_S*D_S-W-f_W*(D_W-A+D_W-S+D_RW+D_W)，

沉积物相：

V_SZ_Sdf_S/dt＝f_W*D_W-S-f_S*(D_S-W+D_RS)，

其中，V_A、V_W和V_S分别为空气相、水相和沉积物相的体积(m³)，E_A和E_W表示污染物向空气和水中的排放速率(mol/h)，C_A和C_W表示污染物平流流入空气和水中的背景浓度(mol/m³)，f_A、f_W和f_S为未知量，分别表示空气相、水相和沉积物相的逸度(pa)；

具体地，步骤二中所述各个环境相中目标污染物的质量变化微分方程组也可以或者表示为：

空气相：df_A/dt＝λ₁*f_W-λ₂*f_A+b₁

水相：df_W/dt＝-λ₃*f_W+λ₄*f_A+λ₅*f_S+b₂

沉积物相：df_S/dt＝λ₆*f_W-λ₇*f_S

其中：λ为未知量f的系数，b为常量；

具体地，步骤二中所述已编入的矩阵、逆矩阵和矩阵乘积运算MMULT函数的语法分别为：

矩阵 $N=\left[\begin{array}{ccc}{\mathrm{\λ}}_1& {\mathrm{\λ}}_2& 0\\ {\mathrm{\λ}}_3& {\mathrm{\λ}}_4& {\mathrm{\λ}}_5\\ {\mathrm{\λ}}_6& 0& {\mathrm{\λ}}_7\end{array}\right]$

$M=\left[\begin{array}{c}{b}_1\\ b_2\\ 0\end{array}\right]$

求矩阵N的逆矩阵的函数表达式：N^-1＝MINVERSE(array)

其中，array表示矩阵N；

矩阵乘积运算MMULT函数表达式： $\left[\begin{array}{c}{f}_A\\ f_W\\ f_S\end{array}\right]=\mathrm{MMULT}(\mathrm{array}1,\mathrm{array}2),$ 其中：array1表示N的逆矩阵，array2表示矩阵M；

具体地，步骤二中将降解反应消减参数、相间质量交换参数和步骤一中污染物的背景排放参数整合到不同的标志单元格中，并在Microsoft Windows环境下实现Crystal Ball软件与Matlab R2009b的连接，整个过程分别为：

标志单元格a＝[b₁ λ₁ λ₂ 0]

标志单元格b＝[b₂ λ₃ λ₄ λ₅]

标志单元格c＝[0 λ₆ 0 λ₇]

标志单元格z＝[Z_A Z_W Z_S]

其中：标志单元格a、b和c分别为矩阵N和M的变形，标志单元格a中的参数分别表示空气相、水相和沉积物相的逸度容量；

实现Crystal Ball软件和Matlab R2009b两种工作环境中数据的交换和同步更新的方法为：

将Crystal Ball工作表中的标志单元格a写入Matlab矩阵a>＝＝MLPutMatrix("a",a)；

将Crystal Ball工作表中的标志单元格b写入Matlab矩阵b>＝＝MLPutMatrix("b",b)；

将Crystal Ball工作表中的标志单元格c写入Matlab矩阵c>＝＝MLPutMatrix("c",c)；

将Crystal Ball工作表中的标志单元格z写入Matlab矩阵c>＝＝MLPutMatrix("z",z)；

将写成字符串形式的命令a传到Matlab中执行>＝＝MLEvalString("save a.txt-ascii a")；

将写成字符串形式的命令b传到Matlab中执行>＝＝MLEvalString("save b.txt-ascii b")；

将写成字符串形式的命令c传到Matlab中执行>＝＝MLEvalString("save c.txt-ascii c")；

将写成字符串形式的命令z传到Matlab中执行>＝＝MLEvalString("save z.txt-ascii z")；

根据研究目的和实际情况确定暴露水平预测方法计算的初始时间、步长和跨度，并以步骤一所获得目标污染物在研究区域的背景排放参数，采用龙格-库塔函数对微分方程组进行计算：

主程序>＝＝MLEvalString("[t,f]＝ode23s(weifen,[t₀,t_s,t_e:],[C_A0C_W0C_S0])")

子程序

在Matlab中执行命令>＝＝MLEvalString("for ii＝1:1:3C(:,ii)＝f(:,ii)*z(ii)")；

将Matlab中的t矩阵放入Crystal Ball工作表中的标志单元格T中

>＝＝MLGetMatrix("t","T")；

将Matlab中的f矩阵放入Crystal Ball工作表中的标志单元格F中

>＝＝MLGetMatrix("f","F")

将Matlab中的c矩阵放入Crystal Ball工作表中的标志单元格C中

>＝＝MLGetMatrix("c","C").

其中，t₀、t_s和t_e分别为模型计算的初始时间，步长和结束时间，C_A0、C_W0和C_S0分别为目标污染物在研究区域中空气相、水相和沉积物相中的背景浓度；

具体地，步骤三中采用软件自带的Monte Carol仿真功能对计算结果进行灵敏度分析和不确定性分析的方法分别为：

1)设定数据表

通过建立数据表将步骤一所获得参数输入到数据中待评估；

2)定义假设的前提

确定单元格中具有波动性的输入参数(除常数外)，并选择服从正态分布；

3)预测结果的确定

确定单元格中目标污染物在各个环境相中的浓度数据是需要预测的指标；

4)选择试验的次数

选择试验的次数，并选择优先运行下的“敏感度分析”；

5)运行模拟

选择开始计算，如果要改变参数重新进行模拟，需要首先重置模拟(点击运行菜单工具栏或者运行菜单下的“重置模拟”按钮)；

6)查看结果

在模拟最后或者运行的过程中，预测窗口会自动显示出目标污染物在地表水环境介质中有机化学品暴露水平预测方法下的不确定性分析和输入参数的灵敏度分析，结果可以复制到工作表中。

本发明是一个结构简单、预测规则透明、易于理解、有效的暴露水平预测方法，并且能够根据研究的目的和要求选择合适级别的模型，对有机化学品的环境暴露和安全管理以及环境生态风险评价具有重要意义。

本发明提供的方法具有如下特点：

本发明适用于描述有机化学品在地表水多介质环境中的分配规律和归趋情况，在环境相的设定上简单且易于理解，并具有较强的实用性。

本发明是建立在表格基础上的暴露水平预测模型，对模型的计算过程和预测规则都是可视化窗口，便于模型的理解和应用。

本发明能够根据研究的目的和要求选择合适级别的模型，不仅可以预测稳态情况下有机化学品在地表水多介质环境中的迁移和转化规律，而且可以得到以时间为基础的暴露水平预测模型的结果，可以较为准确地描述残留污染物在地表水环境体系中长期的浓度变化规律，并可以预测受污染河流的恢复时间，计算区域范围内各环境相中目标污染物的动态质量交换和动态富集情况。

本发明首次将主要应用于商品贸易、金融分析、投资组合分析和工艺过程模拟等领域的Crystal Ball软件与地表水多介质暴露预测模型相结合，直接得到变量灵敏度分析结果，了解影响模型输出结果的关键参数，并通过进一步优化，使预测值与实测值有较好的吻合度，并且采用Monte Carol仿真功能协助分析，给出模型的不确定性信息。

附图说明：

图1为地表水环境系统的概念图；图2为具体实施方式中PCB52在长江南京段30d内的浓度变化图；图3为具体实施方式中空气相中PCB52的灵敏度分析；图4为具体实施方式中水相中PCB52的灵敏度分析；图5为具体实施方式中沉积物相中PCB52的灵敏度分析；图6为具体实施方式中空气相中PCB52的不确定性分析；图7为具体实施方式中水相中PCB52的不确定性分析；图8为具体实施方式中沉积物相中PCB52的不确定性分析；图9为PCB52预测浓度与实测浓度的比较。

具体实施方式

实施例1

本实施例介绍地表水环境介质中有机化学品暴露水平预测的具体方法。

地表水环境介质中有机化学品暴露水平预测方法具体包括以下几步：

步骤一、确定待测地表水环境系统的主环境相及其子相，测量并获得该地表水系统的主要环境参数；确定目标污染物的理化参数及背景排放参数。

步骤二、将所获得的数据输入到该系统中，首先以Mackay多介质逸度模型方法为基础获得各个环境相的逸度容量、降解反应消减参数、相间质量交换参数，并根据步骤一中污染物的背景排放参数，建立以环境相化合物的质量变化＝进入环境相的化合物量-离开环境相的化合物量-降解反应消减的化合物量的方程，该系统将根据已编入的矩阵、逆矩阵和矩阵乘积运算MMULT函数的语法，得到目标污染物在地表水多介质环境稳态系统下各个环境相的逸度和浓度；同时该系统可以将降解反应消减参数、相间质量交换参数和步骤一中污染物的排放参数整合到不同的标志单元格中，并根据研究目的和实际情况确定研究区域在时间上的大小，在Microsoft Windows环境下实现Crystal Ball软件与Matlab R2009b进行连接，实现数据的交换和同步更新，采用龙格-库塔函数对微分方程组进行计算，得到以时间为基础的暴露水平预测模型的逸度，再根据逸度计算得到该地表水环境体系中各个环境相的逸度、浓度和质量分布，可以较为准确地描述残留污染物在地表水环境体系中长期的浓度变化规律，并预测了受污染河流的恢复时间；可计算区域范围内各个境相中目标污染物的动态质量交换和动态富集情况，可为环境管理与决策提供有价值的信息。

步骤三、将地表水环境多介质暴露水平预测方法嵌入到Crystal Ball软件中，并可以根据生成的变量灵敏度结果，得到影响模型输出结果的关键参数，并通过进一步优化，使预测值与实测值有较好的吻合度；采用Monte Carol仿真功能协助分析，给出模型的不确定性信息。

步骤一所述待测地表水环境系统的主环境相及其子相可以划分为空气相、水相和沉积物相，其中空气相包括气相和气溶胶相，水相包括纯水相和悬浮颗粒物，沉积物相包括沉积物固相和水相；

步骤二所述地表水环境介质各主环境相及其子相的逸度容量分别为：

气相：Z_A＝1/(RT)

气溶胶相：Z_Q＝Z_A*6*10⁶/P_S

空气相：Z_TA＝V_Q*Z_Q+V_A*Z_A(V_Q+V_A＝1)

纯水相：Z_W＝S/Ps

悬浮颗粒物：Z_P＝Z_W*ф_P*ρ_P*K_OC

水相：Z_TW＝V_W*Z_W+V_P*Z_P(Vw+Vp＝1)

沉积物固相：Z_S＝Z_W*ф_sed*ρ_sed*K_OC

沉积物相：Z_TS＝Vw*Z_W+V_S*Z_S(Vw+VS＝1)

所述Z为逸度容量(mol·m^-3·pa)，R为气体常数(pa·m³·mol^-1·k^-1)，T为环境温度(k)，Ps为过冷蒸汽压(pa)；V_Q为空气中气溶胶微粒的含量(-)，V_A为空气中气体的含量(-)，S为水中的溶解度(mol/m³)，P_s为饱和蒸汽压(pa)，ф_P为悬浮颗粒物的有机碳含量(-)，ρ_P为悬浮颗粒物的密度(kg/m³)，Koc为有机碳分配系数(-)，Vw为水中纯水的含量(-)，Vp为水中悬浮颗粒物的含量(-)，ф_sed为沉积物的有机碳含量(-)，ρ_sed为沉积物的密度(kg/m³)，Vs为沉积物中固体相的含量(-)；

步骤二所述各个环境相的降解反应消减参数分别为：

空气相：D_RA＝k_RA*V_A*Z_A

$k_{\mathrm{RA}}=\frac{\ln 2}{t_{1/2A}}$

水相：D_RW＝k_RW*V_W*Z_W

$k_{\mathrm{RW}}=\frac{\ln 2}{t_{1/2W}}$

沉积物相：D_RS＝k_RS*V_S*Z_S

$k_{\mathrm{RS}}=\frac{\ln 2}{t_{1/2S}}$

其中，k_RA、k_RW和k_RS分别为目标污染物在空气相、水相和沉积物相中的降解速率(s^-1)，t_1/2A、t_1/2W和t_1/2S为化合物在空气相、水相和沉积物相中的半衰期(s)；

步骤二所述各个环境相间的质量交换参数分别为：

水向空气的挥发：D_W-A＝1/(1/(K_VA*Aw*Z_A)+1/(K_VW*Aw*Z_W))

空气与水的交换：D_A-W＝D_W-A+D_RW+D_QD+D_QW

D_RW＝A_A-W*U_R*Z_W

D_QW＝A_A-W*U_R*Q*V_Q*Z_Q

D_Qv＝A_A-W*U_Q*V_Q*Z_Q

空气的平流流入：D_A＝G_A*Z_A

水的平流流入：D_W＝G_w*Z_w

沉积物与水间的交换：D_S-W＝D_Y+D_RS

D_Y＝1/﹛1/(k_SW*A_S-W*Z_W)+Y₄/B_MW*A_S-W*Z_W)﹞﹜

D_RS＝U_RS*A_S-W*Z_S

水与沉积物间的交换：D_W-S＝D_Y+D_DS

D_DS＝U_DP*A_W-S*Z_P

所述K_VA水-气界面气侧传质系数(m/h)，K_VW为水-气界面中水侧传质系数(m/h)，A_A-W和A_S-W分别被为空气-水界面和沉积物-水界面的面积(m²)，U_R为年平均降雨量(m/h)、Q为清除速率、G_A为空气的平流流入速率(m³/h)，G_w为水径流量(m³/h)，k_sw为底泥-水界面水侧传质系数，Y₄为底泥中扩散路径长度(m)，B_MW为水中分子扩散系数(m²/h)，U_RS为底泥再悬浮速率(m/h)，U_DP为底泥沉降速率(m/h)；

步骤二所述各个环境相中目标污染物的质量变化微分方程组分别为：

空气相：

V_AZ_Adf_A/dt＝E_A+G_A*C_A+f_W*D_W-A-f_A*(D_A-W+D_RA+D_A)，

水相：

V_WZ_Wdf_W/dt＝E_W+G_W*C_W+f_A*D_A-W+f_S*D_S-W-f_W*(D_W-A+D_W-S+D_RW+D_W)，

沉积物相：

V_SZ_Sdf_S/dt＝f_W*D_W-S-f_S*(D_S-W+D_RS)，

其中，V_A、V_W和V_S分别为空气相、水相和沉积物相的体积(m³)，E_A和E_W表示污染物向空气和水中的排放速率(mol/h)，C_A和C_W表示污染物平流流入空气和水中的背景浓度(mol/m³)，f_A、f_W和f_S为未知量，分别表示空气相、水相和沉积物相的逸度(pa)；

步骤二中所述各个环境相中目标污染物的质量变化微分方程组也可以或者表示为：

空气相：df_A/dt＝λ₁*f_W-λ₂*f_A+b₁

水相：df_W/dt＝-λ₃*f_W+λ₄*f_A+λ₅*f_S+b₂

沉积物相：df_S/dt＝λ₆*f_W-λ₇*f_S

其中：λ为未知量f的系数，b为常量；

步骤二中所述已编入的矩阵、逆矩阵和矩阵乘积运算MMULT函数的语法分别为：

矩阵 $N=\left[\begin{array}{ccc}{\mathrm{\λ}}_1& {\mathrm{\λ}}_2& 0\\ {\mathrm{\λ}}_3& {\mathrm{\λ}}_4& {\mathrm{\λ}}_5\\ {\mathrm{\λ}}_6& 0& {\mathrm{\λ}}_7\end{array}\right]$

$M=\left[\begin{array}{c}{b}_1\\ b_2\\ 0\end{array}\right]$

求矩阵N的逆矩阵的函数表达式：N^-1＝MINVERSE(array)

其中，array表示矩阵N；

矩阵乘积运算MMULT函数表达式： $\left[\begin{array}{c}{f}_A\\ f_W\\ f_S\end{array}\right]=\mathrm{MMULT}(\mathrm{array}1,\mathrm{array}2),$ 其中：array1表示N的逆矩阵，array2表示矩阵M；

步骤二中将降解反应消减参数、相间质量交换参数和步骤一中污染物的背景排放参数整合到不同的标志单元格中，并在Microsoft Windows环境下实现CrystalBall软件与Matlab R2009b的连接，整个过程分别为：

标志单元格a＝[b₁ λ₁ λ₂ 0]

标志单元格b＝[b₂ λ₃ λ₄ λ₅]

标志单元格c＝[0 λ₆ 0 λ₇]

标志单元格z＝[Z_A Z_W Z_S]

其中：标志单元格a、b和c分别为矩阵N和M的变形，标志单元格a中的参数分别表示空气相、水相和沉积物相的逸度容量；

实现Crystal Ball软件和Matlab R2009b两种工作环境中数据的交换和同步更新的方法为：

将Crystal Ball工作表中的标志单元格a写入Matlab矩阵a>＝＝MLPutMatrix("a",a)；

将Crystal Ball工作表中的标志单元格b写入Matlab矩阵b>＝＝MLPutMatrix("b",b)；

将Crystal Ball工作表中的标志单元格c写入Matlab矩阵c>＝＝MLPutMatrix("c",c)；

将Crystal Ball工作表中的标志单元格z写入Matlab矩阵c>＝＝MLPutMatrix("z",z)；

将写成字符串形式的命令a传到Matlab中执行>＝＝MLEvalString("save a.txt-ascii a")；

将写成字符串形式的命令b传到Matlab中执行>＝＝MLEvalString("save b.txt-ascii b")；

将写成字符串形式的命令c传到Matlab中执行>＝＝MLEvalString("save c.txt-ascii c")；

将写成字符串形式的命令z传到Matlab中执行>＝＝MLEvalString("save z.txt-ascii z")；

根据研究目的和实际情况确定暴露水平预测方法计算的初始时间、步长和跨度，并以步骤一所获得目标污染物在研究区域的背景排放参数，采用龙格-库塔函数对微分方程组进行计算：

主程序>＝＝MLEvalString("[t,f]＝ode23s(weifen,[t₀,t_s,t_e:],[C_A0C_W0C_S0])")

子程序

在Matlab中执行命令>＝＝MLEvalString("for ii＝1:1:3C(:,ii)＝f(:,ii)*z(ii)")；

将Matlab中的t矩阵放入Crystal Ball工作表中的标志单元格T中

>＝＝MLGetMatrix("t","T")；

将Matlab中的f矩阵放入Crystal Ball工作表中的标志单元格F中

>＝＝MLGetMatrix("f","F")

将Matlab中的c矩阵放入Crystal Ball工作表中的标志单元格C中

>＝＝MLGetMatrix("c","C").

其中，t₀、t_s和t_e分别为模型计算的初始时间，步长和结束时间，C_A0、C_W0和C_S0分别为目标污染物在研究区域中空气相、水相和沉积物相中的背景浓度；

具体地，步骤三中采用软件自带的Monte Carol仿真功能对计算结果进行灵敏度分析和不确定性分析的方法分别为：

1)设定数据表

通过建立数据表将步骤一所获得参数输入到数据中待评估；

2)定义假设的前提

确定单元格中具有波动性的输入参数(除常数外)，并选择服从正态分布；

3)预测结果的确定

确定单元格中目标污染物在各个环境相中的浓度数据是需要预测的指标；

4)选择试验的次数

选择试验的次数，并选择优先运行下的“敏感度分析”；

5)运行模拟

选择开始计算，如果要改变参数重新进行模拟，需要首先重置模拟(点击运行菜单工具栏或者运行菜单下的“重置模拟”按钮)；

6)查看结果

在模拟最后或者运行的过程中，预测窗口会自动显示出目标污染物在地表水环境介质中有机化学品暴露水平预测方法下的不确定性分析和输入参数的灵敏度分析，结果可以复制到工作表中。

实施例2

本实施例为利用本发明进行地表水环境介质中有机化学品暴露水平预测的具体实施案例。

待测水域为长江南京段，确定地表水环境介质的主环境相划分为空气相、水相和沉积物相，其中空气相包括气相和气溶胶相，水相包括纯水相和悬浮颗粒物，沉积物相包括沉积物固相和水相；测量并获得该地表水系统的环境参数，针对环境中检出率和检出浓度较高的持久性有机污染物PCB52进行模拟；确定目标污染物的理化参数及背景排放参数；

计算得到的各主环境相及其子相的逸度容量、降解反应消减参数和质量交换参数，并将降解反应消减参数、相间质量交换参数和PCB52的背景排放参数整合到4个标志单元格中，并实现Crystal Ball软件与Matlab R2009b的连接；

为了模拟PCB52在长江南京段地表水环境下的主要迁移过程，分析其在长江南京段的迁移规律和浓度分布，确定该暴露水平预测方法计算的2014年8月1日，步长为1d，时间跨度为30d，采用龙格-库塔函数对微分方程组进行计算；

采用软件自带的Monte Carol仿真功能对计算结果进行灵敏度分析和不确定性分析。

本发明的一种地表水环境介质中有机化学品暴露水平预测方法已经通过具体的实例进行了描述，本领域技术人员可借鉴本发明内容，适当改变目标污染物、环境相、预测规则等环节来实现相应的其它目的，其相关改变都没有脱离本发明的内容，所有类似的替换和改动对于本领域技术人员来说是显而易见的，都被视为包括在本发明的范围之内。

Claims (1)

1.地表水环境介质中有机化学品暴露水平预测方法，其特征在于：

步骤一、确定待测地表水环境介质的主环境相及其子相，测量并获得该地表水系统的环境参数；确定目标污染物的理化参数及背景排放参数；

步骤二、将所获得的数据输入到模型中，首先以Mackay多介质逸度模型方法为基础获得各个环境相的逸度容量、降解反应消减参数、相间质量交换参数，并根据步骤一中污染物的背景排放参数，建立以“环境相化合物的质量变化＝进入环境相的化合物量—离开环境相的化合物量—降解反应消减的化合物量”的方程，模型将根据已编入的矩阵、逆矩阵和矩阵乘积运算MMULT函数的语法，得到目标污染物在地表水环境多介质稳态系统下各个环境相的逸度和浓度；同时该系统将降解反应消减参数、相间质量交换参数和步骤一中污染物的排放参数整合到不同的标志单元格中，并根据研究目的和实际情况确定研究区域在时间上的大小，在Microsoft Windows环境下实现Crystal Ball软件与Matlab R2009b进行连接，实现数据的交换和同步更新，采用龙格-库塔函数对微分方程组进行计算，得到以时间为基础的暴露预测模型的逸度，再根据逸度计算得到该地表水环境体系中各个环境相的逸度、浓度和质量分布，并预测受污染河流的恢复时间；计算区域范围内各环境相中目标污染物的动态质量交换和动态富集情况；

步骤三、将地表水环境介质中有机化学品暴露水平预测模型嵌入到CrystalBall软件中，并根据生成的变量灵敏度结果，得到影响模型输出结果的关键参数，并通过进一步优化，使预测值与实测值有较好的吻合度；采用Monte Carol仿真功能协助分析，给出模型的不确定性信息；

步骤一所述待测地表水环境系统的主环境相及其子相划分为空气相、水相和沉积物相，其中空气相包括气相和气溶胶相，水相包括纯水相和悬浮颗粒物，沉积物相包括沉积物固相和水相；

步骤二所述地表水环境介质各主环境相及其子相的逸度容量分别为：

气相：Z_A＝1/(RT)

气溶胶相：Z_Q＝Z_A*6*10⁶/P_S

空气相：Z_TA＝V_Q*Z_Q+V_A*Z_A(V_Q+V_A＝1)

纯水相：Z_W＝S/Ps

悬浮颗粒物：Z_P＝Z_W*ф_P*ρ_P*K_OC

水相：Z_TW＝V_W*Z_W+V_P*Z_P(Vw+Vp＝1)

沉积物固相：Z_S＝Z_W*ф_sed*ρ_sed*K_OC

沉积物相：Z_TS＝Vw*Z_W+V_S*Z_S(Vw+VS＝1)

所述Z为逸度容量(mol·m^-3·pa)，R为气体常数(pa·m³·mol^-1·k^-1)，T为环境温度(k)，Ps为过冷蒸汽压(pa)；V_Q为空气中气溶胶微粒的含量(-)，V_A为空气中气体的含量(-)，S为水中的溶解度(mol/m³)，P_s为饱和蒸汽压(pa)，ф_P为悬浮颗粒物的有机碳含量(-)，ρ_P为悬浮颗粒物的密度(kg/m³)，Koc为有机碳分配系数(-)，Vw为水中纯水的含量(-)，Vp为水中悬浮颗粒物的含量(-)，ф_sed为沉积物的有机碳含量(-)，ρ_sed为沉积物的密度(kg/m³)，Vs为沉积物中固体相的含量(-)；

步骤二所述各个环境相的降解反应消减参数分别为：

空气相：D_RA＝k_RA*V_A*Z_A

$k_{\mathrm{RA}}=\frac{\ln 2}{t_{1/2_A}}$

水相：D_RW＝k_RW*V_W*Z_W

$k_{\mathrm{RW}}=\frac{\ln 2}{t_{1/2_W}}$

沉积物相：D_RS＝k_RS*V_S*Z_S

$k_{\mathrm{RS}}=\frac{\ln 2}{t_{1/2_S}}$

其中，k_RA、k_RW和k_RS分别为目标污染物在空气相、水相和沉积物相中的降解速率(s^-1)，t_1/2A、t_1/2W和t_1/2S为化合物在空气相、水相和沉积物相中的半衰期(s)；

步骤二所述各个环境相间的质量交换参数分别为：

水向空气的挥发：D_W-A＝1/(1/(K_VA*Aw*Z_A)+1/(K_VW*Aw*Z_W))

空气与水的交换：D_A-W＝D_W-A+D_RW+D_QD+D_QW

D_RW＝A_A-W*U_R*Z_W

D_QW＝A_A-W*U_R*Q*V_Q*Z_Q

D_Qv＝A_A-W*U_Q*V_Q*Z_Q

空气的平流流入：D_A＝G_A*Z_A

水的平流流入：D_W＝G_w*Z_w

沉积物与水间的交换：D_S-W＝D_Y+D_RS

D_Y＝1/﹛1/(k_SW*A_S-W*Z_W)+Y₄/B_MW*A_S-W*Z_W)﹞﹜

D_RS＝U_RS*A_S-W*Z_S

水与沉积物间的交换：D_W-S＝D_Y+D_DS

D_DS＝U_DP*A_W-S*Z_P

所述K_VA水-气界面气侧传质系数(m/h)，K_VW为水-气界面中水侧传质系数(m/h)，A_A-W和A_S-W分别被为空气-水界面和沉积物-水界面的面积(m²)，U_R为年平均降雨量(m/h)、Q为清除速率、G_A为空气的平流流入速率(m³/h)，G_w为水径流量(m³/h)，k_sw为底泥-水界面水侧传质系数，Y₄为底泥中扩散路径长度(m)，B_MW为水中分子扩散系数(m²/h)，U_RS为底泥再悬浮速率(m/h)，U_DP为底泥沉降速率(m/h)；

步骤二所述各个环境相中目标污染物的质量变化微分方程组分别为：

空气相：

V_AZ_Adf_A/dt＝E_A+G_A*C_A+f_W*D_W-A-f_A*(D_A-W+D_RA+D_A)，

水相：

V_WZ_Wdf_W/dt＝E_W+G_W*C_W+f_A*D_A-W+f_S*D_S-W-f_W*(D_W-A+D_W-S+D_RW+D_W)，

沉积物相：

V_SZ_Sdf_S/dt＝f_W*D_W-S-f_S*(D_S-W+D_RS)，

其中，V_A、V_W和V_S分别为空气相、水相和沉积物相的体积(m³)，E_A和E_W表示污染物向空气和水中的排放速率(mol/h)，C_A和C_W表示污染物平流流入空气和水中的背景浓度(mol/m³)，f_A、f_W和f_S为未知量，分别表示空气相、水相和沉积物相的逸度(pa)；

步骤二中所述各个环境相中目标污染物的质量变化微分方程组或者表示为：

空气相：df_A/dt＝λ₁*f_W-λ₂*f_A+b₁

水相：df_W/dt＝-λ₃*f_W+λ₄*f_A+λ₅*f_S+b₂

沉积物相：df_S/dt＝λ₆*f_W-λ₇*f_S

其中：λ为未知量f的系数，b为常量；

步骤二中所述已编入的矩阵、逆矩阵和矩阵乘积运算MMULT函数的语法分别为：

矩阵 $N=\left[\begin{array}{ccc}{\mathrm{\λ}}_1& {\mathrm{\λ}}_2& 0\\ {\mathrm{\λ}}_3& {\mathrm{\λ}}_4& {\mathrm{\λ}}_2\\ {\mathrm{\λ}}_6& 0& {\mathrm{\λ}}_7\end{array}\right]$

$M=\left[\begin{array}{c}{b}_1\\ b_2\\ 0\end{array}\right]$

求矩阵N的逆矩阵的函数表达式：N^-1＝MINVERSE(array)

其中，array表示矩阵N；

矩阵乘积运算MMULT函数表达式： $\left[\begin{array}{c}{f}_A\\ f_W\\ f_S\end{array}\right]=\mathrm{MMULT}(\mathrm{array}1,\mathrm{array}2),$ 其中：array1表示N的逆矩阵，array2表示矩阵M；

所述步骤二中将降解反应消减参数、相间质量交换参数和步骤一中污染物的背景排放参数整合到不同的标志单元格中，并在Microsoft Windows环境下实现Crystal Ball软件与Matlab R2009b的连接，整个过程分别为：

标志单元格a＝[b₁ λ₁ λ₂ 0]

标志单元格b＝[b₂ λ₃ λ₄ λ₅]

标志单元格c＝[0 λ₆ 0 λ₇]

标志单元格z＝[Z_A Z_W Z_S]

其中：标志单元格a、b和c分别为矩阵N和M的变形，标志单元格a中的参数分别表示空气相、水相和沉积物相的逸度容量；

实现Crystal Ball软件和Matlab R2009b两种工作环境中数据的交换和同步更新的方法为：

将Crystal Ball工作表中的标志单元格a写入Matlab矩阵a>＝＝MLPutMatrix("a",a)；

将Crystal Ball工作表中的标志单元格b写入Matlab矩阵b>＝＝MLPutMatrix("b",b)；

将Crystal Ball工作表中的标志单元格c写入Matlab矩阵c>＝＝MLPutMatrix("c",c)；

将Crystal Ball工作表中的标志单元格z写入Matlab矩阵c>＝＝MLPutMatrix("z",z)；

将写成字符串形式的命令a传到Matlab中执行>＝＝MLEvalString("save a.txt-ascii a")；

将写成字符串形式的命令b传到Matlab中执行>＝＝MLEvalString("save b.txt-ascii b")；

将写成字符串形式的命令c传到Matlab中执行>＝＝MLEvalString("save c.txt-ascii c")；

将写成字符串形式的命令z传到Matlab中执行>＝＝MLEvalString("save z.txt-ascii z")；

根据研究目的和实际情况确定暴露水平预测方法计算的初始时间、步长和跨度，并以步骤一所获得目标污染物在研究区域的背景排放参数，采用龙格-库塔函数对微分方程组进行计算：

主程序>＝＝MLEvalString("[t,f]＝ode23s(weifen,[t₀,t_s,t_e:],[C_A0C_W0C_S0])")

子程序

function dft＝weifen(t,f)

a＝importdata('a.txt')；a＝a(1:4)；

b＝importdata('b.txt')；b＝b(1:4)；

c＝importdata('c.txt')；c＝c(1:4)；

z＝importdata('z.txt')；z＝z(1:3)；

dft＝[a(1)+f(1)*a(2)+f(2)*a(3)+f(3)*a(4)；

b(1)+f(1)*b(2)+f(2)*b(3)+f(3)*b(4)；

c(1)+f(1)*c(2)+f(2)*c(3)+f(3)*c(4)；]；

在Matlab中执行命令>＝＝MLEvalString("for ii＝1:1:3C(:,ii)＝f(:,ii)*z(ii)")；

将Matlab中的t矩阵放入Crystal Ball工作表中的标志单元格T中

>＝＝MLGetMatrix("t","T")；

将Matlab中的f矩阵放入Crystal Ball工作表中的标志单元格F中

>＝＝MLGetMatrix("f","F")

将Matlab中的c矩阵放入Crystal Ball工作表中的标志单元格C中

>＝＝MLGetMatrix("c","C").

其中，t₀、t_s和t_e分别为模型计算的初始时间，步长和结束时间，C_A0、C_W0和C_S0分别为目标污染物在研究区域中空气相、水相和沉积物相中的背景浓度；

所述步骤三中采用软件自带的Monte Carol仿真功能对计算结果进行不确定性分析的方法分别为：

1)设定数据表

通过建立数据表将步骤一所获得参数输入到数据中待评估；

2)定义假设的前提

确定单元格中具有波动性的输入参数(除常数外)，并选择服从正态分布；

3)预测结果的确定

确定单元格中目标污染物在各个环境相中的浓度数据是需要预测的指标；

4)选择试验的次数

选择试验的次数，并选择优先运行下的“敏感度分析”；

5)运行模拟

选择开始计算，如果要改变参数重新进行模拟，需要首先重置模拟(点击运行菜单工具栏或者运行菜单下的“重置模拟”按钮)；

6)查看结果

在模拟最后或者运行的过程中，预测窗口会自动显示出目标污染物在地表水环境介质中有机化学品暴露水平预测方法下的不确定性分析和输入参数的灵敏度分析，结果可以复制到工作表中。