汉江中下游土壤侵蚀及颗粒态非点源磷负荷研究

沈 虹,张万顺,彭 虹

（1.北京大学城市与环境学院,地表过程分析与模拟教育部重点实验室,北京100871;2.武汉大学资源与环境科学学院,武汉430079;3.武汉大学水利水电学院,武汉430072;4.北京大学深圳研究生院,城市人居环境科学与技术重点实验室,广东深圳 518055）

非点源污染研究最早始于20世纪60年代美国和日本所开展的土地利用与湖泊富营养化关系调查。而后,随着计算机、3S技术的发展,以及野外实验、室内分析、水质监测等方法的改进完善,描述地表水文过程、土壤侵蚀过程、污染物迁移转化过程等的数学模型随之出现[1],各种不同类型的非点源模型相继问世并得到广泛应用[2-3]。非点源污染研究也从统计分析向半机理及机理过程转变[4],研究尺度从小流域扩展到中大型流域甚至国家尺度。相对于国外,国内研究起步虽晚但发展迅速,已有的研究案例覆盖了不同尺度、不同类型的非点源,主要的研究关注点包括河网水系的分维影响[5]、负荷的模拟估算及空间分布特征[6]、污染物的迁移转化[7]等。因此,建立定量化模型仍是当前非点源污染研究的重点。

根据负荷产生和迁移过程,非点源负荷可以分为溶解态和颗粒态两种类型[8]。溶解态负荷主要是指水溶性污染物,颗粒态负荷则是指在侵蚀性降雨作用下,被冲刷剥蚀吸附在土壤及泥沙中的污染物,它伴随着泥沙迁移进入河道,最后在河流中解析形成的水体污染[9]。水土流失、土地利用及水肥管理方式是颗粒态非点源负荷形成的重要原因,它不仅携带了大量地表肥沃的土壤,同时还带走了大量累积在地表的磷污染物质[10-11]。因此,有必要对土壤流失与颗粒态负荷的定量关系进行深入研究。目前常用的研究方法是GIS与USLE结合,已广泛应用于黑河[12]、九龙江[13-14]、巢湖[15]等流域。

汉江作为长江的重要支流,其中下游水质状况对于保障沿岸饮水及工农业生产安全具有重要的意义。然而,近年来随着汉江中下游水质的恶化,非点源污染已成为流域水体污染的主要影响因素之一。目前关于汉江中下游水土流失以及颗粒态非点源磷负荷的研究较少[16-17],研究范围也主要集中在流域沿岸,多以行政意义界定流域边界,忽视了污染负荷在迁移转化过程中所依赖的自然地理要素,且方法以综合层次判别分析和经验方法为主,缺少定量的模型化研究。因此,本研究基于GIS平台,采用USLE和负荷输出经验方程对汉江中下游土壤侵蚀和颗粒态非点源负荷进行了定量估算。按照自然地理要素和时空过程一致性原则,以栅格作为基本计算单元,以子流域作为基本的汇总和模拟单元,分别研究了影响土壤侵蚀的五大因素,通过空间叠加方法得到了流域土壤侵蚀状况,并在此基础上估算了非点源负荷及其空间分布特征,以期为流域水土保持和颗粒态非点源防控提供理论依据。

1 流域概况及环境数据

1.1 流域概况

本研究采用的流域范围以自然边界为准,利用DEM基于集水区生成,地跨湖北、河南两省(30°-34°N,100°-114°E),上起丹江口,下至汉江入长江河口,全长616 km,总流域面积6.078 4万km2,涉及湖北省的26个县、市、区和河南省的13个县、市。流域内拥有河南南阳盆地和湖北江汉平原两大农业区,孕育了唐白河、南河、北河、小清河、蛮河、竹皮河、汉北河等河流水系。近年来汉江中下游先后5次发生富营养化现象,水质状况日益堪忧,非点源磷污染成为其重要原因之一。

1.2 基础数据及预处理

本研究中所有的基础数据(表1)均被纳入到统一的坐标系统下:Albers等面积圆锥投影,并转化为100 m×100 mGrid格式。为了体现流域内地理要素的空间异质性,在ArcSWAT中对凹陷填充后的DEM进行了再分割。经过反复对比分析发现,当集水区面积阈值等于30 000 hm2时,生成的水系与实际情况具有最大相似性,此时共划分31个子流域,得到流域总面积为6.0784万km2(图1);同时,在DEM 的基础上进行坡度、坡向分析,将坡度分为五组0～7.5°,7.5°～12.5°,12.5°～ 17.5°,17.5°～ 22.5°,22.5°～ 90°,合并破碎多边形得到坡度图。将坡向按东北、东、东南、南、西南、西、西北、北8个方向进行分组得到坡向图。解译流域2007年Landsat TM影像,经过几何校正、影像拼接、直方图匹配、颜色纠正等,得到疏林地、有林地、水田、旱地、建设用地、水域和未利用地7类土地利用类型;土壤数据来源于1∶400万中国土壤分布图、中国土壤数据库、湖北省和河南省土种志,包括土壤的空间分布、质地组成、有机质含量、TP背景值等,制作土壤磷元素含量和其他气象要素等相关专题图。

表1 基础数据及说明

图1 汉江中下游子流域分布图

2 研究方法

以栅格为非点源的基本计算单元,以子流域为非点源负荷汇总的基本单元,基于GIS、RS技术,采用ULSE方程和负荷输出经验方程对颗粒态非点源负荷量进行估算。

2.1 土壤侵蚀模数估算

USLE(universe soil loss equation)是一种经验统计模型,它由Wischmeier于20世纪60年代提出。模型的基本思想是侵蚀模数大小取决于引起侵蚀的力与土壤抵抗降雨冲刷的能力的组合,计算结果表述为降雨引起的单位面积土壤流失量[18]。目前,将USLE与GIS/RS结合已经成为土壤侵蚀研究的重要方法[19-20]。

式中:A——年均土壤侵蚀量[t/(hm2◦a)];α——单位转换系数,为224.2;R——降雨侵蚀力因子[(MJ◦mm)/(hm2◦h◦a)];K代表土壤可蚀性因子[(t◦hm2◦h)/(hm2◦MJ◦mm)];LS——坡度坡长因子,无量纲;C——经营管理因子,无量纲;P——水土保持措施因子,无量纲。

2.1.1 降雨侵蚀力因子R 降雨侵蚀力旨在评价降雨引起土壤侵蚀的潜在能力,它与降雨量、降雨历时、降雨强度、雨滴大小等因素有关[21]。利用研究区1961-2006年月平均降雨量和年平均降雨量资料,根据Wischmeier经验公式估算R[22],表达式为

式中:Pi——月均降雨量(mm);P——年均降雨量(mm)。

2.1.2 土壤可蚀性因子K 土壤可蚀性因子K是土壤侵蚀性的量度,它反映了不同类型土壤抵抗侵蚀能力的高低,与土壤质地和有机质含量密切相关。目前K估算主要基于小区实验展开,常用的经验估算办法是Williams等人提出的土壤质地、有机质含量测定方法。本研究在采集汉江中下游土壤资料、图件、土壤类型空间分布等资料的基础上,通过查表和参考众多文献研究[23-25],确定了研究区土壤可蚀性因子K值。

表2 汉江中下游土壤可蚀性因子K值

2.1.3 坡长坡度因子LS 坡度坡长因子代表地形因素变化对土壤侵蚀的影响。坡度、坡长与侵蚀作用力呈现正向关系,当坡度和坡长因子越陡或越长时,径流速度将越大,相应的,在坡面汇集的径流所引起的侵蚀作用力也将增强。本研究引用王宁等[26]关于地形因子LS的估算方法,得到坡长坡度因子LS的计算表达式如下:

式中:LS——坡度坡长因子,无量纲;l——坡长(m);θ——坡度(°);m——坡长指数,其值与坡度相关。

2.1.4 经营管理因子C 经营管理因子C是指在相同土壤、坡度和降雨条件下,某种特定土地覆被的土壤流失量与耕作后连续休闲地的土壤流失量比值,可以反映地表覆盖、土地利用以及一般管理措施对土壤侵蚀速度的影响。由于经营管理因子的影响因素非常多,通常采用简化处理,仅考虑最主要影响因素植被覆盖度与经营管理因子C的作用关系[27-28]。资料研究证实,NDVI与植被覆盖具有良好的线性关系[29]。因此本研究将遥感解译的NDVI作为桥梁,利用NDVI推导出植被盖度(公式4),进而建立植被盖度与经营管理因子的经验公式(公式5)。其值介于0～1,当植被盖度为0时,土地裸露,最易造成土壤侵蚀,因此经营管理因子 C最大;当C≥78.3%时,植被盖度较高,可以阻挡降雨动能造成的侵蚀,近似认为此时经营管理最完善,经营管理因子C赋值为0。

式中:NDVI(Normalized difference vegetation index)——植被归一化指数,无量纲;lc——植被盖度,无量纲;C——经营管理因子,无量纲。

2.1.5 水土保持因子P 水土保持因子P指采取水土保持措施后的土壤流失量与顺坡种植时的土地流失量比值,反映了水土保持措施对于土壤侵蚀的控制作用,数值介于0～1。目前,水土保持因子P的求取方法有基于坡度分级赋值和基于不同土地利用赋值等方式。

考虑到研究区的特点,本研究采用土地利用赋值方法。在查阅流域水保措施资料及参考大量文献的基础上[30-31],确定研究区水土保持因子P值。水土保持因子P值越大表明水保措施越差,土壤潜在侵蚀的可能性越大,例如未利用地通常为裸地,可近似认为未采取任何水保措施,因此水土保持因子P值为1;而对于水域,由于不会产生水土流失,因此可视为拥有完善的水保措施,赋值为0。

表3 汉江中下游水土保持因子P值

2.2 颗粒态氮磷负荷估算

参照相关研究确定颗粒态非点源负荷的经验输出模型为[7,13-15,17,21,32]:

式中:Wxp——颗粒态非点源磷负荷[kg/(hm2◦a)];α——单位换算常数,无量纲;Cp——颗粒态磷负荷浓度,来源于中国科学院南京土壤研究所提供的中国土壤数据库、湖北省和河南省土种志;A——土壤流失量[t/(km2◦a)];η——富集系数,无量纲 ;Sd——流域土壤侵蚀迁移到河道的运移比,无量纲。

(1)土壤磷浓度计算。土壤中磷含量是土壤流失的直接补给,它反映了土壤流失的潜在能力。本研究中采用的土壤磷浓度来源于中国土壤数据库、湖北省和河南省土种志(表4)。

表4 土壤中颗粒态非点源磷负荷输出浓度 %

(2)富集比。侵蚀形成的泥沙具有富集养分元素的特性,同雨前表土养分含量相比,侵蚀形成的泥沙养分含量较高,这种现象被称为泥沙富集作用[33],一般用某养分在侵蚀泥沙中的含量与其在被侵蚀土壤中的含量之比表示。本研究在参阅同类研究的基础上[12,17,33],确定研究区富集比为2。

(3)迁移比Sd。颗粒态非点源磷负荷伴随土壤侵蚀、迁移的过程中会出现损耗,只有真正进入河道中的污染负荷才会对流域水质造成影响,并且距离河道的远近直接影响着运移损耗。根据长江水利委员会对长江流域长期定位的研究结果:长江流域的Sd值大约为0.1～0.4[17],据此确定研究区内平均泥沙输移比为0.28,最大值0.4,最小值0.1。利用Arc-GIS空间分析工具计算每个栅格单元到最近河道的距离,按照距离倒数赋值输移比,得到 Sd空间分布图。

3 结果与分析

3.1 土壤侵蚀估算及验证

2007年研究区土壤的平均侵蚀模数为867.55 t/(km2◦a),最高侵蚀模数为子流域 21,高达5 158.75 t/(km2◦a),分布在神农架林区。最低侵蚀模数为子流域28,低至111.42 t/(km2◦a),分布在竹皮河流域附近。按照1997年水利部颁布的《土壤侵蚀分类分级标准》将研究区土壤侵蚀划分为八级,得到侵蚀程度分类图(附图1)。总体上,汉江中下游流域属于轻度侵蚀区,但依然有5.69%的区域侵蚀模数大于2 500 t/(km2◦a),其中,中度侵蚀1.8%,强度侵蚀1.21%,极强度侵蚀1.32%,剧烈侵蚀1.37%,土壤侵蚀的主要分布区以中低山和丘陵地带为主,而土壤侵蚀的最严重区域则以中高山农耕地为主,如神农架林区、房县、保康、谷城部分区域及嵩县、内乡等地带,这与文献[16]的研究结论较为一致。

3.2 非点源负荷空间分布

研究区颗粒态非点源磷负荷为8 496.79 t/a,对比2002年史志华等人[17]在该流域的计算结果8 893.4 t/a,颗粒态非点源磷负荷的总量有所减少,但空间分布具有一致性,支流唐白河、南河、蛮河等非点源磷负荷的入汇量较高。从子流域单位负荷的角度分析,子流域21,16,15,12,22,2,24颗粒态磷单位负荷较高(附图2),其中最高为子流域21,单位负荷高达8.42 kg/(hm2◦a)。具体来说,这些子流域主要分布在三片区域,一是神农架林区、房县、保康、谷城、南漳西北部一带,平均颗粒态磷单位负荷为4.1 kg/(hm2◦a);二是嵩县、内乡一带,平均颗粒态磷单位负荷量为3.0 kg/(hm2◦a);三是钟祥东北部,平均颗粒态磷单位负荷为1.45 kg/(hm2◦a)。在神农架、房县、保康一带虽然有林地面积大,但海拔高,坡度大,平均坡度在22.5°以上,且在经济发展的带动下,将坡地改成耕地种植高山玉米、蔬菜等现象较普遍,因此土壤流失最为严重,随水土流失携带的颗粒态非点源负荷亦最高。在嵩县、内乡一带,其北部平均坡度达22.5°以上,坡耕地亦是造成土壤流失和非点源负荷的最重要原因;而在南部,作为南阳盆地的一部分,磷肥的大量施用是颗粒态负荷增高的重要原因之一。在钟祥东北部,山地丘陵是最主要的地貌类型,土地利用方式多为旱地、疏林地、未利用地等,植被覆盖度较低,导致坡长坡度因子LS、经营管理因子C、水土保持措施因子P较高,因此土壤流失产生的非点源负荷量高,此外,分布在该地区的石灰土,土壤T P含量较高也是部分原因。自钟祥以下的汉江下游河段,颗粒态非点源磷负荷均较低,平均单位磷负荷为0.3 kg/(hm2◦a),见图2。

表5 汉江中下游土壤侵蚀模数分级表

图2 子流域颗粒态非点源单位负荷量和总量关系

表6 各支流非点源磷负荷的汇入情况

在实际负荷迁移转化过程中,颗粒态负荷不是直接通过子流域汇入干流,而是由支流携带进入主干河道,因此将子流域按照主要支流进行汇总分析有利于指导支流的非点源防控工作。在本研究中,将31个子流域按构成分组为9个主要支流和直接入江10种情形,分别为北河、南河、小清河、唐白河、蛮河、双河、利河、竹皮河、汉北河和直接入江。从支流入汇的非点源负荷分析(表6),研究区内颗粒态非点源磷负荷主要来源于支流入汇,直接入江的颗粒态非点源磷负荷仅占16.87%。唐白河由于流域面积大,同时又流经“中州粮仓”南阳盆地,农业磷肥广泛施用使得流失的土壤中含有大量的颗粒态非点源磷负荷,因而成为最大的颗粒态非点源磷负荷贡献者,总量达到3 079.26 t,贡献率达36.24%;南河因流经神农架、保康、房县等土壤流失最严重区,颗粒态非点源磷负荷较高,达到2 454.45 t,贡献率为28.89%,居第二位。这里山地绵延、地势陡峻,且坡地耕作,植被渐被破环,侵蚀流失的负荷量大;蛮河流经南漳区域,磷负荷量也较高,位居第三。北河、汉北河、利河的颗粒态非点源磷负荷贡献率居中,为2.75%～3.72%。小清河、双河、竹皮河因面积小,单位负荷量低,因而对颗粒态氮磷负荷总量的贡献率低,磷的贡献率均小于等于1.23%。

4 结论

(1)汉江中下游2007年土壤的平均侵蚀模数为867.55 t/(km2◦a),总体上属于轻度侵蚀,但依然有5.69%的区域侵蚀模数大于2 500 t/(km2◦a)。流域内由于土壤流失引起的颗粒态非点源磷负荷流失量达8 496.79 t/a,其中南河、唐白河、蛮河等支流的非点源单位负荷量和总负荷量均较高,最高单位负荷量可达8.42 kg/(hm2◦a),与水土流失区严重区基本一致。该区域主要分布在神农架林区、房县、保康等山地丘陵以及南阳盆地等农业发达地区,与坡耕地耕作、农田化肥施用紧密相关,是今后水土流失和非点源防控的关键区,因此在防控非点源工作中,要从颗粒态负荷产生的源头做起,一是要加大森林覆盖,二是要严禁坡地耕作,防控和减少水土流失。

(2)由于部分河口泥沙等相关数据较难获取,因此本研究中无法利用实测数据对估算的颗粒态负荷量进行验证。同时,在本研究中以负荷的经验输出公式代替了负荷复杂的迁移转化机理,在一定程度上会影响估算的精度,因此有待在今后的工作中进一步改进和完善。

[1]程红光,岳勇,杨胜天,等.黄河流域非点源污染负荷估算与分析[J].环境科学学报,2006,26(3):384-391.

[2]Michele Munafòa,Giuliano Cecchib,Fabio Baioccoa,et al.River pollution from non-point sources:a new simplified method of assessment[J].Journal of Environmental Management,2005,77(2):93-98.

[3]Polyakova V,Faresa A,Kubob D,et al.Evaluation of a non-point source pollution model,AnnAGNPS,in a tropical watershed[J].Environment Modelling and Software,2007,22(11):1517-1627.

[4]张秋玲,陈英旭,俞巧钢,等.非点源污染模型研究进展[J].应用生态学报,2007,18(8):1886-1890.

[5]荣琨,陈兴伟,王林,等.流域非点源污染与水系分维关系研究[J].水土保持研究,2009,16(2):128-131.

[6]王晓燕,王晓峰,汪清平,等.北京密云水库小流域非点源污染负荷估算[J].地理科学,2004,24(2):227-231.

[7]李怀恩,沈晋.流域非点源污染模型的建立与应用实例[J].环境科学学报,1997,17(2):141-147.

[8]李继承.嘉陵江流域非点源污染负荷模拟研究[D].重庆:重庆大学,2007.

[9]Bardgett R D,Anderson J M,Behan-Pelletier V,et al.The influence of soil biodiversity on hydrological pathways and the transfer of materials between terrestrial and aquatic ecosystem[J].Ecosystems,2001,4:421-429.

[10]程文娟,史静,夏运生,等.滇池流域农田土壤氮磷流失分析研究[J].水土保持学报,2008,22(5):52-55.

[11]张玉斌,郑粉莉,武敏.土壤侵蚀引起的农业非点源污染研究进展[J].水科学进展,2007,18(1):123-132.

[12]胥彦玲,李怀恩,倪永明,等.基于 USLE的黑河流域非点源污染定量研究[J].西北农林科技大学学报:自然科学版,2006,34(3):138-142.

[13]黄金良,洪华生,张珞平,等.基于GIS和 USLE的九龙江流域土壤侵蚀量预测研究[J].水土保持学报,2004,18(5):75-79.

[14]黄金良,洪华生,张珞平,等.基于GIS的九龙江流域农业非点源氮磷负荷估算研究[J].农业环境科学学报,2004,23(5):866-871.

[15]汪晓辉,巢湖流域非点源N、P污染排放负荷估算及控制研究[D].合肥:合肥工业大学,2006.

[16]史志华,蔡崇法,丁树文,等.基于GIS的汉江中下游土壤侵蚀强度初步分析[J].水土保持学报,2001,15(6):14-17.

[17]史志华,蔡崇法,丁树文,等.基于GIS的汉江中下游农业面源氮磷负荷研究[J].环境科学学报,2002,22(4):473-477.

[18]曾凌云.基于RUSLE模型的喀斯特地区土壤侵蚀研究:以贵州红枫湖流域为例[D].北京:北京大学,2008.

[19]徐清艳,周跃.基于GIS和 USLE的突然侵蚀预测研究:以云南新平大红山铁矿为例[J].水土保持研究,2009,16(4):72-75.

[20]张有全,宫辉力,赵文吉,等.基于GIS和 USLE的密云县土壤侵蚀评价及空间特征研究[J].水土保持研究,2007,14(3):358-364.

[21]杨胜天,程红光,步青松,等.全国土壤侵蚀量估算及其在吸附态氮磷流失匡算中的应用[J].环境科学学报,2006,26(3):366-374

[22]Wischmeier W H,Smith D D.Predicting Rainfall Erosion losses:a Guide to Conservation[M].Washington DC:US Department of Agriculture,Science and Education Administration,Agricultural Handbook,1978.

[23]汪东川.基于ArcGIS的坡面产流产沙过程模拟:以重庆市璧山县为例[D].重庆:西南大学,2005.

[24]邓良基,侯大斌,王昌全,等.四川自然土壤和旱耕地土壤可蚀性特征研究[J].中国水土保持,2003(7):23-25.

[25]张科利,彭文英,杨红丽.中国土壤可蚀性值及其估算[J].土壤学报,2007,44(1):7-13.

[26]王宁,朱颜明,徐崇刚.GIS用于流域径流污染物的量化研究[J].东北师范大学学报:自然科学版,2002,34(2):92-98.

[27]赵琰鑫,张万顺,王艳,等.基于3S技术和USLE的深圳市茜坑水库流域土壤侵蚀强度预测研究[J].亚热带资源与环境学报,2007,2(3):23-28.

[28]蔡崇法,丁树文,史志华,等.应用 USLE模型与地理信息系统IDRISI预测小流域土壤侵蚀量的研究[J].水土保持学报,2000,14(2):19-24.

[29]高龙华.遥感和GIS支持下流域非点源污染模型研究[D].南京:河海大学,2006.

[30]卜兆宏,孙金庄.水土流失定量遥感方法及应用研究[J].土壤学报,1997,34(3):235-245.

[31]许月卿,邵晓梅.基于GIS和RUSLE的土壤侵蚀量计算:以贵州省猫跳河流域为例[J].北京林业大学学报,2006,28(4):67-71.

[32]施为光.四川省清平水库流域非点源污染负荷计算[J].重庆环境科学,2002,22(2):33-36.

[33]薛素玲.基于GIS的黑河流域非点源氮磷模拟[D].西安:西安理工大学,2006.