澜沧江中下游流域土壤侵蚀时空演变特征

陈晨晨, 武 谦, 张占友, 董 卫, 程 瑶, 王 金

(1.河北工程大学 河北省智慧水利重点实验室, 河北 邯郸 056038; 2.河北工程大学 水利水电学院,河北 邯郸 056038; 3.邯郸市漳滏河灌溉供水管理处, 河北 邯郸 056001; 4.北京湜沅科技有限公司, 北京 100070)

土壤侵蚀是指在水力、风力、冻融或重力等外部地质营力作用下，对土壤及其母质进行破坏、剥蚀、搬运和沉积的全部过程[1]。目前，它已成为我国重大的生态环境问题。研究发现，土壤侵蚀会导致土地退化、土壤肥力下降、泥石流和滑坡等自然灾害，也会严重影响河流的自然泄洪能力和梯级水电大坝的建设和利用[2]。土壤侵蚀不仅对土壤功能产生严重影响，而且因其产生的大量泥沙进入河流生态系统中，成为水库沉积物重金属污染的主要来源，同时还导致下游水库蓄水减少和水轮机磨损，给水库水质和水电站的功能带来了很大的负面影响[3-5]。山地丘陵区由于地形起伏、地质构造复杂，自然灾害较为频繁，受强降雨、工程建设及城镇化造成的植被破坏等因素的影响，发生的土壤侵蚀尤为严重[6]，尤其是云南地势起伏，高山丘陵居多，澜沧江中下游流域生态环境脆弱，流域分布众多的断裂带，地震、崩塌、滑坡、泥石流等地质灾害频繁发生，25°以上的坡耕地面积达1 500 km2，水土流失严重。因此澜沧江流域云南段—澜沧江中下游流域土壤侵蚀时空演变特征的研究对于流域内水土保持规划及梯级水电开发与运行具有重要意义。

目前，国内关于地势险峻流域内土壤侵蚀的计算主要是依据通用的土壤流失方程(USLE)模型和改进的土壤流失方程(RUSLE)模型，利用遥感影像数据和地理信息系统(GIS)技术运行模型。吴芳等[2]使用USLE模型和GIS技术探究湄公河流域云南段土壤侵蚀空间特征及其优先治理区。Chuenchum等[7-8]利用RUSLE模型和GIS技术，从空间分布和泥沙淤积两个方面估算澜沧江流域年均土壤侵蚀，评估2030年和2040年未来情景下澜沧江流域土壤侵蚀的空间分布及其产沙趋势。以上学者在澜沧江流域的研究结果证明ULSE和RUSLE模型在澜沧江流域计算土壤侵蚀的可靠性。中国西南地区喀斯特流域地势险峻，高山、深谷居多，邻近澜沧江流域，受亚热带季风气候的影响，Chen等[9]利用现有土壤侵蚀数据验证RUSLE模型在地势险峻流域计算土壤侵蚀的准确性。众多研究表明，RUSLE模型加上GIS技术是用来评价地势险峻流域尺度上土壤侵蚀及其时空分布特征的有效工具[10-12]。关于澜沧江流域土壤侵蚀空间分布的研究目前主要集中在土壤侵蚀强度和空间分布趋势上，很少有研究探索控制澜沧江流域土壤侵蚀的主要因素。因此，本次研究计划采用RUSLE-GIS工具，同时引入一种基于统计学的机器学习算法——随机森林算法，朱青等[13]在赣江上游流域使用随机森林算法系统地分析土壤侵蚀因子在不同子流域间的重要程度，该算法已被应用于准确掌握影响流域土壤侵蚀的主要因素，对于提升水土保持工作的精准化具有现实指导意义。本次研究计划在探明澜沧江中下游流域土壤侵蚀时空分布特征的基础之上，对降雨侵蚀因子、地形因子、植被覆盖管理因子和土壤可蚀性因子进行重要性分析，以期针对澜沧江中下游流域土壤侵蚀问题提出合理的水土防治措施和对策，为流域管理者实施针对性的水土保持措施提供参考依据。并为接下来探究土壤侵蚀过程与梯级水库重金属迁移转化机制的影响提供数据支撑。

1 研究区概况

澜沧江发源于中国青海省唐古拉山脉北部，东临中缅国际边界，北临中老国际边界。本文以澜沧江中下游流域为研究区域。研究区域地形由南到北、由中低山到宽谷、由高山到深谷，地势起伏较大，海拔为478～4 209 m，加之西南季风的影响使得流域内气温和降雨量呈现自下游向上游递减的趋势。但是局部地区会因为地形影响，出现山体迎风坡海拔越高，降水量越多的“地形雨”；山体背风坡海拔越低，降雨量越少的“雨影区”现象[14]。本次研究区域雨影区较多，经纬度在98°—101°E，22°—27°N。研究区域面积为78 200 km2，约占澜沧江流域面积的50%。

2 研究方法与数据来源

2.1 改进的土壤流失方程(RUSLE)

RUSLE模型是一种经验的土壤侵蚀模型，已被公认为计算土地平均土壤侵蚀风险的标准方法。该模型与GIS和遥感技术相结合，是目前最流行的土壤侵蚀模数估算模型[7-8]。在ArcGIS栅格计算器中将各因子相乘就可以得到澜沧江中下游流域段的土壤侵蚀模数：

A=R×K×LS×C×P

(1)

式中：A为土壤侵蚀模数[t/(hm2·a)]；R为降雨侵蚀因子[(MJ·mm)/(hm2·h·a)]；K为土壤可蚀性因子[(t·hm2·h)/(MJ·hm2·mm)]；LS为地形因子；C为植被覆盖管理因子；P为水土保持措施因子；LS，C和P均为无量纲单位。

2.1.1 降雨侵蚀因子(R)R因子代表由降雨引起的潜在土壤侵蚀能力，是评价流域内土壤侵蚀状况的主要动力指标。本文采用Wischmeier月尺度计算公式[15]，公式(2)近几年被广泛应用于计算澜沧江流域R因子：

(2)

式中:Pi为月降雨量(mm)。

2.1.2 土壤可蚀性因子(K)K因子是指在某一特定土壤类型的单位上，每一个侵蚀指数单位所测量的土壤流失率，它受土壤的各种物理、化学和矿物学性质的影响。本文在考虑土壤质地和有机碳的基础上，利用公式(3)[16]计算K因子：

(3)

式中:OM为有机质,土壤中有机质含量可以用土壤中总的有机碳换算得到,我国目前沿用的“Van Bemmelen,1.724”,土壤有机质=土壤有机碳×1.724[17];M为原始粒径组分(Clay,Silt,Sand)的乘积;S和P分别为土壤的结构和种类。

2.1.3 地形因子(LS) 在RUSLE模型中，地形对土壤侵蚀的影响是用LS因子来表示的，LS因子是坡长因子(L)和坡度因子(S)的乘积。本研究计划采用公式(4)—(5)[18]计算LS因子：

(4)

(5)

式中:λ为坡长(m)，是从坡面径流的起点开始,垂直于等高线沿坡向下直到坡度减缓到足以发生沉积的水平距离,根据坡长和坡度间的负相关关系[19],来进一步确定LS因子的大小(表1)。θ为坡度;m为坡长指数。

表1 坡长-坡度关系

2.1.4 土地覆盖管理因子(C)C因子是指种植模式和管理措施对土壤侵蚀产生的影响，也是指在特定条件下种植的土地与平整、连续休耕造成的土壤侵蚀损失之比。本研究采用归一化差异植被指数(NDVI)进行尺度变换来近似C因子[18,20]，相关公式(6)—(7)如下：

(6)

(7)

式中：AVHRR2为近红外波段影像图；AVHRR1为红波段影像图。研究发现[20]α，β值可分别取2，1。

2.1.5 水土保持措施因子(P)P因子代表实施土壤侵蚀保护措施过程中的成效，并被定义为特定耕作方式下的土壤流失与上下坡耕作方式下相应损失的比值。已有研究中列举大量P值取值范围[7,13,21]，根据不同的土地利用类型，将澜沧江中下游流域P因子进行赋值，见表2。

表2 澜沧江中下游流域段P因子值

利用GIS技术载入改进的土壤流失方程，计算出澜沧江中下游流域的R，K，LS，C和P因子，部分结果见图1。

图1 澜沧江中下游流域土壤侵蚀因子分布

2.2 随机森林算法(Random Forest)

随机森林算法(RF)是Breiman[22]在2001年提出的一种组成式的有监督学习算法，通过对样本单元和变量进行抽样，生成大量决策树，从而对决策树依次进行分类，来建立评价指标的重要性所需要的分类树。变量的相对重要性是通过划分变量节点异质性的减少量并对所有分类树进行平均得到的。节点异质性是由简单二元分类树规则的基尼系数来定义。RF对多元共线性不敏感，其结果对缺失和非平衡的数据有较好的容忍度，且无需对数据进行归一化处理。相对重要性评估是RF算法的一个重要特征，土壤侵蚀因子的相对重要程度越大，表明对土壤侵蚀的影响越大[13,23]。本研究基于R软件中的RandomForest程序包，利用回归算法来计算土壤侵蚀因子的相对重要程度值。

2.3 数据来源

为了提供探究澜沧江中下游流域土壤侵蚀与梯级深大水库重金属迁移转化机制影响的支撑数据，根据澜沧江深大水库(小湾和糯扎渡水库)投产建成时间选择本次研究年份。数据年份为2005年、2010年、2015年，主要数据来源分别为：澜沧江中下游流域20个雨量站逐日降雨量和逐年降雨量数据，分别来源于中国气象科学数据共享服务网(https:∥data.cma.cn/)和中国科学院资源环境数据中心(https:∥www.resdc.cn/)；土壤数据来源于中国科学院资源环境数据中心和世界土壤数据库(HWSD)中的中国土壤数据集；数字高程模型(DEM)来源于地理空间数据云(https:∥ www.gscloud.cn/)，30 m空间分辨率；归一化植被指数(NDVI)来源于地理空间数据云中Landsat 4-5，7，8遥感影像，下载范围按照澜沧江流域经纬度(北纬21°30′—32°40′，东经94°40′—101°50′)下载。该影像提供7个波段的影像图，我们选取近红外和红波段影像来计算NDVI；土地利用类型，来源于中国科学院资源环境数据中心提供的中国土地利用遥感监测1 km栅格数据。以上遥感数据均经过统一化处理，利用ArcGIS裁剪提取研究区域，定义投影为WGS_1984_UTM_Zone_47N，数据重采样像元大小为30 m×30 m。

3 结果与分析

3.1 澜沧江中下游流域土壤侵蚀时间变化特征

本文选取澜沧江中下游流域的土壤侵蚀数据资料进行分析。运用ArcGIS软件，载入RUSLE模型计算出澜沧江中下游流域2005—2015年平均土壤侵蚀量(图2)。结果显示2005—2015年澜沧江中下游流域多年平均土壤侵蚀模数为2.52×103t/(km2·a)，初入中风险侵蚀级别。由于该研究区域土壤侵蚀方向的研究较少，缺少验证本次研究结果可靠性的数据，且本次研究与众多学者在澜沧江流域平均土壤侵蚀模数结果相比[2,7-8]，存在一定差异，但是整体空间分布趋势基本一致。2005年土壤侵蚀模数的变化范围在0～1.89万 t/(km2·a)，平均土壤侵蚀模数为2.37×103t/(km2·a)；2010年土壤侵蚀模数的变化范围在0～0.99万 t/(km2·a)，平均土壤侵蚀模数为2.51×103t/(km2·a)，相比于2005年平均土壤侵蚀模数增加了约5.98%；2015年土壤侵蚀模数的变化范围在0～1.58万 t/(km2·a)，平均土壤侵蚀模数为2.67×103t/(km2·a)，相比于2010年平均土壤侵蚀模数增加了约5.68%。2005—2015年澜沧江中下游流域平均土壤侵蚀模数呈现缓慢增加的趋势。

图2 研究区域及子流域土壤侵蚀时空分布特征

本研究按照目前沿用的土壤侵蚀风险标准，见表3，对澜沧江中下游土壤侵蚀风险进行划分，并计算相应风险区域的侵蚀面积。从表3中可以看出，澜沧江中下游流域2005—2015年较低风险土壤侵蚀面积占比分别为8.55%，7.72%，6.56%，低风险区域面积占比分别为48.37%，40.38%，35.91%，较低和低风险区域呈现逐个时间段减小的趋势；中风险土壤侵蚀面积占比分别为30.09%，42.82%，49.37%，随着时间的变化呈现递增的趋势；较高和高风险侵蚀区域多年变化趋势较小，但是始终呈现递减的趋势。通过澜沧江中下游流域土壤侵蚀整体时间变化规律，我们可以发现，2005—2010年土壤侵蚀空间分布并不均匀，具有很强的空间异质性，但是自2010年以后，澜沧江中下游流域土壤侵蚀空间分布的异质性变得很弱。研究区域土壤侵蚀强度变化趋势由低风险和高风险向中风险侵蚀强度转移，并且中风险侵蚀范围在2010年以后大幅度增加。

表3 土壤侵蚀风险划分标准与侵蚀面积

利用ArcGIS的区域统计功能计算出澜沧江中下游流域土壤侵蚀风险等级面积转移矩阵，见表4，2005—2010年流域内侵蚀等级转移矩阵主要集中在低和中风险，约占侵蚀面积的86%。2010—2015年流域内侵蚀等级转移矩阵也主要集中在低和中风险，约占侵蚀面积的89%。但是与2005—2010年相比，2010—2015年低风险面积转移矩阵减小了约11%，中风险面积转移矩阵增加了约17%，这表明2005—2015年，研究区域土壤侵蚀风险呈现一种良性发展的趋势。

表4 土壤侵蚀风险等级面积转移矩阵 %

3.2 澜沧江中下游流域土壤侵蚀空间变化特征

本次研究利用QGIS(Quantum GIS)软件载入SWAT(Soil and Water Assessment Tool)模型，根据研究区域澜沧江中下游流域地形和水系分布进行子流域划分，经过水文计算后共得10个子流域S1—S10，子流域S1—S5，S6—S10分别位于澜沧江中游和下游流域。根据土壤侵蚀风险标准，利用ArcGIS对中下游子流域的平均土壤侵蚀模数进行划分，得到子流域侵蚀风险的分布特征与侵蚀模数统计结果(图2—3)。

从图2—3中可以看出，中下游子流域平均土壤侵蚀模数范围为1 208.4～3 870.4 t/(km2·a)，基本处于较低风险以上和中风险侵蚀以下的区域。2005—2015年研究区域中游子流域平均土壤侵蚀模数差异性较大，其中2005年、2010年基本处于低风险侵蚀级别，2015年属于中风险侵蚀级别，且2010年中游子流域S1—S3的平均土壤侵蚀模数是2005年的一倍左右，平均土壤侵蚀模数随时间呈现递增的趋势。2005—2015年研究区域下游子流域土壤侵蚀模数差异性较小，且侵蚀风险基本处于中风险侵蚀级别，平均土壤侵蚀模数随时间呈现波动递减的趋势。2005年中下游子流域S1，S2，S8和S10从低风险侵蚀级别到2010年变成了中风险侵蚀级别，下游子流域S9从中风险变成了低风险侵蚀级别。2010年中游子流域S3—S5从低风险侵蚀级别到2015年变成了中风险侵蚀级别。由表5可以看出，2005—2010年中游子流域S1—S2，下游子流域S8—S10平均高程增加，平均降雨量和植被覆盖率减小；下游子流域S9平均高程增加，平均降雨量和植被覆盖率大幅度减小。中下游子流域S1和S9平均海拔均在2 000 m以上，平均降雨量在1 000 mm以下；中下游子流域S2和S10平均海拔在1 500 m以下，平均降雨量在1 500 mm以上。2010—2015年中游子流域S3—S5平均高程和降雨量增加，S3，S4平均植被覆盖率增加、平均海拔在1 500 m以下，平均降雨量接近1 500 mm，S5减小。中游子流域S2—S4和下游子流域S6与S10平均降雨量较高。研究表明澜沧江中下游流域在西南季风的影响下，迎风坡地形较多且相对较高，潮湿气流在爬升过程中会产生大量的降雨，易形成以“地形雨”为主的降雨区域；下游子流域S5和S9的平均降雨量较低，背风坡较多且相对较低，易形成以“雨影区”为主的少量降雨区域。从整体上来看，澜沧江中下游流域受西南季风和地形的影响，降水呈现一种自下游向上游递减的趋势，局部受地形影响存在一定差异，且2005—2015年该流域植被覆盖率逐渐下降。因此本次研究发现澜沧江中下游流域地势的起伏、降雨量空间分布的不一、植被覆盖率的逐年降低是影响该流域土壤侵蚀时空分布特征的主要原因。

图3 子流域2005-2015年平均土壤侵蚀模数统计

表5 子流域地理环境状况

3.3 土壤侵蚀因子重要性分析

本次研究基于R语言中的随机森林模型，计算出澜沧江中下游流域2005—2015年土壤侵蚀因子指标的Gini减小平均值，然后将其转换为各个侵蚀因子指标的相对重要程度并绘制成条形柱状图，结果见图4。由图4可以发现，影响澜沧江中下游流域土壤侵蚀的主要侵蚀因子是地形因子(LS)和植被覆盖管理因子(C)，两者的重要程度之和在60%以上。基于联合信息熵方法，吴芳等[2]发现土地利用类型、坡度和海拔是影响湄公河流域土壤侵蚀的主要因素；陈龙等[24]研究指出澜沧江流域土壤侵蚀与坡度和海拔密切相关；姚华荣等[25]发现澜沧江流域土地利用变化对土壤侵蚀变化具有强烈影响。以上研究结果与本次利用随机森林算法研究结果基本一致。2010年植被覆盖管理因子(C)的相对重要程度降低，地形因子(LS)对澜沧江中下游土壤侵蚀影响程度增大，相对重要程度高达45%，这与2010年流域海拔的增加以及植被覆盖率的减小有着很大关系。植被覆盖管理因子(C)在2005—2015年都在很大程度上影响着澜沧江中下游流域的土壤侵蚀，这也体现了流域内整体的相似性[26]。此外水土保持措施因子(P)、降雨侵蚀因子(R)和土壤可蚀性因子(K)对土壤侵蚀的重要程度偏低，均未超过20%，表明这3种因子在整体上的差异性相对较小，各子流域间的变化起伏不大。通过随机森林计算发现植被覆盖管理因子(C)和地形因子(C)是澜沧江中下游流域土壤侵蚀治理应重点关注对象。因此，流域相关管理部门在治理土壤侵蚀的过程中，应调整森林中乔木树冠遮蔽地面的程度“郁闭度”，因为达到一定郁闭度的林草植被有保护土壤不被侵蚀的作用，还能增加土壤的肥沃程度。尽管郁闭度越高，保持水土的能力越强，但是李宗勋等[27]发现马尾松林郁闭度过高或过低不利于水土保持。因此对于减缓澜沧江中下游流域土壤侵蚀的不同种类的森林郁闭度还有待进一步研究。同时应禁止基本建设时对土地不合理的利用、毁林毁草、陡坡开荒，这种严重破坏了地面植被和稳定地形的行为。实施适地适树措施，大力开展防护林和封山育林，对旱地实施坡耕地退耕改梯田平种，耕作、轮作和培肥为一体化的水保措施，以降低土壤侵蚀强度和土壤侵蚀量。在出现大量“地形雨”和“雨影区”的子流域引入遥感分析技术，以便探测哪些地形会引起迎风/背风坡效应[28]，来进一步改善降雨分布不均匀导致的土壤侵蚀加剧状况。

图4 澜沧江中下游流域土壤侵蚀因子的相对重要程度

4 结 论

本文通过调查澜沧江中下游流域2005—2015年的地形资料、降雨资料、土壤类型、土地利用类型和植被覆盖率，利用ArcGIS载入改进的土壤流失方程模型，估算出2005年、2010年、2015年的平均土壤侵蚀量，研究发现流域2005—2015年多年平均土壤侵蚀属于中度侵蚀风险级别，中、低风险侵蚀面积高达80%以上，随着时间的变化，该流域土壤侵蚀空间分布特征呈现高、低侵蚀风险范围收缩，中风险侵蚀范围增大的趋势。同时引入当前热门机器算法之随机森林算法研究流域内土壤侵蚀因子的重要程度，进一步精准识别影响澜沧江中下游流域土壤侵蚀的主要因素是植被覆盖管理因子(C)和地形因子(LS)，地形因子(LS)对流域土壤侵蚀的影响程度最大，植被覆盖管理因子(C)次之。因此相关水土保持部门在制定水保措施时，应考虑引入遥感技术(RS)和地理信息系统(GIS)来实时监测较高风险侵蚀地区的地形和植被覆盖情况，一方面实施合理的工程施工进行微地形生态修复，来进一步改善由于地形引起的“地形雨”和“雨影区”现象，有利于减小流域内土壤侵蚀时空分布的差异性；另一方面改善土地利用类型的使用情况，注重植被覆盖度的提高，寻找合适的森林“郁闭度”，禁止人类活动对森林生态系统造成进一步破坏。为流域管理者针对性治理澜沧江中下游流域土壤侵蚀提供科学的依据。同时本次研究结果为探索澜沧江中下游流域土壤侵蚀时空分布与梯级水库重金属时空分布的关系提供了基础条件。