有机污染物环境逸度模型构建及软件化

李俊明，王斌，衡利苹，余刚

(1.北京航空航天大学化学学院，北京 100191 ；2.清华大学环境学院，环境前沿技术北京实验室，新兴有机污染物控制北京市重点实验室，北京 100084;3.清华苏州环境创新研究院，江苏 苏州 215163)

0 前言

自然环境具有复杂性和动态性，仅利用环境监测手段，耗时耗力，且难以系统反映污染物在环境中迁移与转化过程[1-2]。通过建立环境归宿模型，可以推算有机污染物的环境污染和归趋[3-5]。在众多的有机污染物环境归宿模型中，Donald Mackay教授创建的多介质环境逸度模型(Multimedia environmental fugacity model)有较好的可靠性[6-9]。不论是对研究人员还是对决策者来说，一个能模拟有机污染物环境污染过程和归趋的可靠工具显得尤为重要。

随着计算机技术的迅猛发展，目前国外有较为成熟的逸度模型模拟软件且已被作为一项研究工具应用，例如Simple Box[10]、Chem CAN[11]、GREAT-ER[12]、LOTOS-EUROS[13]等，而国内少有这类逸度模型软件工具，且绝大多数已有软件中没有集成参数灵敏度分析和模型不确定分析。以往我们使用逸度模型预测有机物归趋状况或进行结果分析时，常常利用MATLAB等软件编程计算[14-16]，但这样的做法不方便计算程序的移植，并且对于非专业人员来说直接编程求解较为困难。现希望借助计算机软件平台，在Ⅲ级逸度模型的基础上，建立可靠的有机逸度模型集成软件，并集成参数灵敏性和模型的不确定性分析，为有机污染物环境归趋分析、风险评估和预警提供支撑。

1 软件构建理论

1.1 质量平衡模型

基于三级逸度模型，在非平衡、稳态、流动的条件下，将逸度容量作为浓度与逸度的纽带，结合质量平衡模型[17]，可求解得到有机污染物在环境各介质中的浓度。式(1)—(4)分别为空气相、水体相、土壤相、底泥相的质量平衡方程。

(1)

(2)

(3)

(4)

式中：Ei(i=1,2,3,4)——有机污染物向介质i的排放速率，mol/h；GAi(i=1,2)——该有机污染物随介质i平流输入速率，m3/h；cAi(i=1,2)——对应介质i平流输入浓度，mol/m3；fi(i=1,2,3,4)——介质i中有机污染物的逸度；Dji(i,j=1,2,3,4)——有机污染物由j介质传输到介质i的传输通量系数, mol/(h·Pa)；DRi(i=1,2,3,4)、DAi(i=1,2,4)分别为在介质i中有机污染物因为反应和对流的消减系数，mol/(h·Pa)。以上i，j= 1，2，3，4，分别代表空气相、水体相、土壤相、底泥相。

通过求解各相逸度，通过式(5)可进一步求解有机物在各相的浓度。

(5)

式中：c——浓度，mol/m3；Z——逸度容量，mol/(m3·Pa)；f——逸度，Pa。

1.2 灵敏度计算

输入参数的可获得性和准确性是保证模型计算结果准确的重要因素，灵敏度用于评估某个参数对计算结果的影响程度大小，能够帮助用户将最大精力用于高灵敏度参数获取上，进而优化模型计算结果，提高预测值的准确度。灵敏度CS的计算公式如下：

(6)

式中：Rx——将某个输入参数改变为其x倍后得到的计算结果值，ng/L。

本研究中，输入参数的灵敏度数值的绝对值>0.2时，视为该参数对模拟计算结果影响大；反之视为对结果影响小。

1.3 不确定性分析

不确定性分析是模型输出结果变异性评估的重要手段之一，本软件不确定性分析的实现依赖于蒙特卡洛模拟，其具体方法是根据参数的概率分布或对参数设定一个合适的概率分布，生成大量服从该分布的随机数，然后将这些随机数组输入到软件进行计算，从而得到一组计算结果进行概率统计分析。本软件通过对该序列进行排序，最终序列的上下分位数值可得到不确定度范围[18]。求解过程见图1。在不能获得参数真实分布的情况下，为简化处理过程，通常可以设定参数服从正态分布、对数正态分布或三角分布等。

图1 不确定度求解过程

另外，为对比各物质的不确定度大小，引入相对不确定度，计算公式如下：

(7)

2 软件系统设计及实现

2.1 需求分析

2.1.1 用户需求分析

对于决策者来说，通过逸度模型定性分析某种污染物在该地区的分布情况，掌握污染物的归趋，辅助制定当地环境综合管理方案；对于研究人员来说通过逸度模型分析某种有机物环境污染和归趋情况，为其风险评估提供支撑。另外，技术人员还需在得到模拟结果后进行灵敏度和不确定度分析，确定模拟结果的可靠性。

2.1.2 功能需求分析

2.1.2.1 参数输入

作为软件的输入端，包括有机物物化参数如摩尔质量、密度、水溶解度、蒸气压、半衰期等；区域参数如区域面积、环境温度、介质密度、介质体积分数等；传质及排放参数如降雨速率、干沉降速率、排放量等。

作为软件的核心功能，通过相关参数输入，在三级逸度模型基础上，计算有机物空气相、水体相、土壤相、底泥相的浓度。该模块所得计算结果在输出窗口上显示。

2.1.2.3 结果分析

在模拟计算完成后，可对模拟结果进行分析，包括灵敏度分析和不确定度分析部分。灵敏度分析帮助用户分析各输入参数的重要程度；不确定度分析能够表示因为参数变异性和不确定性，以及模型本身的不确定性引起的计算结果的不确定性或变化范围。

2.1.2.4 帮助

主要包括软件的使用方法及注意事项。

2.2 开发平台

软件集成平台见表1。由表1可见，逸度模型软件的构建共有4层，第一层的用户界面使用了Pyside 2函数库进行设计。第二层是逸度模型算法的逻辑层，使用Python作为主要的实现工具，由于在计算中有大量矩阵需处理，因此选择NumPy函数库帮助快速矩阵运算，而Matplotlib可方便实现数据的可视化。第三层是操作系统层，为Microsoft Windows 10。第四层为系统提供计算的硬件层。

表1 软件集成平台

2.3 系统功能实现

2.3.1 技术路线

在双侧向测井响应中，低角度裂缝的电阻率值在致密高电阻率背景下明显降低，曲线形状尖锐，深、浅侧向测井值一般呈负差异。裂缝的张开度与电导率成正比。当张开度增大时，低角度裂缝的深、浅侧向电阻率都下降，幅度差值增大，但它们之间的幅度差远小于高角度裂缝。在给定条件下，对于油气层，当侵入深度大于深侧向探测深度时，为低阻负差异；对于浅层侵入，无论流体性质如何，均为高阻正差异。

利用Python语言和NumPy函数库完成有机污染物归趋模拟的基本功能，结合Matplotlib函数库完成系统的灵敏度分析和不确定度分析功能，最后由Pyside2函数库完成软件的图形界面制作，将各功能模块组装，最终完成逸度模型软件化。技术路线见图2。

图2 技术路线

2.3.2 程序结构

系统的程序结构见图3。

图3 程序结构

由图3可见， 通过输入参数，软件将模拟计算有机物在该地区的归趋情况，用户可根据不确定度分析及灵敏度分析结果调整参数，并进一步优化模拟值，并与观测值进行比较，直到达到验证要求为止。

2.3.3 图形界面

本系统为多文档界面(MDI)，能够同时显示多个文档，每个文档都显示在其自己的窗口。系统总体分为三类窗口，分别为主窗口，输入窗口和输出窗口。其中输入窗口包括有机污染物物化参数子窗口，区域参数子窗口，传质和输入参数子窗口；输出窗口包括模拟结果子窗口，灵敏度分析子窗口和不确定度分析子窗口。各窗口之间的功能及联系见图4。

图4 各窗口之间功能及联系

3 应用实例及验证

近年来，药品和个人护理产品(PPCPs)作为一类新污染物，已成为全球关注的重点问题，我国已经开始开展新污染物治理行动[19]。随着北京市城市化发展，PPCPs已成为地表水中普遍存在的污染物，对生态系统和人体健康造成潜在风险[20]，因此选取北京地区地表水中的几种PPCPs进行实例模拟以验证软件的可靠性。

3.1 参数输入

3.1.1 模拟区域概况

北京市地处华北平原北部，总面积16 370 km2，水域面积212 km2。北京市属于温带季风气候，夏季高温多雨，冬季寒冷干燥，四季分明，年平均温度11.8 ℃，年均降雨量644 mm。

3.1.2 软件主要输入参数

所模拟的有机物包括苯扎贝特、吲哚美辛、双氯芬酸、避蚊胺4种典型的PPCPs，软件的主要输入参数包括4种有机物的理化性质参数与排放速率以及传质参数[21-22]，具体见表2、表3。

表2 4种PPCPs的理化性质参数及其排放速率

表3 主要传质参数 m/h

3.2 模拟结果

将苯扎贝特、吲哚美辛、双氯芬酸、避蚊胺的物化参数、北京区域参数以及相关的传质与输入参数录入到软件中，经由软件模拟得到计算结果。水相中PPCPs的软件模拟浓度和实测浓度[22]见表4。

表4 水相中PPCPs的软件模拟和实测浓度 ng/L

为方便评估模型的准确性，模拟值和实测值的对比见图5。由图5可见，4种物质的模拟浓度和实测浓度较为接近，模拟值与实测值之差均在0.15个对数单位内(取绝对值)，表明模拟值与测量值吻合性好，模型拟合良好，软件预测结果较为可信。

图5 北京地区PPCPs在水相中的软件模拟值与实测值对比分析

3.3 灵敏度分析结果

利用软件的灵敏度分析功能，筛选出对输出结果影响较大的主要输入参数，见图6。由图6可见，除辛醇-水分配系数对避蚊胺的模拟结果影响较小外，各物质在水体相中模拟结果影响较大的输入参数大致相同，分别为区域面积、区域水占比、水体排放速率、水体深度、水中的半衰期、辛醇-水分配系数以及水体停留时间。其中，区域面积、区域水占比、水体排放速率对模拟结果的影响结果最大，在模型软件应用过程中应优先确定这些参数的准确性。

图6 水体相中各有机污染物灵敏度分析

3.4 不确定度分析结果

利用软件的不确定分析功能，可以获得污染物环境介质浓度概率分布，并求解模拟结果的相对不确定度。各有机污染物的不确定度分析见表5。由表5可见，4种物质的相对不确定度接近。其中，苯扎贝特的相对不确定度较高，说明其环境浓度变异性相对较大。

表5 PPCPs的不确定度分析 ng/L

3 结论

基于三级逸度模型理论设计的软件实现了预测有机物在多介质环境中分布的基本功能，并将模型计算、灵敏度分析、不确定度分析等功能集成到一套软件中，设计了友好的用户界面和可视化结果可供专业或非专业人员使用。利用本软件对北京地区地表水中的PPCPs进行了浓度预测，结果显示软件模拟值与实测值较接近，说明了软件具有可靠性。模型在进行不确定度分析时，应尽量采用参数的真实概率分布，以期得到更准确可靠的结果，为污染物风险评估提供支撑。