北沙河上游流域潜在非点源污染风险时空变化分析

2020-05-02时迪迪张守红

环境科学研究 2020年4期

时迪迪，张守红*，王红

1.北京林业大学水土保持学院,北京 100083 2.北京市昌平区水务局,北京 102200

点源、非点源污染及城市雨水径流合流制溢流污染等导致水环境恶化[1]，严重影响生态环境和区域可持续发展.近年来，我国的点源污染逐步得到有效控制，而非点源污染随着雨水径流等流入水体中，成为水环境质量下降的主要原因[2-4].非点源污染受降雨、地形、土壤性质、水文过程及人类活动等因素影响，时空变化规律复杂.土地利用变化是人类活动的集中体现，使下垫面的土壤成分、水质特征、水文过程等发生变化，进而对非点源污染产生的影响重大[5-7].所以，土地利用的合理布局是解决非点源污染的重要手段[8].

由于驱动力的可变性和人类活动的异质性，非点源污染具有随机性、广泛性、滞后性和潜伏性等特征，导致其治理难度相对较大[9-11].非点源污染治理的关键在于污染风险的识别和空间解析，以及风险分区和分区治理[12-13].已有研究中常采用模型模拟法[14]、径流分割法[15]、水质水量相关法[16]等方法进行非点源污染风险识别和空间解析，其中模型模拟法的应用最为广泛[17].

非点源污染风险的识别和空间解析模型可分为机理模型和经验模型两类[18-19].机理模型也称物理模型，根据污染物的产生和运移转化机理构建数学模型，模拟污染物的输入、迁移转化和输出，代表性模型主要有SWAT模型[20-21]、AGNPS模型[22-23]、ANSWERS模型[23-24]等.机理模型操作复杂，同时需要输入大量空间数据，在缺乏监测数据地区的应用性相对较差[25].经验模型也称黑箱模型[26]，即综合考虑可能对污染物产生潜在影响的因素[27]，通过统计或因果分析建立经验关系模拟流域非点源污染负荷[28-29].代表性的经验模型有输出系数模型[30]和农业非点源污染指数(APPI)模型[31-32]等.经验模型未充分考虑污染物的迁移转化复杂过程，对数据精度要求低，易操作，可应用于数据资料缺乏地区[29].

意大利学者Munafo等[33]提出的PNPI (potential non-point pollution index,潜在非点源污染指数)模型通过简化污染物迁移过程模拟分析区域非点源污染风险.该模型不需要高度详细的输入数据，且计算较为简便，近年来已逐渐被广泛应用于非点源污染研究中.如Cecchi等[34]应用PNPI模型评估了意大利维泰博省非点源污染对地表水生态环境的压力，并解析了污染风险的空间分布；李婷等[35]将该模型应用于北京市密云区曹家路小流域非点源污染风险分区；李婷等[36]对该模型进行了土地利用类型扩展、土地利用分类比对和土壤渗透性分级3个方面改进，得到潜在污染指数(potential pollution index,PPI)，并在此基础上建立了河段潜在污染指数法.然而，上述研究在土地利用指标、径流指标和距离指标三者的权重赋值方面多依赖专家评价，决策结果具有较大的主观随意性.客观赋权法根据各指标在评价单位中实际数据的离散程度形成最终结果，客观性较强，主要包括主成分分析法、离差最大化方法、离差决策法及均方差决策法等.主成分分析法和离差最大化方法计算繁琐，不便于推广使用；离差决策法反映随机变量的离散程度不如均方差准确，且对属性值的无量纲化有一定的要求;均方差决策法作为一种重要的客观赋权方法，克服了上述方法的不足，该方法不依赖于主观判断，对属性值的无量纲化无限制要求，适用于任何场合，且具有概念清楚和计算简便等特点，在各领域得到了广泛应用[37].如李博翀等[38]通过构建城市地价质量指标体系，运用该方法和协调度评价模型等，对河北省县域尺度城市地价质量进行测度和空间差异分析；邢璐平等[39]通过构建耕地集约利用水平评价指标体系，采用该方法确定指标权重，对江苏省2000—2015年耕地集约利用水平的时空变化规律进行探索分析；田俊峰等[40]运用该方法和耦合协调度模型，对东北三省2003—2015年城市土地利用综合效益及耦合协调关系进行分析研究.因此，可望通过采用均方差决策法改进PNPI模型基于专家评价法的不足.

该研究采用均方差决策法替代专家打分法以改进该PNPI模型主观赋权的不足，并以北沙河上游小流域为例，基于1980年、1990年、2000年、2010年和2017年的5期土地利用，以及土壤类型和DEM等基础数据，借助ArcGIS软件和Fortran语言，分析研究区近40年土地利用类型变化对非点源污染风险时空变化的影响，划分非点源污染风险等级分区，提出非点源污染防治措施建议，以期为流域水环境治理和水生态修复提供科学支撑与参考.

1 研究区概况及数据来源

1.1 研究区概况

北沙河位于北京市昌平区西北部，由虎峪沟、关沟、狡猊沟、兴隆沟、白洋城沟、柏峪沟、高崖口沟汇合而成，流经马池口镇、流村镇、南口地区等人口密集区，是北运河上游温榆河水系的重要组成之一[41].近年来，北沙河流域内经济快速发展，人口迅速增长，对水生态环境造成极大压力.北沙河流域是北运河的重要源头区之一，其水环境污染风险对下游北京城区和通州副中心的水环境和水生态安全具有重要影响.北沙河上游流域总面积114.4 km2(见图1)，地形和土壤类型分布见图2.研究区地形复杂，海拔范围在30～1 403 m之间，整体上呈西北向东南倾斜，由山区逐渐过渡到平原.土壤类型以褐土为主，占研究区总面积的90%以上，潮土、水稻土主要分布在东南角，另有少量区域为粗骨土.

1.2 数据来源与数据处理

该研究基础数据包括研究区边界图、行政区划图、DEM、土地利用类型图、土壤类型图等，数据来源及说明见表1.鉴于栅格和矢量数据来源不同，采用几何校正将所有数据统一为WGS84坐标系，UTM投影，以便后续分析计算.

图1 研究区地理位置与水系Fig.1 Location and streams of the study area

图2 研究区高程和土壤类型Fig.2 DEM and soil types of the study area

表1 基础数据及其来源Table 1 Sources of basic data

2 研究方法

该研究利用改进的PNPI模型客观评估流域内不同土地单元对受纳水体造成的潜在非点源污染风险，评估流程见图3.借助GIS和Fortran语言，将研究区划分为30 m×30 m的网格单元.每个单元格的潜在污染风险可表达为土地利用指标(land cover index,LCI)、径流指标(run off index,ROI)、距离指标(distance index,DI)三者的函数[33-34]，对以上3个指标归一化后进行加权计算，得到每个单元格的PNPI值(PNPI=a×LCI+b×ROI+c×DI，其中a、b、c分别为土地利用指标、径流指标和距离指标的权重).

图3 PNPI的计算路线Fig.3 Flowchart of the PNPI model

2.1 土地利用指标计算

土地利用指标表示不同土地利用类型对受纳水体造成的潜在污染风险.该指标由专家们根据自己的专业知识对不同土地利用类型可能产生的污染量进行判断打分(0～10分)得出.参与打分的专家人数越多，所涉及的研究领域越广泛，土地利用标准差越小，专家打分法的结果越可信.为了使结果更加科学全面，PNPI模型中土地利用指标的赋值咨询了包括生物学家、自然科学家、生态学家和环境工程师等具有不同专业背景的专家，获得了土地利用指标平均值及其标准差[33-34](见表2).

表2 各土地利用指标平均值及其标准差[33-34]Table 2 Mean and standard deviation of LCI values[33-34]

2.2 径流指标计算

径流指标表示地形、土地利用类型和土壤类型等对污染物流动性产生的阻碍和过滤等作用，反映径流流速和流量对污染物迁移扩散的影响.每个网格单元的径流指标值为沿着该网格单元到河网水系之间的水力路径上每个单元径流系数的平均值[33-34].每个网格单元的径流系数由该单元的土壤类型、土地利用类型和坡度决定.该研究参照Cecchi等[34]提供的径流系数(见表3)给每个网格单元的初步径流系数赋值，并进行坡度校正.其中，土壤渗透性分级参考车振海[42]的四级渗透分级法(见表4)，坡度校正参考表5.

表3 不同土壤类型和土地利用类型的径流系数[34]Table 3 Runoff coefficients for different types of soil and land use[34]

注：A、B、C、D表示土壤渗透等级.

表4 土壤渗透分级[42]Table 4 Classification of soil permeability[42]

表5 坡度修正系数[33]Table 5 Slope correction coefficients[33]

2.3 距离指标计算

距离指标表示土地单元到受纳水体的水力距离，即各单元格到达地表水体沿水流路径的长度.距离越长，污染物到达受纳水体的过程中流失越多，排入水体时造成的污染风险越少.距离指标(DIi)计算公式为

DIi=exp(-[Di]×k)

(1)

式中：[Di]为第i个单元格到河网的水力距离，用单元格数表示；k为常量，取值为 0.090 533[30].

2.4 均方差决策法

(2)

(3)

(4)

(5)

式中，xij为指标集Gj中的指标值，Zij为标准化后的指标值，xmax和xmin分别为Gj指标集中的最大值和最小值，E(Gj)为第j项评价指标的平均值，σ(Gj)为第j项评价指标的均方差，W(Gj)为第j项评价指标的权重.

3 结果与分析

3.1 土地利用类型时空变化特征解析

北沙河上游流域的西北山区主要以林地为主，东南平原区主要为旱地、园地、建设用地、农村居民地、城镇用地(见图4).由表6可见：研究区土地利用类型变化最大的是城镇用地和旱地，城镇用地面积在1980—2017年增加了18.5%，旱地面积减少了15.9%；园地面积先增后减，在2000年达到最大值，为17.58 km2，到2017年降至3.88 km2；其余土地利用类型面积变化不明显.

表6 1980—2017年研究区土地利用类型面积及占比
Table 6 Areas and percentage of different land use types from 1980 to 2017

3.2 潜在非点源污染风险空间分布特征解析

3.2.1土地利用指标空间分布

以2017年土地利用数据为基础，讨论潜在非点源污染风险空间分布特征.研究区土地利用指标值差异较大，最低值为0.140(水库坑塘)，最高值为8.220(城镇用地)，整体呈现东南高、西北低的空间特征(见图5).东南平原区村庄城镇集中、人口建筑密集，以种植和养殖等污染排放风险较严重的生产方式为主，对周围水体造成的污染风险较大，因此土地利用指标值较高.而西北山区地势陡峭、人口稀少，土地利用类型以林地为主，对周围水体造成的污染风险较小.

图5 研究区土地利用指标值空间分布Fig.5 Spatial distribution of the land use indicator value

3.2.2径流指标空间分布

径流指标的空间分布特征主要受土壤渗透性、土地利用类型和坡度等因素影响(见图6).研究区土壤类型包括褐土、潮土、粗骨土和水稻土，这4种土壤的渗透系数在0.05～0.40 m/d之间，渗透等级均为C (见表4)，因此研究区内土壤类型对径流指标值的空间分布特征影响较小.各土地利用类型的初步径流系数差异较大，水库坑塘作为污染物的受纳水体，初步径流系数为0，而城镇用地不透水面积比例高，导致污染物在迁移路径中渗漏较小，到达受纳水体的污染量较多，初步径流系数为0.90.研究区地形复杂，尤其是西北部山区，坡度较大，径流汇流速度快，导致迁移过程中渗漏较小，对初步径流系数的修正系数较高(见表5)，因此径流指标值较高区主要分布于西北部山区.地势平坦的东南部平原区，因坡度较小，初步径流系数的修正系数较低，径流指数值较低.

图6 研究区径流指标值空间分布Fig.6 Spatial distribution of the runoff indicator value

研究区径流指数值较高区集中在西北部山区，较低值集中在东南部平原区.径流指数值与坡度相关性较强，主要因为土壤类型在研究区范围内变化不大，对径流指数值的相对影响较小.另外，虽然土地利用对初步径流系数有较强的相关性，但是因研究区内坡度变化较大(见图2)，对初步径流系数的修正影响强烈(见表5)，导致土地利用类型与径流指标值的相关性被弱化.

3.2.3距离指标空间分布

距离指标表示污染源迁移路径上的降解，与到达受纳水体的水力距离呈负相关.距离指标值越高，距离河网越近，表示迁移路径上的降解量越小，污染风险越高.图7中距离指标值与河网及其他受纳水体的空间分布表现出明显的相关性，距离指标高值区主要分布于受纳水体及其附近，较低值分布于距离水体较远的区域.

图7 研究区距离指标值空间分布Fig.7 Spatial distribution of the distance indicator value

3.2.4潜在非点源污染风险空间分布特征

通过对土地利用指标、径流指标和距离指标进行加权计算，得到研究区潜在非点源污染风险空间分布图.利用自然断点法对潜在非点源污染风险进行等级划分，可依次得到风险极低、低、中等、高、极高区(见图8)和各风险等级区的面积占比(见表7).研究区潜在非点源污染风险极低、低、中等、高和极高区的面积占比分别为31.8%、20.1%、17.3%、20.0%和10.9%.其中风险极低区占比最大，主要土地利用类型为林地；风险极高区占比最小，主要土地利用类型为旱地、农村居民地、城镇用地.由于均方差决策法确定的土地利用指标、径流指标和距离指标的权重分别为48.5%、19.1%和32.4%(见表8)，说明土地利用指标对非点源污染风险相对贡献率最高.

图8 基于均方差决策法潜在非点源污染风险空间分布特征Fig.8 Spatial distribution characteristics of potential non-point source pollution risks based on the mean-square deviation decision method

表7 研究区潜在非点源污染风险等级及其面积占比
Table 7 The classification of potential non-piont pollution risks

风险等级潜在非点源污染风险值范围面积∕km2占比∕%极低0.190～0.30636.3931.8低0.306～0.46322.9520.1中等0.463～0.62519.7417.2高0.625～0.77022.8220.0极高0.770～0.97212.5110.9

研究区非点源污染风险分布呈现出东南高、西北低的空间特征，且高风险区集中分布在平原区河网两侧.主要原因是：①东南部为平原区且靠近河网下游，人口建筑等集中分布于河网两侧，不透水面积占比高；②东南部主要为种植业和养殖业，产污高且降解少；③西北部主要为山区林地，人类活动少，产生的污染风险较少，且大部分区域远离水体，污染风险相对较低.

3.3 潜在非点源污染风险时间变化特征解析

3.3.1指标权重变化

经均方差决策法计算，各年份土地利用指标、径流指标、距离指标的权重及其平均值见表8.由于各土地单元距受纳水体的远近、坡度及土壤类型均不变，土地利用类型变化是导致各指标权重变化的主要原因.1908—2017年，土地利用指标权重呈微弱减少趋势(见表8)，主要原因在于土地利用类型呈现均一化的趋势，尤其是2000年以后研究区城镇化的快速推进，城镇用地面积增加，导致研究区土地利用数值的离散程度减弱.这表明相较于专家评价法，均方差决策法可以体现指标权重随土地利用类型变化而变化的优势.

表8 1980—2017年土地利用指标、径流指标、距离指标的权重Table 8 Weights of land use indicator,runoff indicator and distance indicator from 1980 to 2017

3.3.2潜在非点源污染风险时间变化特征

图9展示了研究区在1980—2017年潜在非点源污染风险时间变化特征.对不同时期极高风险区统计得到各土地利用类型占极高风险区面积比例的变化特征如表9所示.由表9可见：1980年和1990年，极高风险区的主要土地利用类型为旱地(66.5%～66.9%)、园地(20.9%～21.1%)和农村居民地(11.6%～12.1%)；随着城镇化发展，2000年旱地面积占比由31.9%减至19.2%，园地面积占比由8.4%增至15.4%，农村居民地面积占比由6.6%增至8.1%，建设用地面积占比则由0.3%增至3.9%，极高风险区的主要土地利用类型变化为旱地(37.5%)、园地(30.8%)、农村居民地(16.1%)和建设用地(14.4%)；随着城镇进一步扩展，2010年城镇用地面积大量增加，极高风险区的主要土地利用类型转变为城镇用地(56.2%)、旱地(19.1%)、农村居民地(14.4%)和建设用地(6.4%)；2017年城镇化进程达到顶峰，极高风险区的主要土地利用类型为城镇用地(64.2%)、旱地(13.7%)、农村居民地(9.5%)和建设用地(9.0%).

图9 1980—2017年基于均方差决策法潜在非点源污染风险时间变化特征Fig.9 Temporal change characteristics of potential non-point source pollution risks based on the mean-square deviation decision method from 1980 to 2017

表9 1980—2017年潜在非点源污染极高风险区中各土地利用类型面积占比
Table 9 The area proportion of land use types in the extremely high risk areas from 1980 to 2017

4 讨论

该研究采用均方差决策法计算得出1980—2017年土地利用指标、径流指标、距离指标的权重平均值分别为0.49、0.18、0.33(见表8).与Cecchi等[34-35]通过专家评价法得出的土地利用指标、径流指标、距离指标的权重值(依次分别为0.48、0.26、0.26)相比，均方差决策法确定的径流指标权重值相对较小，距离指标权重值较大.其主要原因可能在于，研究区内的土壤类型对径流指标值的空间差异影响较小.坡度和土地利用类型的空间差异虽然较明显，对初步渗透系数的空间差异影响也较大，但因坡度变化较大区域集中于西北部山区，而土地利用类型变化较剧烈区域集中于东南部平原区，坡度修正过程对土地利用类型决定的初步渗透系数进行一定程度的均一化，导致径流系数的离散程度相对较低.因此，均方差决策法所得结果中研究区西北山区河网部分为中等风险区(见图9)，与专家评价法所得结果中的极低风险区(见图10)明显不同.

PNPI模型通过简化污染物迁移过程模拟分析区域非点源污染风险[29]，因此，该模型不能详细模拟和区分具体污染物(如磷、氮、重金属等)的迁移转化过程.但该模型计算简便，不需要高度详细的输入数据，且能给出符合实际的流域潜在非点源污染风险分布特征，故而被广泛应用于评估地表水生态环境的压力和解析污染风险的空间分布[34-36]等方面.水污染源解析法主要有输出系数法、源特征因子法、源解析方法等，其中一些方法对基础数据的要求较高，而我国相关基础数据较为缺乏或者不具备系统性和完整性，监测网络建设滞后，一些主要污染物的源成分谱数据库还未能有效的建立.输出系数法需要数据资料最少，但要开展大量的现场观测验证或修正工作输出系数[44]，点面源污染切割较难，模型的适用性验证有待开展.

图10 1980—2017年基于专家评价法潜在非点源污染风险时间变化特征Fig.10 Temporal change characteristics of potential non-point source pollution risks based on expert evaluation method from 1980 to 2017

PNPI模型计算过程(见图3)中，土地利用指标、径流指标和距离指标的计算并非完全独立，如土地利用指标和径流指标都受土地利用类型影响.然而，二者考虑的是土地利用的不同特征，土地利用指标主要由人类活动(如耕种区的农药和化肥应用)决定，体现的是污染的产生；而径流指标主要与不同土地利用类型的下垫面(某土地单元将降雨转化为径流的效率以及从该单元到排水网络的路径)物理特征有关，体现的是污染的运移.由于土地利用指标主要由专家们根据专业知识对不同土地利用类型可能产生的污染量进行判断打分得出，依然存在主观判断误差.非点源污染产生主要包括人畜粪便排放、化肥农药的使用、大气沉降及生物固氮等方式[45-46]，可以通过计算这些污染源的产生量化土地利用指标.

由北沙河上游潜在非点源污染风险空间分布可见，风险较高区域主要集中于东南部平原区.东南平原区集中了大量的畜禽养殖污染、农林业化肥农药污染、城乡居民点和道路地表径流污染等，应为非点源污染控制优先区或流域综合治理的重点治理区域.畜禽养殖污染风险可以通过调整养殖规模和结构、推行生态养殖理念[47]、粪污工程化处理和粪污水资源化利用[48]来削减.针对农林种植业产生的化肥农药污染，应采用优化种植结构、指导农地和果园科学施肥施药、生物肥药替代传统肥药、实施绿色防控等技术措施降低农业和种植业面源污染风险[49].而对于城乡居民点和道路地表径流污染，可采用植被过滤带[50]、生物滞留池[51]、下凹式绿地[52]、透水铺装[53]、植草沟和人工湿地[51]等低影响开发措施降低其污染风险.