基于GIS和MUSLE的东北沟小流域非点源污染关键区识别

刘 楠, 谢永生,, 索改弟, 景民晓, 陈 磊

(1.西北农林科技大学 资源环境学院, 陕西 杨凌 712100; 2.中国科学院 水利部 水土保持研究所, 陕西 杨凌 712100)

基于GIS和MUSLE的东北沟小流域非点源污染关键区识别

刘 楠1, 谢永生1,2, 索改弟1, 景民晓1, 陈 磊2

(1.西北农林科技大学 资源环境学院, 陕西 杨凌 712100; 2.中国科学院 水利部 水土保持研究所, 陕西 杨凌 712100)

以河北省平泉县东北沟小流域DEM、土地利用类型和土壤调查数据为基础，在地理信息系统(GIS)的支持下，运用修正的通用土壤流失方程(MUSLE)对研究区非点源污染进行关键区识别及分级。结果表明：关键区按所占面积从大到小依次为低风险关键区、高风险关键区和中风险关键区，且各区主要沿水域分布。关键区中农业用地占到了流域农业用地总面积的37.3%，而随着风险等级的提高，关键区中农业用地所占比例显著增加。研究区非点源污染的主要来源是沿水域分布的农业用地，与实际情况吻合，研究结果能为研究区非点源污染治理提供科学指导和理论支持。将GIS与MUSLE模型相结合能够快速识别非点源污染关键区，该方法科学有效，具有很强的适用性。

非点源污染; 关键区; GIS; MUSLE

国内外大量研究表明，在点源污染得到基本控制后，伴随水土流失而产生的非点源污染成为流域水环境恶化的主要原因[1]。据统计，全世界30%～50%的地表已经受到非点源污染物的影响[2]，而在我国，继点源污染后，来自农业的非点源污染已经成为最重要的环境污染源[3]。但是由于非点源污染具有以下特征：污染发生随机性；机理过程复杂性；排放途径不确定性；污染负荷的时空差异性，因此对其的监测、模拟与控制变得十分困难[4]。而通常少数景观单元输出的污染物往往占了整个流域污染负荷的大部分[5]，这些地区就成为非点源污染的关键源区。有研究显示，磷的关键源区只占整个流域的小部分，磷流失量的90%来源于10%的区域，且主要分布于靠近河流的区域[6]，因此对关键区的识别和治理变得尤为重要。但传统的非点源污染监测方法需要耗费大量的人力、物力和财力，且难以获得科学的流域范围污染评价结果[7]，在我国环境问题突出和环境保护投入有限的背景下，如何识别非点源污染关键区，有针对性地制定控制与管理措施，从而提高污染治理的成效成为迫切需要解决的课题[8]。

地理信息系统(GIS)具有强大的空间数据管理能力，并具有空间数据采集、可视化和全面的空间分析功能。这些优点契合了非点源污染监测和模拟过程中的难点，为非点源污染关键区识别提供了现实的工具和方法，且具有费用低、周期短、简便易行等特点，因此GIS和非点源污染的有机结合成为了必然。

1 研究区概况

东北沟小流域位于河北省平泉县西部，地理坐标为118°3′39″—118°37′9″E，41°3′7″—41°6′19″N，属燕山山脉中低区，流域总面积19.4 km2，海拔593～1 179 m。多年平均降雨量为540.8 mm，气候为温带大陆性季风气候，因此降水季节性变化明显，5—8月降雨量占全年降水量的75%左右，且降雨多以暴雨形式出现，降雨强度大、侵蚀力强，为非点源污染的发生创造了条件。流域内土壤类型以褐土、棕壤和石质土为主，共占流域面积的86%，其中褐土占到了流域面积的38%。流域地势由东北和西北向西南倾斜，地貌类型属于土石山区，同时它也属于流入潘家口水库的滦河水系，是京津地区重要的水源地。

2 研究方法

土壤侵蚀和非点源污染是密不可分的。Maas等[9]发现，土壤流失严重的地区通常也是农业非点源污染发生的关键区域。Schauble[10]认为，土壤侵蚀不仅包括泥沙颗粒的运移，而且还包括随之发生的营养物和污染物的运移。这两个过程都是通过地表径流同时发生的，具有不可分离性。所以，侵蚀模型可以用来识别非点源污染的关键区。

2.1 模型选择

本研究采用Sivertun等[11]提出的MUSLE模型，它是在USLE的基础上改进而来的。表达式如下：

P=K·S·W·U

(1)

式中：P——发生非点源污染的风险；K——土壤因子；S——坡度因子；W——河道因子；U——土地利用因子。和USLE相比，该模型没有给出土壤侵蚀和污染物的定量信息，但可以直观分辨出对水质产生影响的高风险区域。该模型具有以下几个优点：(1) 快速识别关键区，能够为后续的深入研究提供基础；(2) 省时省力，对硬件和操作人员没有特殊要求，在一般的商业GIS软件上即可运行；(3) 易于更新，当其中一个因子发生变化时，并不影响其它因子，整个模型可以方便进行重新计算，所以它也能够对区域内的污染风险进行时空变化分析。

2.2 MUSLE各因子的计算方法

2.2.1 土壤K因子 土壤K因子即土壤可蚀性因子，即土壤遭受侵蚀的敏感程度，是指单位侵蚀力在标准小区上所产生的土壤流失量。在同等条件下，土壤K值越大，土壤受侵蚀的潜在危险性越大。土壤可蚀性可以通过小区测定直接得到，但这种方法非常繁琐，而且要求十分苛刻。所以本研究采用EPIC[12]模型计算土壤因子K

K={0.2+0.3exp[-0.0256Sd(1-Si/100)]}×[Si/(Cl+Si)]0.3×{1-0.25C/[C+ exp(3.72-2.95C)]}×{1-0.7Sn/[Sn+ exp(-5.51+22.9Sn)]}

(2)

式中：Sd——砂粒质量分数；Si——粉粒质量分数；Cl——黏粒质量分数；C——有机碳质量分数。Sn=1-Sd/100。

经过公式计算求得不同土壤类型的K值：褐土为0.237 (t·h)/(MJ·mm)，棕壤为0.172 (t·h)/(MJ·mm)，石质土为0.172 (t·h)/(MJ·mm)，新积土为0.154 (t·h)/(MJ·mm)，粗骨土为0.302 (t·h)/(MJ·mm)，黄绵土为0.332 (t·h)/(MJ·mm)。

2.2.2 地形因子 在模型的早期应用中，只考虑了坡度因子S，而忽略了坡长因子L对土壤侵蚀的影响，这是不科学的。因为在土壤侵蚀的过程中，当坡度相同时，坡长越长，径流速度越大，汇流的径流越大，侵蚀力越强。因此，Mitasova等[13]对该模型进行了改进，用LS因子即地形因子，取代了S因子，使模型更加完善，其计算公式如下：

LS(r)=(m+1)[A(r)/a0]m[sinb(r)/b0]n

(3)

式中：A——上坡每单位等高线宽度的贡献面积(m2)；b——坡度(°)；m，n——常数，分别取0.6和1.3；a0，b0——标准小区的坡长和坡度，即a0=22.1 m，b0=5.16°=0.09。

2.2.3 河道因子 距离河道的远近是污染物质扩散的一个非常重要的因子，Sivertun等[14]在河道因子的计算中建立了一个权重函数(公式4)，根据每个栅格距离河道的实际距离赋予栅格权重值，从而可以更加精确地表达距离河道不同距离的泥沙或污染负荷进入河道的比例。

(4)

2.2.4 土地利用因子 不同的土地利用方式决定了不同营养物质的含量及发生土壤侵蚀的风险，McElroy等[15]综合考虑了各种因素，经过验证，为不同土地利用类型确定了不同的权重函数。本研究参考了这个标准，各土地利用方式的权重值分别为：水体为0，林地为0.005，草地为0.010，居民点0.030，果园为0.075，耕地为0.100。

3 模型运算

MUSLE模型适用性很强，可以在一般的商业GIS软件上运行，如ArcVIEW或ArcGIS等。本研究使用ArcGIS 9.3对模型进行运算，使用的功能包括水文分析、重分类、数据格式转换、栅格计算器等。计算流程如下：

计算过程所需数据包括土壤类型图、DEM、水系图和土地利用图。所有因子图均采用Albers等积圆锥投影，栅格大小为5 m×5 m。由于MUSLE方程推荐在较小地块上进行运算，因此本研究符合这一要求，并能够较为详细地识别非点源污染关键区。

(1) 土壤因子图。在研究区土壤类型矢量文件的属性表中新建属性列，命名为栅格值。将计算得到的K值按照类别输入到这个属性列中。然后运用ArcToolbox中的转换功能，以栅格值属性列为输出域，将土壤类型矢量图转换为土壤K值栅格图。

(2) 地形因子图。在ArcGIS 9.3下使用空间分析模块，对小流域DEM填洼，然后生成坡度图，命名为Slope；填洼后进行水流方向和汇流累积量计算，将汇流累积量图中所有河流及水体赋值为0，这样可以避免产生不符合实际的高流水累积值，把处理过的流水累积图命名为Accumulation。然后在栅格计算器中运行以下表达式，得到地形因子图。

Pow([Accumulation]×5/22.1,0.6)×Pow(sin([slope]×0.01745)/0.09,1.3)×1.6

(3) 土地利用因子图。对东北沟小流域土地利用矢量图进行修改，同样在属性表中新建属性列，命名为栅格值，按照上文确定的权重值分别输入到栅格值属性列中，并以其为输出域，将土地利用矢量图转换为土地利用栅格图。

(4) 河道因子图。在ArcGIS 9.3中应用空间分析功能，以小流域DEM为源数据，以小流域矢量河道为目标，计算每个栅格距离水体的距离，命名为watercourse。然后在生成的栅格上利用以下表达式进行计算，使各栅格得到合理的权重值。

0.6/{exp[0.002×(watercourse)]-0.4}

4 结果与分析

运用以上方法得出各因子图，如附图1和图1所示。使用ArcMap中的Raster Calculator对4个因子图进行相乘合并，成图的栅格大小为5 m×5 m，栅格值范围为0～22.34，代表了污染风险的大小。根据标准差对成图进行重分类，可将流域分为4个类型区：安全区、低风险区、中风险区和高风险区。其中低风险区、中风险区和高风险区即为非点源污染关键区。分级标准如表1所示。

图1 地形因子和河道因子栅格图

污染风险分区栅格值范围安全区0～0.095(低于平均值)低风险区0.095～0.257(高于平均值0～1个标准差)中风险区0.257～0.419(高于平均值1～2个标准差)高风险区0.419～22.339(大于2个标准差)

在ArcMap中给出了各类栅格的统计信息，通过计算栅格数可以得出各分区的面积及各区所占比例。如表2所示。

表2 风险区统计

从表中可以看出，非点源污染关键区的面积占到流域总面积的4.2%；另经计算，关键区中农业用地占流域农业用地总面积的37.3%，而关键区中非农业用地仅占流域非农业用地总面积的1.5%；关键区中农业用地面积占到了66.1%，中高风险关键区中农业用地面积占到了86.3%，高风险关键区中农业用地面积所占比例则高达92.4%。这说明农业非点源污染是小流域非点源污染的主要来源，且随着风险等级的提高，农业用地逐渐占据主导地位。

通过分析东北沟小流域非点源污染风险分级图(附图2)可以看出，污染关键区主要沿水域分布，中上游均有集中分布区，这是因为：(1) 流域内地形东北、西北高，向西南逐渐降低，高海拔的污染关键区主要由陡坡农业用地构成，而低海拔区域虽然地形因子的影响相对较小，但受到其它因素的影响发生污染的可能性也很大。(2) 从土壤类型来看，黄绵土和粗骨土沿河流分布，发生土壤侵蚀的风险较高，因此这些区域也成为了非点源污染的关键区。(3) 从土地利用方式来看，流域内农业用地多分布在水域附近，且密集程度高，施用的农药化肥随着水土流失的发生对水体造成污染，因此这些区域具有较高的非点源污染风险。

5 结 论

(1) 流域中非点源污染关键区面积为0.82 km2，占流域总面积的4.2%，且主要沿水域分布。(2) 流域中非点源污染的主要来源是农业用地，且随着风险等级的提高，农业用地所占比例显著增加。(3) 沿水域分布的黄绵土和粗骨土发生土壤侵蚀的风险较高，这对污染物的运移起到了非常积极的作用，因此这些区域也是非点源污染的关键区。

针对以上问题，应从源头和过程两个方面采取相应措施：(1) 鼓励农民科学合理配备化肥和农药施用量，推广高效、低毒、低残留农药。(2) 积极实施农田水土保持工程，减少农田水土流失，抑制污染物的运移过程。(3) 在非点源污染关键区与水体间建立植被缓冲区等措施，对污染物进行阻截、吸收和转化。

[1] 张维理,吴淑霞,冀宏杰,等.中国农业面源污染形势估计及控制对策.Ⅰ:21世纪初期中国农业面源污染的形势估计[J].中国农业科学,2004,37(7):1008-1017.

[2] Pimental D. World soil erosion and conservation[M]. Cambridge: Cambridge University Press,1993.

[3] 胥彦玲,王苏舰,李怀恩.土地覆盖变化对流域非点源污染的影响研究:以黑河流域为例[J].水土保持研究,2010,17(3):251-253.

[4] 周健.试论农业非点源污染的危害[J].农业环境保护,1990,9(1):22-25.

[5]Gburek W J, Sharpley A N. Hydrologic controls on phosphorus loss from upland agricultural watersheds[J]. Journal of Environmental Quality,1998,27(2):267-277.

[6] 张淑荣,陈利顶,傅伯杰,等.农业区非点源污染潜在危险性评价:以于桥水库流域磷流失为例[J].第四纪研究,2003,23(3):262-269.

[7] 杨菁荟,张万昌.SWAT模型及其在水环境非点源污染中的应用研究进展[J].水土保持研究,2009,16(5):260-266.

[8] 周军,张颖,高凤杰,等.阿什河流域农业非点源污染优先控制区域识别[J].水土保持研究,2013,20(3):180-184.

[9] Maas R P, Smolen M D, Dressing S A. Selecting critical areas for nonpoint-source pollution control[J]. Journal of Soil and Water Conservation,1985,40(1):68-71.

[10] Schauble T. Erosionsprognosen mit GIS und EDV-Ein Vergleich verschiedener Bewertungskonzepte am Beispiel einer Gaulandschaft[D]. Geographisches Institut, Universitat Tubingen,Germany,1999.

[11] Sivertun A K E, Reinelt L E, Castensson R. A GIS method to aid in non-point source critical area analysis[J]. International Journal of Geographical Information System,1988,2(4):365-378.

[12] Williams J R, Renard K G, Dyke P T. EPIC: A new method for assessing erosion′s effect on soil productivity[J]. Journal of Soil and Water Conservation,1983,38(5):381-383.

[13] Mitasova H, Hofierka J, Zlocha M, et al. Modelling topographic potential for erosion and deposition using GIS[J]. International Journal of Geographical Information Systems,1996,10(5):629-641.

[14] Sivertun Å, Prange L. Non-point source critical area analysis in the Gisselö watershed using GIS[J]. Environmental Modelling & Software,2003,18(10):887-898.

[15] McElroy A D. Loading functions for assessment of water pollution from nonpoint sources[M]. US Environmental Protection Agency, Office of Research and Development,1976.

IdentificationofCriticalNon-pointSourcePollutionAreasinDongbeigouWatershedBasedonGISandMUSLE

LIU Nan1, XIE Yong-sheng1,2, SUO Gai-di1, JING Min-xiao1, CHEN Lei2

(1.NorthwestA&FUniversity,Yangling,Shaanxi712100,China; 2.InstituteofSoilandWaterConservation,ChineseAcademyofScienceandMinistryofWaterResources,Yangling,Shaanxi712100,China)

Based on DEM and the survey data of land use and soil of Dongbeigou watershed, this paper combined GIS and MUSLE model to identify the critical non-point source pollution areas, and then graded the areas. The results showed that: (1) in descending order, the critical areas were ranged as low risk area, high risk area and moderate risk area. The critical areas mainly distributed along the river; (2) among the critical areas, the agricultural land took 37.3 percent of total agricultural land in the watershed, the proportion of the agricultural land increased significantly as the risk level upgraded, the agricultural lands distributing along the river were the main sources of the pollution, which was in accord with actual situation. The result can provide the basis for further pollution control. Combined with GIS technology, MUSLE model can rapidly identify the critical non-point source pollution areas. This method is scientific and reasonable, and it is of strong adaptability.

non-point pollution; critical areas; GIS; MUSLE

2013-06-10

：2013-09-02

国家科技支撑计划课题“农田水土保持工程与耕作关键技术研究”(2011BAD31B01);水利部公益性行业科研专项“工程开挖面与堆积体水土流失测算技术研究”(201201048);“风力作用下扰动地表侵蚀预报关键技术研究”(201201047)

刘楠(1989—),男,河北石家庄人,硕士研究生,主要从事GIS在环境监测与评价中的应用。E-mail:freestyleconan@nwsuaf.edu.cn

谢永生(1960—),男,河南开封人,研究员,主要从事环境评价、水土保持及土地资源等方面研究。Email:ysxie@ms.iswc.ac.cn

X171.5;TP79

：A

：1005-3409(2014)02-0123-04