土壤侵蚀不仅使土地资源遭到破坏，生态环境恶化，而且给区域社会及生态环境的可持续发展带来潜在威胁^[1-3]；国内外对土壤侵蚀的研究以水蚀研究较多，且多以USLE(universal soil loss equation)模型为评价模型^[4]；近几年风蚀研究逐渐增加，土壤风蚀模型RWEQ (revised wind erosion equation)可以定量评价土壤风蚀量动态变化^[5]，是水土保持生态评价的重要工具，已在新疆、甘肃^[3]、青海^[6]、内蒙古^[7-8]、河北^[1]、天津^[9]等地应用。对于土壤侵蚀的驱动因素研究，国内外多对降雨侵蚀力^{[2, 10-12]}、土壤特征^[13-14]和植被覆盖^[15-16]单个驱动因素进行研究，或者分别对水蚀或风蚀进行驱动因素分析，很少同时对风蚀和水蚀进行综合驱动因素分析。土壤侵蚀和土壤、地形、气象、植被等因子相关^{[1, 15]}，但在区域尺度上，土壤类型和地形特征在短时间尺度内基本保持稳态，体现土壤侵蚀的空间分布特征，而气象因子和植被覆被的动态变化是区域土壤侵蚀时空变化的主导因素^{[1, 11, 15-16]}；因此，分别研究静态、动态因子对土壤侵蚀的影响，能够更全面地探索区域内土壤侵蚀规律。

坝上高原与华北北部山区水土流失严重、水资源短缺，这些潜在威胁着该区及周边地区的生态安全和可持续发展。近几十年来，由于气候变化、植被修复工程及人类活动的驱动^[16-17]，该区生态系统发生了相应变化，而对更能体现区域差异的地理特征，如土壤、地形地貌等的研究欠缺；因此，分析土壤侵蚀的空间差异与时空变化对相应驱动因素的响应，对于目前留存和新出现的问题具有重要意义。选取该区1998—2015年气象数据、NDVI(normalized difference vegetation index)数据集及地理相关数据，应用USLE、RWEQ模型和趋势分析，以求反映土壤侵蚀及其驱动因素的时空特征、变化趋势。研究包括：1)计算1982—2015年研究区土壤侵蚀时空分布特征及变化趋势，分析研究区不同区域的土壤侵蚀程度差异；2)研究土壤侵蚀变化的相对静态因子土壤、地形的空间区域差异特征；3)分析土壤侵蚀变化的相对动态因子气象、植被覆盖度(fraction vegetation cover, FVC)在时间和空间上的动态变化。研究结果可为制订适应综合环境的土地资源保护措施提供科学依据。

1 研究区概况

坝上高原和华北北部山区地处我国农牧交错带中部，北起坝上高原，西起黄土高原北部，向东至华北平原，南至太行山。地理位置位于E 110°47′~119°25′，N 38°49′~42°37之间，面积约16.84万km²。行政区域包括北京、河北、内蒙古、山西、天津的66个县(区、旗)。海拔范围从东南到西北介于-2~2 982 m之间，是我国第2、第3阶地生态交错带的典型多海拔地区。研究区横跨半湿润半干旱气候带，以温带大陆性气候和温带季风气候为主，四季分明，年降水量在327~788 mm之间，年平均气温为-0.8~13.3 ℃。土壤类型以褐土、棕壤、栗钙土、栗褐土和潮土为主。土地利用类型、植被类型受土壤、地形和气候条件的影响，结构空间各异，可分为农田、林地、草地、城市用地、水域和未利用地6种类型。作为一个区域面积较大的统一整体，很难理解时空变异的生态地理条件，因此，依据行政区域的完整性，以及地形、降水、温度、土壤植被类型等物理地理条件的相似性，我们将研究区划分为3个一级区、8个二级区(图 1)。

2 材料与方法 2.1 数据来源与处理

我们在资源环境数据云平台(http://www.resdc.cn)收集了模型所需的植被、地理数据，包括DEM数据、土壤属性数据、植被类型数据、土地利用类型数据集(2000—2015年)、长时间序列(1998—2015)NDVI数据集。利用DEM数据计算坡长坡度因子和地表糙度，土壤属性数据计算土壤可蚀性和土壤结皮值。利用多年土地利用类型数据赋予其相应的P值，以NDVI数据集来计算植被覆盖度。可从中国气象科学数据共享服务网平台(http://data.cma.cn)获得研究区域内及附近的49个气象站的气象数据集，研究区内气象站位置如图 1所示，主要包括降水、温度、相对湿度、风速、太阳辐射和日照时间等气候数据。利用这些日(月)气象数据计算降雨侵蚀力、气象因子，并应用空间插值法将点数据转换为栅格数据。在ArcGIS中应用模型计算参数因子时，所使用数据的坐标系统、空间分辨率均保持一致。

2.2 研究方法

1) 采用USLE模型评估土壤水蚀量^[4]

$ {Q_{\rm{U}}} = RKLSCP。$

(1)

式中：Q_U为土壤水蚀量即侵蚀模数，t/(hm²·a)；R为降雨侵蚀力因子^[18]，MJ·mm/(hm²·h·a)；K为土壤可蚀性因子^[13]，t·hm²·h/(hm²·MJ·mm)；L为坡长因子；S为坡度因子^{[4, 19-20]}；C为植被覆盖与生物措施因子^[17]；P为工程措施因子^[20]；L、S、C和P量纲均为1。

2) 采用RWEQ模型评估土壤风蚀量^{[1, 3, 5-9]}

$ {Q_{\rm{R}}} = \frac{{2Z}}{{{S^2}}}{Q_{\rm{m}}}{{\rm{e}}^{ - {{\left( {\frac{2}{s}} \right)}^2}}}; $

(2)

$ {Q_{\rm{m}}} = 109.8\left( {{\rm{WF}} \cdot {\rm{EF}} \cdot {K^\prime } \cdot {\rm{SCF}} \cdot {\rm{VC}}} \right); $

(3)

$ S = 150.71{\left( {{\rm{WF}} \cdot {\rm{EF}} \cdot {\rm{SCF}} \cdot {K^\prime } \cdot {\rm{VC}}} \right)^{ - 0.3711}}。$

(4)

式中：Q_R为土壤风蚀量即侵蚀模数，t/(hm²·a)；Z为下风向最大风蚀出现距离，m；S为土壤转移量达到最大转移量63.2%的地块长度，m；Q_m为最大土壤转移量，kg/m；WF为气象因子，kg/m；EF为土壤可蚀性成分；K′为地表糙度因子; SCF为土壤结皮因子；VC为植被因子；WF、EF、K′、SCF和VC量纲均为1，取值范围为0~1。

3) 采用一元线性回归方程的斜率^[21]来分析区域尺度上各因子的时空变化，计算方法为

$ {\rm{ Slope }} = \frac{{n\sum\limits_{i = 1}^n i {X_i} - \left( {\sum\limits_{i = 1}^n {{X_i}} } \right)}}{{n\sum\limits_{i = 1}^n {{i^2}} - {{\left( {\sum\limits_{i = 1}^n i } \right)}^2}}}。$

(5)

式中：X_i为第i年因变量; i为年份，1为1998年的第1年，2为1999年的下一年，以此类推，直到n(n=18)，n为1998—2015年研究期间的累计年份。斜率Slope为在一定时间序列期间因变量的变化趋势。一般来说，如果坡度值为正值，因变量的变化呈现总体上的增加或恢复趋势；反之，如果为负值，则表示整体上的下降或退化趋势。

3 结果与分析 3.1 土壤侵蚀特征 3.1.1 空间分布

图 2a为年平均水力侵蚀模数(1998—2015年)的空间分布，整体格局由东南向西北方向呈下降趋势，水力侵蚀多发生在B区，其空间格局随海拔呈梯度式改变，山谷平原的侵蚀强度属于微度，丘陵沟壑区属于轻度，高山区属于中度和强烈^[22]，其中BT区(352.52 t/(km²·a))的高山岩石、丘陵沟壑、农田地带受到更为严重的水力侵蚀；而BY区受到较低强度的水力侵蚀。图 2b为年平均风力侵蚀模数(1998—2015年)的空间分布，整体由西北向东南方向呈下降趋势，全区的侵蚀强度以微度为主^[22]。风力侵蚀多发生在A区，其相对较大的侵蚀模数主要分布在A3区的沙地和稀疏草原，A2及BT区的农地，A区与BY区过渡带的丘陵沟壑区林地；相比之下，东部风蚀程度较低。

3.1.2 时空变化

图 3a为研究区年平均水力侵蚀模数的年际变化(1998—2015年)，平均值为164.99 t/(km²·a)，最大值(2000年)和最小值(2014年)分别为273.94和86.27 t/(km²·a)，整体呈下降趋势。山地丘陵区的水力侵蚀模数大部分呈下降趋势(88.8%)(图 2c)，其中BY区的降低幅度最大(图 3c)；仅有少数地区的水蚀情况加剧(11.2%)，其中BT区加剧程度最大。图 3b为研究区年平均风力侵蚀模数的年际变化(1998—2015年)，平均值为8.49 t/(km²·a)，最大值(2010年)和最小值(2003年)分别为17.49和2.60 t/(km²·a)，整体呈轻微下降趋势；其中在2001—2003年期间侵蚀模数急剧下降，但在2008—2010年期间逐渐上升。在研究期期间，整个研究区50.1%的地区风力侵蚀模数呈现出下降趋势(图 2d)，其中A1区的降低幅度最大(图 3d)；而16.6%的地区情况有所加剧，其中BY区的加剧幅度最大；没有明显变化的占30.3%。

3.2 相对静态因素 3.2.1 地形特征

从各类型区的地形特征统计值和空间分布来看，研究区的高程范围在-2~2 982 m之间(图 4a)，组成以1 200~1 600 m为主(35.4%)，其次为800~1 200 m(29.0%)；其中A1~A3区在各类型区中形成最大高程值区，而C区在各类型区中高程均值最小。研究区的坡度因子值介于0~7.82之间(图 4b)，组成以0~0.59为主(71.3%)，其次为0.96~1.84(12.0%)；其中BT区(0.75)在各类型区中达最大值，而A2和C区在各类型区中值最小。研究区的单元格出口坡长因子范围介于0~2.62之间(图 4c)，组成以0~1.22为主(61.5%)，而＞2.36的区域仅占1.1%(BT区)；坡长因子与坡度因子值域分布状况类似，高山区达最大值，而高原、平原区则较小。从不同类型区的组成来看，A区的高程组成多以1 200~1 600 m为主，B区的高程组成以多海拔梯度组成(图 5a)。全区的坡度因子组成均以0~0.59为主，其中A和C区所占本区域的比例更是超过80%，显示出其地面和缓的状态；B区坡度因子的组成呈多梯度式变化(图 5b)。各区的单元格出口坡长因子组成均以0~1.22为主，所占本区域的比例均超过50%，而地势起伏较大的高山区(BT区)在坡长因子值较高的分级所占比例较高(图 5c)。

3.2.2 土壤特征

研究区地处半干旱半湿润地区，土壤发育程度相对较高，土壤类型相对丰富，以钙层土、半淋溶土、淋溶土、半水成土为主，占到了全区面积的89.33%(图 4d)。由图 4和图 5可知，水蚀中土壤可蚀性K值介于0.005 8~0.0124 t·hm²·h/(hm²·MJ·mm)之间，均值为0.008 3 t·hm²·h/(hm²·MJ·mm)，风蚀中土壤可蚀性EF值介于0.08~0.60之间，均值为0.38，土壤可蚀性较大的土壤类型均为褐土、潮土、栗褐土；而风蚀中土壤结皮因子SCF值介于0.10~0.94之间，均值为0.29，SCF值较大的土壤类型为栗褐土、栗钙土。

从不同区域的土壤类型组成特性来看，研究区东南部(BY和C区)的土壤类型以土壤黏粒和有机质含量较高的褐土、棕壤、潮土为主，由于水分条件相对较好，土壤发育良好，土壤可蚀性明显偏高；西部(A4、A5和BT区)的土壤类型变化较大，这主要与区域土壤的空间变化及过渡有关，其土壤类型以黏粒、粉粒含量、有机质含量较高的褐土、栗褐土为主，土壤更容易被侵蚀，土壤可蚀性也相对较高，SCF明显偏高；北部(A1~A3区)的土壤类型以沙粉粒含量、碳酸钙含量较高的栗钙土为主，栗钙土发育的气候条件明显偏干燥寒冷，地表物质风化过程缓慢，土壤淋溶作用弱，土壤黏粒含量和有机质含量偏低，土壤可蚀性偏低，SCF较高(图 5e)。

3.3 相对动态因素 3.3.1 气候特征

图 6a为年平均降雨侵蚀力空间分布图，其由东南向西北方向呈递减的趋势，变化范围在566.01~3 037.10 MJ·mm/(hm²·h·a)之间；其中C区及BY区形成高值区，而A区及BT区的降雨侵蚀力值均在1 100 MJ·mm/(hm²·h·a)以下。从图 3a中可知，1998—2015年逐年降雨侵蚀力总体呈上升趋势，表明该地区面临发生水蚀的潜在压力有所增加。进一步分析降雨侵蚀力的空间变化趋势(图 3e和6b)，B区和C区降雨侵蚀力的倾向率均为正值，呈现出不同程度的上升趋势，其中，C区的上升趋势最为明显；而A区多为负值，呈下降趋势，仅A5区出现上升趋势。相比较而言，B和C区由降雨引起的土壤侵蚀风险存在或增加的可能性更显著。

图 6c为年平均气象因子空间分布图，其由西北向东南方向呈递减的趋势，变化范围在0.33~20.15 kg/m之间；其中A1和A2区形成高值区，而A5、B和C区气象因子值均在7.00 kg/m以下。从图 3b中可知，1998—2015年逐年气象因子总体呈下降趋势。进一步分析气象因子的空间变化趋势(图 3f和6d)，各地区气象因子的倾向率多为负值，呈现出不同程度的下降趋势，其中，A1区的下降趋势最为明显；而BY区为正值，呈上升趋势。相比较而言，A区由风力引起的土壤侵蚀风险有所减缓，而BY区由风力引起的土壤侵蚀风险存在的可能性有所增加。

3.3.2 植被特征

研究区总体呈西北方区域地势高、干燥少雨、蒸发量大，东南方区域地势相对较低、温暖湿润、降雨相对充沛。由于复杂的地理气候条件，研究区年平均植被覆盖度的空间分布异质性较高(图 6e)，总体上由东南向西北逐步递减(0.01~0.96)；其中东南部山区及平原的植被覆盖度较高(FVC＞0.65)，而西北部高原区的植被覆盖度较低(FVC＜0.50)。研究区的植被类型结构及特征主要包括栽培种植作物、草本植物、乔灌木3类，其中栽培种植作物分布范围广泛，以单作及耐寒栽培(FVC=0.55)、单作及落叶果园(FVC=0.67)、复作(FVC=0.73)为主；草地资源在研究区西北部较丰富，草原主要分布在A区(FVC=0.49)，草甸多分布在A区(FVC=0.65)及B区南部(FVC=0.57)；林地资源在研究区东南部较丰富，针叶林、阔叶林主要分布在BY区(FVC=0.94)，其次为BT区的灌木林(FVC=0.77)。图 3h为研究区年平均植被覆盖度的年际变化(1998—2015年)整体呈显著上升趋势。如图 3g和6f所示，经过植被修复后，大部分区域植被覆盖度呈上升趋势，以A4~A5区的南部、BY区北部、BY区与BT区的相邻地带最为明显，但A区中部的植被覆盖度有所下降。

4 讨论

由USLE和RWEQ模型估计的土壤侵蚀模数是基于自然变化和人类活动驱动的，从另一个角度也可看作相对动态因子(如气候条件、植被覆盖等)和相对静态因子(如地形、土壤特征等)的驱动。土壤侵蚀的空间分布主要由相对静态因子土壤、地形控制，不同的地形形成不同的气候条件，从而形成了不同的土壤类型，孕育生长出不同的植被类型，形成地域性特征。而土壤侵蚀的时空变化主要由相对动态因子植被、气候驱动，气候变化对土壤侵蚀起降低或加剧作用，且具有较强的空间异质性，体现为降雨侵蚀力或风力因子时，作为反映降雨、风速等对土壤侵蚀影响的潜在能力，其时空分布的定量研究是进行土壤侵蚀评价的基础；植被覆盖度的时空变化主要受气候变化和植被恢复工程的共同影响，所形成的新的植被结构和特征，可有效改善生态环境、减少土壤流失，是水土保持的重要因素。

在1998—2015年期间，植被覆盖度呈明显上升趋势，植被改善区域超过总面积的80%，在降雨侵蚀力增加的趋势下，水蚀能以显著趋势降低，说明植被修复工程对缓解土壤侵蚀发挥了重要作用, 这主要是由于植被覆盖可以通过截留降雨和保持水分等方式有效地降低降雨侵蚀力、保护土壤表面、减少泥沙排放，从而控制土壤侵蚀。高原区地势平缓，地表粗糙度相对较低，土壤覆被度低且较易裸露，植被覆盖的增加可以有效固定风沙、提高地表粗糙度，减缓风速，京津风沙源区的植被修复工程对土壤风蚀的作用也不可忽视^{[1, 15-16]}。

土壤侵蚀加剧的区域往往是由于地势起伏、坡度坡长因子相对较高和土壤可蚀性相对较高的原因，在降雨侵蚀力的驱动下极易产生水蚀，如BT区；而BY区北部山地丘陵地带风蚀的加剧不仅仅由于本没有足够的土壤条件、地形特征为植被提供充分的生长环境，再加之逐渐增加的植被覆盖，土壤由此易消耗水分而变得易侵蚀，同时亦易形成温度差，加大了风力，从而形成新的风沙源。综上，在风蚀中土壤性质和风速起着重要作用，对于水蚀则为植被覆盖和降雨侵蚀力，相对动态因子中植被的恢复能够有效的降低水土流失，而对于目前留存和新出现的土壤侵蚀加剧状况，被忽视的相对静态因子给出了有效的解释。

5 结论

1) 研究区东南部以水力侵蚀为主，在1998—2015年期间其侵蚀模数显著下降(Slope=-12.21 t/(km²·a))，研究区88.8%的区域呈下降趋势，主要分布在BY区；研究区西北部以风力侵蚀为主，虽自2009年其侵蚀模数有所上升，但研究区的整体状况呈下降趋势(Slope=-0.026 t/(km²·a))，仅少部分区域(16.6%)有所加剧。

2) 相对静态因素中，土壤、地形是控制土壤侵蚀空间分布的主要因子。研究区A区海拔高、地势平缓，其中A1~A3区坡度因子较低，地面缓和，粗糙度低，土壤类型以栗钙土为主，黏粒和有机质含量低，易产生风力侵蚀；B区地势起伏，海拔、坡度坡长因子呈多梯度分布变化的特征，主要土壤类型以褐土、栗褐土为主，土壤可蚀性相对较高，易产生水力侵蚀。

3) 相对动态因素中，降雨侵蚀力、气象(风力)因子是土壤侵蚀产生的主要驱动因素，具有较强的时空差异特征。降雨侵蚀力主要分布在研究区东部，1998—2015年总体呈上升趋势(Slope=7.01 MJ·mm/(hm²·h·a))，主要增长区域分布在研究区南部，气象因子仅在BY区西北部具有上升趋势，但整体不具有明显趋势；在1998—2015年期间植被覆盖度呈显著上升趋势，年均增速达0.005 8/a，研究区超过80%的区域植被覆盖度增加，仅少部分区域出现加剧状况，主要分布在A区中部。