1 引言(Introduction)

VOCs是形成臭氧(O₃)和细颗粒物(PM_2.5)污染的重要前体物(王伶瑞等, 2020), 会造成雾霾及光化学污染等大气污染问题.同时VOCs具有毒性和刺激性, 严重威胁着人们的健康.据统计2015年全国338个地级及以上城市中, 16%的城市O₃日最大8 h平均浓度第90百分位数超过国家二级标准, 超标天数中以O₃为首要污染物的占17%, 成为影响我国环境空气质量的又一重要污染物(宁淼等, 2016).与传统的各行政区域各自为政单独治理相比, 区域联防联控以大气环境功能区为单位, 让该单位的省市从该功能区的整体需求出发, 共同规划实施大气污染防治计划, 协调安排, 相互监督, 相互协调, 最终实现控制复杂的空气污染并改善区域空气质量, 共享治理成果并塑造整体区域优势的目的(李馨, 2017).区域VOCs污染控制因素的偏差, 范围的扩大, 不确定性因素更多以及更加精细和复杂的控制方法是当前治理工作中最棘手的难题(刘磊等, 2019).在空气污染控制中, 人们提出了多种实用的控制方法.其中, 基于优化数学方法的理论控制策略一直是污染控制的重要研究内容(赫洁, 2012).

目前, 国内学者对大气污染的减排方案做了大量研究, 但针对VOCs单个污染物的研究仍聚焦在分行业分技术的基础研究方面, VOCs的科学联防联控研究还有待深入.陈文泰等(2013)介绍了基于二产物和基于挥发性分级两种常用的参数化估算方法, 并总结分析二次有机气溶胶的估算结果, 同时还指出未来的研究方向是如何准确量化VOCs的转化过程.我国VOCs排放源多, 成分复杂, 识别VOCs排放源特征可以通过排放源成分谱(莫梓伟等, 2014), 典型的排放源VOCs成分谱组分特征中, 不完全燃烧产物(烯烃和醛酮类)是机动车尾气的重要组分, 烯烃在我国油品挥发成分谱中的比例较高, 含氧VOCs在溶剂使用源和燃烧排放源中不容忽视, 化学组分复杂是工业排放的重要特征.从2016年开始我国对VOCs进行总量控制(张嘉妮等, 2018), 这也对VOCs排放数据的检测提出了更高的要求, 王鸣等(2019)基于8次中国城市大气VOCs外场监测, 通过挖掘VOCs浓度、组成和化学活性的内在规律, 对VOCs监测数据质量进行评估来保证数据的质量.而牵制VOCs减排总量最大的影响因素是减排产生的成本运行费用, 羌宁等(2020)从排放浓度、排放风量及连续还是间歇排放等污染源排放特征角度确定出12种典型的VOCs排放情景, 针对这些典型情景估算了目前主流减排方案的运行费用, 并以减排单位质量VOCs的运行费用为费用效能指标进行了技术经济比较, 以期为依据VOCs排放特征选择适宜的减排方案提供参考.

现在的VOCs污染控制模型大多是单目标优化模型, 主要以减排总量最大或者减排费用最低为目标.但是在实际的操作中, 一味的地讲究控制减排总量势必造成当地经济发展水平下降, 政府需兼顾考虑减排总量和治理成本的共同影响.针对上述问题, 本文在已有的研究基础上提出构建VOCs污染控制多目标优化模型, 以便能依据需要快速自动生成VOCs联防联控多级最优减排方案.对于任意给定的计算区域和时间段, 确定在多级减排总量约束下的VOCs联防联控最优减排方案.

2 区域VOCs联防联控多目标优化模型建立(Establishment of a multi-objective optimization model for regional VOCs joint prevention and control) 2.1 目标函数

区域联防联控以大气环境功能区为单位, 让该单位的省市县从该功能区的整体需求出发, 共同规划实施大气污染防治计划.建立VOCs联防联控多目标优化模型首先选取一片区域, 将该区域划分为m个子区域.子区域i的人为排放源按照排放类别划分为n_i个排放源单元.设定当前时间为T, 时间T的前a期为T-a, 其中a=1, 2, …, T-1, 时间T的后b期为T+b.假设在时间t、子区域i、排放源单元j的VOCs减排总量为X_{i, j}^t, 其中i=1, 2, …, m；j=1, 2, …, n_i；t=T-a, T-a+1, …, T+b-1, T+b.时间T的前a期VOCs的减排总量为X_{i, j}^T-a, 数据为已知.时间T的后b期VOCs的减排总量为X_{i, j}^T+b, 数据待求.设置区域挥发性有机物减排总量最大, 减排成本及减排经济影响最低为函数目标.

① VOCs减排总量最大的目标函数

(1)

式中, X_{i, j}^t为在时间t子区域i排放源单元j的VOCs减排总量.

② VOCs减排措施成本最低的目标函数

(2)

式中, k_{i, j}^t为在时间t子区域i排放源单元j单位VOCs减排成本系数, y_{i, j}为在时间t子区域i排放源单元j的减排措施导致的减排率.减排总量越大, 剩余排放量减排难度越大, 即减排成本系数也会增大.

③ VOCs减排限产产生的经济影响最小的目标函数

(3)

式中, x_{i, j}为在时间t子区域i排放源单元j的限产导致的减排率, g_{i, j}为单位减排量对应的货币化经济损失.

2.2 约束条件

① VOCs排放源O₃生成潜势约束

VOCs臭氧生成潜势与区域VOCs各组分排放量及各源谱物种反应活性有关.

(4)

式中, E_{j, l}^t为在时间t的排放源j中组分l的排放总量；X_j^t为在时间t、排放源j的VOCs排放量；P_{j, l}为排放源j中l组分占VOCs排放总量的比例因子；不同排放源的成分谱参考美国环保署(USEPA)的Speciate Database数据库.

(5)

式中, OPF_{maxj, l}^t为在时间t、排放源单元j的VOCs的最大臭氧生成潜势值, Q_{j, l}为在时间t排放源单元j的组分l的实际排放量, MIR_{j, l}^t(Carter, 2009)为在时间t排放源单元j的组分l的臭氧最大反应活性生成系数.

② VOCs减排总量非负约束

(6)

③ VOCs减排总量政策约束

(7)

式中, P_E为从时期T-a到时期T+b的VOCs总减排控制目标.

④ VOCs的减排总量应大于上一年减排总量

(8)

其中, t=T, …, T+b.

⑤ VOCs减排成本约束

(9)

其中, t=T, …, T+b.

3 模型求解(Model solving) 3.1 多目标优化模型

本研究在VOCs减排问题上, 综合考虑了年减排目标、减排成本及减排限产产生的经济影响的多目标动态优化问题描述见式(10).式(10)是复杂的多目标优化模型, 为此本文提出MOEA/DAMA算法对该模型求解.

(10)

3.2 MOEA/DAMA算法

近年来, 学者关于多目标优化做了许多工作.Srinivas (1995)于1995年提出非支配排序遗传算法, 它采用的是非支配分层方法, 这样可以使更优良的个体有更大的几率遗传到下一代, 适应度共享策略则使得准Pareto面上的优良个体分布的更加均匀, 在保持了群体多样性的同时克服了超级个体的过度繁殖, 防止早熟收敛.Deb(2002)于2002年提出一种带精英策略的非支配排序的多目标进化算法(NAGA-Ⅱ), 解决了NSGA的主要缺陷, 实现快速、准确的搜索性能.Deb(2013)提出了处理多个优化目标的进化优化算法NSGA-Ⅲ, 其主要思路是在NSGA-Ⅱ的基础上, 引入参考点机制, 对于那些非支配并且接近参考点的种群个体进行保留.张青富(2007)于2007年提出了基于分解的多目标进化算法(MOEA/D), 它将一个多目标优化问题分解成若干个标量优化子问题, 并同时对其进行优化, 通过利用与每一个子问题相邻的子问题的优化信息来优化它本身, 计算复杂度更低.Al(2011)于2011年提出强化Pareto进化算法(SPEA-Ⅱ), 该算法将非支配解存储在第二个连续更新的种群中, 然后根据支配它的非支配点的数量来评估其个体的固定性, 使用pareto优势关系保留种群多样性, 并结合聚类程序来减少非支配解的集合而不会破坏其特征.Maquez(2013)于2013年提出R2-MOGA和R2-MODE算法, 算法进一步采用种群分解策略和R2指标选择策略进行多样性管理, 在处理不同形状的PF时具有良好的性能.Wang(2015)于2015年提出多种群合作的多目标差分进化算法(CMODE), 该算法采用不用的差分变异算子和交叉策略对精英种群进行迭代, 以此来加快算法的收敛速度.

基于自适应变异率的分解多目标进化算法(Decomposition Multi-objective Evolutionary Algorithm based on Adaptive Mutation Rate, 简称MOEA/DAMA), 是在张青富教授MOEA/D算法的基础上将算法固定的变异率p_m这一参数重新设计本文提出的一种全新算法.MOEA/DAMA算法将一个多目标优化问题转换为多个标量子问题, 且每一个子问题由一个均匀分布的权重向量构成.对于每生成一个新解, 则基于聚合函数对该子问题附近的解进行替换.该算法通过分析领域问题来优化方法, 这样可以避免在优化过程中的解陷入局部最优, 所以在保持Pareto最优解的解的分布性方面具有很大优势.

同时MOEA/DAMA重新设计了变异率p_m, 使得变异率p_m在个体适应度较大时的变化程度加快, 而在个体适应度值稍小时, 其变异率p_m的变化程度减缓, 变异率p_m的表达式为：

(11)

式中, m₁, m₂为变异率调节系数, 针对不同的问题保证其变异率一直处于合理的范围内, 从而提高算法的精确性；n_p为支配个体的p个体的个数, N为种群大小, x_num为决策变量的个数, .

和传统的MOEA/D算法一样, MOEA/DAMA算法主要包括种群初始化、交叉、变异、选择4项操作.所不同的是在种群初始化设定变异率参数p_m时, 引入了自适应变异算子来调节不同支配个体的变异率.

(1) MOEA/DAMA算法的执行步骤如下：

① 种群初始化.G=0, 根据式(12)产生N个初代种群构成初始化种群, 根据式(13)产生初始参考点z.

(12)

(13)

式中, r为个体编号, r=1, 2, …, N；w=1, 2, …, x_num；X_{1, n1}_min, X_{2, n2}_min, …, X_{m, nm}_min和X_{1, n1}_max, X_{2, n2}_max, …, X_{m, nm}_max分别为指定决策变量X_{i, j}^t约束的上限与下限；rand(0, 1)为0~1之间的随机数.z_k为初始化的理想点, 称为初始参考点, k=1, 2, …, f_num, f_num为目标函数个数.

② 交叉变异.从初始种群权重向量索引中随机挑选两个序列对应的种群个体作为父本1、父本2, 对父本进行模拟二进制交叉变异获得新的个体.交叉变异之前设定交叉算子p_c为固定常数, 变异算子p_m为式(11).

(14)

(15)

式中, Y_{r, w}^t为经过交叉变异获得的新的个体, U_{r, w, p1}^t, U_{r, w, p2}^t为随机挑选的父本1、父本2, γ_{r, w}、θ分别为交叉变异参数.

③ 选择更新种群.获得的新解Y与初始非支配种群个体进行比较, 获得更优的目标函数f的参考点z, 对选择的更小的个体进行式(18)切比雪夫聚合策略见图(1)选择.

(16)

(17)

④ 终止.迭代更新至最大代数停止, 输出Pareto最优解集X.

(2) MOEA/DAMA算法的分解策略

切比雪夫聚合方法(Tchebycheff Approach)：

(18)

z^*为参考点, 计算方法见式(12), z={z₁, z₂, …, z_{f_num}}, *代表最好.也就是在各个个体中, 各个目标函数所获得的最好的目标变量作为参考点向量中的一个变量.等式左边代表的是MOEA/DAMA中的适应度值, 其含义是在已知z^*和各个权重向量的情况下, 各个个体中适应度最小的值.每个个体都有一个权重变量, 所以λ_w代表每个个体中对应第w个目标函数的权重.f_k(x)－z代表第i个目标变量到历史最优目标变量个体的距离.距离越短, 证明解越好.切比雪夫聚合过程见图 1.

3.3 VOCs最优减排模型求解过程 3.3.1 MOEA/DAMA算法参数设置

模型参数取值方法见表 1.

3.3.2 VOCs多目标优化模型求解

利用3.2节MOEA/DAMA优化算法流程对VOCs多目标优化模型求解的步骤如下：

1   Step1：程序初始化

2   Step1.1：按式(12)计算方式获得在可行空间均匀随机产生初始种群, 即U_{1, w, 0}^t, U_{2, w, 0}^t, …, U_{N, w, 0}^t, w=1, 2, …, x_num；

3    Step1.2：计算任意两个权重向量间的欧式距离, 同时查找每权向量最近的T个权重向量的索引, 令B(v)={v₁, v₂, …, v_T}, λ^v1, λ^v2, …, λ^vT是λ^v最近的T个权重向量, v=1, 2, …, N；

4    Step1.3：按照式(13)初始化参考点z=(z₁, z₂, …, z_k)^T, k=1, 2, …, f_num, z_k=min{f_k(U_{1, w}^t), f_k(U_{2, w}^t), …, f_k(U_{N, w}^t)}

5    Step1.4：检验初始化种群个体U_{1, w, 0}^t, U_{2, w, 0}^t, …, U_{N, w, 0}^t是否为非支配个体；

6   Step2：迭代更新For G=1, 2, …, max_gen, Do

7    Step2.1：父本个体选择：从B(v)中随机选取两个序号p、q, 构成两个父本U_{r, w, p}^t和U_{r, w, q}^t；

8    Step2.2：交叉：对父本U_{r, w, p}^t和U_{r, w, q}^t按照式(14)进行模拟二进制交叉操作获得子代新的个体U_{1, w, G+1}^t, U_{2, w, G+1}^t, …, U_{N, w, G+1}^t；

9    Step2.3：变异：按照式(15)将新的子代个体U_{r, w, G+1}^t进行变异操作, 变异产生新的个体Y_{r, w, G+1}^t；

10    Step2.3：按照式(16)选择更新种群：若z_k<f_k(Y_{r, w, G+1}^t), 则z_k=f_k(Y_{r, w, G+1}^t)；

11    Step2.4：按照式(17)更新领域解：若g^te(Yλ^w, z)≤g^te(Xλ^w, z), w∈B(w), 则X_{r, w, G+1}^t=z_k对应的X_{r, w, G}^t；

12    Step2.5：更新X_G+1

13    End for

14   Step3：终止条件：如果迭代更新达到最大迭代代数, 则停止并输出X, 否则转Step2.

4 关联区域VOCs多目标最优减排方案实例研究(A Case Study on the Multi-Objective Optimal Emission Reduction Plan for VOCs in Related Area) 4.1 研究区域与对象

以陕西某城市群为研究对象, 将该区域划分为6个子区域, 将每个子区域VOCs的人为排放源划分为化石燃料燃烧源、生物质燃烧、工业过程源、溶剂使用源、移动源、储存与运输源6类.在基准年每个子区域各类排放源VOCs实际减排量参考《挥发性有机物无组织排放控制标准》(GB 37822—2019)中最低控制效率标准计算结果见表 2.

4.2 数据获取与指标体系构建

采用排放因子法计算VOCs排放量(陆秋琴等, 2018).

(19)

式中, A为第i种排放源的原料消耗量或产品生产量, 即为活动水平；EF为原料类型或产品类型的排放因子.

具体的活动水平数据包括宏观统计数据和微观调研数据, 其中, 宏观统计数据来源于当地省市统计年鉴、统计部门提供的其他统计数据；微观调研数据来源于项目调研人员的走访、问卷调查及政府相关部门的统计数据; 少量无法获取的活动水平数据通过文献调研、公式推算和估算的方法获取.排放因子参考近年来大气挥发性有机物排放清单编制技术指南(2014)和国内类似的挥发性有机物排放清单研究及国内相关标准, 同时结合该区域具体情况确定具体的排放因子.计算每个子区域各类排放源VOCs实际排放量见表 3.

每个子区域各类排放源VOCs的减排成本系数的计算参考Gains-China模型(Kanada et al., 2013)计算的单位年减排成本方法来粗略计算, 其单位年减排成本包括减排设备投资(与工厂安装、维护等相关的初始成本)、固定运营成本(维护和运营减排设备所需的总成本)及可变运营成本(通过减排设备实际使用过程总产生的成本, 包括劳动力需求、运行所需的额外能源设备等), 粗略计算每个子区域各种类排放源减排成本系数见表 4.

每个子区域各排放源减排经济损失系数的计算采用CGE模型(刘昌新等, 2013), 通过CGE模型研究环境类问题的工作较多, 但近几年才有专家运用到VOCs问题中来, 借鉴CGE模型在二氧化碳减排问题中资源-经济-环境均衡的运用, 将VOCs排放源分为化石燃料燃烧源、生物质燃烧、工业过程源、溶剂使用源、移动源、储存与运输源6个部门, 同时部门进行排放源细化构建CGE模型, 计算各区域排放源减排经济损失系数见表 5.

4.3 结果分析与模型对比

将表 2~表 7中数据归一化处理后, 分别代入到多目标减排模型中, 同时运用MOEA/DAMA算法进行优化运算.MOEA/DAMA参数依据上述表 1设定的参数输入.算法采用Matlab编程, 为了测试算法的有效性, 将测试结果与MOEA/D及NAGA-Ⅱ算法的优化结果进行比较, 在结果比较之前先对数据进行处理, 即在对f₁的数据进行处理, 对NAGA-Ⅱ算法最后一代输出的Pareto解集进行升序处理, MOEA/D与MOEA/DAMA算法的Pareto解集保持不变.运算结果如图 2所示.

经过250次迭代获得了300组Pareto最优解集, 由图 2可以明显看出MOEA/DAMA算法的Pareto最优解的收敛性和覆盖品质均比NSGA-Ⅱ更好, 在解的覆盖品质上稍优于MOEA/D.

由图 3明显看出MOEA/DAMA规避了MOEA/D算法固定变异率而过早的收敛的问题, 同时收敛性能比MOEA/D和NSGA-Ⅱ更佳.

由图 4可以看出减排总量和减排措施成本呈负相关, 与减排产生的经济影响呈负相关从Parrto最优解集的概念角度考虑, 在对多目标问题进行优化时, 各个目标本身存在矛盾, 即Pareto最优解集都是最优解.融合决策者偏好的主观赋权法将决策者根据自己偏好指定的信心值分为3种类型：成本偏好、发展偏好、环境偏好.

3个偏好场景：①成本偏好型.此在这种情况下, 优先考虑成本, 并且减排的经济影响目标的优先级低于成本目标.②发展偏好型.此情景下, 行业的发展优先受到保护, 减排的同时尽可能先减少经济的影响, 最后才考虑成本.③均衡偏好型.成本目标与经济影响目标偏好度相等.

对每一类决策偏好均赋予主观权重, 本案例在每一类型决策偏好均设置一组权重来选取一组具有代表性的解, 根据其决策者偏好, 在300组Pareto最优解集中选取的3组最优解见表 8.

选取的3组解分别代表了3种不同决策者偏好的解, 每一组解均为最优解, 且均可达到VOCs区域联防联控的目的.现陕西某城市群即将举办国家级赛事, 而雾霾问题仍未有效解决, 急需将VOCs减排到可承受的最低值, 可采用序号3的减排方案, 即偏好权重为(0.67、0.25、0.08)的环境偏好型减排方案, 其年VOCs减排总量为58720 t, 具体到每个区域每个排放源的减排方案如表 9所示.

5 结论(Conclusions)

本文构建了一种新的区域联防联控VOCs最优减排优化模型, 并提出了基于自适应变异率的分解多目标进化算法(MOEA/DAMA), 来寻求区域VOCs减排模型的优化策略.本文以某地区VOCs减排为研究对象, 代入至区域VOCs联防联控多目标优化模型, 采用MOEA/DAMA算法对模型进行求解, 验证了模型和算法的可行性, 并与MOEA/D及NSGA-Ⅱ算法的优化结果进行比较, 结果证明MOEA/D在求解多目标优化问题上更加优于MOEA/D算法和NSGA-Ⅱ算法.最后根据决策者偏好在Pareto最优解集中选择了3个具有代表性的解以满足各类决策者需求, 同时为满足案例的需求根据满意度准则选择了一组解来验证模型算法的可行性、实用性、高效性.但是基于决策偏好的决策方式仍然存在大量主观判断, 且没有解决提前将决策偏好权重设置好, 从而在求解过程中直接输出决策者需求的解的问题, 后续可将决策偏好理论融入到VOCs联防联控多目标优化模型及MOEA/DAMA求解算法中, 从而更加高效、精确地输出区域VOCs最优减排方案.