小流域非点源污染模拟与生态修复影响评价

吴敬东, 杨胜天, 叶芝菡, 胡晓静, 张耀方

(1.北京市水科学技术研究院， 北京100048； 2.北京师范大学 地理学与遥感科学学院， 北京 100875)

非点源污染或面源污染(non-point source pollution)是指在降雨径流(包括灌溉)的淋溶和冲刷作用下，大气、地面和土壤中的污染物进入江河、湖泊、水库和海洋等水体而造成的水污染[1]。目前，大量研究证明，非点源污染是导致水环境恶化的主要原因之一。据美国、日本等国的科学研究报道，即使点源污染得到全面控制，江河的水质达标率也仅为65%，湖泊的水质达标率为42%，海域水质达标率为78%。在我国，云南省滇池、天津市于桥水库、安徽省巢湖、江苏省太湖和辽宁省大伙房水库等水域，其非点源污染比例均超过点源污染，非点源污染已上升为威胁人类饮用水的主要原因[2]。

当前，对密云库上游非点源污染产生、排放和迁移转换机理研究不是很深入，对流域不同区域种植业、养殖业和农村生活污水产生、排放系数及源强等重要参数还不掌握，对河湖水库的水污染变化规律也不是很清楚，对各种非点源污染控制措施的效果效益研究也较少。因此，在点源控制技术相对成熟之后，对非点源污染机理及防控技术和措施的研究，成为当前非点源污染研究的热点。密云水库是北京城市供水的唯一大型地表水源，也是南水北调工程通水后余水存蓄库，保护好密云水库水源对北京城市供水系统的保障、调节和储备安全有重大意义。由于人为活动的影响，密云水库流域的生态环境发生了急剧的变化，入库径流污染物浓度有增加趋势，水库水质安全压力大增。本文通过开展水源地非点源污染研究，摸清非点源污染主要来源、传输途径及污染贡献，旨在为保护密云水库水源地、开展流域山水林田湖一体化修复和生态环境建设提供参考。

1 研究区概况

选择北京密云水库上游蛇鱼川小流域作为研究区。蛇鱼川小流域位于北京市密云水库西北侧，东经116°47′—116°57′，北纬40°35′—40°39′，属于白河水系蛇鱼川河流域，蛇鱼川河直接注入密云水库，小流域地跨密云水库一级(西湾子村)和二级(黄峪口村)水源保护区。蛇鱼川小流域面积25.66 km2。蛇鱼川小流域属燕山山脉，地势北高南低，海拔高度为145～940 m。全流域坡度变化于0°～79°，坡度大于35°的面积占到全流域面积的56%，陡坡在流域内普遍存在。土壤以褐土为主，包括淋溶褐土和普通褐土。土壤质地较粗，多为砂壤土，颗粒松散，黏粒含量低，保水性较差，容易产生土壤侵蚀。流域属暖温带大陆性半湿润半干旱气候，四季分明，冬季干燥寒冷，夏季炎热，雨量集中且多暴雨，春季干旱多风。流域多年平均降水量为652 mm，75%集中在6—9月份。1月份平均气温-6.6 ℃，7月份平均气温24.8 ℃。无霜期为176 d，年日照时数为2 762 h，最大冻土深为85 cm。小流域植被类型有人工的针叶林、落叶阔叶林、灌丛、灌草从、草丛、水生植物以及经济林、农田等植被。

2 研究方法

本研究从流域土地利用类型出发，利用EcoHAT-LCM模型和传统大尺度非点源污染模型相结合的模式，分别计算各流域场次降雨径流过程和非点源污染负荷，利用EcoHAT-LCM模型计算流域水文过程，再利用大尺度非点源污染模型，分别模拟流域多年平均、5 a一遇和10 a一遇暴雨产生的溶解态非点源污染负荷和吸附态非点源污染负荷，设置情景分析，评价小流域生态建设方案对密云水库流域水文过程和污染防控效果影响。

2.1 EcoHAT-LCM模型

EcoHAT-LCM模型是刘昌明等[3-8]在陕西、甘肃、四川、青海、西藏和新疆等地区做了大量坡地降雨实验，通过实验途径探讨了暴雨径流的超渗产流理论，并在这一理论的基础上提出的暴雨径流模型。

产流模型是LCM暴雨径流模型的核心内容，其通过LCM入渗公式计算单位时段内入渗量，并耦合水量平衡方程，以参数优化选择方法率定模型不确定性参数，再求算各个水文分量。对于地表及地表以下的产流过程，LCM模型采用出流系数法计算壤中流和地下径流。

LCM模型中产流模型将整个过程分为地表层、土壤层和地下水层3步进行：首先计算土壤的下渗能力，采用经验公式计算：

f=R·Pr

(1)

R=0.878 1×ln(r)+1.342 2

(2)

式中：f——下渗量/mm；P——降雨量/mm；R，r——经验系数，与土地利用/覆盖和土壤水含量有关，是计算产流过程的关键。

地表径流量的计算基于水量平衡方程，即地表径流量等于降水量减去下渗量：

Qd=P-f=P-R·Pr

(3)

计算壤中流量时，考虑壤中流与土壤湿度和降水入渗量成正比，因此，经验公式为：

Ql=La·(Ws/Wsm)·f

(4)

式中：La——壤中流系数；Ws——非饱和土壤含水量(mm)；Wsm——土壤最大蓄水量(mm)；Ql——壤中流量(mm)。

地下水补给量与基流的计算同样基于水量平衡方程，计算公式为：

REC=Rc·(Ws/Wsm)·(f-Ql)

(5)

Qb=Kb·(GWs+REC)

(6)

式中：REC——地下水补给量(mm)；Rc——地下水补给系数；Kb——基流系数； GWs——地下水储量(mm)。

坡面汇流采用等流时线法。基于DEM数据，进行子流域划分，结合流域等流时线分布图，划分子流域内等流时线空间分布，以等流时线区间为单元，计算流域内坡面汇流过程。子流域等流时线计算公式为：

Tsub_i=Ti-Tmin_i

(7)

式中：Tsub_i——调整后第i个子流域内某点汇流到子流域出口时长；Ti——第i个子流域内某点汇流到全流域出口时间；Tmin_i——第i个子流域内汇流到全流域出口最短时间。

河道汇流采用马斯京根法。基于ArcGIS子流域属性计算结果，结合子流域上下游拓扑关系，应用马斯京根方程对子流域间河道汇流过程进行模拟计算，计算公式为：

Qout,2=C1·Qin,2+C2·Qin,1+C3·Qout,1

(8)

式中：Qin,1，Qin,2——河道时段初和时段末入流量(m3/s)；Qout,1，Qout,2——河道时段初和时段末出流量(m3/s)；C1，C2，C3——调节系数。

2.2 大尺度非点源污染模型

大尺度非点源污染模型是郝芳华、杨胜天、程红光等[9-10]在充分借鉴统计性经验模型和机理性过程模型优势的基础上，结合我国在非点源污染调查工作中的实际情况构建的，将非点源污染负荷产生和运移过程分成溶解态污染负荷和吸附态污染负荷分别进行：在溶解态污染模型中，按城市径流、农村生活、畜禽养殖和农田径流4种类型进行模拟；在吸附态污染模型中，先调用土壤侵蚀模型，估算出研究区内的土壤侵蚀量，然后借助构建的表层土壤氮磷含量数据库，完成水蚀环境下吸附态氮磷的负荷量估算。

2.2.1 模型原理 从影响非点源污染产生要素来看，污染形成主要取决于两个方面[11-15]： ①地表污染物的量级大小及存在形态，即源强要素； ②降水冲刷产流过程的影响程度，即产流要素。前者是流域社会经济特征的综合反映，包括居民的生活方式(如污水、垃圾排放)、城市发展水平(如排水管道铺设、垃圾处理)、农业活动(如化肥、农药的施用)和畜牧业发展情况等；后者则体现出区域的自然环境特征，包括气候条件(如降水频次、雨强、雨量)、土壤类型、地形特征和植被盖度等，这两方面因素是模型构建过程中主要考虑的因素。其次，从污染物产生、输移和转化规律看，非点源污染物有两类：一类为溶解态污染物，另一类为吸附态污染物。溶解态污染物即表明该物质具有水溶性，能随地表径流一起发生迁移，整个过程受水循环控制，主要来自城市径流、农田径流、农村生活和畜禽养殖等人类活动；吸附态污染物即污染物通过附着于颗粒体(泥沙)实现迁移运动，其过程同水土流失密切相关。

2.2.2 溶解态模型 溶解态污染负荷采用二元结构模型，认为影响非点源产生和分布的因子主要是两个方面，即区域的自然环境特征和社会经济活动。这两个方面的影响又分别归纳为土壤类型、土地利用、地形、降雨以及开发强度等。在上述因素中，地形、降雨和土壤结构属于自然因子，而土地利用和开发强度属于社会因子。综合考虑社会、自然二元因子的影响，按农田、农村居民点、畜禽养殖(大牲畜、小牲畜)几种类型分别计算。

(9)

式中：i——非点源污染类型，分别是城市径流、农村居民点、畜禽养殖和农田径流4个类型；n——非点源污染类型数；C——单位面积非点源污染负荷(t/km2)；Qi——单位面积非点源污染源强(t/km2)；ε——径流系数；ε0——标准径流系数，反映不透水硬化地面情况；k——地面冲刷系数；R——标准雨强(mm/d)；t——降雨历时(d)；Ni——自然因子修正系数，表示自然因子修正，如在非点源污染产生过程中需要进一步修正的坡度、植被覆盖和土壤等因子；Si——社会因子修正系数，表示社会发展程度对非点源污染源强的削弱程度。

2.2.3 吸附态模型 土壤侵蚀与非点源污染是一对密不可分的共生现象，在农业非点源污染中表现的尤其明显。从本质上来说，土壤侵蚀产生的泥沙本身就是一种特殊的非点源污染物，会增加水体浊度，破坏水生生物的栖息生存环境，而且泥沙中往往携带了大量的有机物、磷酸盐、铵离子等，从而造成受纳水体污染。非点源污染类型泥沙的产生受到众多外界因素的影响，污染源与受纳水体间的空间距离、地形坡度与土壤糙度、地表植物的滞留作用、积水和低洼区等均会影响进入水体的输沙量。吸附态非点源污染负荷受土壤侵蚀状况的控制，其计算模型为：

Ca=X·Qa·η

(10)

式中：Ca——吸附态非点源污染负荷(t/km2)；X——土壤侵蚀量(t/km2)；Qa——吸附态非点源污染源强(%)；η——吸附态非点源污染富集系数。

2.2.4 土壤侵蚀模型 本研究采用修正的MUSLE方程计算土壤侵蚀量，计算公式为：

Sed=11.8·(RS·qpeak·Apixel)0.56·

KUSLE·CUSLE·PUSLE·LSUSLE·CFRG

(11)

式中：Sed——土壤流失量(t)；RS——地表径流量(mm)；qpeak——洪峰径流(m3/s)；Apixel——栅格单元面积(hm2)；KUSLE——土壤可蚀性因子；CUSLE——植被覆盖和作物管理因子；PUSLE——保持措施因子；LSUSLE——地形因子； CFRG——粗碎屑因子。

洪峰径流计算公式为：

(12)

式中：αtc——汇流时段内日降雨比例，对于日尺度模型取值1；tconc——栅格汇流时间，对于日尺度模型取值24； 3.6——单位转换因子。

土地利用信息提取采用2017年GF-1卫星多光谱数据，数据空间分辨率为10 m；同时应用无人机机载多光谱传感器获取研究区内敏感区域下垫面数据，数据空间分辨率为1 m，再结合地面测量，按照3级分类标准，进行土地利用信息提取，获取整个研究区的土地利用与土地覆盖信息。

3 结果与分析

3.1 小流域土地利用特征

根据土地利用特征分析结果进行蛇鱼川小流域土地利用类型专题图制作，如附图3和表1所示，小流域内林地面积最大，占流域面积的93%，其中天然林地2 056.49 hm2，占流域面积的80%，以板栗树和核桃树为主的经济林(属其他林地)占13%。林地面积与旱地和草地面积合计，流域植被覆盖面积达到2 451.07 hm2，占全流域的96%。

表1 蛇鱼川小流域土地利用类型面积统计

3.2 小流域EcoHAT-LCM模型计算

应用EcoHat-LCM模型分别模拟蛇鱼川小流域多年平均、5 a一遇和10 a一遇的洪水径流过程线。由于流域内无雨量站，无法获取历史雨量资料，选取离蛇鱼川小流域最近的张家坟雨量站来代表当地的雨量状况。

根据《北京水文手册》，获得张家坟雨量站多年平均、5 a一遇、10 a一遇的24 h暴雨量，不同频率场次暴雨量详见表2，同时雨量累积曲线采用《北京水文手册》中确定的北京山区24 h降雨雨型。

表2 蛇鱼川小流域不同频率场次暴雨量

基于ASTER-GDEM数据进行研究区河网提取，子流域划分和等流时线计算[16]，并对各子流域内部的等流时线进行调整，将等流时线划分到各子流域中，实现对子流域内部坡面单元的划分，为坡面汇流过程提供空间数据支持。

由于北京地区近年来降水持续偏少，气候偏旱，因此模型土壤入渗系数选择主要参考张亦弛[17]的黄河多沙粗沙区LCM模型分布式构建与应用硕士论文结果来分析蛇鱼川小流域不同频次设计暴雨的洪峰流量，入渗参数选择如表3所示。

表3 偏旱状况下的土壤入渗系数查找表

通过对设计暴雨数据进行IDW(inverse distance weight)空间插值处理，获取小时尺度面雨量数据，结合DEM流域汇流属性数据[18]，分布式构建确定LCM模型对降雨空间异质性的响应机制，基于土地利用、土壤质地等栅格数据对LCM模型下渗系数实现空间离散化处理，最后应用分布式LCM模型对蛇鱼川小流域多年平均、5 a和10 a这3种设计频率暴雨分别进行径流模拟。偏旱状况下蛇鱼川小流域在多年平均、5 a一遇和10 a一遇场次暴雨下的洪峰流量，分别为34.07，76.45和131.80 m3/s，与吴敬东[19]用WMS模型的模拟结果近似，分别为39.82，81.43和132.89 m3/s。

3.3 小流域非点源污染负荷计算

应用大尺度非点源污染模型结合LCM模型的模拟结果，分别模拟蛇鱼川小流域多年平均、5 a一遇和10 a一遇暴雨产生的溶解态非点源污染负荷和吸附态非点源污染负荷。首先模拟得到流域土壤侵蚀状况，不同频率场次暴雨下研究区土壤侵蚀的空间分布如图1所示。多年平均暴雨下，蛇鱼川小流域大部分坡面区域土壤侵蚀属微度侵蚀，轻度侵蚀主要分布在山区旱地。5 a一遇暴雨下，蛇鱼川小流域大部分坡面区域土壤侵蚀仍属微度侵蚀，但在北部和西部部分山区，出现了轻度侵蚀。10 a一遇暴雨下，蛇鱼川小流域大部分坡面区域以轻度侵蚀为主，其次是微度侵蚀，山区旱地为中度侵蚀或强度侵蚀。

对研究区在不同频率场次暴雨下的土壤侵蚀量进行统计(详见表4)。多年平均暴雨下全流域平均侵蚀模数为5.99 t/km2，微度侵蚀占全流域面积的98.3%，轻度侵蚀占1.7%；5 a一遇暴雨下全流域平均侵蚀模数为13.86 t/km2，微度侵蚀占80.7%，轻度侵蚀占19.3%，中度侵蚀占0.1%；10 a一遇暴雨下全流域平均侵蚀模数为24.96 t/km2，微度侵蚀占42%，轻度侵蚀占57.5%，中度侵蚀占0.4%，强度侵蚀占0.1%。

总体上，蛇鱼川小流域由于植被覆盖率高，占全流域的96%，土壤侵蚀以微度侵蚀为主，但是随着暴雨量的增加，侵蚀级别向轻度侵蚀方向移动，当暴雨量为10 a一遇时，将出现强度侵蚀。不同频率场次暴雨产生的溶解态非点源污染负荷和吸附态非点源污染负荷的空间分布分别如图2—4所示。流域内溶解态N，P负荷最高的为经济林地和居民点，其次为旱地，这3种土地利用类型受人类活动影响，施用大量化肥或产生大量垃圾、生活污水，导致溶解态非点源N，P负荷高。吸附态N，P负荷与土壤流失量密切相关，由于流域内土壤侵蚀主要发生在林地，通过林地水土流失带走的吸附态N，P负荷相对较高。

表4 不同频率场次暴雨下蛇鱼川小流域土壤侵蚀统计情况

图2 多年平均暴雨下产生的溶解态非点源污染负荷

对研究区在不同频率场次暴雨下产生的非点源N，P污染负荷进行统计，并根据等流时线计算不同子流域污染负荷进入河道的总量(详见表5)。蛇鱼川小流域由于植被覆盖率高，土壤侵蚀程度较轻，大部分坡面区域土壤侵蚀属轻微侵蚀，导致流域总体吸附态污染负荷较其他区域偏小。以板栗树为主的经济林，对氮肥需求量大，导致流域溶解态N负荷较其他区域偏高。此外，不同频次暴雨对吸附态N，P负荷的影响大于对溶解态N，P负荷的影响，如表6所示，随着场次雨量的增加，吸附态N、P负荷占污染负荷总量的比重在逐渐增加，相应地溶解态N，P负荷占污染负荷总量的比重在逐渐减少。

图3 5 a一遇暴雨下产生的溶解态非点源污染负荷

图4 10 a一遇暴雨下产生的溶解态非点源污染负荷

表5 不同频率场次暴雨下产生的非点源N，P污染负荷统计情况

表6不同频率场次暴雨下溶解态和吸附态N、P污染负荷所占比例%

暴雨频率溶解态N吸附态N溶解态P吸附态P多年平均435736645 a一遇3367287210 a一遇26742179

3.4 结果验证

夏立忠等[20]对密云水库流域的氮素流失进行了监测，密云水库流域氮素流失通量为12.7 kg/hm2。吴敬东[19]对水土流失、农业生产、畜禽饲养和村庄生活等污染源负荷的计算，得到蛇鱼川小流域年总氮流失量为91.21 t, 总磷31.28 t，折全流域平均氮素流失通量17 kg/hm2，磷素流失通量5.20 kg/hm2，农业生产的溶解态氮流失通量6.20 kg/hm2，生活污水的溶解态氮流失通量和溶解态磷流失通量分别为25.49和2.61 kg/hm2。

由于缺乏场次暴雨产污负荷的实测数据，统计研究区在不同频率场次暴雨下单位面积产生的非点源N，P污染负荷，按照降水量占多年平均降水量的比重对场次暴雨的产污负荷进行折算，与夏立忠[20]和吴敬东[19]等人的研究结果进行对比，结论较为一致(表7)。

表7 蛇鱼川场次暴雨产生的非点源污染负荷验证

3.5 情景分析

情景设置是情景分析中的核心内容，直接决定情景分析的科学性以及由此提出建议、措施的合理性等。根据非点源污染负荷估算结果可知，农村居民点和农业生产活动为蛇鱼川小流域非点源污染的主要来源。在此基础上，本研究以研究区社会、经济现状为基础，分别针对农村居民点和农业生产活动进行情景设置。其中基准情景(S0)采用研究区在现状条件下3种频率场次暴雨产生的非点源污染负荷的模拟结果。

(1) 情景1。农村居民点设置0.2的垃圾处理率和0.1的垃圾入网率，其他条件不变下，设置农村居民点具有小城镇级别的垃圾处理率0.2和垃圾入网率0.1，通过模拟计算，农村生活的N，P污染负荷能够削减28%(表3—7)。

(2) 情景2。经济林施肥量削减为原施肥水平的50%，其他条件不变下，将经济林的施肥量变为原施肥水平的1/2，通过模拟计算，农业生产活动的N，P污染负荷能够削减43.9%和48.9%(表8)。

表8 不同情景条件下蛇鱼川场次暴雨产生的非点源污染负荷模拟结果

注：3种情景设置的TN，TP负荷单位为t。

4 结 论

(1) 偏旱状况下蛇鱼川小流域在多年平均、5 a一遇和10 a一遇场次暴雨下的洪峰流量分别为34.07，76.45，131.80 m3/s。

(2) 蛇鱼川小流域植被覆盖率较高，土壤侵蚀以微度侵蚀为主，但是随着暴雨量的增加，侵蚀级别向轻度侵蚀方向移动，当暴雨量为10 a一遇时，将出现强度侵蚀。蛇鱼川流域由于以板栗树为主的经济林，对氮肥需求量大，导致流域溶解态N负荷较高。不同频次暴雨对吸附态N，P负荷的影响大于对溶解态N，P负荷的影响，随着场次雨量的增加，吸附态N，P负荷占污染负荷总量的比重在逐渐增加，相应地溶解态N，P负荷占污染负荷总量的比重在逐渐减少。

(3) 小流域内居民点具有小城镇级别的垃圾处理率和垃圾入网率时，农村生活N，P污染负荷可以削减28%；减少经济林的现有施肥量，蛇鱼川小流域农业生产活动的N，P污染负荷能够削减43.9%和48.9%。