基于累积风险评价与多目标优化分配的区域水污染物总量控制

闵雪峰，俞阳*，赵锐，2，刘佳宁

（1.西南交通大学地球科学与环境工程学院，成都 611756；2.西南交通大学高速铁路运营安全空间信息技术国家地方联合共建工程实验室，成都 611756）

城乡区域是典型的自然-经济-社会复合系统，人口增长和社会经济活动加剧给水环境质量造成一定压力[1]，水污染带来的环境风险已成为制约区域经济、社会与环境协同发展的主要因素之一[2]。水污染物总量控制是改善水环境质量的有效途径，实施总量管理方案的关键在于准确识别出水污染高风险重点管控区域[3]。为减轻和改善人类活动对水环境的不利影响，开展区域水环境累积风险评价及水污染物总量控制系统研究至关重要。

水环境风险评价是系统制定水污染控制方案的依据[4]，包括环境风险评价指标体系的构建和评价指标权重的确定两个方面。压力-状态-响应（Pressure-State-Response，PSR）模型可有效体现经济发展与环境状态之间的交互影响，被广泛应用于环境系统风险评价指标体系构建[5-6]。谢蓉蓉等[7]引入风险源-风险受体-风险响应指标体系，构建了水环境综合风险评价模型。施恩等[8]基于累积性水环境风险来源、影响范围、受影响的相关者、管理方法、问题关注点等特性，构建了物理风险-监管风险-声誉风险一体化水环境风险评估体系，识别了风险影响因素及风险管控水平。由于人类活动和自然因素引起的水环境累积风险受风险源特征、水文水质、区域人口、经济发展等多因素共同影响，且各因素间存在空间异质性[9]，所以需对区域进行控制单元划分以实现区域水环境精准管控。刘思瑶等[10]构建了基于水文、水质和水生态多要素指标驱动的控制单元风险评价指标体系，并识别了四川南河小流域各空间单元的风险管控优先级。然而，当前的风险评价研究主要关注控制区污染物排放量及水质指标[11]，忽视了风险受体的资源要素、水污染利益关联主体的风险感知程度及应急响应能力，尚未从环境风险防范化解潜能角度评价水环境压力。在评价指标权重系数确定方面，常用方法包括主观赋权法和客观赋权法[12]，其中：主观法能够反映决策者的偏好，如层次分析法、专家评价法等；客观赋权法则仅根据指标本身的特征确定其重要性，如熵值法、变异系数法等。在实际应用中，主观法获得的权重反映出决策者对不同指标的重视程度，但具有较强不确定性；客观方法虽然能有效避免主观因素影响，但权重值会与真实值产生差距。现有研究常采用单一方法分析权重[13]，一定程度上忽视了定性、定量指标的客观性，限制了对环境污染事故可能造成的人体健康风险的准确评价。

污染负荷分配是以水环境容量或目标控制总量为基础，通过制定决策方案将污染负荷分配至各区域或污染源[14]。国外学者对水污染负荷分配的研究多以经济收益为基础，实现污染负荷减排最大化或/和污染治理费用最小化[15]。我国水污染总量控制是分层逐级开展的，分配过程充分考虑经济最优、公平、可持续和技术可行性等原则[16]，存在多个相互冲突、不可通约的目标，属于多目标优化问题（Multi-objective Optimal Problems，MOPs）[17]。2012 年，Zhang 等[18]提出的多目标污染负荷分配方法首次实现了流域、水功能区和排污口多对象融合，确定了公平和高效的分配方案。随后，Meng 等[19]开发了兼顾公平与效率的两阶段随机规划模型，并将其应用于四个主要排污断面（工业、市政、畜牧业和农业）的污染负荷分配。Wu等[20]基于公平性原则建立了面向流域和区域两层利益主体的多指标基尼系数污染负荷分配框架。张璇等[16]针对环境保护税法与水功能区水质达标实施背景下的水污染负荷分配问题，采用双层多目标优化模型为各管理层提供了分配方案。近年来，为了寻求MOPs 的Pareto 最优解，学者们提出了多目标进化算法（Multi-ObjectiveEvolutionary Algorithms，MOEAs），并将其广泛应用于污染控制及水资源规划和管理[21]。其中，Bi等[22]提出的基于淘汰算子的改进非支配排序遗传算法Ⅲ（Improved NSGA-Ⅲalgorithm based on Elimination Operator，NSGA-Ⅲ-EO）能有效处理约束MOPs，显著保持解集的多样性和收敛性，同时能克服维数灾难和计算成本高的问题。本研究将其应用于求解构建的多目标水污染负荷分配模型。

目前，国内水环境治理成效显著，但水污染物排放量依然偏高[23]，其中四川省2017 年化学需氧量（CODCr）、氨氮（NH3-N）和总磷（TP）排放量分别高达116.36 万、5.91 万t 和1.58 万t[24]，污染物总量控制工作依然面临严峻形势。此外，现有区域水环境管理研究尚未在识别水污染重点管控单元的基础上开展污染负荷优化分配研究[25]，使得环境管理部门缺乏高效可行的水污染防治与管理决策体系，难以切实指导区域水环境管理策略的制定。鉴于此，为便于区域水环境精准管控，本研究选择四川省典型乡镇行政单位为最小控制单元，将环境风险防范潜力纳入评价指标体系，采用主客观赋权耦合分析法计算权重系数，通过累积综合风险评价确定各控制单元的风险等级，识别重点风险源；结合经验系数类比法及产排污系数法获得不同类型点源、面源污染物入河量，基于“十三五”规划的污染总量控制目标，构建多目标污染负荷优化分配模型，得出科学合理的优化分配方案，为制定区域水环境质量改善策略提供科学参考。

1 材料与方法

1.1 研究区概况

地处长江上游青衣江流域的四川洪雅县（102°49′～103°32′E，29°24′～30°00′N）属于典型的农业与城市污染源叠加区域，辖区覆盖青衣江58.8 km 长的河段，属于中亚热带湿润气候区，年平均气温16.9 ℃，年平均降雨量1 759.8 mm，产业结构依靠第二产业带动第一、第三产业协同发展。研究以洪川镇、余坪镇、将军镇等12 个乡镇为基本控制单元，区域内经济高速发展导致污染问题频发，青衣江干流部分水质指标不稳定，难以满足水质功能要求。

考虑到断面水质基本为2017 年之前监测的数据，因此本研究选用2017 年洪雅县行政区划进行水污染总量分配。

1.2 数据来源

本研究所用数据包括该地区的社会基本情况数据、社会经济发展数据、污染物排放量及水质数据、政府监管数据等，这些数据主要来源于第二次全国污染源普查数据、2011—2019 年青衣江木城镇断面水质监测数据、《2020 年眉山统计年鉴》等。参数赋值主要参考《眉山市城市总体规划（2009—2020）》《全国污染源普查城镇生活源产排污系数手册》《眉山市水污染防治目标责任书》以及洪雅县人民政府的公开报道和文件。

1.3 区域水环境综合风险评价

1.3.1 区域水环境累积风险评价指标体系

本研究在明确区域环境风险受体自然要素、利益关联主体风险感知及应急响应能力的基础上，以水环境累积风险为目标层，建立了包含压力（风险受体自然要素）、状态（风险源体排污现状）、响应（风险受体风险防范潜力）3 个准则层的水环境累积风险评估指标体系（图1）。

图1 水环境累积风险综合评价指标体系Figure 1 Index system for comprehensive evaluation of cumulative water environment risk

压力指标描述人口发展和经济发展等人类活动或自然因素给环境带来的胁迫或影响，包括人口增长和社会经济发展指标，其中经济发展方向和产业结构设置决定了污水排放量和污染负荷类别。状态指标描述人类活动对水环境质量的影响，包含不同类型污废水排放量、畜禽养殖量、耕地面积和CODCr、NH3-N、TP 的排放量。响应指标反映各关联主体感知到环境风险时，能够通过防治污染和修复环境来持续改善水环境质量的潜力。在第一、第二、第三产业并存的城乡污染叠加影响区域，居民、政府、企业是3 个重要利益主体，因此，可选取居民环境保护意识、政府污水处理效率、环境监测能力以及应对突发污染安全事故的应急响应能力和企业的风险防控能力代表各主体的环境风险防范潜力。

1.3.2 指标风险等级划分

为比较各控制单元的相对风险等级，进行分级、分区精准管理，在水环境累积风险等级评分标准中，将所有控制单元的综合风险评分阈值（Rmin,Rmax）划分为极低风险、低风险、中风险和高风险4 个等级[7]（表1）。

表1 累积风险等级评分标准Table 1 Criteria for the cumulative risk evaluation

控制单元中各单项指标风险分级主要依据国家标准或地方标准，以及风险评价相关的研究[26]，包括《国民经济行业分类》（GB/T 4754—2017）、“关于印发《上市公司环保核查行业分类管理名录》的通知”（环办函〔2008〕373 号）、《畜禽养殖业污染物排放标准》（GB 18596—2001）、《环境影响评价技术导则地表水环境》（HJ 2.3—2018）、《地表水环境质量标准》（GB 3838—2002）等，各项指标风险分级见表2。

表2 区域水环境累积风险评价指标分级Table 2 Index classification of cumulative risk assessment for regional water environment

1.3.3 综合风险评分法

综合风险评价首先通过数学理论模型确定各指标权重，对控制单元内各指标的评分与相应权重相乘并求和，得到综合风险评分：

式中：Ri为控制单元i的综合风险值；ωj表示指标j的权重；rj表示各指标的风险评分；p为评价指标个数。

本研究采用耦合模糊层次分析法和变异系数法确定各指标权重，具体计算公式如下：

式中：uj是指标j由模糊层次分析法获得的权重；vj是指标j由变异系数法确定的权重。

（1）模糊层次分析法

模糊层次分析法是定性与定量相结合的多因素决策分析方法，根据多专家对评价指标的两两比较，形成一个由重要系数组成的判断矩阵，可应用于不精确数值的定量和定性权重分析，具体计算步骤参照相关研究[27]。

（2）变异系数法

变异系数利用各指标所包含的信息，通过比较数据的标准差与均值来反映各评价指标的相对重要性，避免专家在确定各指标权重时的主观性，同时可削弱极值指标对评价结果的影响，具体计算步骤参照相关研究[28]。

1.4 多目标水污染负荷优化分配

1.4.1 多目标优化分配模型构建

污染负荷分配应遵循公平与效率原则[14]，其中公平原则体现为各控制单元的综合环境基尼系数和削减率不均度最小，效率原则体现为经济产值最大。各目标函数及约束条件如下：

（1）环境基尼系数最小

均衡考虑资源、经济、社会发展对污染物排放量的影响，选取人口、GDP、土地面积3 个指标作为评价水污染物总量分配方案公平性的资源因素[29]，其综合环境基尼系数（G）最小可表示为：

式中：m为资源因素的个数，即m=3；wk表示资源因素k对污染物排放量的贡献，可采用信息熵计算，具体公式如下：

式中：Wi,0表示控制单元i中的污染物现状入河量；Zi,k表示控制单元i中评价因素k的指标值；Pi,k表示控制单元i中单位评价因素k的污染物入河量；Ni,k为控制单元i中Pi,k在各评价因素中所占比例；Ek为评价因素k单位污染物入河量的信息熵。

基尼系数Gk的计算方法包括梯形面积法和积分法，对于资源划分份数较少的评价模型通常采用梯形面积法进行计算[16]：

式中：Gk为资源因素k的基尼系数；Xi,k为控制单元i中指标k占其总量百分比的累计值；Yi为控制单元i中污染物分配量占其总量百分比的累计值；n为控制单元个数，n=12。当i=1时，设Xi-1,k、Yi-1等于0。

（2）削减率不均度最小

式中：σ为各控制单元削减率的方差；pi表示控制单元i的污染物削减率；pˉ表示削减率pi的均值；xi表示控制单元i的污染负荷分配量；xi,0表示控制单元i的污染物现状入河量。

（3）经济产值最大

式中：Ci表示控制单元i中单位污染负荷产生的国内生产总值；GDPi表示控制单元i的国内生产总值，亿元。

（4）约束条件

污染负荷分配模型应满足分配总量不超过允许入河总量、削减率不超过有效范围、环境基尼系数不高于现状值的条件，即：

式中：XT为满足削减目标的年排放总量。

“十三五”规划要求研究区总削减率为8.83%，考虑各行政控制单元的经济与环境管理效率差异以及各控制单元的实际排放情况和纳污能力，设置各控制单元削减率上限、下限分别为20%、1%[29]：

基于公平原则，优化分配方案的综合环境基尼系数应不大于分配前的基尼系数值，即：

式中：G为优化分配方案的基尼系数；G0为现状综合基尼系数。

1.4.2 水污染负荷分配模型求解与方案优选

本研究采用NSGA-Ⅲ-EO 算法求解多目标水污染负荷分配模型，算法详细过程参照相关研究[22]。设置种群数量800，迭代次数60，可获得一组Pareto最优解。随后分别以环境基尼系数最小（情景一）、削减率不均度最小（情景二）、经济产值最大（情景三）和综合最优（情景四）四种方案开展情景分析。

实际应用中通常只需一组优化分配方案用于指导区域水污染管控。本研究采用模糊集理论对Pareto解集中各最优解的目标函数值进行满意度评价，确定综合最优分配方案。目标函数最小化的模糊隶属度函数[30]为：

式中：un表示第n个目标函数fn的隶属度值；fn为Pareto 最优解中第n个目标函数的函数值；分别表示第n个目标函数值fn的最小值和最大值。

对于Pareto 解集中的每个最优解，采用式（16）求出该解对应各目标函数的隶属度，根据式（17）求出该解的综合满意度，Pareto 解集中综合满意度最大值对应的解即为最优分配方案。

式中：uN为第N个Pareto 解的满意度值；N为Pareto 解集中最优解的个数。

2 结果与讨论

2.1 控制单元水环境累积风险评价

分别采用模糊层次分析法和变异系数法计算各指标权重，进而得到各指标耦合权重系数（表3）。

表3 水环境风险评价指标权重系数Table 3 Weighting coefficients of the water environment risk evaluation index

由表3 可知，在指标层中，人口发展，经济发展，水环境质量中的污废水排放量、养殖量和耕地面积指标，以及居民响应指标中的生活污水处理率均为定量指标[7]，其耦合权重系数与变异系数法的分析结果更接近，保障了各指标之间的权重公平。居民响应指标中的公众环保意识和政府响应指标为定性指标，其耦合权重系数与模糊层次分析法的结果更接近，表明主客观方法的耦合使指标权重系数的确定更具合理性和科学性。

基于各指标风险等级评分和相应权重，计算研究区12 个控制单元累积综合风险评分。根据表1 风险等级评分标准，确定高风险和中风险控制单元分别为2 个和3 个，低风险和极低风险控制单元分别为3 个和4个，各控制单元的风险分布及评分值见图2。

如图2 所示，位于研究区青衣江流域上游的槽渔滩镇控制单元综合风险最低，下游洪川镇和余坪镇风险等级较高，其各评价指标的风险值如图3 所示。分析可知，洪川镇的高环境风险是由工业废水排放中高负荷的CODCr、NH3-N 和TP 造成的，余坪镇累积环境风险较高是因为规模化畜禽养殖和工业生产产生了大量的CODCr、NH3-N和TP负荷。由此可见，各控制单元的环境风险不仅受内部水污染物排放量影响，还受流域内水污染累积效应影响，这与尚彦辰等[26]分析辽河流域水环境累积风险时得到的结论一致。此外，高风险控制区内居民、政府和企业的风险防范潜力有待提升，后期应将每个控制单元的环境风险控制纳入区域风险统筹管理，提升利益相关者的系统性风险防控能力[31]。

图2 控制单元综合风险等级空间分布图Figure 2 Spatial distribution of the comprehensive risk evaluation result for the control units

图3 高风险控制单元的评价指标值Figure 3 Evaluation index value of high risk control units

2.2 多目标水污染负荷分配方案情景分析

《四川省水污染防治目标责任书》要求2020 年青衣江流域水质目标为：TP≤0.33 mg·L-1，其他指标达到《地表水环境质量标准》（GB 3838—2002）Ⅱ类水质。由洪雅县2017年重点监控断面水质可知，CODCr和TP浓度明显低于目标限值，NH3-N 浓度存在超标风险[32]。王玉云等[33]在青衣江流域水环境生态补偿研究中也将NH3-N 作为测算流域生态补偿的水质因子，因此本研究选用NH3-N 为代表性污染物控制因子进行水污染负荷分配。根据当地统计年鉴和《第二次全国污染源普查公报》数据，基于经验系数类比法及产排污系数法获得基准年各控制单元不同类型点源、面源NH3-N 入河量，汇总获得研究区NH3-N 入河总量[34]。《眉山市环境保护“十三五”规划》要求2020年洪雅县NH3-N 削减率达到13.9%，假定“十三五”期间NH3-N 削减率逐年相等，可得2017 年NH3-N 允许入河总量为622.05 t。由NSGA-Ⅲ-EO算法获得的所有分配方案的Pareto前沿面见图4。

图4 NSGA-Ⅲ-EO 获得分配方案的Pareto前沿Figure 4 Pareto front of allocation schemes obtained by NSGA-Ⅲ-EO

根据多目标优化分配模型，计算得到4 种情景下洪雅县各控制单元NH3-N 总削减率分别为16.4%、10.3%、9.1%和16.3%，进一步得到不同情景下各控制单元削减方案（图5）。由图5 可知，NH3-N 入河总量中洪川镇贡献度最高，为27.3%；其次为余坪镇、将军镇和东岳镇，分别占15.7%、11.2%和10.8%；其余8个控制单元的占比均小于10%，并且七里坪镇的占比最低，仅为2.0%。分析可知，情景三追求经济产值最大时，洪雅县NH3-N 总削减率最低，其中NH3-N 入河量贡献度最高的洪川镇的削减率仅为3.2%，表明区域仅追求经济利益时忽视了环境影响，可能会增加区域环境风险。情景二各控制单元均匀削减污染物排放量时，虽然各控制单元的削减率接近，但洪雅县总削减率较低，原因在于低排放控制单元减污能力有限，同时未充分利用高排放控制单元的控污能力[14]。在4种情景中，综合最优方案（情景四）与环境基尼系数最小方案（情景一）的总削减率较高，且环境基尼系数较小，表明考虑各控制单元的人口、GDP、土地面积等资源要素可较大限度地实现区域水污染物总量削减[16]。此外，情景四中各控制单元削减率方差为0.002 9，小于情景一的方差值0.004 6，表明综合最优方案下各控制单元之间削减率差异更小，各控制单元中的相关利益主体更易接受，方案有效性和可执行性更强[14]。

图5 各控制单元的NH3-N入河量占比及不同情景分配方案Figure 5 The allocation scheme of NH3-N of control units under different scenarios and proportion of discharge

基于上述分析，选择综合最优分配方案作为区域水污染总量控制的指导方案，即NH3-N 总削减量为111.15 t，整体削减率为16.3%。此时洪川镇控制单元的NH3-N 削减率最高，为19.7%，其次为余坪镇，削减率为18.9%。程一鑫等[29]研究发现区域污染物分配方案与现状排放量之间不具有一致性，主要原因在于其仅考虑区域人口、资源、经济和水污染物承受能力等客观因素，忽视了各乡镇排污量的经济产值效率及削减率公平性。本研究中优化后的控制单元NH3-N 削减比例排序与综合累积风险评价等级基本一致，即污染物入河量贡献越大，综合环境风险等级越高，污染负荷削减比例相对较高，证实了分配方案的公平性、效率性和合理性。

3 结论

本研究结合环境累积风险评价与水污染负荷多目标分配模型，提出了基于高风险优先管控区识别的区域水污染系统管控方法，对四川洪雅县的NH3-N总量进行分配，结论如下：

（1）研究区存在2 个高风险和3 个中风险控制单元，高风险控制单元为洪川镇和余坪镇。其中洪川镇控制单元的主要风险源为工业废水排放，余坪镇的风险源是规模化畜禽养殖及工业废水排放。

（2）针对高风险控制单元，应加强重点行业、企业整治和排污口监测，加强对畜禽散养的监管力度，探索建立分散养殖粪污收集、贮存、处理与利用体系。针对中风险控制单元，应加快完善乡镇污水管网建设和污水处理厂建设，提高生活污水处理率，同时提升公众环保意识，提高政府对突发污染事件的应急响应能力。

（3）在综合最优分配情景中，洪雅县NH3-N 总削减量为111.15 t，削减率为16.3%，其中洪川镇和余坪镇控制单元的削减率分别为19.7%和18.9%。本研究在实现各控制单元公平、高效削减水污染物的同时，获得了较高的区域总削减率，说明分配方案科学合理，可切实指导区域水污染物总量控制方案的制定。