退耕还林(草)等工程改变了黄土高原的土地利用格局，导致土壤侵蚀发生相应变化。我国学者针对黄土高原小流域土地利用与土壤侵蚀的关系，利用流域与径流小区监测、核素示踪等方法开展大量研究^[1]。然而，流域监测仅能提供小流域总体流失量，无法反映土壤侵蚀的空间分布；径流小区监测与核素示踪可反映不同土地利用的土壤侵蚀状况，但往往无法直接用于分析较大空间尺度的土壤侵蚀分布格局^[2]。

通用土壤流失方程(universal soil loss equation, USLE)、修正通用土壤流失方程(revised universal soil loss equation, RUSLE)和中国土壤流失方程(Chinese soil loss equation, CSLE)等土壤侵蚀模型的大力发展为探讨土地利用变化对小流域乃至更大尺度的流域和区域土壤侵蚀的影响提供了有力的科学支撑^[3-4]。USLE与RUSLE虽可有效反映小流域土壤侵蚀的时空格局，但其多用于缓坡地，在坡度整体较大的黄土高原应用较少。此外，这2个模型所涉及的作物覆盖与管理因子主要基于美国的水土保持措施开展研究与计算，与我国的相关措施无法一一对应。相比之下，Liu等^[5]结合我国的地形特征，改进陡坡地形因子，并充分考虑我国生物、工程和耕作措施对土壤侵蚀的影响而提出的CSLE模型更适合黄土高原地区的土壤侵蚀评价。

六道沟小流域位于毛乌素沙地与黄土丘陵区的过渡地带，是典型的水蚀风蚀交错区，长期以来受到国内外学术界的高度关注。自20世纪80年代以来，众多学者在此开展大量定位观测试验，多集中在土地利用动态变化分析及其对土壤理化性质的影响方面，兼有少数小区和沟道尺度上土地利用与土壤侵蚀关系的探讨^[6-7]。笔者拟结合航片、遥感影像和实地调查，分析六道沟小流域1990—2017年近30年来的土地利用格局变化，并利用CSLE模型对其土壤侵蚀空间分布及演变开展定量评价，以期为黄土高原水蚀风蚀交错区水土流失防治及水土保持规划提供科学依据。

六道沟小流域(E 110°21′~110°23′，N 38°46′~ 38°51′)位于陕西省神木市以西14 km(图 1)，为典型的盖沙黄土丘陵沟壑地貌。小流域面积689.2 hm²，主沟道长4.2 km，自南向北流入窟野河的一级支流三道河。气候类型为中温带半干旱气候，夏秋季节暴雨集中，冬春季节干旱多风。

图 1(Fig. 1)

图 1 六道沟小流域地理位置 Fig. 1 Geographic location of Liudaogou catchment

小流域的地面组成物质主要为第四纪黄土沉积物^[8]，多为壤土和砂壤土，土壤结构疏松，易遭受侵蚀。自20世纪50年代起，六道沟小流域便开展水土流失治理，但方式单一、规模较小，以零星修建水浇地、沟台地为主。直至“七五”期间(1986—1990年)，小流域被纳入黄河中游的重点治理流域，工程措施与生物措施双管齐下，才有了较大进展。2002年，退耕还林(草)工程开始在神木市全面实施，小流域的土地利用类型与水土保持状况进一步发生变化。

选取1990、1995、2002、2010和2017年5个年份分别进行土地利用调查。其中：1990年数据结合野外实地考察和航片判读资料转绘获得^[9]；1995和2002年数据是在1 :5万彩红外航片基础上，经过实地调查获得^[10]；2010年数据以World View-1卫星遥感影像为基础，开展实地调查获得^[7]；2017年数据是在Google Earth影像上开展野外实地调查的结果。

为定量描述小流域的土地利用动态变化，分别计算各土地利用的单一动态度和综合动态度^[11]。单一动态度K_i为某一时期Δt土地利用类型i的面积变化比例(%)：

式中：ΔU_i为Δt土地利用类型i变化的面积，hm²；U_{i, t0}、U_{i, tf}分别为Δt开始时t₀和结束时t_f土地利用i的面积，hm²。综合动态度K_c反映小流域内各土地利用在某一时期Δt的总体变化速率：

式中：ΔU_ij为Δt土地利用i转变为土地利用j的面积，hm²；n为小流域内的土地利用类型总数；S_T为小流域总面积，689.2 hm²。

小流域各时期的土壤侵蚀量采用CSLE^[5]计算：

式中：A为土壤侵蚀模数，t/(hm² ·a)；R为降雨侵蚀力因子，MJ ·mm/(hm² ·h ·a)；K为土壤可蚀性因子，t ·hm² ·h/(hm² ·MJ ·mm)；L、S分别为坡长、坡度因子，量纲为1；B、E、T分别为生物、工程、耕作措施因子，量纲为1。考虑降水的季节变化，R因子基于陕西省神木市1957—2014年的日降雨数据，利用半月降雨侵蚀力方法计算^[12]。K因子源自第一次全国水利普查水土保持专项调查中的土壤可蚀性因子估算成果^[13]，在此基础上利用黄土高原的径流小区资料进行修正，再根据六道沟小流域边界裁切得到对应数据。L因子基于小流域1 :1万地形图，采用Foster等^[14]提出的分段坡长因子及刘宝元等^[15]对坡长指数m的研究成果进行计算：

式中：L_i为第i个栅格的坡长因子；λ_out、λ_in分别为栅格出口及入口的坡长，m；m为坡长指数；θ为坡度，(°)。S因子基于小流域1 :1万地形图，采用McCool等^[16]与Liu等^[17]的公式进行计算：

图 2示出小流域2017年的生物、工程与耕作措施类型及分布。与降水类似，黄土高原的植被覆盖也存在明显的季节变化，二者的耦合关系是土壤侵蚀发生和发展的关键。六道沟小流域的生物措施主要为林草生态恢复与人工造林(图 2(a))，林地、草地和园地的B因子也在半月尺度上进行计算，基础数据源为30 m分辨率的Landsat数据。由于云层覆盖等原因，某些时段的Landsat影像无法提供有效信息，故采用250 m分辨率的MODIS植被指数产品(MOD13Q1)进行插补，获得各年24个半月的盖度图像。结合土地利用数据，进一步计算对应的土壤流失比率，然后以半月降雨侵蚀力占年侵蚀力的比例作为权重，加权计算得到B因子^[12]。旱地、农村居民点、独立工矿用地、水域及水利设施用地和其他土地(裸岩和沙地)直接依次赋值为1、0.025、1、0和0。

图 2(Fig. 2)

图 2 2017年六道沟小流域水土保持措施 Fig. 2 Soil and water conservation practices in Liudaogou catchment in 2017

E、T因子分别参照第一次全国水利普查水力侵蚀普查成果和Guo等^[18]对中国作物覆盖与管理因子的研究成果赋值得到。小流域内的工程措施主要为土坎水平梯田、坡式梯田和鱼鳞坑(图 2(b))，E因子分别赋值为0.084、0.414和0.181；耕作措施主要为反映作物覆盖影响的轮作(图 2(c))，T因子为0.496。

随着退耕还林(草)工程的推进，六道沟小流域的土地利用格局发生了明显变化(图 3)。1990、1995和2002年，小流域皆以草地为主，旱地和灌木林地次之，3种土地利用的平均占比分别为43.3%、28.2%和21.4%。2002年以来，草地进一步扩大；旱地面积急剧减少，主要集中在主沟和支沟沟底水热条件较好的区域。至2010年，草地占比高达60.9%，旱地比例缩减至7.7%。然而，2017年，草地仅占小流域总面积的33.2%，以柠条和沙柳为优势种的灌木林地成为主要的土地利用类型，占比40.8%，占据了大量山坡。

图 3(Fig. 3)

图 3 1990—2017年各时期六道沟小流域土地利用分布 Fig. 3 Spatial distribution of land use from 1990 to 2017 in Liudaogou catchment

按调查年份依次划分出1990—1995、1995—2002、2002—2010、2010—2017年4个时段，分别计算各土地利用的单一动态度。结果表明，近30年来，小流域旱地持续减少，灌木林地、农村居民点和独立工矿用地持续增加(表 1)。其中：旱地以2002—2010年间减幅最大，年平均减少8.3%；灌木林地除2002—2010年间小幅减少外，在其余时段均增加。农村居民点和独立工矿用地的扩张也主要集中在2002—2010和2010—2017年2个时段。由于小流域下游煤矿资源的大规模开采及附近住宅用地的扩张，独立工矿用地年均增幅分别为23.0%和12.7%，农村居民点分别为11.3%和7.9%。有林地面积主要在2010—2017年间大幅增加，年均增幅高达36.5%。

表 1(Tab. 1)

表 1 1990—2017年各时段六道沟小流域各土地利用单一动态度 Tab. 1 Dynamic indices of each land use during different periods from 1990 to 2017 in Liudaogou catchment  %/a 土地利用类型Land use type 1990—1995 1995—2002 2002—2010 2010—2017 旱地Dryland -3.1 -2.9 -8.3 -5.0 果园Orchard 0 -9.6 -12.5 — 有林地Forestland 1.0 0 -0.3 36.5 灌木林地Shrubland 5.8 1.3 -1.2 12.1 其他林地Other woodland -5.5 0 12.7 -8.1 草地Grassland -0.3 1.3 4.1 -6.5 农村居民点Rural residential land -0.8 0 11.3 7.9 独立工矿用地Industrial and mining land — 0 23.0 12.7 水域及水利设施用地Water and water conservancy facilities land 12.7 0 0 -4.3 其他土地Other land -1.4 0 -6.6 -14.3     注：“—”表示无对应时段的动态度结果，2010年果园面积和1990年独立工矿用地面积都为0 hm2。Notes: “—” indicates no results of dynamic indices in corresponding periods, as the areas of orchard in 2010 and industrial and mining land in 1990 are both equal to 0 hm2.

表 1 1990—2017年各时段六道沟小流域各土地利用单一动态度 Tab. 1 Dynamic indices of each land use during different periods from 1990 to 2017 in Liudaogou catchment  %/a

综合考虑各土地利用的变化速率，发现其以2002年为节点呈先缓慢后急剧的变化趋势。2002年之前，仅有水域及水利设施用地在1990—1995年间年均变化率超过10%；2002年之后，果园、有林地、灌木林地、其他林地、农村居民点、独立工矿用地及其他土地的年均变化率在对应的1个或2个时段中超过10%。计算综合土地利用动态度也发现，1990—2002年，土地利用年均变化速率为1.7 %，2002—2017年则增至3.9%。这一趋势与神木市2002年开始实施退耕还林(草)工程相契合。该工程极大地减少了旱地面积，增加了林草地面积，加快了土地利用变化。同时，煤炭资源的大规模开采导致独立工矿用地及农村居民点面积的相应增加。

结合土地利用的总体变化趋势，以2002年为节点，分别计算1990—2002、2002—2017年六道沟小流域的土地利用转移矩阵(表 2和3)。由于流域面积较小，转移矩阵采用一级土地利用分类，即旱地划入耕地，果园划入园地，有林地、灌木林地、其他林地纳入林地，农村居民点与独立工矿用地纳入居民点及工矿用地，草地、水域及水利设施用地、其他土地保持不变。

表 2(Tab. 2)

表 2 1990—2002年六道沟小流域土地利用转移矩阵 Tab. 2 Land use transition matrices from 1990 to 2002 for Liudaogou catchment hm2 耕地 Cropland 园地 Garden land 林地 Woodland 草地 Grassland 居民点及工矿 用地Residential and industrial land 水域及水利设施用地 Water and water conservancy facilities land 其他土地 Other land 1990年总计 Total in 1990 耕地Cropland 137.4 0 31.1 58.8 3.1 0.9 0 231.2 园地Garden land 0.8 0.4 0 0 0 0 0 1.2 林地Woodland 0.7 0 146.7 4.0 1.6 0 0 153.0 草地Grassland 16.7 0 21.9 252.5 1.1 0 0 292.2 居民点及工矿用地 Residential and industrial land 0.8 0 0 0 3.2 0 0 3.9 水域及水利设施用地 Water and water conservancy facilities land 0 0 0 0 0 1.3 0 1.3 其他土地Other land 0 0 0.4 0 0 0 5.9 6.3 2002年总计Total in 2002 156.4 0.4 200.1 315.3 9.0 2.2 5.9 689.2

表 2 1990—2002年六道沟小流域土地利用转移矩阵 Tab. 2 Land use transition matrices from 1990 to 2002 for Liudaogou catchment hm2

表 3(Tab. 3)

表 3 2002—2017年六道沟小流域土地利用转移矩阵 Tab. 3 Land use transition matrices from 2002 to 2017 for Liudaogou catchment  hm2 耕地 Cropland 园地 Garden land 林地 Woodland 草地 Grassland 居民点及工矿用地 Residential and industrial land 水域及水利设施用地 Water and water conservancy facilities land 2002年总计 Total in 2002 耕地Cropland 29.2 4.7 47.0 69.4 5.9 0.1 156.4 园地Garden land 0.2 0 0.1 0.1 0 0 0.4 林地Woodland 0.8 2.3 141.8 47.8 7.4 0 200.1 草地Grassland 3.7 3.6 181.3 107.2 19.2 0.2 315.3 居民点及工矿用地 Residential and industrial land 0.1 0 1.3 1.5 6.1 0 9.0 水域及水利设施用地 Water and water conservancy facilities land 0.1 0 0.6 0.4 0 1.2 2.2 其他土地Other land 0 0 2.7 2.7 0.5 0 5.9 2017年总计Total in 2017 34.2 10.6 374.7 229.0 39.1 1.5 689.2

表 3 2002—2017年六道沟小流域土地利用转移矩阵 Tab. 3 Land use transition matrices from 2002 to 2017 for Liudaogou catchment  hm2

1990—2002年间，小流域共发生土地利用转移141.9 hm²，占总面积的20.6%(表 2)。对应的土地流转主要发生在耕地、林地和草地之间。耕地至草地的净转出面积为42.1 hm²，占1990年耕地面积的18.2%。同时，30.4 hm²的耕地与17.9 hm²的草地净转出为林地。

2002—2017年，土地利用转移更为剧烈，转移总面积增至403.7 hm²，占小流域总面积的58.6%(表 3)。耕地净减少122.2 hm²，主要转出为草地和林地，主要与当地的退耕还林(草)政策有关。据2017年调研，当地人均仅可保留1亩(1 hm²=15亩)耕地，且须同时种植1亩柠条(灌木)。草地面积的减少主要源自其向林地的流转。2002—2017年，133.5 hm²草地净转出为林地，占2002年草地面积的42.3%。此外，随着当地煤炭资源的开发，大量耕地、林地和草地转变为居民点及工矿用地。2002年之前，小流域内小煤矿散布，面积较小；2010年左右，六道沟村办煤矿、蛇圪垯煤矿等整合扩建成腾伟矿业有限公司，占据了原煤矿附近的草地、林地和部分主沟耕地。而居民点面积的增加主要为紧邻煤矿新建的蛇圪垯新村。

利用CSLE计算得到1990、1995、2002、2010和2017年六道沟小流域的土壤侵蚀模数，分别为27.8、25.9、15.7、12.5和7.2 t/(hm² ·a)，呈逐渐降低趋势，与小流域所在窟野河流域的把口站——温家川站的监测数据趋势一致。由于CSLE反映的是多年平均情况^[5]，因此对温家川实测数据进行平均处理，即以上述年份为中点，前后各推及2年计算其5年平均输沙量，2017年只推及之前2年的数据。结果表明，1988—1992、1993—1997、2000—2004、2008—2012和2015—2017年，窟野河的年均输沙模数分别为21.8、16.4、2.1、0.3和0.1 t/(hm² ·a)^[19]，量级基本与六道沟小流域的CSLE模拟结果相当。黄土高原各流域在水土保持治理前的多年平均泥沙输移比约为1，因此输沙模数与侵蚀模数接近；2000年左右大规模开展退耕还林(草)工程以来，泥沙输移比明显<1，输沙模数小于侵蚀模数^[20]，这可能是造成六道沟小流域CSLE模拟结果与窟野河实测数据异同的主要原因。可见，CSLE较好地反映了六道沟小流域的土壤侵蚀强度及趋势。

根据SL 190—2007《土壤侵蚀分类分级标准》^[21]，对小流域各年份的土壤侵蚀强度进行分级(图 4)。1990—2017年，小流域侵蚀皆以微度侵蚀为主，且随时间总体呈扩张趋势(图 4)，主要发生在淤地坝及坡度和缓的塬地和峁顶，面积占比从1990和1995年的约50%逐渐增加至2017年的84.3%。中度侵蚀面积总体减少，从1990和1995年的占比约15 %降低至2017的2.9%。1990和1995年，强烈、极强烈和剧烈侵蚀主要发生在陡坡及支沟与主沟沟壁，3种强度共占小流域总面积的17%左右。2002年后，耕地大面积转变为草地和林地，强烈及以上侵蚀多集中在主沟和支沟沟壁，三者占比总和不超过7%。这一趋势与马瞳宇等^[7]利用淤地坝泥沙淤积量及粒径分布得到的结果一致。

图 4(Fig. 4)

图 4 1990—2017年各时期六道沟小流域土壤侵蚀强度分级分布 Fig. 4 Spatial distribution of soil erosion intensity classes from 1990 to 2017 in Liudaogou catchment

比较各土地利用的土壤侵蚀模数及强度分级(图 5)：独立工矿用地最严重，年均土壤侵蚀模数为45.1~66.0 t/(hm² ·a)；旱地次之，年均土壤侵蚀模数为21.5~35.7 t/(hm² ·a)；草地和灌木林地最弱，分别为3.8~25.8 t/(hm² ·a)和4.1~25.5 t/(hm² ·a)。这与佘冬立^[22]在该小流域利用径流小区监测得到的结果接近。

图 5(Fig. 5)

1990年独立工矿用地面积为0。 The area of industrial and mining land was 0 in 1990. 图 5 1990—2017年六道沟小流域不同土地利用的土壤侵蚀强度分级占比 Fig. 5 Proportion of each soil erosion intensity class for 4 different land uses from 1990 to 2017 in Liudaogou catchment

煤矿开采造成地表裸露、土地变形和塌陷，大量疏松的弃土弃渣堆积于沟道。1990年，小流域独立工矿用地未成规模，调查面积为0。1995和2002年独立工矿用地的土壤侵蚀强度接近，中度、强烈和极强烈侵蚀面积占比总和超过65%。2010和2017年强烈及以上侵蚀面积明显减少，仅占约35%，由于扩张的独立工矿用地主要分布在下游主沟道西侧的缓坡(图 3)，L与S因子均降低，二者乘积相比1995和2002年减少25.9%。

1990—2017年间，旱地均以微度侵蚀为主。1990、1995和2002年，旱地年均土壤侵蚀模数约33 t/(hm² ·a)，主要分布在上游及主沟道东侧陡坡(图 3)，虽有少量梯田，但因年久失修，水土保持作用差。随着科学合理的退耕，2010年，旱地多集中在主沟和支沟沟底水热条件较好的区域，L因子总体降低，平均侵蚀模数降至21.5 t/(hm² ·a)。2017年，小流域上游西南侧草地和灌木林地被重新开垦且无任何水保措施，使旱地平均土壤侵蚀模数增至30.6 t/(hm² ·a)。

草地和灌木林地土壤扰动小，侵蚀程度总体较弱。1990—2017年间，二者皆以微度侵蚀为主，但占比呈增大趋势，从1990和1995年的约42%增至2017年的约90%(图 5(c)和(d))。1990和1995年，草地和灌木林地盖度皆在36%以下，平均侵蚀模数在21 t/(hm² ·a)以上。伴随着退耕还林(草)工程，2002与2010年草地和灌木林地盖度增加，B因子降低约50.0%，相应土壤侵蚀模数均减少至约11 t/(hm² ·a)。2017年植被盖度最高，均超过65%，B因子仅为0.04左右，土壤侵蚀模数分别低至3.8和4.1 t/(hm² ·a)。

1) 1990—2017年，小流域土地利用格局发生明显变化，旱地持续减少，灌木林地、农村居民点和独立工矿用地持续增加。以2002年为节点，各土地利用的变化速率呈先缓慢后急剧趋势。

2) 近30年来小流域土壤侵蚀总体呈减弱趋势。年均土壤侵蚀模数从1990年的27.8 t/(hm² ·a)减少至2017年的7.2 t/(hm² ·a)，减幅74.1%。微度侵蚀面积不断增加，强烈、极强烈和剧烈侵蚀面积不断减少，2017年三者面积占比总和低至2.6%。

3) 各土地利用类型中，独立工矿用地侵蚀最严重，但随时间推移侵蚀模数逐渐减少，主要源于L、S因子的减小。旱地次之，侵蚀强度随随退耕总体递减，但部分陡坡的重新开垦使侵蚀有所增强。草地和灌木林地侵蚀强度最小，随着退耕还林(草)工程的推进，植被盖度增加，侵蚀强度减弱。