基于结构化植被指数的延河流域土壤侵蚀时空动态分析

刘欢欢, 刚成诚, 温仲明,, 陈同德, 刘 悦, 陈 印

(1.西北农林科技大学 草业与草原学院, 陕西 杨凌 712100; 2.西北农林科技大学 水土保持研究所,陕西 杨凌 712100; 3.中国科学院 水利部 水土保持研究所, 陕西 杨凌 712100)

土壤侵蚀是一个全球性的环境问题，不仅降低土壤质量、破坏土地资源，而且淤积河道、加剧洪涝灾害，对农业生产和生态环境构成了极大的威胁，严重制约了社会经济的可持续发展[1-2]。据估算，全球土壤侵蚀的面积约为1.64×107km2，约占陆地总面积的10.95%[1]。水利部2018年全国水土流失动态监测结果显示，2018年全国发生土壤侵蚀的面积约为2.73×106km2，占国土总面积的28.6%[3]。我国每年因土壤侵蚀造成的经济损失高达百亿元以上[1,4]。因此，科学计算流域土壤侵蚀动态变化是深入认识土壤侵蚀发生过程与机理的重要步骤，对于评估区域水土保持措施有效性及合理调整具有重要的理论与现实意义。

土壤侵蚀模型的发展为计算不同时空尺度的土壤侵蚀提供了有力的技术手段[5-7]。目前在我国应用比较广泛的土壤侵蚀模型包括通用土壤流失方程(Universal Soil Loss Equation, USLE)、修正的通用土壤流失方程(Revised Universal Soil Loss Equation, RUSLE)和中国土壤流失方程(Chinese Soil Loss Equation, CSLE)。如王涛等[8]利用RUSLE模型计算了2000—2014年陕北无定河流域土壤侵蚀的时空变化；Teng等则利用RUSLE模型预测了未来青藏高原的土壤侵蚀[9]；USLE被用于计算祁连山[10]、湄公河流域[11]及长株潭城市群等[12]地区的土壤侵蚀和水土流失特征。CSLE也被用于不同县域、省(市)区及流域尺度的土壤侵蚀研究[13-15]。

黄土高原由于其特殊地理位置及气候条件，加之长期以来不合理的人类活动，使其成为中国乃至全球土壤侵蚀最为强烈的地区之一，生态环境极其脆弱[16-17]。自上个世纪末，国家先后实施了“退耕还林还草”等一系列生态工程，黄土高原的生态环境得到了有效的改善，植被覆盖度增加，多年平均入黄泥沙量从16亿t锐减至2～3亿t[18-19]，可见植被在抵抗土壤侵蚀中的重要作用。现有土壤侵蚀模型中一般采用给不同土地利用类型赋值或者基于植被覆盖度来表征植被在水土保持中的作用，即植被覆盖管理因子(C)[6-7,15,20]。但由于植被覆盖度仅能体现植被的投影盖度，即植被的水平结构信息，无法体现植被的层次结构，即垂直结构信息，因此其不能完全真实的反映植被的水土保持功能。据此，雷婉宁和温仲明等在野外大量观测数据的基础上，提出了结构化植被指数(Cs)的概念[21]。与传统的归一化植被指数(NDVI)相比，Cs在高植被覆盖区，能够有效避免植被指数饱和现象，而且能够将不同层次的植被信息相耦合，以表达植被群落更为完整的水土保持功能[21-22]。

为进一步了解植被结构对土壤侵蚀的影响，本研究分别采用Cs和NDVI计算植被覆盖管理因子，在此基础上，利用RUSLE模型模拟2000—2018年延河流域逐年土壤侵蚀量，分析其时空动态特征，比较NDVI和Cs在流域水土流失评价中的可靠性。研究结果一方面为土壤侵蚀模型的植被因子参数优化提供科学依据，另一方面为黄土高原植被结构优化调控和管理提供科学参考。

1 材料与方法

1.1 研究区概况

延河是黄河中游的一级支流，发源于靖边县境内，由西北向东南流经志丹、安塞和延安，在延长县南河沟乡凉水岸附近汇入黄河，全长286.9 km[7]。主要支流包括杏子河，西川、潘龙川和南川等，整个流域面积为7 725 km2，位于陕北黄土高原中部丘陵沟壑区(36°21′—37°19′N, 108°38′—110°29′E)，地势西北高东南低(图1)，地形破碎复杂，黄绵土占流域总面积的85%以上，抗侵蚀能力差[7,23]。

该流域属于温带大陆性半干旱气候区，年均温约为9℃，年降水量500 mm，大部分地区干旱少雨，降雨季节分配不均，夏季多暴雨，降水变化率大。植被类型以次生植被和中、低覆盖度草地为主，天然林集中部分在流域的西南部。常见的植物种从南向北有辽东栎、刺槐、油松。延安至安塞之间为柠条和白羊草，安塞以北为百里香、长芒草等[23-24]。

图1 延河流域地形

1.2 数据来源

本研究所用的数据包括：(1) DEM数据，来源于中国科学院资源环境科学与数据中心，空间分辨率为30 m；(2) 气象数据：流域内及周边气象站数据来源于中国气候数据共享网；(3) 土地利用数据：来源于中国科学院资源环境科学与数据中心，空间分辨率为30 m；(4) 遥感影像数据：2000—2018年延河流域6—9，11—2月的Landsat 7 ETM和Landsat 8 OLI遥感影像(影像编号127034，127035，126035)，来源于美国USGS网站，利用ENVI软件进行缺失条带插补、辐射定标和大气校正等操作，提取遥感指数应用于RUSLE模型中。

1.3 土壤侵蚀模型

本研究采用修正的通用土壤流失方程RUSLE来计算延河流域的土壤侵蚀模数：

A=R·K·LS·C·P

(1)

式中:A为年均土壤侵蚀模数,是指单位面积单位时间和空间的平均土壤流失量〔t/(km2·a)〕;R为降雨侵蚀力因子〔MJ·mm/(km2·h·a)〕;K为土壤可蚀性因子〔t·hm2·h /(km2·MJ·mm)〕;L和S分别为坡长和坡度因子;C为植被覆盖与管理因子;P为水土保持措施因子。

1.3.1 降雨侵蚀力因子R利用逐月雨量估算延河流域的降雨侵蚀力，其计算方法为：

(2)

R=αFβ

(3)

式中:pi为逐月降雨量(mm);p为年均降水量(mm);n为产生侵蚀性降雨(10 mm)的月份数[25];F为修正参数;R为年均降雨侵蚀力〔MJ·mm/(km2·h·a)〕;α和β为模型参数,依据现有研究[7,26],其值分别取0.183 3,1.995 7。将延河流域内及周边气象站的月降水数据利用薄板样条插值方法插值为空间分辨率为30 m的空间图,再利用上述公式计算出降雨侵蚀力的空间分布。

1.3.2 土壤可蚀性因子K利用EPIC模型计算延河流域不同土壤类型的K值[27-28]，其计算公式为：

K={0.2+0.3exp〔0.0256SAN(1-SIL/100)〕}×

(4)

式中:K为土壤可蚀性因子;SAN,SIL和CLA分别是砂粒、粉粒和黏粒的有机碳含量(%),其中SN1=1-SAN/100。

1.3.3 坡度和坡长因子LS在流域尺度，坡度和坡长可以用过DEM提取。本研究采用刘宝元等在黄土高原建立的坡度坡长因子计算方法[29-30]：

L=(λ/22.1)α

(5)

(6)

α=β/(1+β)

(7)

β=(sinθ/0.089)/〔3.0(sinθ)0.8+0.56〕

(8)

式中:L为坡长因子;λ为由DEM提取的坡长(m);22.1为标准小区坡长;α为坡度坡长指数;S为坡度因子；θ为由DEM提取的坡度值。

1.3.4 植被覆盖管理因子C目前，已有大量关于植被指数和C因子的关系式，本研究采用下列公式进行计算：

(9)

式中:α和β的取值分别为2,1。由于水体和建筑用地的NDVI较低,接近于0,但他们的土壤侵蚀强度很小,因此将水体和建筑用地的C值赋值为0[7]。

NDVI的计算公式为：

(10)

式中:NIR和R分别为Landsat 7,8影像的近红外和红光波段。

在计算流域尺度Cs时，土壤调整植被指数(MSAVI)能够有效的去除土壤背景值的影响，用来反映植被绿色覆盖状况。在提取枯枝落叶层信息时，归一化耕作指数(NDTI)与Cs的相关性要高于其他植被黄度指数[21-22]。因此，利用夏季和冬季的遥感影像提取MSAVI和NDTI来综合表征Cs，其公式为：

(11)

(12)

Cs=0.175MSAVI+0.522NDTI

(13)

式中:NDTI为归一化耕作指数;SWIR1和SWIR2分别为Landsat 7,8影像的短波红外波段1,2;MSAVI为修正的土壤调整植被指数。关于Cs的验证过程详见参考文献[21-22]。

1.3.5 水土保持措施因子P依据现有研究[7,20]，采用给不同土地利用类型赋值的方法确定P值。水体、建筑用地、林地和草地一般未采取水土保持措施，因此赋值为1。对于耕地而言，通常坡度越大，水土保持措施的作用越重要，因此，依据表1按照坡度范围对耕地进行赋值。

2 结果与分析

2.1 2000-2018年延河流域土壤侵蚀强度的时间动态

2000—2018年，基于Cs的延河流域平均土壤侵蚀模数为8 354.62 t/(km2·a)，高于基于NDVI得到的6 421.13 t/(km2·a)(图2)。基于Cs和NDVI的土壤侵蚀模数的变化趋势相同，在2000—2006年下降明显，之后趋于平稳。最高值均出现在2001年，分别为13 916.34，15 867.67 t/(km2·a)，之后逐渐下降，在2006年分别下降为4 092.22，5 962.21 t/(km2·a)。2014年土壤侵蚀模数达到最低值，之后又呈现上升的趋势。基于Cs和NDVI的2018年土壤侵蚀模数分别为7 102.62，5 053.78 t/(km2·a)，与2000年相比，分别下降了44.15%和56.41%。

表1 不同坡度范围耕地的P值

依据水利部颁发的《土壤侵蚀分类标准》，将土壤侵蚀强度分为微度侵蚀〔0～500 t/(km2·a)〕、轻度侵蚀〔500～2 500 t/(km2·a)〕、中度侵蚀〔2 500～5 000 t/(km2·a)〕、强烈侵蚀〔5 000～8 000 t/(km2·a)〕、极强烈侵蚀〔800～15 000 t/(km2·a)〕和剧烈侵蚀〔>15 000 t/(km2·a)〕。统计得到2000—2018年延河流域不同土壤侵蚀等级面积百分比(图3)。2000—2018年，基于NDVI的土壤侵蚀强度为强烈、极强烈和剧烈侵蚀的面积分别下降了29.66%，44.48%和62.75%，而基于Cs的土壤侵蚀强度为强烈、极强烈和剧烈侵蚀的面积分别下降了12.17%，28.43%和51.21%。轻度和微度侵蚀的面积都呈现上升的趋势。基于Cs和NDVI的结果表明，在2000—2018年，延河流域的土壤侵蚀强度明显下降。在2018年，基于Cs的强烈、极强烈和剧烈侵蚀的面积占流域总面积的36.84%，而基于NDVI的三者面积占比为27.22%。

图2 2000-2018年延河流域土壤侵蚀模数

图3 基于NDVI(A)和Cs(B)计算的2000-2018年延河流域不同侵蚀强度比例变化

2.2 2000-2018年延河流域土壤侵蚀的空间动态

2000年、2009年、2018年延河流域土壤侵蚀的空间分布见图4。在2000年，土壤侵蚀模数大于5 000 t/(km2·a)的区域广泛分布于延河流域的上中下游地区，轻度和微度侵蚀的区域主要分布在延河沿岸地带及流域的西南部；在2009年，轻度和微度侵蚀的面积增加，主要分布在延河流域的中游地区，而且基于NDVI的增加面积大于基于Cs的面积；在2018年，延河流域中游和下游地区微度和轻度侵蚀的面积进一步增加，土壤侵蚀模数大于5 000 t/(km2·a)的区域主要集中在延河流域的上游地区。

根据2.1中的研究结果，为进一步揭示土壤侵蚀强度变化的空间分布特征，分析了延河流域2000—2007年和2007—2018年土壤侵蚀强度空间转化(图5)。基于NDVI的2000—2007年，土壤侵蚀强度等级不变的面积占流域总面积的36.13%，主要分布在延河的中游和沿岸地区；土壤侵蚀等级降低的区域则占62.42%，其中下降一级的面积占34.17%，其广泛分布于整个流域；在2007—2018年，土壤侵蚀强度等级不变的面积增加到45.20%，侵蚀等级下降的面积降低为35.38%，主要分布在流域的中游和下游，侵蚀强度增加的区域主要分布在流域的上游。在2000—2018年，土壤侵蚀强度下降的面积占流域总面积的62.45%。

在2000—2007年，基于Cs的延河流域内土壤侵蚀强度无变化的区域占流域总面积的54.84%，侵蚀强度下降的面积占42.14%；在2007—2018年，流域内土壤侵蚀强度无变化和下降的面积分别下降到46.44%和36.71%。从整个研究时间范围来看，土壤侵蚀强度以下降趋势为主，面积占53.70%，其主要分布在延河流域的中游和下游地区。在2000—2018年，土壤侵蚀强度下降地区的面积占流域总面积的53.70%，低于基于NDVI的结果。总体而言，基于Cs和NDVI的结果均表明延河流域土壤侵蚀状况得到明显改善。

图4 基于NDVI(A, B, C)和CS(D, E, F)的延河流域2000年、2009年、2018年土壤侵蚀模数变化

图5 基于NDVI(A, B, C)和Cs(D,E,F)的延河流域2000-2007年、2007-2018年、2000-2018年土壤侵蚀强度的变化

3 讨 论

植被和降雨是影响土壤侵蚀动态的两个主要因素。本研究发现，在2000—2018年，延河流域的NDVI和Cs均呈现上升的趋势，增加速率分别为0.012/a和0.014/a，尤其是在2000—2006年，其增加速率均高达0.05/a(图6)。自20世纪90年代初开始，延河流域开展了大规模的退耕还林还草、封山育林等植被恢复工程，造林种草的实施力度增强，截至2006年，延河流域林草植被种植面积高达2 370.05 km2[31-32]。因此，延河流域土壤侵蚀模数下降的主要原因是黄土高原大范围开展的植被恢复措施，这与其他研究的结果一致[7,33]。但由于某些年份的降雨量较大，导致土壤侵蚀模数仍较高，如2011年、2013年、2017年，降雨侵蚀力分别达到了2 371.14，2 950.10，2 378.16 MJ·mm/(km2·h·a)，其对应的土壤侵蚀模数高于6 000 t/(km2·a)。尽管如此，土壤侵蚀强度仍低于2001年、2003年，进一步说明了生态工程措施的有效性。

不同土壤侵蚀模型及因子的不同计算方法均会导致结果的差异，如李天宏等利用RUSLE模型计算延河流域2006年的土壤侵蚀模数约为3 000 t/(km2·a)[7]，而谢红霞等利用CSLE模型计算的同年结果为5 009 t/(km2·a)[23]。本研究基于NDVI和Cs得到的2006年土壤侵蚀模数分别为4 092.22，5 962.21 t/(km2·a)。本研究中，水土保持措施因子P采用的是给不同土地利用类型及不同坡度耕地赋值的方法，而未采用Lufafa等[34]利用坡度计算P的方法，该方法在黄土高原亦被广泛采用[8,24,35]。该方法仅按照坡度计算p值，不区分土地利用类型，即同一坡度的农田和自然植被有相同的水土保持措施因子。若本文采用此方法计算P，进而得到的土壤侵蚀模数将低于现研究结果。此外，本研究未将淤堤坝和梯田等工程措施计算在内，也会对结果造成一定的误差。

图6 2000-2018年延河流域NDVI，Cs及年降雨侵蚀力变化

植被覆盖管理因子的准确估算对于土壤侵蚀预测和植被水土保持效益评价具有重要的意义。乔木林中若无灌木或草本植被，林下植被结构的不完整可导致植被群落的水土保持功能下降，其水土流失可能仍然较为严重。灌草层和枯枝落叶层能够有效的拦蓄径流，减少雨滴动能，增加土壤入渗，从而降低土壤侵蚀，在控制土壤侵蚀过程中发挥着重要作用[36]。传统的NDVI指数由于只考虑了植被盖度，无法反映林下植被结构状况，因此，会存在低估土壤侵蚀强度的现象。而结构化植被指数Cs不仅考虑了植物群落的水平结构特征，还考虑了包括植被乔木层、灌木层、草本层及枯枝落叶层不同垂直结构特征在水土保持作用中的差异，并且较好了克服了NDVI在高植被覆盖区易饱和的现象[21]。基于遥感影像，以MSAVI为绿度指数和NDTI为黄度指数提取的Cs能够在流域尺度更好的反映植被群落的水土保持作用与效益。本文研究结果表明2018年延河流域上游地区的土壤侵蚀模数仍较高，且Cs的土壤侵蚀强度要大于基于NDVI的结果，说明这些地区未来水土保持工作仍需进一步加强或者调整。

4 结 论

(1) 在2000—2018年，基于Cs和NDVI得到的延河流域土壤侵蚀模数分别下降了44.15%和56.41%；基于Cs计算的强烈侵蚀、极强烈和剧烈侵蚀的面积分别下降了12.17%，28.43%和51.21%，而基于NDVI计算的比例分别下降了29.66%，44.48%和62.75%；

(2) 在2000—2007年，土壤侵蚀强度下降的面积占流域总面积的42.14%，广泛分布于整个延河流域；在2007—2018年，侵蚀强度持续下降的面积占流域总面积的36.71%，主要分布在延河的中游和下游；在整个研究时间范围内，延河流域土壤侵蚀强度下降的面积占流域总面积的53.70%；

(3) Cs的变化趋势表明生态工程的实施促进了植被的恢复，有效降低了土壤侵蚀强度。延河流域的上游土壤侵蚀强度仍较高，未来应进一步加强实施水土保持措施。

(4) 基于Cs的植被覆盖管理因子不仅能够反映生态工程实施以来植被的改善，而且由于其包含植被垂直结构特征，能够更好的描述植被在控制土壤侵蚀中的作用，为区域尺度水土保持措施中植被结构优化和调控提供了参考。未来可考虑利用无人机遥感或更高时空分辨率的遥感影像在不同的区域尺度开展进一步研究。