余子贤,钱 瑶,李家兵,李小梅,唐立娜,*
1 福建师范大学环境科学与工程学院,福州 350007 2 中国科学院城市环境研究所,城市环境与健康重点实验室,福建省流域生态重点实验室,厦门 361021 3 数字福建环境监测物联网实验室,福州 350117
随着点源污染逐步得到有效控制,非点源污染导致水环境问题日益突出[1]。非点源污染与点源污染相对应,是指溶解的或固体污染物从非特定的地点,在降水和径流冲刷作用下,通过径流过程而汇入受纳水体(如河流、湖泊、水库、海湾等),引起的水体污染[2]。我国学者对中国五大流域、五大湖泊进行非点源污染的调查发现非点源污染已经成为我国水质恶化、水体富营养化的主要原因[3—6]。城市化进程导致的土地利用变化对小流域的非点源污染有着显著影响。
非点源污染的研究经过几十年的探索与发展,当下的非点源污染研究主要涵盖了非点源污染负荷量核算和风险识别的指标体系定性评价两个方面[7]。污染负荷量化研究广泛结合了数学模型,如通用土壤流失方程(RUSLE)、化学污染物径流负荷和流失模型(CREAM)和小流域尺度评价农业非点源污染的模型(AGNPS)等。历经不断发展,美国农业部于1994年推出的SWAT模型成为了目前国内外应用最广泛的模型之一[8]。采用定量模型的研究虽然可以估计相对准确的污染物流出,但这些污染负荷参数由于其空间代表性有限,难以对区域或流域尺度进行表征且模型本身使用最大的困难是对大量参数的验证,盲目的使用模型会直接影响数值模拟的精度[9—10]。雷能忠等[11]认为使用各种复杂手段获取某些地区精确污染负荷不仅十分困难,也不是污染控制所必须的,为了提高非点源污染的治理成效,识别流域内污染关键源区(CSAs: Critical Source Areas),从而使控制与管理措施更具针对性,已经被公认为是减轻非点源污染危害的关键技术并且得到了广泛应用[12]。目前识别关键的非点源污染区域的方法主要包括磷指数法(PI)[13]、潜在非点源污染指数(PNPI)和潜力指数法(APPI)三种。此外,还有研究运用修正的通用土壤流失方程(MUSLE)对非点源污染关键区进行识别[14]。张丽等[15]采用半定量、经验性的流域尺度磷流失分级方案结合坡度、高程和至河流距离等因素确定了非点源磷输入关键区域。
经济社会高速发展背景下剧烈的城市化过程极大程度的改变了土地利用,土地利用变化可以进一步影响流域的营养物富集程度、水环境容量和非点源污染的发生。在“山水林田湖草是共同生命体”的系统思想指导下,要求我们在生态环境治理中更加注重统筹兼顾。在对非点源污染进行治理时要寻根寻源,找到问题的关键所在。因此,研究土地利用变化同流域非点源污染的关系是明确流域非点源污染程度和来源并进行流域污染物控制的关键。
小流域是城市化进程中受土地利用变化影响强烈的区域,往往是非点源污染发生的重点风险区域。而半城市化的小流域由于受限于土地利用配置的快速转变以及研究区域范围小而导致难以布设全面的水文观测站点来获取庞大的数据量进行复杂模型的参数调控和模拟,产生了无法对半城市化小流域区域进行非点源污染进行详尽模拟的困难[16—17]。陈利顶等[18]最先阐明了异质景观中非点源污染的动态变化,提出以研究“源”“汇”景观空间分布格局来分析非点源污染形成的影响,基于“源-汇”理论提出了景观空间负荷对比指数,量化了景观空间与非点源污染之间的关系[19]。刘芳等[20]探讨了长江上游流域景观空间格局和非点源污染之间的定量关系,研究说明景观空间负荷对比指数对非点源污染负荷响应显著。王金亮等[21]使用“源-汇”理论对小流域非点源污染风险格局进行划分与评价。“源-汇”理论可以作为缺资料半城市化小流域开展非点源污染空间风险评价的有用方法。
后溪小流域地处典型半城市化区域,其近十年来经历着厦门岛外快速城市化进程下带来的非点源污染发生风险。小流域上游是厦门市重要饮用水源地石兜-板头水库,下游是景观水体杏林湾水库[22]并与入海河口相连,可以被视作厦门市经济社会可持续发展的重要的水环境资源。因此,后溪小流域也是应用基于生命共同体理念开展“山水林田湖草沙海城”系统治理的典型区域。
本文基于采用“源-汇”理论构建非点源污染风险评价指数,并选取典型半城市影响下的后溪小流域为研究实例开展非点源污染风险评价和关键源区识别。为“十四五”推动生态文明实现新进步和“山水林田湖草沙海城”系统治理提供可靠依据。
后溪流域位于厦门市集美区,流域总面积192.11 km2,地处北纬24°34′02′′—24°45′48″,东经117°55′14″—118°06′52″,后溪发源于戴云山脉与博平岭山脉交界的老寮仓山西麓,是厦门市第二大河流。地处沿海丘陵地带,地势自西北向东南倾斜,地势起伏较大,呈丘陵和山地、平原的梯状分布(图1)。流域属于南亚热带海洋性季风气候,全年平均气温20.6 ℃左右,多年平均降水量1206 mm,主要集中在4月—10月。主要土壤类型为红壤、赤红壤和黄壤。
图1 研究区地理位置Fig.1 The geographical position of research area
至2020年,流域所在地区城镇化率已达86.99%,2010—2020年流域范围内土地利用变化率为33.43%。流域水质受人为干扰大,总氮(TN)和总磷(TP)存在不同程度的超标[23]。当前,流域正面临着生态退化、生物多样性丧失及非点源污染导致水环境恶化问题。
研究数据需求包括土地利用、地形、水系分布。(1)基于Google Earth遥感影像(分辨率1.19 m),分别生成2010年、2015年和2020年后溪流域土地利用与覆被图,在使用易康(eCognition 9.0)监督分类的预处理基础上配合高分辨率影像在ArcGIS中进行目视解意修正,共将土地利用类型划分为6类,包括林地、草地、耕地、建设用地、园地和水体。(2)基于DEM(分辨率30 m),以及提取的遥感影像数字化水系,运用ArcGIS 10.6软件,基于Agree算法对DEM修正,通过Hydrology模块提取流域边界、河道信息、坡度等地形或水文数据,数据来源于国家测绘地理信息局(http://www.sbsm.gov.cn/)。(3)使用Arcgis 10.6中的渔网工具将后溪流域进行网格化的划分,共将研究区域以30m×30m的网格尺度划分为2307个研究单元,以网格中心点计算距河最近距离。各网格的属性表赋值有土地利用数据、平均坡度、平均高程、距河距离和景观格局指数等数据。
本研究采用网格污染指数评价非点源污染发生风险,基于陈利顶提出“源-汇”理论的指导思想及其发展的景观空间负荷对比指数[19]为参考设计了以网格为研究单元的网格污染指数,分别计算网格单元非点源氮污染和磷污染的网格单元负荷值,最后得到研究区网格单元污染指数,以此为依据对关键源区进行识别。计算步骤如下:
(1)研究区网格单元污染物负荷指数(Grid cell Pollution Load, GPL)
(1)
(2)
GPLNP=GPLN+GPLP
(3)
式中,GPLN、GPLP和GPLNP分别是总氮、总磷、氮磷总体的污染负荷指数;m是“源”景观类型数,n是“汇”景观类型数;WiN是第i类“源”景观排放总氮的权重,WiP第i类“源”景观排放总磷的权重、WjN是第j类“汇”景观吸纳或截留总氮的权重,WjP是第j类“汇”景观吸纳或截留总磷的权重。Pi是第i类“源”景观在网格中的面积比例,Pj是第j类“汇”景观在网格中的面积比例。各类景观类型的面积比例是基于遥感影像使用易康(eCognition 9.0)监督分类并目视解意修正后统计得到的结果。
(2)地理要素数据的标准化处理
由于网格单元非点源污染发生风险不仅受自身氮、磷负荷值的影响,同时还由网格单元的高程(Elevation)、距河距离(Distance)和坡度(Slop)共同决定。故对上述三地理要素进行标准化处理,考虑其产生的影响。
(4)
式中,X指地理要素,包括高程(E),距离(D)和坡度(S)。由于高程和离水系的距离是负向指标,坡度是正向指标,因此坡度计算的分子为X-Xmin。
(3)研究区网格单元的网格污染指数(Grid cell Pollution Index, GPI)
(5)
式中,GPI是在GPL的基础上多考虑了网格单元坡度、高程和离河流水体距离3个地理要素对网格单元非点源污染物排放的影响。
(4)非点源污染发生风险等级划分
统计后溪流域研究区共2307个风险识别单元中“源”、“汇”景观类型所占据该单元的面积比重,根据历年不同土地利用类型的氮、磷污染物输出排放和拦截吸收权重计算系数确定污染物负荷,最后结合地理要素产生的影响计算研究区网格单元的网格污染指数(GPI)来识别研究区非点源污染风险的分布情况。通过计算,若GPI值小于0时,表明该研究单元发生非点源污染风险较低,此区域是以“汇”景观类型起主导作用的区域;GPI值大于0时,表明该研究单元存在一定的非点源污染发生风险,则此区域是以“源”景观类型起主导作用的区域。进一步地将不同GPI值在Arcgis中进行重分类操作,GPI对应的污染风险等级如表1。
表1 非点源污染发生风险等级
不同景观类型对非点源污染的作用大小不同,为了避免专家打分法导致的主观性误差并且更加贴合研究区的实际情况,本文基于全国污染源普查手册、统计年鉴等资料来确定污染物输出权重。
其中,农田污染输出权重按照农业污染普查手册中“地表径流-南方湿润平原区-平地-水田-双季稻”区域的标准农田流失系数作为依据计算得出,总氮流失量(TNL标准)为13.99 kg hm-2a-1,总氮流失量(TPL标准)为1.15 kg hm-2a-1。其中标准农田的总施氮量为360.88 kg/hm2,总施磷量为109.19 kg/hm2,降水在1000—1200 mm之间。根据历年厦门经济特区年鉴中区域实际的总氮施肥量、总磷施肥量和年降水量,确定耕地的氮施肥修正系数(N修正系数)、磷施肥修正系数(P修正系数)和降水修正系数(R修正系数),从而计算流失量(表2)。计算研究区历年耕地总氮流失量(TNL)、总磷(TPL)流失量公式如下:
TNL=TNL标准×N修正系数×R修正系数
(6)
TPL=TPL标准×P修正系数×R修正系数
(7)
表2 耕地氮磷营养物质年流失量
园地氮磷排放量依据农业污染源普查手册中监测类型为地表径流、地形为缓坡区、土地利用方式为旱地、种植模式为园地的标准农田的模式,其总氮流失量(TNL标准)为19.95 kg hm-2a-1,总磷流失量(TPL标准)为1.60 kg hm-2a-1,标准农田中氮磷施肥量分别为418.33 kg/hm2、193.34 kg/hm2,结合历年厦门经济特区年鉴统中区域实际的总氮施肥量、总磷施肥量、降水量。确定园地的氮磷施肥修正系数和降水修正系数从而计算流失量(表3)。计算研究区历年园地总氮流失量(TNL)、总磷(TPL)流失量如下:
TNL=TNL标准×N修正系数×R修正系数
(8)
TPL=TPL标准×P修正系数×R修正系数
(9)
表3 园地氮磷营养物质年流失量
建设用地包括了城镇用地和农村居民点的氮、磷输出系数参考黄金良等[24]在九龙江流域农业非点源氮磷负荷估算的研究结果,确定城镇用地总氮、总磷的流失量为12.09 kg hm-2a-1、1.14 kg hm-2a-1。农村居民点总氮、总磷的流失量为17.74 kg hm-2a-1、2.10 kg hm-2a-1。结合历年厦门经济特区年鉴统中城镇用地和农村居民点占比,确定其修正系数分别为0.75和0.25。计算研究区建设用地总氮(TN)、总磷(TP)流失量如下:
TN=12.09 kg hm-2a-1×0.75+17.74 kg hm-2a-1×0.25=13.50 kg hm-2a-1
TP=1.14 kg hm-2a-1×0.75+2.10 kg hm-2a-1×0.25=1.38 kg hm-2a-1
裸地的总氮(TN)、总磷(TP)流失量参考黄宁等[25]在九龙江流域的研究,确定总氮(RTNR)、总磷(RTPR)流失量分别为8.09 kg hm-2a-1和1.90 kg hm-2a-1。
詹书侠等[26]在中亚热带丘陵红壤区森林演替典型阶段土壤氮磷有效性的研究中发现不同林地利用类型土壤氮、磷截留值分别在0.15—1.43g/kg 、0.19—0.54 g/kg。后溪流域地处南亚热带常绿阔叶林带,根据研究结果确定氮、磷截留系数分别为1.43 g/kg 、0.54 g/kg。王卫霞等[27]对我国南亚热带几种人工生态系统生态系统氮储量研究发现南亚热带常绿阔叶林土壤层氮储存量为16.48 t/hm2,研究区所处赤红壤区的土壤全磷含量为0.28 kg/m2[28]。计算研究区林地总氮(RTNR)、总磷(RTPR)平均吸收能力分别为:
RTNR=1.43 g/kg×16.48 t/hm2=23.57 kg/hm2
RTPR=0.54 g/kg×0.28 kg/m2×10=1.51 kg/hm2
由于不同河流对氮、磷的净化功能与温度、水位、面积等有关[29],同一河流在不同时期差异较大,难以用统一标准衡量污染物去除量。此处参考许芬等[30]研究中依据河流与标准农田中耕地对氮磷排放的相对关系从而确定水体对氮、磷吸收权重分别为0.01和0.03。
综合上述各景观类型的污染物流失量,最终得到不同景观类型对主要污染物氮、磷的污染排放权重,见表4。
表4 各景观类型污染物输出权重
通过对遥感影像解译得到后溪流域研究区2010年、2015年和2020年的土地利用图(图2)。综合历年影像的解译结果来看,在研究区的南部主要以“源”景观集中分布,土地利用类型主要是以建设用地和耕地为主导,而北部地区以“汇”景观集中分布,土地利用类型主要为林地,研究区域的“源”、“汇”景观出现明显景观的异质性,总体上呈现出南“源”北“汇”的分布格局。
图2 后溪流域2010—2020年土地利用Fig.2 Land use in Houxi Basin from 2010 to 2020
图3 历年土地利用类型占比Fig.3 Proportion of land use types in past years
统计了历年土地利用类型面积占比(图3)。结果表明,林地在三年中始终为主导的土地利用类型,历年面积分别为:112.90 km2、106.98 km2、111.92 km2,基本保持不变。建设用地为次主导类型,呈增加趋势,其在2010—2015年间增幅明显,这也是厦门岛外城市化扩张最剧烈时期,截止至2020年共增加了16.06 km2,现有建设用地面积53.19 km2。耕地面积则呈现逐年下降的趋势,十年中耕地面积共减少了13.64 km2,仅余6.68 km2。裸地、水体和园地的面积变化程度不大。
为了明确十年中土地利用互相转移情况,计算了2010—2020年后溪流域土地利用转移矩阵(表5)。其中耕地转出面积占比最大,达71.61%,主要向建设用地和林地发生转移,共转移6.21 km2和7.30 km2。林地转出面积最大,达25.50 km2,其中向建设用地转移10.39 km2。建设用地为转入面积最大的土地利用类型,十年中共有25.50 km2其他土地利用类型转化为建设用地。裸地土地利用变化率最高,有95.65%的裸地是由其他土地利用转入。水体和园地变化程度不大。
对研究区的2307个风险识别单元进行网格污染指数(GPI)值计算来识别非点源污染关键源区,对计算结果重分类后分别得到2010年、2015年和2020年三年的非点源污染发生风险分布图(图4)。
“汇”景观区始终是面积占比最大的区域(图5),现有“汇”景观区面积131.54 km2。历年占比分别为64.89%、64.93%和68.44%,平均GPI值为-0.80,-0.75和-0.74。
关键源区是非点源污染可能发生的区域,现有关键源区面积共60.65 km2,历年占比为35.11%、35.07%和31.56%,平均GPI值为0.25、0.27和0.26;由发生风险大小细分关键源区为以下三个区域:
低风险区是非点源污染发生风险区(GPI>0)的主要构成部分,是非点源污染发生风险较低的区域。现有低风险区面积32.16 km2。历年占比分别为18.55%、16.65%和16.73%,平均GPI值为0.13,0.13和0.14。
中风险区现有面积为23.16 km2。历年占比分别为14.04%,15.35%和12.05%,平均GPI值为0.35,0.35和0.36。
表5 2010—2020年后溪流域土地利用面积转移矩阵/km2
图4 后溪流域历年非点源污染发生风险Fig.4 Historical risk of non-point source pollution in Houxi Basin
图5 后溪流域历年非点源污染发生风险区面积占比 Fig.5 Percentage of non-point source pollution risk area in Houxi Basin in past years
高风险区是“源”景观污染负荷显著大于“汇”景观消纳作用的区域,存在较大的非点源污染发生风险,在非点源污染发生风险区占比最少,现有面积5.33 km2。历年占比分别为2.51%,3.08%和2.78%。平均GPI值为0.58,0.57和0.56。
风险区分布同样呈现出明显的空间分异性,北部以林地为主要土地利用类型的典型“汇”景观构成了对非点源污染具有吸纳截留作用的“汇”景观区域,而南部的非点源污染发生区主要是由耕地、建设用地等非点源污染输出强烈的“源”景观构成。
计算得到了2010—2020年各风险区面积转移矩阵(表6)。十年中,由非点源污染发生风险区(GPI>0)向“汇”景观区域(GPI<0)转移面积为21.97 km2,占研究区面积的11.43%。相反,由“汇”景观区域共转出15.14 km2土地面积变为非点源污染发生风险区,占研究区面积的7.88%。十年中,非点源污染风险区的土地利用转移率均大于50%,高风险区土地利用转移明显,共有99.80%的土地面积发生转移。这表明非点源污染的关键源区尤其是高风险区是随土地利用的变化转移呈现出明显的变化。
表6 2010—2020年后溪流域非点源污染风险等级面积转移矩阵/km2
不同风险区的土地利用类型构成在时间序列上呈现出明显变化(图6)。整体上,林地为“汇”景观区域的主导土地利用类型,非点源污染发生风险区的主导土地利用类型为建设用地。
图6 2010—2020年各风险区景观占比Fig.6 Landscape proportion of each risk area from 2010 to 2020
“汇”景观区在2010—2020年均主要由林地构成,占比分别为80.62%、76.81%和77.40%,其次由建设用地构成,占比分别为7.53%、8.99%和9.78%。较高的林地占比使得该区域具有较强的截留和吸纳非点源污染物的能力,发生非点源污染风险的概率很低。
低风险区在2010—2020年中主要由建设用地构成,占比分别为42.17%、48.06%和53.44%,其次由林地构成,占比分别为27.08%、27.30%和27.71%。低风险区的“源”景观主要以非点源污染物排放输出较弱的建设用地和裸地构成,相对应的由林地和水体构成“汇”景观占据较少的面积,这使得该区域网格单元的GPI值大于0,但又由于该区域“源”景观非点源污染物输出能力有限且受到“汇”景观对非点源污染吸收和截留作用下的一定中和,具有较低程度非点源污染发生风险。
中风险区在2010年主要由建设用地和耕地构成,占比分别为45.19%和30.95%。在2015年和2020年则以建设用地为主导,占比为67.94%和73.03%。中风险区中“源”景观的占比进一步扩大,“汇”景观对污染物吸纳截留作用已经明显小于产生的污染物负荷,具有相对于低风险区更高的非点源污染发生风险。
高风险区在2010年由耕地类型为主导,占比为70.31%。在2015年和2020年以建设用地占据主要类型,占比分别为86.29%和81.70%。典型的“源”景观耕地和大面积的建设用地均产生了较高的污染物负荷,而此区域“汇”景观占据的面积很少,几乎起不到对“源”景观的制约作用。高风险区是典型的非点源污染极可能发生区域且是可能是非点源污染较为严重的区域。
(1)非点源污染发生风险区分布与土地利用分布规律趋同。北部以林地为主要类型“汇”景观区域对应了绝大部分“汇”景观区域。以耕地为典型的“源”景观,对应了非点源污染发生的高风险区。中低风险区分布在中南部地区与建设用地分布相近。这说明景观类型及其自身的排放和吸纳污染物的水平很大程度上决定了所在网格单元发生非点源污染的风险水平。非点源污染发生风险区中“汇”景观区域面积上升,中、低风险区面积减少但其风险值均升高。这是由于“汇”景观占比减少,“源”景观占比增加导致。后溪历年非点源污染发生平均GPI值分别为0.24,0.25和0.27,呈上升趋势。可见,风险值的大小由各景观单元自身的污染排放水平以及所处地理环境共同决定,在2010—2020年间耕地的面积不断减少以及其污染权重的不断下降,这是导致高风险区面积减少及其风险值下降的主要原因。在城市化的持续影响下,建设用地面积不断扩张,大量的其他土地利用类型向建设用地发生转移,这使得建设用地在险区中占比不断上升,进而导致中、低风险区的平均GPI升高。从时空尺度上分析GPI同景观占比的变化,可以发现城市化很大程度上对非点源污染发生风险的关键源区产生了影响。
(2)景观类型的分布、占比及自身排污能力决定了区域非点源污染风险的发生,城市化影响下土地利用结构的复杂性势必要求我们区分出不同情境来优化调整“源-汇”配置结构对小流域非点源污染进行调控。基于“源-汇”理论对非点源污染调控的思路是控制污染物在景观单元内达到收支平衡或在进入水体之前通过拦截达到平衡,实质上就要求我们从源头和过程两个角度出发对源、汇景观的配置进行考虑。
从整体景观格局调控、关键区景观组合方式或景观类型的转换和关键区单一“源”的局部调控三个尺度制定出对应的高、中和低三种强度的调控对策,以此满足流域地区对调控强度不同需求。高强度调控——整体景观格局调控:从宏观尺度上,要让“源”、“汇”景观在离河距离、坡度和高程等地理要素方面适宜性分布,从整体上构建合理的流域“源-汇”景观空间的布局模式;中强度调控——景观组合方式或景观类型的转换:从中尺度上,基于“源”、“汇”景观组合产出非点源污染较少的原理,构建“源-汇”景观合理组合模式即将高污染负荷“源”景观的周边相关景观转为低污染负荷“源”景观类型或“汇”景观类型;低强度调控——单一“源”景观的局部调控:从微观尺度上,在总体不改变“源”景观类型的前提下,以增“汇”减“源”为目标,在“源”景观中增补或镶嵌若干斑块的“汇”景观。
后续开展政策制定可从上述三种调控强度入手,探讨如何合理地调整土地利用配置模式从而达到“源”“汇”景观的收支平衡进而构建较低的非点源污染发生风险区,达到控制非点源污染发生风险的目的。
本文在“源”“汇”景观格局理论指导和考虑研究区历史资料及地理要素的前提下构建非点源污染风险评价指数,并以半城市化的厦门市后溪小流域为例,使用GIS技术进行非点源污染发生风险的关键区域识别,同时对风险区构成、转移情况进行分析。后溪小流域目前“汇”景观区域占流域面积的67.86%,非点源污染发生风险区占32.14%。风险区非点源污染发生风险值处于低风险水平但呈上升趋势。“汇”景观区域集中分布在研究区北部,主要以林地为主要构成。非点源污染发生风险集中在后溪流域南部,建设用地是各风险区的主要构成,并呈上升趋势。研究说明城市化的进程会进一步对风险区结构产生影响,同时对如何进行关键源区调控进行了讨论,能够为城市化影响进程下的中小尺度流域的饮用水源地管理规划提供借鉴。
虽然对厦门市后溪流域进行了非点源污染风险的定性评价,但由于流域内降雨、土壤属性、气候条件相对统一,因此未对此类均质性因素的影响进行考虑。若将此方法应用于复杂环境状况的研究区域,不同地区的降水、土壤类型等也会对非点源污染的分布产生影响,后续研究应进一步探讨除现有因素外多因素对非点源污染风险的影响;文中各污染系数主要参考已有的研究来进行确定,并不一定完全适用于后溪小流域,后续应对研究区各景观类型的典型区域进行详尽的试验来确定各景观单元的详细数值来进一步得到更为准确的结果;其次由于仅对三期的影响进行了解译,仅能对当前的非点源污染风险形势进行估计,提出的调控政策建议也需要后续研究可对未来不同的土地利用情境下非点源污染风险进行模拟,估算政策有效性并择优进行选择,从而使得政策调控能够先行落地,集中资源对可能潜在的非点源污染进行调控。