东江是珠江流域三大水系之一, 发源于江西省寻乌县, 进入广东省后流经河源、惠州、东莞等城市。东江流域水系发达, 水资源丰富, 流域上游至下游土地利用类型由森林区、农田区向工业、制造区和居民密集住宅区过渡^[1]。同时, 东江是广州、深圳、东莞等沿岸经济最发达城市的饮用水源, 承担着向香港供水的任务, 总供水人口达4×10⁷人^[2]。然而, 近年来随着经济的快速发展与人口的大规模集聚, 东江流域产生了一系列环境问题, 其中突发性水污染事件风险尤为突出^[3], 据统计, 2006—2015年, 生态环境部调度的东江流域突发水污染事故约26起。突发性水污染事件具有大量污染物突然集中排放的特点, 严重影响城镇居民的生活生产^[4-6], 对流域生态环境和人体健康构成潜在危害^[7]。因此, 研究东江流域突发水污染风险具有重要的理论和现实意义。

区域水污染风险分区评估是环境风险评估中的重要组成部分, 是环境风险管理的基础和有效工具, 通过评估可识别水污染高风险的热点区域, 进而为决策者提供更具针对性的环境风险管理措施^[8-10]。目前, 国内外许多学者在区域水环境风险分区评估等领域做了大量的研究^[11-18], 如邢永健等^[18]采用集对分析法构建了基于风险场的区域突发性水污染风险评估方法, 该方法充分考虑了水体扩散的流域性特征；肖瑶等^[19]以水功能区控制单元为基础, 采用压力-状态-响应(Pressure-State-Response)环境分析模型构建环境风险源危险性指标体系, 进行流域风险源识别及危险性等级评价, 并对不同等级区提出相应的管理方案。薛鹏丽等^[10]在环境风险系统理论的指导下, 借鉴“自上而下”和“自下而上”传统区划方法对上海市突发环境风险进行区划研究, 针对上海市布局型环境风险和不同风险区提供了相应的管理措施。综上, 区域环境风险分区研究是识别突发水污染高风险热点区域的重要工具, 但现有研究大多基于县级或流域以上行政单元, 以栅格为单元的微观尺度研究相对较少, 评估结果不能较好地反映区域内风险的差异性。

基于此, 本研究以2015年为基准年, 利用环境统计数据、DEM数据、水质监测断面数据和基础地理数据, 引入地理信息系统空间分析法、区域生长法, 综合考虑了水系流向、水系级别及水质等因素, 以1 km×1 km网格为基本单元, 对东江流域开展突发水污染风险分区, 分析其区域内部风险分布特征, 以期为东江流域生态环境管理提供基础数据和理论支撑。

1 数据与方法 1.1 研究区概况

东江是珠江的主要河流之一, 地理范围为113.32°—115.53°E, 22.32°—25.11°N, 流域面积约32749 km²(图 1)。其主要支流有浰水、新丰江、秋香江、公庄河、西枝江、淡水河、石马河和增江等。各支流流经市县见图 1所示, 主要产业从流域上游以农业、农产品加工业和采矿业为主逐渐过渡到流域下游以大型产业、制造加工业、对外贸易为主^[20]。

1.2 数据来源

(1) DEM数据。来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/), 空间分辨率为1 km×1 km。

(2) 中国1:100万河流分布数据。来源于国家基础地理信息中心(http://www.webmap.cn/main.do?method=index)。

(3) 中国1:100万湖泊水库分布数据。来源于国家基础地理信息中心(http://www.webmap.cn/main.do?method=index)。

(4) 2015年国家水质监测断面数据。来源于中国环境监测总站(http://123.127.175.45:8082/)。

(5) 2015年环境统计数据。来源于中国环境监测总站。

(6) 2006—2015年生态环境部调度的突发水污染事故。来源于生态环境部。

1.3 数据预处理

(1) 2015年环统数据Q值获取

Q值指风险源环境风险物质最大存在量与临界量的比值^[21]。本研究以生态环境部2010年在全国范围内开展的“重点行业企业环境风险及化学品调查”数据及江苏省、天津滨海新区、乌鲁木齐等其他区域调查数据(企业名称、经纬度、Q值)为基础, 筛选与2015年环境统计数据中相对应的企业；然后, 以这些企业的Q值为基础, 基于行业类别代码、企业规模等参数, 推导其他企业的Q值, 即其他企业Q值取与其具有相同企业规模、相同行业已有企业Q值的平均值。

(2) 中国1:100万河流水质功能区划分

本研究基于国家水质断面数据将中国1:100万河流进行分段, 分段河流水质由相对应的国家监测断面水质确定, 即国家监测断面上游分段河流水质取该监测断面水质。

(3) 中国1:100万湖泊水库水质功能区划分

本研究基于国家水质监测断面数据确定湖泊水库水质, 如一个湖泊(或水库)有多个监测断面, 取水质最好的监测断面水质作为该湖泊(或水库)的水质。

1.4 方法

本研究参考《行政区域突发环境事件风险评估推荐方法(环办应急[2018]9号)》^[22]及基于风险场的区域突发性网格化环境风险评估方法^{[18, 23-24]}, 开展东江流域突发水污染风险分区评估, 主要有两个方面的改进：(1)引入区域生长法来具体确定风险源真实的影响范围(通过高程搜索计算)；(2)根据不同河流湖库的级别和水体功能区涉及的不同区域的敏感性确定水环境风险受体易损性指数。

1.4.1 网格划分

本研究基于东江流域区划数据, 利用ArcGIS的create fishnet功能将东江流域划分为1km×1km的网格, 并对其进行编号。

1.4.2 网格突发水环境风险场强度计算

(1) 水环境风险场强度计算

水环境风险主要通过水系(或流域)扩散, 本方法采用线性递减函数构建水环境风险场强度计算模型, 假设最大影响范围为10 km(可根据评估区域地理水文特征适当调整)^[24]。区域内某一个网格的水环境风险场强度可表示为：

(1)

式中, E_{x, y}为某一个网格的水环境风险强度；Q_i为第i个风险源环境风险物质最大存在量与临界量的比值；P_{x, y}为风险场在某一个网格出现的概率, 一般可取10^-6/a^{[18, 23-24]}(可根据评估区域风险源特征适当调整)；l_i为网格中心点与风险源的距离, 单位为km；n为风险源的个数。

由于水具有一定的流动性, 一般会从高处往低处流, 风险源所处的网格的高程应高于周边网格的高程, 才能对周边的网格造成影响。因而, 本文首先将每个风险源的影响范围限定在风险源周边10km范围内, 然后采用区域生长法来确定风险源真实的影响范围, 即从距离风险源最近的第一圈网格开始, 比较这些网格的高程值是否小于风险源的高程值, 如果小于, 则风险源周边的这个网格被视为受影响的网格, 这样依次由近至远搜索, 一旦遇到大于风险源的高程值的网格或达到了10km之外的范围, 则停止搜索。

为了便于各个网格水环境风险场强度的比较, 本方法对各个网格的水环境风险场强度进行标准化处理, 公式如下：

(2)

式中, E_{x, y}为某一个网格的水环境风险场强度；E_max为整个研究区网格的最大水环境风险场强度；E_min为整个研究区网格的最小水环境风险场强度。

1.4.3 网格环境风险受体受损性

(1) 水环境风险受体易损性计算

水环境风险受体易损性指数V_{x, y}主要参考Cao等^[24]相关研究成果, 根据不同河流的级别以及一级水体功能区涉及的不同区域的敏感性确定, 具体方法见表 1所示。

1.4.4 网格水环境风险值

网格环境风险值计算参考《行政区域突发环境事件风险评估推荐方法(环办应急[2018]9号)》^[22], 具体公式如下：

(3)

式中, R_{x, y}为网格水环境风险值；E_{x, y}为网格水环境风险场强；V_{x, y}为网格水环境风险受体分值。

1.4.5 网格水环境风险分区

本研究基于《行政区域突发环境事件风险评估推荐方法(环办应急[2018]9号)》^[22], 根据网格环境风险值的大小, 将环境风险划分为四个等级：高风险(R > 80)、较高风险(60 < R ≤ 80)、中风险(30 < R ≤ 60)、低风险(R ≤ 30)。

1.4.6 子流域提取

研究针对东江干流及主要的一级支流进行子流域提取。基于ArcGIS水文模块, 以一级支流与干流的汇流点为泄流点, 提取干流和各一级支流的不嵌套子流域(图 1)。

2 结果分析 2.1 水污染风险场强及风险受体易损性分布

东江流域突发水污染风险场强及风险受体易损性分布如图 2所示。从风险场强分布来看, 流域上游地区风险场强相对较低, 流域中下游地区风险场强较大, 风险场强较高区域主要集中分布在浰水流域西部、新丰江流域东部、秋香江流域东北部、公庄水流域西部、淡水河流域、石马河流域, 其中, 浰水流域西部风险场强较高, 主要是由于其周边选矿及采矿厂较多所致；新丰江流域风险场强较高, 主要是由于其周边水泥及钢铁企业较多所致；秋香江流域风险场强较高, 主要是由于其周边水泥及制药企业较多所致；公庄水流域西部风险场强较高, 主要是由于其周边石灰厂、纺织厂较多所致；淡水河流域、石马河流域风险场强较高, 主要是由于其周边化工厂、制造加工厂较多所致；从风险受体分布来看, 受体较为敏感的区域主要分布在东江干流两岸。

2.2 风险分区评估结果与分析 2.2.1 风险分区评估结果

通过突发水污染风险评估方法, 得到东江流域突发水污染风险评估结果, 如图 3所示。高风险区主要集中分布在龙川县、博罗县、东莞市、深圳市等地, 这主要是由于这些地区风险企业较多, 且距离东江干流较近。高风险区面积为696.7 km², 占评估区总面积的1.99%；较高风险区面积为1458.4 km², 占比4.17%；中风险区面积为2762.0 km², 占比7.90%；低风险区面积为30031.7 km², 占比85.93%(表 2)。

从各子流域平均风险值来看, 石马河流域风险最高, 其次为淡水河流域、公庄水流域(图 4)；从各子流域风险区面积来看, 石马河流域高风险区域面积最大, 其次为公庄水流域、淡水河流域(表 3)；从各子流域风险区面积所占比例来看, 石马河流域高风险区域面积占比最大, 其次为公庄水流域、淡水河流域(图 5)；石马河流域较高风险区域面积占比最大, 其次为淡水河流域、公庄水流域(图 5)。

2.2.2 高风险热点区域分析

评估结果显示, 在东江流域各子流域中, 石马河流域、淡水河流域、公庄水流域风险相对较高(图 4), 通过进一步分析, 这主要是由于这些地区企业个数较多, 总计890个, 占东江流域全部企业的23.9%, 且每100 km²企业个数远超过东江流域的平均值, 石马河流域约是平均值的4.5倍, 淡水河流域约是平均值的3.8倍(表 4)。

从沿东江干流两岸分布来看, 高风险区面积沿东江上游至下游呈逐渐增加趋势, 主要集中分布在龙川县、博罗县、东莞市、深圳市等地, 通过分析, 东江干流1km范围内风险企业个数有1048家, 10km范围内风险企业个数有5389家, 约占总数的8.7%；从沿东江主要支流两岸分布来看, 高风险区主要分布在惠州市、深圳市等地(图 3)。

2.2.3 水污染风险与水污染事件数对比分析

根据2006—2015年生态环境部调度的突发水污染事故空间分布显示(图 6), 突发水污染事故主要发生在东江流域下游, 以龙川县、惠州市、东莞市、深圳市居多, 这与本文估算的东江流域风险空间分布结果基本一致, 表明本研究计算的结果具有一定的可靠性。

3 讨论 3.1 水污染风险空间分布原因分析

从突发水污染风险空间分布来看, 石马河流域、淡水河流域风险较高, 这与段锦^[3]、刘琦^[2]等研究成果基本一致；由于突发水污染风险主要受风险源强度及风险受体敏感性两个因素影响, 进一步分析可发现, 石马河流域、淡水河流域风险较高, 主要是由于风险源场强较大所致(图 2), 其中, 石马河流域单位面积风险企业数约是平均值的4.5倍, 淡水河流域约是平均值的3.8倍；从行业类型来看, 淡水河流经深圳市、惠州市, 而石马河流经深圳市、东莞市, 以上城市均属于广东省经济发达城市, 以化学原料和化学制品制造业、制造加工业及对外贸易为支柱产业, 这些行业均为高风险行业、涉水环境风险物质使用、存储、转运量大, 同时, 取水量和排水量大, 易由生产安全事故或非法排污等引发突发水污染事件。

3.2 东江流域重大以上水污染事件与水污染风险对比

从2006—2015年生态环境部调度的重大以上突发水污染事件来看, 东江流域共发生2起, 发生在浰水流域中部、秋香江流域东北部, 与本文估算的子流域网格平均风险值有所偏差, 这主要是由于本文主要是通过建立风险源和受体之间的空间关联客观评估可能发生的所有突发水污染风险, 并未考虑各地环境风险防控与管理对风险的抑制和消减作用。对于管理因素抑制和消减作用, 薛鹏丽等^[10]、肖瑶等^[19]、毕军等^[25]人已有一些理论研究, 环境风险系统应考虑风险源危害性、传输途径、风险受体脆弱性等3个方面, 而受体脆弱性应包括暴露受体的敏感性和适应力两个因素, 适应力主要包括人口密度、经济密度、事故预警、应急能力等指标, 暴露受体敏感性越强, 脆弱性就越大, 而适应力越强, 则脆弱性越低。然而, 目前有关适应力因素的数据尚无法客观准确的量化到1km×1km的空间网格。

从流域突发水污染事件来看, 流域中上游水污染事件总数相对较少, 但重大水污染事件数较多, 流域下游水污染事件数较多, 但重大水污染事件数较少, 这可能与地区经济发展水平、科研能力、事故预警、应急能力有关, 下游地区经济发达, 科研投入、应急管理与事故预警投入较多, 应急能力较高, 使得下游地区虽然具有较高的水污染事件发生风险, 但是却鲜有重大水污染事件的发生。

需要指出的是, 由于不确定性是风险的一个本质特征, 东江流域下游等地高风险的分布特征仍可以从总体上反映区域的突发水污染风险水平, 即具有较高的发生可能性, 并在风险管控和应急措施失效的情况产生重大的后果；因此, 从应急管理上, 不仅仅需加强东江流域中上游应急能力建设, 流域下游也需继续加强。

3.3 对水环境管理的意义

与之前研究相比, 本研究考虑水的流向问题, 避免了以往研究中水环境风险影响范围过大的问题^{[8, 17]}；同时研究引入了DEM数据(数字高程模型), 并考虑到水从高处往低处流, 对水环境风险场强影响范围加以处理, 提高了评估的精度。此外, 研究在水环境风险受体赋值的过程中, 充分考虑了水环境风险受体级别及水质的影响, 相对于以往的研究, 提高了评估结果的精度。

本研究参考人口数据空间化的方法^[26], 并将该方法扩展到社会经济统计数据。统计数据的格网化打破了统计数据以行政界线为单元进行研究的限制, 使研究的空间尺度缩小到格网范围, 消除了行政界线变动带来的影响, 并可以反映区域内部的空间差异性, 提高研究结果的可信度和说服力。该方法能更好地帮助决策者辨别突发水环境风险热点区域, 为环境风险管理提供参考, 同时也可为其他类似风险评估提供一个可行的解决思路。

4 结论与展望

本研究基于2015年环统数据、DEM数据、国家水质监测断面数据和基础地理数据, 开展了东江流域突发水污染风险评估, 获得了如下主要结论：

(1) 高风险区面积为696.7 km², 占评估区总面积的1.99%；较高风险区面积为1458.4 km², 占比4.17%；中风险区面积为2762.0 km², 占比7.90%；低风险区面积为30031.7 km², 占比85.93%。

(2) 从各子流域平均风险值来看, 石马河流域风险最高, 其次为淡水河流域、公庄水流域。

(3) 从沿东江干流两岸分布来看, 高风险区面积沿东江上游至下游呈逐渐增加趋势, 主要集中分布在龙川县、博罗县、东莞市、深圳市等地；从沿东江主要支流两岸分布来看, 高风险区主要分布在惠州市、深圳市等地。

此外, 由于数据缺失的原因, 本研究仅仅评估了固定源的影响, 并未考虑危化品运输(陆运、水运)的影响, 且未采用水环境风险物质排放量来计算突发水污染风险场强, 而是利用Q值, 导致结果具有一定的不确定性, 这也是下一步需解决的问题。尽管本文存在一些局限, 但并不影响展示东江流域突发水污染风险的空间分布总体特征。