平原河网区降雨径流污染负荷测算——以太湖流域望虞河西岸为例*

庞琰瑾，袁增伟

(南京大学环境学院污染控制与资源化研究国家重点实验室，南京 210023)

伴随我国城镇化和工业化的快速推进，不透水地面和大气污染物沉降量快速增加，这不仅加剧了城市内涝[1]，而且造成降雨径流污染负荷激增. 有学者指出，在城市化程度较高的地区，若不重视面源污染控制尤其是降雨径流污染，即使点源污染得到全面控制，受纳水体的水质仍无法得到明显改善[2-3]. 在“河长制”先行的区域，污染治理全面推行精准化，逐渐实现从粗放治理向精准治污转变. 而精细量化降雨径流污染负荷是流域尺度实现面源精准治污全过程控制的重要前提.

传统的降雨径流污染负荷测算往往采用输出系数法，输出系数指单位时间单位面积负荷量，常采用多年平均值，这种方法虽然简便，但未能充分考虑污染负荷产生的时空异质性[4-5]. 随着对降雨径流污染研究的精度要求不断提高，各水文模型被广泛应用于降雨径流污染研究中[6-8]，如SWMM、SWAT、AnnAGNPS等模型. 但这些模型涉及参数过多、资料不易获取且计算过程复杂[9-10].

SCS-CN(soil conservation service-curve number)是被广泛应用于降雨径流污染研究的经验水文模型之一[11-13]. 该方法基于土壤的前期湿度条件(AMC)、土壤类型和土地利用类型特征综合确定了CN取值(变化在0～100之间)，并通过经验公式进一步推算出土壤的最大可能入渗量(S). 最后，基于降雨量(P)、初损值(Ia)和最大可能入渗量(S)之间的数学关系，即可计算得到降雨径流深(Q). 一般情况下，根据降雨事件前5天的累积降雨量P5，将AMC划分为3种等级，即干旱条件(AMCⅠ)CN1、一般条件(AMCⅡ)CN2和湿润条件(AMCⅢ)CN3. 该模型具有计算过程简单、所需参数少、操作简便、资料易于获取等特点，适用于资料缺乏且下垫面复杂的地区[14].

然而，早期的SCS-CN方法没有明确区分前期湿度条件且AMC分级方式会造成初损值突变. 针对这些问题，学者们相继提出了AMC的定量化修正方法，如Mishra等对前期土壤湿度做了定量化修正，构建了Mishra-Singh(MS)模型[15];Geetha等将AMC修正方法与前期湿度含量(AMA)修正方法进行了对比[16]. Li等利用修正后的MS模型测算了沈阳市2015年不同用地类型的降雨径流深[17]. 另外，SCS-CN方法将初损值(Ia)与最大可能入渗量(S)的比值定义为初损率(λ)，通常取标准值0.2，但该取值未考虑其与研究区地理和气候因素的相关性[18]，因此，邓景成等[19-20]分别对我国黄土区和三峡地区的降雨径流初损率做了本地化率定. 总体来看，前期基于SCS-CN模型的降雨径流研究多是为了预测城市雨洪，只关注了地表径流量，并未关注径流造成的水污染物迁移传输，但这些改进和修正模型为提高地表径流量的计算精度提供了有力支持[21]，也为精准测算降雨径流污染负荷提供了方法保障.

望虞河西岸综合示范区位于太湖流域，是正在经历快速城市化的地区，又是“引江济太”工程的清水廊道流经区. 该地区水污染较为严重，其水环境治理对于保障太湖水质具有重要意义. 因此，本研究通过在研究区内开展不同土地利用类型的降雨观测实验，修正SCS-CN模型中的初损率，并基于土地利用类型遥感解译和降雨径流污染物浓度测定，精细刻画望虞河西岸河网区的降雨径流污染负荷时空分布格局. 研究结果为流域面源污染负荷测算和水污染防治方案制定提供了更为详细具体的决策支撑.

1 数据与方法

1.1 研究区概况

研究区(图1)隶属于江苏省无锡市，属亚热带季风气候，降水丰沛，年降水量在1100 mm左右[22]. 研究区占地面积486 km2，位于望虞河西岸南部，以新吴区和锡山区为主，另有约8%面积属于梁溪区. 新吴区是高新产业和现代工业区，锡山区是工农业共同发展的城村复合区，梁溪区为现代商住区.

图1 研究区内主要水系及行政区划

1.2 污染负荷测算

1.2.1 降雨径流深测算 降雨径流深采用SCS-CN模型进行计算[23]，模型假设降雨实际入渗量与其潜在最大入渗量之比等于降雨径流深与潜在最大径流深之比. 其模拟方程为：

(1)

式中，Q为降雨径流深(mm);P为降雨量(mm);S为最大可能入渗量(mm);λ为初损率，是初损值Ia与S的比值，无量纲. 在硬质地表上，由于水分蒸发很快，土壤的前期湿度条件这一因素对径流产生的影响不明显，因此降雨径流深通过公式(1)计算. 然而，对于绿地和农田这两种下垫面类型，本研究采用对前期土壤湿度做了定量化修正的Mishra-Singh(MS)模型[15]，引入M参数，降雨径流深通过公式(3)计算.

(2)

(3)

式中，P5为降雨事件前5 d的累积降雨量，用来表征前期土壤湿度.

1.2.2 污染负荷测算 研究区产生的降雨径流污染，可由“降雨径流深”、“污染物浓度”与“土地利用类型面积”这3个空间化图层数据在ArcGIS软件中用栅格计算器相乘叠加获得，即：

(4)

式中，Wj为降雨径流第j类污染物总负荷;Qi为第i种土地利用类型的降雨径流深，采用本地化修正后的SCS-CN模型计算;Cij为第i种土地利用类型上第j类污染物的浓度，采用间隔采样的监测数据计算获得降雨事件平均浓度EMC;Si为第i种土地利用类型面积，通过遥感解译获得;n为流域内土地利用类型数.

1.3 数据收集与处理

1.3.1 土地利用类型解译 土地利用类型划分参考《土地利用现状分类》国家标准(GB/T 21010-2019)，研究区的土地利用类型分为不透水地表、透水地表两大类，具体包含9种，即商服用地、工矿仓储用地、城镇住宅用地、农村住宅用地、道路、林地、农田、草地、未利用地. 选取拍摄于2017年的Landsat 8卫星影像图和GeoEye-1高分辨率商业遥感影像图作为遥感数据源，用ENVI 5.2软件进行遥感解译. 对监督分类得到的土地利用类型图进行人工目视解译调整，最终得到5 m×5 m高分辨率的土地利用类型图(图2). 解译结果显示，农田在各土地利用类型中占比最大，其次是林地，城镇居民用地、工矿仓储用地、商服用地和道路等硬质地表约占总面积的40%. 从各行政区来看，梁溪区近40%的面积为城镇住宅用地，而锡山区50%以上的面积为农田和林地，新吴区农田、林地、工矿仓储用地分别占该区总面积的20%左右. 考虑到透水和不透水两大类下垫面类型在产流原理和产流特征上有较大差异，同时植被类型和人类活动(如耕作)也会影响产流特征，因此，本研究在研究产流特征时将土地利用类型划分为绿地(包括草地和林地)、农田、硬质地面(包括商服用地、工矿仓储用地、城镇住宅用地、农村住宅用地、道路和未利用地)3类;在研究产污特征时，考虑10种不同的土地利用类型(道路又细分为城市道路和农村道路，两者在遥感解译时难以区分，后续计算时按城镇住宅和农村住宅的比例加权平均处理).

图2 2017年研究区土地利用类型

1.3.2 日降雨量 本研究收集了研究区内10个降雨观测站的日降雨量数据,涵盖了2017年1月1日至2019年12月31日，并使用泰森多边形插值的方法进行空间插值(图3).

图3 降雨径流实验空间点位分布及降雨站的泰森多边形分布

1.3.3 初损率修正 为实现初损率λ的本地化修正，本研究于2018年7月22日至2019年11月30日分别开展了3次天然降雨观测实验和9次模拟降雨实验，测定3种下垫面类型的降雨初损值. 其中，模拟降雨实验采用NLJY-10型便携式人工模拟降雨系统，模拟大、中、小3种强度的降雨(小雨难以产流)，径流场长2.5 m、宽2 m、高3 m，其自带雨量感测装置，可以输出每20 s的降雨强度. 实验监测点位的空间分布如图3所示. 不同下垫面类型的降雨初损值见表1. 通过与手册值进行对比，同时考虑到研究区全年降雨以中雨为主，因此采用模拟降雨实验中雨条件下的初损值进行后续的参数本地化修正. 在此基础上，本研究对3类下垫面的初损率从0开始做步长为0.05的迭代计算，并将初损率和实测的初损值代入SCS-CN模型，计算出相应的降雨径流深. 由降雨径流深除以降雨量计算出径流系数，并与《室外排水设计规范》(GB 50014-2006)给出的参考径流系数、文献中的径流系数[25-26]进行对比，从而推定可靠性较高的初损率取值.

表1 不同下垫面类型的降雨初损值

1.3.4 污染物浓度测定 在降雨径流观测实验中，同时采集径流样品测定其中的污染物浓度. 具体步骤是: 对于各土地利用类型的观测点，自产流起前半小时每间隔5 min采集一次样品，后半小时每15 min采集一次样品. 采样完成后，将各采样点收集到的8～9瓶样品于24 h内送至实验室进行监测分析. COD、总磷(TP)、总氮(TN)、氨氮(NH3-N)4种污染物的浓度监测方法分别采用最新国标方法，即重铬酸钾法(HJ 828-2017)、钼酸铵分光光度法(GB/T 11893-1989)、碱性过硫酸钾消解紫外分光光度法(HJ 535-2009)和纳氏试剂分光光度法(HJ 535-2009). 各土地利用类型上的EMC取值与文献值[25-35]的对比如图4所示.

图4 四种污染物EMC与文献值对比

2 结果与讨论

2.1 地表径流深

研究表明，研究区域绿地和农田的初损率值为0.3，硬质地表的初损率值为0.9. 这一参数修正结果对于城市化程度较高的平原河网区具有一定的参考价值. 在硬质地表上，从降雨向地表径流转化这一过程的主要损失在于产流前的初损部分，因此初损值本身的大小对后续计算影响较大，而硬质地表上产流较快，在获取初损值时容易产生较大误差，未来可以适当增加实验次数，提高参数修正可靠性.

2017-2019年研究区内各降雨观测站在3类下垫面上的径流深计算结果如表2. 从表2可以看出: 1)各站产生的径流深与降雨深之间存在显著的正相关性，即随着各站降雨量的减少而递减，且各观测站径流深的差异程度与降雨深差异基本相同;2)各用地类型上产生的径流深呈现如下特征: 硬质地表>绿地>农田.

表2 各降雨观测站三类下垫面上的径流深

2.2 地表径流量

用径流深与各降雨观测站各类下垫面的面积相乘获得2017-2019年各站在三类下垫面上的径流量(表3). 从表3可以看出: 硬质地表产生的径流量远远大于绿地和农田产生的径流量，这表明城市化增加了不透水地表面积，城市降雨径流和洪涝发生风险也随之增加.

表3 各降雨观测站三类下垫面上的径流量

表4 研究区2017-2019年污染负荷量

2.3 降雨径流污染物负荷

根据上述获得的各用地类型EMC值、地表径流量以及用地面积，核算出研究区2017-2019年径流污染负荷，如表4所示.

从图5中可以看出，单位面积TN、NH3-N、COD负荷的高值均分布于硬质地表，其中单位面积TN和COD负荷最高值在农村住宅用地，原因是农村住宅用地属于硬质地表，雨水下渗量少，产生的地表径流较多，并且农村住宅用地位于农村，周围环境中存在大量有机肥、化肥、农药等物质，导致污染物浓度较高;而单位面积TP负荷最高值在未利用地，原因是未利用地产生的地表径流较多，并且未利用地上的TP浓度偏高(图4)，这说明未利用地是TP污染的重要来源之一，应加强对未利用地的规划利用;此外，草地的四种污染物单位面积负荷均最低，林地次之，在城市建设中可以适当增加草地、林地的面积，或将硬质化的建设用地与草地、林地合理配置，减少污染负荷的产生.

图5 研究区2017-2019年各用地类型的单位面积污染负荷

考虑到污染负荷总量往往与占地面积紧密相关，并不能反映各行政区污染程度的高低，因此本研究采用单位面积污染负荷量来识别重点污染管控区域. 从图6中可以看出，锡山区单位面积TP、TN、NH3-N负荷最高，新吴区次之，梁溪区最低，尤其是单位面积TP负荷，锡山区有大面积农田覆盖，而农田径流的TP浓度更高，导致锡山区单位面积TP负荷更高;新吴区工业较发达，其单位面积COD负荷最高，锡山区次之，梁溪区最低. 由此看出，研究区径流水质管控需重点关注锡山区和新吴区.

图6 各行政区的单位面积污染负荷

2.4 可靠性分析

将本研究测算出的单位面积污染负荷与我国其他区域相比较[25,31-32,35-36]，如表5所示. 可以看出，本研究测算出的各污染物单位面积负荷值均落在已有研究的区间范围内. 不同区域的各类污染物单位面积负荷相差较大，主要是因为各研究区域地理位置、气候条件、土地利用类型、区域发展阶段的差异，相比于辽河流域、潮河流域、白河流域、十五里河流域、抚仙湖流域、海河干流流域，望虞河西岸综合示范区城市化程度更高，硬质化地表面积大.

表5 不同研究区域单位面积径流污染负荷比较

3 结论

1)研究区绿地和农田的建议初损率取值为0.3，硬质地表的建议初损率取值为0.9，修正后的SCS-CN模型能更好地模拟降雨径流深.

2)径流深与降雨深之间存在显著的正相关性，随着各站点降雨量减少，径流深也呈现递减趋势;各用地类型上的径流深表现为: 硬质地表>绿地>农田;研究区径流深最大的地方集中在梁溪区，该区城市化程度最高，分布着大量硬质地表，雨洪管理需重点关注.

3)硬质地表产生的径流量远远大于绿地和农田产生的径流量. 城市化增加了不透水地表面积，城市洪涝发生风险也随之增加. 这一结论对于高度城市化的太湖流域具有普遍性.

4)污染负荷与降雨深之间存在显著的正相关性，研究区降雨径流中TP、TN、NH3-N、COD 4类污染物的负荷量分别从2017年的190、1359、445和16041 t减少到2019年的118、949、314和11250 t. 单位面积TP、TN、NH3-N负荷锡山区最高，新吴区次之，梁溪区最低;单位面积COD负荷新吴区最高，锡山区次之，梁溪区最低. 从用地类型上来看，单位面积TN和COD负荷的最高值在农村住宅用地，单位面积TP负荷最高值在未利用地. 此外，草地的4种污染物单位面积负荷均最低，林地次之.