随着全球经济快速发展, 水环境污染已成为一个亟需解决的全球性问题, 受到学者和管理者关注^[1].水体污染按污染形式可以分成点源污染和非点源(面源)污染, 近年来, 随着点源污染管控措施逐步完善, 非点源污染已成为全球水环境防治的最大挑战之一^[2~4].相较于点源污染, 非点源污染具有来源广泛、暴发具有随机性、空间差异明显和监测困难等特点, 使得非点源污染防治难度远远高于点源污染^{[5, 6]}.目前, 国内外非点源污染研究以非点源污染负荷定量估算和非点源污染风险定性评估为主^{[7, 8]}, 但有研究认为, 相对于非点源污染负荷的精确估算, 如何准确识别污染输出高风险区域并制定相应防治措施才是非点源污染研究的首要工作^{[9, 10]}.

输出风险模型作为非点源污染输出风险评估常用方法之一, 以输出系数法为基础, 核心思想是根据土地利用类型空间分布特征实现区域非点源污染输出风险评估^{[11, 12]}.该方法不考虑非点源污染形成及迁移转化等复杂物理过程, 模型具有结构简单、基础数据易获取及参数较少等特点, 且模型对时空尺度不敏感, 可移植性好^{[13, 14]}.目前, 基于输出风险模型的非点源污染研究, 大多已考虑降水和地形时空差异影响, 并通过引入新的模型参数, 对经典的输出风险模型进行改进, 模型模拟精度得到了较大改善^{[11, 15~17]}.但通过调研, 发现目前研究存在以下3点不足之处: ①土地利用数据空间分辨率和精度不足, 大多数研究以30 m分辨率Landsat系列卫星影像为基础, 采用单一分类器提取土地利用数据; ②基于气象站点插值的降水数据精度不足, 尤其在站点数据稀少区域, 这使得降水因子计算存在不确定性; ③目前非点源污染风险研究以年尺度为主, 但实际上, 大多数区域年内降水时空分布差异较大, 非点源污染风险也存在较大差异.

水库作为人类赖以生存的资源, 在防洪、灌溉和城市用水等方面起到重要作用^[18].潘家口水库是“引滦入津”工程主要水库, 近年来由于受上游农业和城市排污等影响, 水库水质持续恶化^{[19, 20]}.本文以潘家口水库流域为研究区, 首先以高空间分辨率(16 m)国产GF-6 WFV影像作为土地利用分类基础数据, 采用多分类器组合方法, 实现研究区土地利用高精度提取; 其次以高时空分辨率(时间分辨率0.5 h, 空间分辨率0.1°)的卫星反演降水产品GPM(global precipitation measurement mission, GPM)作为流域降水基础数据, 通过PSO-BP(particle swarm optimization- back propagation, PSO-BP)降尺度方法, 得到研究区高精度降水数据; 最后构建了考虑流域降水和地形差异的输出风险模型, 区别于传统研究, 在月尺度上评估潘家口水库流域非点源污染输出风险, 并提出针对性的污染防控措施, 以期为流域非点源污染研究和水环境保护提供科学参考.

潘家口水库流域(115°~119°E, 40°~43°N)位于海河之北, 滦河中上游, 跨内蒙古自治区和河北省, 主要面积在承德市、张家口市、赤峰市和锡林郭勒盟, 流域控制面积约33 400 km², 分别占滦河流域和海河流域面积的75%和10%.流域海拔147~2 236 m, 由西北向东南递减, 流域地处亚热带季风区, 夏季炎热多雨, 冬季寒冷干燥, 年平均气温9.6℃, 年蒸发量400 mm, 流域多年平均降水量540 mm.年内降水集中, 7月和8月降水量占年降水量的50%以上, 且流域降水空间分布不均.研究区地理位置及海拔如图 1所示.

图 1

Fig. 1

图 1 研究区地理位置及高程 Fig. 1 Location and elevation of the study area

传统输出风险模型适用于降水均匀和地势平坦区域, 而随着模型的广泛应用, 降水和地形因子已成为输出风险模型不可忽视的模型参数^[21].因此, 本研究综合考虑潘家口水库流域降水时空差异以及地形空间差异, 引入动态降水因子和地形因子对模型进行改进, 构建输出风险评估模型, 识别了潘家口水库流域非点源污染输出风险时空分布特征.本研究技术路线如图 2所示.

图 2

Fig. 2

图 2 技术路线示意 Fig. 2 Technique flow chart

本研究所需数据包括潘家口水库流域DEM、18景GF-6 WFV影像、卫星反演降水产品GPM、降水站点数据和NDVI(normalized difference vegetation index, NDVI), 所有数据时间均为2018年.数据处理包括: ①子流域是输出风险模型的最小计算单元, 以DEM为基础数据, 通过ArcGIS 10.6软件进行填洼、流向计算、汇流量计算、河网水系提取和河网分级等处理, 将研究区划分为105个子流域.②土地利用分类是输出风险模型的基础, 以18景GF-6 WFV影像为土地利用提取基础数据, 采用多分类器组合方法^[22], 将研究区土地利用分为水体、草地、耕地、林地、建设用地和未利用地这6类.③降水作为非点源污染产输的主要驱动力, 无疑是最重要的参数之一, 采用文献[23]的卫星反演降水产品GPM空间降尺度PSO-BP模型, 得到研究区高精度降水数据.④土壤侵蚀是非点源污染产输的主要途径之一, 一般认为土壤侵蚀强度与坡度正相关, 即坡度越大, 发生侵蚀的可能性越高^[24].以DEM为基础, 利用ArcGIS 10.6软件中Spatial Analysis工具计算研究区坡度, 作为模型地形因子计算的基础数据.

输出系数法是输出风险模型的基础, 输出系数法最早由Frink提出, 经Johnes等改进, 形成经典输出系数模型^[25].模型不涉及非点源污染复杂的产输机制, 参数少、计算简单、模拟精度高, 且对监测数据要求不高, 所以在大尺度非点源污染研究中应用广泛^{[13, 26]}.非点源污染风险评估需要明确的阈值, 有研究表明, 相对于其他土地利用类型, 林地产污少, 且对污染物的迁移有一定抑制, 所以一般认为林地是非点源污染的“汇”, 其他土地利用类型是非点源污染的“源”^{[27, 28]}.本文在非点源污染风险计算中, 以林地输出系数作为风险计算阈值, 大于则认为存在输出风险, 反之则不存在输出风险^{[13, 29]}.考虑降水动态因子和地形因子的输出风险模型计算公式为:

(1)

(2)

式中, L_k为第k个子流域内污染物j带权重的平均输出系数; j为污染物类型, 本研究为总氮(total nitrogen, TN)和总磷(total phosphorus, TP); i为流域内土地利用类型; E_ij为污染物j在第i类土地利用类型的输出系数; A_ik为第k个子流域内第i种土地利用类型的面积百分比; α_k和β_k分别为第k个子流域内降水因子和地形因子; G_k为第k个子流域非点源污染输出概率; E_h表示非点源污染输出计算阈值, 即林地的最大输出系数.

降水时空差异对非点源污染形成与迁移影响显著, 本文动态降水因子(α_k)计算需要从时间差异(α_k^t)和空间差异(α_k^s)两个方面进行考虑^{[30, 31]}, 计算公式如式(3)所示:

(3)

式中, α_k为第k个子流域的某月降水因子; α_k^t和α_k^s分别为第k个子流域降水因子的时间差异因子和空间差异因子; R_k^month和R_k^month分别为第k个子流域的某月降水侵蚀力和月平均降水侵蚀力; R^month为全流域某月降雨侵蚀力.降水侵蚀力计算参考徐丽等^[32]降水侵蚀力简易计算方法, 计算公式如式(4)所示:

(4)

式中, R为降水侵蚀力; P为月降水量.计算潘家口水库流域2018年降水影响因子如图 3所示.从中可以看出, 从4~9月, 降水因子整体呈现先增后减趋势, 且降水因子空间差异显著, 其中7月和8月降水因子相对较高.

图 3

Fig. 3

图 3 2018年潘家口水库流域降水因子 Fig. 3 Factors of precipitation in Panjiakou Reservoir basin in 2018

有研究表明^[33], 地形与非点源污染产输存在显著关系, 并采用坡度因子来表示非点源污染由于下垫面起伏而产生的空间差异, 计算公式如式(5)所示:

(5)

式中, β_k为第k个子流域的地形因子; S_k^b为第k个子流域的平均坡度; S^b为全流域平均坡度; b为常数, 参考相关非点源污染研究中地形因子计算方法^{[34, 35]}, 取b=0.610 4.计算潘家口水库流域地形因子如图 4所示.从中可以看出, 潘家口水库流域地形因子取值在0.41~1.50之间, 流域上游位于内蒙古坝上高原, 地势平坦, 所以坡度因子相对较小; 中下游以丘陵和山地为主, 地貌复杂, 坡度因子相对较大.

图 4

Fig. 4

图 4 地形因子 Fig. 4 Topographic factor

输出系数确定是构建非点源输出风险模型的关键步骤, 国内外输出系数确定方法包括野外监测、统计分析和文献查阅.由于研究区监测数据不足, 本文输出系数采用文献查阅方法确定, 结合文献[36]以及国内学者在相似区域的研究成果^{[21, 37~40]}, 确定非点源污染输出系数(表 1).

表 1 (Table 1)

表 1 潘家口水库流域非点源污染输出系数 Table 1 Export coefficients of non-point source pollution in Panjiakou Reservoir basin 项目 土地利用类型 耕地 林地 草地 建设用地 未利用土地 TN 29 2.38 10 11 14.9 TP 0.9 0.15 0.2 0.24 0.51

表 1 潘家口水库流域非点源污染输出系数 Table 1 Export coefficients of non-point source pollution in Panjiakou Reservoir basin

目前, 非点源污染输出风险并没有明确的等级划分, 本文参考方广玲等^[41]对非点源污染输出风险等级划分, 以0.2为间隔, 将TN和TP污染输出风险概率划分为5个等级, 分别为低风险(0~0.2)、较低风险(0.2~0.4)、中风险(0.4~0.6)、较高风险(0.6~0.8)和高风险(0.8~1.0).国内外学者对非点源污染输出风险研究多集中在年尺度, 较少对月尺度污染输出风险进行评估, 故本文在年尺度和月尺度分别计算了潘家口水库流域非点源污染输出风险, 分析研究区非点源污染输出风险时空分布特征.

2018年潘家口TN和TP的非点源污染输出风险空间分布如图 5所示.总体而言, 两种污染物的输出风险均在流域的上游和下游较高, 中游较低.结合行政区划分和非点源污染输出风险面积统计, 从图 5(a)可以看出, TN污染输出高风险和较高风险区分别占流域总面积的23.9%和46.7%, 主要分布在流域上游多伦县、围场满族蒙古族自治县西部、沽源县东北部、丰宁满族自治县东部和流域下游的滦平县、鹰手营子矿区、双滦区、宽城满族自治县、兴隆县东北部, 这些区域土地利用方式以产污量较高的耕地和城市为主, 所以非点源污染输出风险相对较高.中风险区主要集中在流域中游, 约占流域总面积的21.9%, 主要分布在围场满族自治县西南部、隆化县、承德县和平泉县西部, 这些区域土地利用方式以林地为主, 村落和耕地分散在整个区域, 虽然林地对污染物有一定拦截能力, 但由于山区坡度较大, 地表径流对污染物迁移能力较高.较低风险和低风险区面积约占流域总面积的7.5%, 主要集中在流域中游的围场满族自治县中部、隆化县东部和承德县东北部, 这些区域林地占比高, 且多为天然林地, 不需要施肥, 污染物来源较少, 非点源污染风险相应也就较低.

图 5

Fig. 5

图 5 2018年潘家口水库流域TN和TP污染输出风险 Fig. 5 TN and TP pollution output risk of Panjiakou Reservoir basin in 2018

从图 5(b)可以看出, TP污染输出无高风险区.较高风险区和中风险区分别占流域面积的29.4%和38.0%, 主要集中在土地利用类型以耕地为主的上游、中游西部和下游城市区, 包括正蓝旗北部、太仆寺旗和沽源县东部、多伦县、围场满族蒙古族自治县西部、丰宁满族自治县、隆化县中部、滦平县东部、兴隆县和宽城满族自治县等.较低风险区主要集中在流域的中游和下游东部区域, 包括围场满族蒙古族自治县中部、隆化县、承德县西部和平泉县西部等, 约占流域总面积的26.9%.低风险区和无风险区分别占流域总面积的5.2%和0.5%, 集中分布在流域中游东南部, 包括隆化县东部和承德县北部.

潘家口水库流域2018年4~9月TN和TP污染输出风险时空分布如图 6所示, 污染输出风险等级面积占比统计结果如图 7所示.结合图 6和图 7可以看出, 4月TN和TP污染输出的无风险区域均超过全流域面积的一半, 分别达到60.7%和77.9%; 低风险区分别为9.7%和4.2%; 较低风险区和中风险区面积相当, 分别约为9.1%和8.7%; TN污染输出的较高风险区约为11.6%, 主要位于流域下游; 两类污染物输出均无高风险区.5月降水量较小, 污染输出风险相对较低, 流域内污染输出均无较高风险区和高风险区, TP污染输出仅存在无风险区和低风险区.其中, TN和TP污染输出的无风险区占比分别为68.4%和94.8%, 低风险区占比分别为5.9%和5.2%; TN污染输出的较低风险区和中风险区分别为6.6%和19.1%.6月流域汛期开始, 污染输出风险明显高于4月和5月, 但仍不存在高风险区域.其中, TN和TP污染输出的无风险区域面积明显减少, 分别为23.4%和55.0%; 低风险区、较低风险区、中风险区面积占比分别为21.2%、16.3%、9.3%和14.9%、19.3%、10.8%, TN污染输出较高风险区约为29.7%.

图 6

Fig. 6

(a)TN污染输出风险等级; (b)TP污染输出风险等级 图 6 潘家口水库流域TN和TP污染输出风险时空分布 Fig. 6 Temporal and spatial distribution of TN and TP pollution output risk of Panjiakou Reservoir basin

图 7

Fig. 7

图 7 流域非点源污染输出风险等级面积占比 Fig. 7 Proportion of area of non-point source pollution output risk grade in watershed

潘家口水库流域2018年7月和8月总降水量达到年降水量的50%以上, 非点源污染输出风险相应也较高^[38].其中, 7月污染输出风险等级均在较低风险以上, TN和TP污染输出的较低风险区、中风险区、较高风险区、高风险区面积占比分别为0.9%、12.6%、40.0%、46.4%和6.0%、36.3%、46.8%、10.9%.流域TN污染输出以较高风险和高风险为主, 总面积占比达到86.4%, TP污染输出以中风险和较高风险为主, 约占全流域的83.1%.8月降水量最大, 相应的污染风险也达到最高, 其中TN和TP污染输出的中风险、较高风险和高风险区面积占比分别为3.1%、41.7%、55.2%和23.7%、67.4%、7.8%, TP污染输出的较低风险区约为1.1%.9月是流域汛期的最后一个月, 降水量相对于7月和8月明显下降, 污染输出风险也相应减小, TN和TP污染输出的无风险区、低风险区、较低风险区、中风险区面积占比分别为17.4%、25.4%、19.1%、10.1%和62.8%、7.9%、23.0%、6.2%, TN污染输出的较高风险区约为28.0%; 其中TN各等级污染输出风险面积占比平均为20%, 而TP污染输出以无风险区为主, 较低风险区次之.

结合图 6和图 7分析, 总体而言, 从4~9月TN和TP污染输出风险呈现先增大后减小, 趋势与年内降水量变化一致, 这与张婷等^[42]利用SWAT模型的模拟结果一致, 呈现出“汛期高, 非汛期低”的特点.对比分析不同类型污染物输出风险, 发现研究区TN污染输出风险均比TP高, 这与杨金凤等^[43]在海河流域的研究结果一致.分析不同等级污染风险区内主要的土地利用类型, 发现以林地为主要土地利用类型的中游, 非点源污染输出风险相对较低, 而在以耕地和建设用地为主的上游和下游, 非点源污染输出风险较高.

充分考虑流域降水和地形的时空差异, 运用改进的输出风险模型, 评估潘家口水库流域2018年非点源污染输出风险时空分布特征, 并结合流域降水分布和农业活动, 提出非点源污染防控措施.模型模拟过程中土地利用、降水等数据均来源于高分辨率卫星遥感数据, 能够准确得到流域非点源污染风险时空分布特征, 合理划分污染风险区.从结果来看, 流域内以耕地和建设用地为主的区域污染风险更高, 以林地为主的区域风险相对较低, 这与大多数流域非点源污染风险评估研究的结论一致^{[13, 27, 40, 43]}.结合潘家口水库流域已有的非点源污染负荷研究, 汛期的TN和TP污染负荷输出占到全年的75%以上^[42], 相应地在汛期污染输出风险也就越高, 这与本文结论中7月和8月污染风险较高一致.

非点源污染受到诸多因素影响, 其防控是牵扯到流域综合治理的系统工程, 需要根据实际情况, 在污染源和污染迁移方面同时进行防控.结合本文和已有潘家口水库流域非点源污染研究^{[37, 38, 44]}, 提出以下非点源污染防控建议: ①潘家口水库流域上游非点源污染以农业面源污染为主, 污染物来自于农药和化肥流失, 所以合理科学的农业管理是非点源污染防治关键.其中, 农业化肥流失主要受施肥量、时间和方式等影响, 农作物对化肥的吸收能力一般情况下保持不变, 所以合理的施肥量可以一定程度上减少污染物的产生量; 研究区年内降水集中在6~9月, 雨季前的施肥在汛期流失严重, 根据已有研究^[45], 若将6月底施肥改在月初施用, 污染物总量将会下降10%; 而深耕施肥相较于播撒施肥方式, 也能有效降低化肥等污染物的流失.②合理景观布局, 可以一定程度对污染物起到截留作用, 在与土地利用规划不冲突的情况下, 增加林地和草地等汇景观面积, 可以有效降低流域非点源污染^[9].有研究表明, 河岸的植被缓冲带可以有效地降低各种污染物进入水体, 达到非点源污染防治目的^[46].所以针对流域非点源污染防治, 在允许的情况下, 合理布局“源”景观和“汇”景观, 可提高景观单元内对污染物的消纳和降解能力.考虑到研究区河道周围分布着大量农田, 应避免在河道附近种植化肥施用量较大的农作物, 并结合付婧等^[47]对国内外学者关于缓冲带最佳宽度的研究, 推荐设置宽度为5 m的植被缓冲带, 减少近河农田污染物进入水体.

(1) 潘家口水库流域地形复杂, 降水时空差异显著, 传统输出风险模型适用性较差, 引入降水因子和地形因子对输出风险模型进行改进.研究结果表明, 考虑降水和地形的输出风险模型对潘家口水库流域非点源污染输出风险评估更合理.

(2) 潘家口水库流域非点源污染输出风险时空差异显著.结合土地利用空间分布分析, 以耕地为主的上游和以建设用地为主的下游, 非点源污染风险较高; 而以林地为主的流域中游, 非点源污染风险相对较低.从时间上分析, 4~9月非点源污染风险先增后减, 在汛期达到最高, 与潘家口水库流域年内降水变化趋势相符.

(3) 分析潘家口水库流域非点源污染输出风险时空分布特征, 可以发现, 流域农业非点源污染贡献较高, 且在汛期输出风险更高, 所以应制定科学的农业管理措施, 减少非点源污染物产量.尽可能减少在河道附近种植高污染输出农作物, 并设置植被缓冲带.