基于GIS的贡山县土壤侵蚀动态变化研究

李益敏, 袁 静, 王东驰, 刘心知

(1.云南大学 地球科学学院, 昆明 650500; 2.云南省地理研究所, 昆明 650500)

土壤侵蚀是土体在特定时间和空间条件下的迁移过程，是世界上最重要的土地退化问题[1]。土壤侵蚀是最活跃、最敏感的生态灾害之一，它可能会导致土壤肥力和质量下降、生态环境恶化、粮食减产、库塘湖泊淤积等一系列危害，已成为制约人类生存和发展的重要因素之一[2]。随着RS和GIS技术的发展，土壤侵蚀研究向定量化空间化的方向发展，并呈现动态变化研究应用趋势[3-4]。国内外的学者积极地开展了土壤侵蚀时空演变研究以反映区域土壤侵蚀在地理空间上的变化趋势。范建容等[5]基于GIS空间分析功能,开展了四川省李子溪流域土壤侵蚀动态变化研究,探讨了土壤侵蚀动态变化与土地利用变化的关系、侵蚀强度动态变化与侵蚀量变化的关系。信忠保等[6]利用GIS技术分析了黄土高原近50 a来土壤侵蚀的时空变化特征，揭示了降水减少和人类活动是黄土高原输沙强度减弱的重要原因。基于GIS，RS和RUSLE 模型，Alexakis等[7]分析了塞浦路斯Yialias流域土壤侵蚀速率的变化。史彦林[8]以GIS为平台，利用土壤侵蚀模型分析了饶河县土壤侵蚀的变化情况。Teng等[9]利用遥感影像的可见光—近红外光谱特征得到了土壤可蚀性因子K，并基于土壤侵蚀模型(RUSLE)对澳大利亚土壤侵蚀的时空演变规律进行研究。目前，土壤侵蚀动态变化研究已经成为水土保持领域的研究热点之一[10]，然而，针对高山峡谷地区的研究相对较为缺乏[11]。贡山县属于典型的高山峡谷地貌，生态环境脆弱，自然灾害频发，土壤瘠薄，加之气候变化的影响，导致贡山县水土流失情况较为严重[12]，因此，对贡山县土壤侵蚀的时空演变规律进行分析是必要的，可为该区域的水土保持工作和可持续发展提供重要的理论依据和数据支持。基于此，本文利用RUSLE模型，开展贡山县1996—2018年土壤侵蚀强度的时空动态变化特征研究，旨在为贡山县今后水土流失治理和生态环境建设提供科学依据。

1 研究区概况

贡山独龙族怒族自治县(下文称贡山县)隶属云南省怒江傈僳族自治州，位于云南省西北部，地理坐标为东经98°08′—98°56′，北纬27°29′—28°23′，下辖2镇3乡，东邻德钦县、维西傈僳族自治县，西部和西南部与缅甸接壤，南连福贡县，北接西藏自治区察隅县。全县面积4 506 km2，国境线长达172.08 km，贡山县属北亚热带季风气候，年平均气温14.8℃，雨量充沛，年降雨量为2 700～4 700 mm。贡山县土壤和植被的垂直地带性明显，土壤类型多样，多为自然土和黑泥土，还有部分棕泥土及黑泡土，植被类型复杂多样，森林覆盖率为77.2%，森林面积达3 478 km2，贡山县地势呈“三山夹两江”高山峡谷地貌，地质构造复杂，地形起伏大，最高海拔5 128 m，最低海拔1 170 m，海拔高差达3 958 m。全县共有大小河流42条，碧罗雪山、高黎贡山和怒江、独龙江纵贯贡山县境内，地势险陡，天然落差大，水能资源极为丰富，年径流量577亿m3，总径流量达84.5亿m3。

2 试验材料与方法

2.1 模型选择

美国农业部于1997年在通用水土流失方程USLE基础上改进的修正通用土壤流失方程RSULE模型，因其需要参数少，形式简单灵活，参数获取方便并且评价土壤侵蚀的各种影响因素较为全面而得到世界各地的普遍使用。与USLE模型相比，RUSLE模型在众多方面得到了改进，例如：填补了之前数据的空白，数据来源不再单一，能够模拟不同的系统等[13]。RUSLE模型在我国西南山区已有较长时间应用，适用性已得到有效验证[14]，因此本文选取RUSLE模型进行研究区土壤侵蚀的量化，公式为：

A=R·K·L·S·C·P

(1)

式中:A为由于降雨及其径流作用于坡面使其出现细沟或细沟间发生侵蚀所产生的土壤年平均侵蚀量[t/(hm2·a)];R为降雨侵蚀力因子[(MJ·mm)/(hm2·h·a)];K为土壤可蚀性因子[(t·hm2·h)/(hm2·MJ·mm)];LS为坡长和坡度因子,无量纲；C为植被覆盖和管理因子,无量纲；P为水土保持措施因子，无量纲

2.2 数据源及预处理

论文分析所用数据包括1996年、2006年、2012年、2018年日降雨量数据(来源于云南省气象台)；土壤属性和土壤类型数据通过HWSD查询砂粒含量、粉粒含量、黏粒含量和有机碳含量得到；由地理空间数据云平台下载得到DEM数据和贡山县1996年、2006年、2012年、2018年4个时期的Landsat TM遥感影像，对获取的DEM数据根据贡山县矢量边界进行裁剪，Landsat TM遥感影像空间分辨率为30 m×30 m,影像时间为10—12月，每景影像的含云量均小于5%，数据质量较好，并利用ENVI 5.3软件对遥感影像进行辐射定标、FLAASH大气校正、图像增强、图像镶嵌和裁剪等处理。NDVI由遥感影像反演获得，土地利用类型数据由遥感影像目视解译提取得到。

2.3 RUSLE各因子的确定

2.3.1 降雨侵蚀力因子 降雨产生土壤侵蚀主要是通过雨滴溅落对土壤颗粒的冲击和地表径流的搬运作用而产生的,它是引起土壤侵蚀最重要的外部驱动因素。降雨侵蚀力因子(R)是评估区域降雨引起土壤分离和搬运的动力指标[15]。章文波等[16]通过大量的对比分析发现基于日雨量计算多年平均降雨侵蚀力的精度最高。因此，本文根据日降雨量数据采用半月逐日雨量模型计算降雨侵蚀力，模型表达式如下：

(2)

(3)

α=21.586β-7.1891

(4)

式中:Rj为第j个半月时段的降雨侵蚀力值[(MJ·mm)/(hm2·h·a)];Pi为半月时段内第i天的侵蚀性日降雨量(mm),模型中要求Pi≥12，否则Pi值计入0;K为该半月时段内的天数；α，β均为模型待定参数。Pd12为日降雨量≥12 mm时的日均降雨量；Py12为日降雨量≥12 mm时的年均降雨量。虎雄岗等[17]通过对5种空间数据插值方法进行对比与分析发现，协同克里金法考虑了地形高程对降雨量的影响,因而更适合贡山县高山峡谷地区降雨数据的空间插值，具有更好的插值效果。因此，本文将各雨量站点收集到的贡山县日降雨数据代入公式(2)计算出各站点的降雨侵蚀力R,然后使用协同克里金插值法得到贡山县1996年、2006年、2012年、2018年4个时期降水侵蚀力R的空间分布图(图1)。

图1 1996-2018年贡山县降雨侵蚀力因子空间分布

2.3.2 土壤可蚀性因子 土壤可侵蚀因子(K)可以反映土壤被降雨侵蚀力剥蚀和搬运作用的敏感程度，也可以评价土壤是否易受侵蚀营力的破坏，是土壤性质中的一个重要因子[18]。土壤可侵蚀因子K值越大，表明土壤更易受到侵蚀[19]。K值的估算方法较多,最常用的有通过天然小区直接测定、作物模型(EPIC模型)、中值粒径法。本文因研究区域较小，且已获取到研究区的砂砾、粉粒、黏粒和有机碳的含量，故采用被广泛用来进行K值计算的Williams等[20]在EPIC模型中的估算方法,该方法由蔡崇法等[21]在小流域进行应用研究，其公式为:

(5)

式中：Wd为土壤砂粒含量(%)；Wi为土壤粉粒含量(%)；Wt为土壤黏粒含量(%)；Wc为土壤有机碳含量(%)，其中：

(6)

通过公式(5)—(6)计算不同土壤类型的K值，其单位为(t·hm2·h)/(MJ·hm2·mm),然后利用GIS空间属性赋值功能给各种土壤类型赋予相应的K值，得到研究区土壤可蚀性因子K的空间分布图(图2A)。

2.3.3 坡长坡度因子 地形因子(LS)包含斜坡长度因子L和坡度因子S，是影响土壤侵蚀强度的重要因素，其中坡度的影响最大。我国西南土石山区山高坡陡、地形复杂，大部分区域存在25°以上的坡地，甚至有些坡耕地的坡度已达到45°[22-23]。因此根据贡山县坡度大于25°的区域占比面积较大的特点，借鉴刘斌涛等[24]提出的适用于西南土石山区土壤流失方程坡度因子计算公式：

(7)

式中：S为坡度因子；θ为坡度(°)。

采用Wischmeier等[25]提出的经验公式估算坡长因子L值，计算公式为：

(8)

式中:λ为水平投影坡长;22.13为标准小区的坡长;m为可变的坡长因子指数,不同的坡度选取不同的m值，其中θ≤1°时m=0.2，1°<θ≤3°时m=0.3，3°<θ≤5°时m=0.4；θ>5°时m=0.5。由公式(7)—(8)得到LS分布图(图2B)。

图2 贡山县土壤可侵蚀因子和坡长坡度因子空间分布

2.3.4 植被覆盖与管理因子 植被覆盖与管理因子(C)是指在一定植被覆盖和管理措施条件下土壤流失量与同等条件下适时翻耕、连续休闲对照地上土壤流失量之比[26],属于无量纲数,其值为0～1。植物根系对土壤有固着作用，因此，植被覆盖度越高的地区，越不容易发生土壤侵蚀[27]。目前主要有3种方法来确定C值：小区试验法、经验法和基于植被覆盖度的遥感数据定量估算法。本文采用蔡崇法等[21]的算法来计算植被覆盖与管理因子C,其计算公式如下：

(9)

(10)

式中：C为植被覆盖因子；f为植被覆盖度；NDVI为归一化植被指数。

2.3.5 水土保持措施因子 水土保持措施因子(P)是指在一定水土保持措施的作用下产生的土壤侵蚀量与对应未采取保持措施的土壤侵蚀量之比，其值为0～1，0代表基本不发生土壤侵蚀；1表示未采取水土保持措施或措施完全失效，土壤侵蚀十分严重。本文基于彭建[28]、张岩[29]等对山区土壤侵蚀研究基础上,参照其他学者的研究成果[30-31],对研究区不同土地利用类型的水土保持措施因子P进行赋值(表1)。在ArcGIS中计算得到研究区CP因子的空间分布图(图3)。

表1 不同土地利用类型的P值

图3 1996-2018年贡山县植被覆盖与水土保持措施因子(CP)空间分布

3 结果与分析

3.1 土壤侵蚀强度的时间变化

运用ArcGIS栅格空间分析中的地图代数功能对上述各因子计算结果进行叠加运算，得到贡山县的土壤侵蚀模数分布图，并根据我国《土壤侵蚀分类分级标准》(SL190—2007)[32]进行重分类制图，得到贡山县1996年、2006年、2012年、2018年土壤侵蚀模数图(图4)。由图4可知，4个时期大部分区域为轻度侵蚀，土壤侵蚀模数为500～2 500 t/(km2·a)，其余部分以微度侵蚀为主，为0～500 t/(km2·a)，侵蚀模数很小，为无明显侵蚀区域;强度侵蚀、极强度侵蚀和剧烈侵蚀占比非常小。通过GIS统计分析，获得1996年、2006年、2012年和2018年贡山县土壤侵蚀模数的平均数值和植被覆盖率(图5)。由图5A可知，研究区1996年、2006年、2012年和2018年年均土壤侵蚀模数分别为1 627.39,750.80，831.82,788.25 t /(km2·a)。总体来说，1996—2018年贡山县年平均土壤侵蚀模数下降趋势明显，土壤侵蚀状况得到很大的改善。

植被覆盖率基本上都是从同一时期的TM影像数据中获得的，因此，可以客观地反映一些问题。由图5A—B可知，1996—2006年，土壤侵蚀模数明显下降而植被覆盖率增加了4.63%，这与云南省2002年全面启动退耕还林等政策有关；2006—2012年土壤侵蚀模数小幅度上升，同时，植被覆盖率有微弱下降的趋势；2012—2018年土壤侵蚀模数下降，植被覆盖率却大幅度增加8.7%，这与云南省自2012年起实施陡坡地生态治理项目以及2014年启动新一轮退耕还林等政策有关。比较图5A—B,发现土壤侵蚀模数与植被覆盖率有关，植被覆盖率越高的地区，越不容易发生土壤侵蚀。

图4 1996-2018年贡山县土壤侵蚀模数

图5 1996-2018年贡山县土壤侵蚀模数和植被覆盖率

对研究区土壤侵蚀模数>500 t/(km2·a)的区域进行面积统计，得到4个时期不同侵蚀强度的侵蚀面积与百分比的变化情况(表2—3)。结果显示：1996—2018年22 a间土壤侵蚀面积减少了25.74%，总体呈现下降趋势，说明土壤侵蚀状况明显趋于好转。并且轻度侵蚀呈明显下降趋势，面积由1996年的2 475.84 km2(占比58.98%)减少到2018年的1 562.92 km2(占比37.23%),降幅达到36.87%，强度侵蚀、极强烈侵蚀和剧烈侵蚀总体上也均呈现下降趋势。

3.2 土壤侵蚀强度的空间动态变化

从4期土壤侵蚀分布图(图4)可知，微度和轻度侵蚀呈现连续块状分布的特征，贡山县每个乡镇均有大面积分布，中度侵蚀主要分布于独龙江乡和茨开镇，强度及以上土壤侵蚀主要发生于怒江和独龙江沿河两岸，对比观察土壤可蚀性因子K的分布图，可发现该区域K值也较大，因此，该区域土壤类型易发生土壤侵蚀，再加上降雨量充足，使得该区域土壤侵蚀较为严重，是水土保持治理的重点区域。

表2 贡山县1996-2018年土壤侵蚀面积变化 km2

表3 贡山县1996-2018年土壤侵蚀面积百分比变化 %

分别对1996—2018年4个时期各侵蚀等级的面积数据进行统计，得到研究区土壤侵蚀强度转移矩阵(表4—6)：3个时段内分别有75.81%，95.69%，67.09%左右的区域侵蚀强度未发生变化，1996—2006年期间有18.63%区域侵蚀减轻，6.11%区域发生恶化，2012—2018年有26.41%的区域向低侵蚀等级转移，6.50%的区域向高侵蚀等级转移，土壤侵蚀发生恶化的面积明显低于侵蚀状况得到改善的区域，说明研究区的水土保持措施得当，贡山县整体的侵蚀状况已逐步得到控制。但是2006—2012年不仅没有区域侵蚀减轻，反而有4.31%区域发生恶化，说明土壤侵蚀有轻微恶化的趋势。从总的变化趋势来看，3个时间段内微度、轻度和中度侵蚀向强度及以上侵蚀等级的转化率非常低，均在1.91%以下，说明研究区近22 a来土壤侵蚀并不严重。

表4 1996年、2006年贡山县土壤侵蚀强度转移矩阵 km2

表5 2006年、2012年贡山县土壤侵蚀强度转移矩阵 km2

表6 2012年、2018年贡山县土壤侵蚀强度转移矩阵 km2

3.3 不同坡度土壤侵蚀强度的空间分布

地形是影响土壤侵蚀的重要因素之一，其中坡度的大小对土壤侵蚀的影响最大[33]。将坡度划分为0°～15°,15°～25°,25°～35°,35°～50°,>50°共5个等级,通过ArcGIS空间叠置分析功能得到不同坡度条件下各土壤侵蚀等级的分布情况(表7)，结果表明：贡山县不同等级的土壤侵蚀在不同坡度范围的分布特征差异明显。4个时期的微度侵蚀均主要分布在15°～25°和25°～35°；轻度侵蚀主要分布在25°～35°和35°～50°；中度侵蚀和强烈侵蚀主要分布在25°～35°和35°～50°；极强烈侵蚀主要分布在35°～50°和>50°；剧烈侵蚀主要分布在25°～35°和35°～50°。研究区22 a来土壤侵蚀强度与坡度大小有明显的正相关关系，整体表现为坡度越大，越容易被剥蚀，土壤侵蚀等级越高，特别是在坡度大于25°时，土壤侵蚀等级呈明显加速趋势，土壤侵蚀面积增加。

表7 贡山县不同年份不同坡度土壤侵蚀面积

4 结 论

(1) 在整个研究期内，贡山县土壤侵蚀模数经历了先明显下降后稍有上升，然后再下降的变化过程，1996—2018年22 a间土壤侵蚀模数下降趋势明显，降幅为51.56%，植被覆盖率上升12.9%，总体上土壤侵蚀状况趋于好转，尤其是强度、极强度和剧烈侵蚀区面积下降明显。这与云南省大力实施退耕还林和生态修复等政策紧密相关。

(2) 贡山县4个时期的土壤侵蚀以微度和轻度侵蚀为主，均占研究区侵蚀总面积89.4%以上，土壤侵蚀较为严重的区域是怒江和独龙江沿河两岸，主要发生强度及以上的土壤侵蚀，今后应将怒江和独龙江沿河两岸作为水土保持治理的重点区域，改善生态环境，控制水土流失。

(3) 通过土壤侵蚀强度等级转移矩阵可知，微度和轻度侵蚀变化率较小，并且高侵蚀等级向低侵蚀等级转换的比率均高于低侵蚀等级向高侵蚀等级转换的比率，说明自2002年以来贡山县全面启动的退耕还林政策实施效果较好，侵蚀区域得到了治理，生态环境有了明显的改善。

(4) 坡度在地形因子中对土壤侵蚀影响最大，坡度越陡，水土流失便越严重，并且当坡度大于25°时，土壤侵蚀等级呈明显的加速趋势，表明在土壤侵蚀过程中存在25°的临界坡度值。坡度在25°～50°的区域侵蚀强度大，侵蚀等级高，因此，今后贡山县应特别加强该区域的水土流失防治措施，这将大大减少土壤侵蚀面积和土壤侵蚀总量。