2000—2020年黄河流域土壤侵蚀及其驱动因素分析

尹作堂, 常 军

(山东师范大学 地理与环境学院, 山东 济南 250358)

土壤侵蚀是指土壤及土壤母质等成分在水、风力和冻融等外力作用下被破坏、运输和沉积的过程,这是近年来粮食安全和生态安全的主要威胁[1,2]。因此,探索土壤侵蚀和驱动因素的时空特征,并采取有效措施保护土壤资源和修复生态环境至关重要[3,4]。当前,越来越多的研究开始关注土壤侵蚀的空间异质性和尺度效应,如考虑行政边界变化对土壤侵蚀影响因素造成的不确定性,以便将研究结果应用于水土保持规划与治理[1,3,5,6]。如Guo等[3]分析了城市尺度、县级尺度与乡镇尺度上的京津冀地区土壤侵蚀空间分布特征与主要控制因素的变化情况;Liu等[6]分析了县级尺度下,珠江三角洲水土保持与城市化水平、平均高程和归一化植被指数(NDVI)地理加权回归系数的空间分布。这些研究提出了更有针对性的水土保持治理与土地利用政策。

基于修正通用土壤流失方程(RUSLE)[7]进行土壤侵蚀分析是当前常用的研究方法之一。在中国,RUSLE模型已在青藏高原[8,9]、黄土高原[10-13]、喀斯特地貌[14,15]和其他地区[16-18]得到了广泛应用,并取得了良好的评估效果。地理探测器是探测地理现象空间分异性,揭示地理现象与因变量关系的一种空间统计学方法[19],因其能够量化不同驱动因素及其交互作用对土壤侵蚀的实质性贡献,在最近几年被广泛应用于土壤侵蚀驱动因素的研究[3-5,15,20]。但相关学者在使用地理探测器进行土壤侵蚀驱动因素分析时,主要依靠王劲峰等[19]提出的方法及先验经验进行人为设定,忽视了数据离散化方法与分类数可能会对探测结果造成的影响[21,22],故本研究选用参数最优地理探测器(optimal parameters-based geographical detector, OPGD)[22]对数据离散方法与间断数量进行最优化设置,以提高地理探测精度。

作为土壤侵蚀最严重的区域之一,1990—2015年黄河年均土壤流失量2.0×109t[23],造成土地资源退化、土壤养分流失、河道淤积,给黄河流域的生态保护和高质量发展造成巨大阻力。近年来,大量学者对黄河流域及相关地区开展了土壤侵蚀及其驱动因素分析[4,23-25],但关于黄河流域土壤侵蚀驱动因素的空间分布还有待深入研究。鉴于此,本文基于RUSLE模型与热点分析,探讨2000年、2010年和2020年黄河流域及县域尺度下的土壤侵蚀状况,并基于参数最优地理探测器评估黄河流域与县域尺度下的土壤侵蚀驱动因素的解释力,分析县域尺度下土壤侵蚀驱动因素的时空分布,以期为黄河流域水土保持及高质量发展提供有效建议。

1 研究区概况、方法及数据

1.1 研究区

黄河流域地处95°02′ E ～ 119°43′ E,31°28′N～ 41°33′ N之间,整体地势西高东低,流域内海拔高差达6 200 m,地形起伏差异明显,自西向东形成横跨青藏高原、河套平原、鄂尔多斯高原、黄土高原和黄淮海平原的三级阶梯(图1)。区域内气候差异明显,季节差异大,且降水时空分布不均,夏季降水量约占总降水量的70%,气温日较差和年较差较大;土壤包括高山土、干旱土、半淋溶土等多种类型,自然植被包括高寒草甸、草原、落叶林等。流域内地理特征复杂,土壤侵蚀以水力侵蚀营力为主,同时存在风力、冻融和重力等侵蚀营力[26]。近年来,黄河流域因水资源短缺、生态脆弱和水沙关系不协调等问题受到广泛关注。

图1 研究区概况图Fig.1 Overview map of the study area 注:此图基于国家自然资源部标准地图服务网站审图号 为GS2020(4619)的标准地图制作,底图无修改。

1.2 研究方法

1.2.1修正通用土壤流失方程

修正通用土壤流失方程(RUSLE)为:

A=R·K·LS·C·P

(1)

式中:A为年土壤侵蚀模数(t/(hm2·a));R为降雨侵蚀力((MJ·mm)/(hm2·h·a));K为土壤可蚀性因子((t·h)/(MJ·mm));LS为坡度因子S与坡长因子L(无量纲);C为地表植被覆盖管理因子(无量纲);P为水土保持措施因子(无量纲)。

其中,基于章文波等[27]提出的日雨量估算侵蚀力方法,计算流域各气象站点年降雨侵蚀力,通过Kriging空间插值处理,得到2000年、2010年和2020年黄河流域降雨侵蚀力因子;基于泛第三极(20国)土壤可蚀性因子(K)数据集[28],并在陈同德[26]、李天宏[29]、张科利[30,31]等人研究基础上做出修正,K因子取值范围为0.017 2～0.086 9 (t·h)/(MJ·mm);基于泛第三极(20国)坡度坡长因子(LS)数据集[32]提取黄河流域LS因子数据集,LS因子取值范围为0～75.59;C因子以蔡崇法等[33]提出的基于植被覆盖度的C因子估算方法进行估算;基于RS与GIS提取法,参考文献[4]、[23]、[34]并结合研究区实际情况,根据土地利用类型(LUCC)数据与坡度数据对P因子进行赋值(表1)。

表1 不同土地利用类型的P值Tab.1 P values of different land use types

1.2.2热点分析

热点分析是一种空间自相关方法,可用于识别土壤侵蚀在空间上的高值(热点)与低值(冷点)的聚类情况[35]。基于ArcGIS Pro热点分析(Getis-OrdGi*)工具,统计识别具有统计显著性的土壤侵蚀热点与冷点,计算公式为:

(2)

(3)

式中:Gi*为空间关联指数;E(Gi*)为Gi*(d)的期望值;Var(Gi*)为变异系数;Wij(d)为空间权重矩阵;Xj为第j级土壤侵蚀强度等级的面积;n为土壤侵蚀强度等级数;Z(Gi*)为Gi*的标准化Z值。标准化Z值为正且值大,表示土壤侵蚀高值的空间聚类(热点),Z值为负且值小,则表示土壤侵蚀低值的空间聚类(冷点)。

1.2.3参数最优地理探测器(OPGD)

地理探测器是探测驱动因素的一种统计学方法,包括分异及因子探测、交互作用探测、风险区探测与生态探测[19]。本研究使用参数最优地理探测器R语言程序包GD[22]来分析单个因子与因子交互作用对土壤侵蚀的解释力,其中探测结果q值表征解释力的大小。具体设置如下:以土壤侵蚀模数作为因变量,以LUCC(X1)、地貌类型(X2)、土壤类型(X3)、年平均降雨量(X4)、坡度(X5)、海拔(X6)、植被覆盖度(X7)、人口密度(X8)为驱动因素,选用五种连续变量离散化算法(自然断点法、标准差法、分位数法、等间隔法和几何间隔法),间断数量设置为3～8类,进行土壤侵蚀驱动因素分析。

1.3 数据来源

主要数据来源:①黄河流域矢量面数据、地貌类型数据(1 km×1 km)、土壤类型数据(1 km×1 km)和土地利用类型数据(1 km×1 km)来源于中国科学院资源环境科学与数据中心(https://www.resdc.cn/);②ASTER GDEM 数据(30 m×30 m)来源于中国科学院地理空间数据云平台(http://www.gscloud.cn/);③日降雨资料(站点数:113)来源于中国气象数据网(https://data.cma.cn/);④NDVI数据(250 m×250 m)是在美国国家航空航天局提供的MOD13Q1数据产品基础上,采用最大值合成法生成;⑤人口密度数据(1 km×1 km)来源于WorldPop project官方网站(https://www.worldpop.org/)。其中栅格数据在经过精度分析、纠偏等处理后,重采样为1 km分辨率。

2 结果分析

2.1 土壤侵蚀特征分析

依据水利部颁布的《土壤侵蚀分类分级标准》(SL 190—2007)[36],将土壤侵蚀强度分为微度(0 ～1 000 t/(hm2·a))、轻度(1 000 ～ 2 500 t/(hm2·a))、中度(2 500 ～ 5 000 t/(hm2·a))、强烈(5 000 ～ 8 000 t/(hm2·a))、极强烈(8 000 ～ 15 000 t/(hm2·a))和剧烈(> 15 000 t/(hm2·a))六个级别。图2中(a)～(c)分别2000年、2010年和2020年黄河流域以栅格像元为单元划分的土壤侵蚀强度图,图2中(d)～(f)分别为2000年、2010年和2020年黄河流域以县域为单元划分的土壤侵蚀强度图。

图2 土壤侵蚀强度图Fig.2 Map of soil erosion intensity 注:此图基于国家自然资源部标准地图服务网站审图号为GS2020(4619)的标准地图制作,底图无修改。

2000—2020年,黄河流域土壤侵蚀状况总体呈现好转趋势,土壤侵蚀强度减弱区域约1.9×105km2,占流域总面积的24%,土壤侵蚀强度增强区域约6.5×104km2,占流域总面积的8%。2000年、2010年和2020年,黄河流域平均土壤侵蚀模数分别为1 994、1 860和1 384 t/(km2·a),水土流失面积分别为2.80×105、2.62×105和2.14×105km2。侵蚀强度较高的区域集中于黄土高原地区,并且呈现东北—西南的条带状走向。但在此期间,河套平原北部的阴山山脉地区,土壤侵蚀强度具有显著的增强趋势。

由图2(d)～(f)可知,县域尺度下不同区县土壤侵蚀的空间差异性明显。2000年,土壤侵蚀强度呈现极强烈的县分别为陕西省子洲县、子长县、延川县、清涧县、延长县,山西省石楼县和甘肃省会宁县,各县平均土壤侵蚀模数分别为8 572、10 160、10 812、11 449、9 502、9 056和8 651 t/(hm2·a);2010年,土壤侵蚀强度呈现极强烈的县为甘肃省白银区、平川区、甘谷县、环县和庆城县,各县平均土壤侵蚀模数分别为8 907、8 622、9 104、9 873和8 045 t/(hm2·a);2020年,土壤侵蚀强度呈现极强烈的县为宁夏回族自治区大武口区、惠农区和内蒙古自治区乌拉特后旗(流域内部分区域),平均土壤侵蚀模数分别为14 575、8 329和13 288 t/(hm2·a)。县域尺度下,土壤侵蚀热点分析结果(图3)与土壤侵蚀强度呈现相似的时空特征,侵蚀聚集于黄土高原地区,宁夏回族自治区大武口区、惠农区和内蒙古自治区乌拉特后旗逐渐成为新的侵蚀聚集中心。

2.2 侵蚀驱动因素分析

土壤侵蚀受自然因素(降水、坡度、土壤类型等)和人类活动(LUCC、人口密度、植被覆盖度)的共同影响,因此确定自然因素与人类活动对土壤侵蚀的解释力对制定后续的水土保持措施尤为重要。黄河流域土壤侵蚀驱动因素地理探测结果如图4所示。

总体来说,各因子对土壤侵蚀的解释力均有限,其中植被覆盖度因子对流域土壤侵蚀的解释力最强(2000年q值最大,为0.160 5),海拔、土壤类型与年平均降雨次之,另外,两个人为因素(LUCC和人口密度)的解释力均很有限。因子交互作用对土壤侵蚀的解释力略强于单一因子,但最多也仅在2000年坡度与植被覆盖度交互作用时,解释了约32%的土壤侵蚀。植被覆盖度作为主要的控制因素,与其他因素交互作用时q值相对较大。在时间上,土壤侵蚀因子探测的结果总体呈减小趋势,植被覆盖度q值从0.160 5减小至0.125 5,人口密度q值从0.096 9减小至0.009 9,LUCC的q值从0.064 5减小至0.047 0。因子交互作用探测的结果变化趋势相似,但植被覆盖度与地貌类型交互时的q值呈增加趋势,植被覆盖度与LUCC交互时的q值呈先增加后减少的趋势。驱动因素的变化趋势表明,黄河流域土壤侵蚀状况在好转的同时,驱动机理变得更加复杂。

图3 热点分析结果图Fig.3 Map of hot spot analysis result 注:此图基于国家自然资源部标准地图服务网站审图号 为GS2020(4619)的标准地图制作,底图无修改。

在县域尺度下,对各县土壤侵蚀驱动因素因子探测q值求平均,结果如表2所示。在以县域为单元划分黄河流域进行土壤侵蚀驱动因素分析后,土壤侵蚀的驱动机理更加明显,平均q值大大超过以流域为单元开展因子探测的q值,且LUCC这一人为因素成为第二重要的驱动因素。这表明在一定程度内,自然因素对土壤侵蚀的影响能力随尺度减小而逐渐减小,而人为因素对土壤侵蚀的影响随尺度减小而增大,人为因素在局部区域发挥了更大作用。这一发现与Guo等[3]、Li等[5]的结论相似。

图4 黄河流域地理探测结果图Fig.4 Map of geographical exploration results in the Yellow River Basin

表2 因子探测结果平均值

县域尺度下,各县土壤侵蚀驱动因素因子探测结果表明(图5中,空值主要因地理探测结果中的sig大于0.05,无法通过显著性检验造成),黄河流域各县土壤侵蚀驱动因素的空间分异明显,但总体上无明显分布规律。在空间上,LUCC、地貌类型、土壤类型、降水、海拔和坡度在宁夏回族自治区北部部分区县(平罗县、贺兰县、西夏县等)对土壤侵蚀的解释力较强;地貌类型与海拔在内蒙古自治区部分区县(土默特右旗、土默特左旗、东河区和赛罕区)与山西省部分区县(离石区、中阳县和汾阳市等)对土壤侵蚀的解释力较强;人口密度在绝大多数区县对土壤侵蚀的解释力有限。此外,植被覆盖度与其他因素相比,展现出不同的空间分布格局,其在山西省、陕西省和甘肃省对土壤侵蚀的解释力较强,而这些区县的土壤侵蚀强度减弱最为明显,表明植被恢复能有效减少土壤侵蚀。

除植被覆盖度外,其他因素的q值在时间上总体呈现减小趋势,这与流域土壤侵蚀驱动因素的时间变化特征相似。为探究植被覆盖度因子q值增强的原因,将植被覆盖度因子探测结果与植被覆盖度数据进行叠加分析,发现2000—2020年,黄河流域植被覆盖度从49%上升至59%,而植被覆盖度因子在植被覆盖度较高的区县对土壤侵蚀的解释力更强。又因在RUSLE模型中,当植被覆盖度大于78.3%时C因子被赋值为0,说明当研究区植被覆盖度高于特定值时,植被控制土壤侵蚀的益处趋于稳定,这与Teng等[37]的研究结论相似。

图5 各县因子探测结果图Fig.5 Map of factor detection results for each county 注:此图基于国家自然资源部标准地图服务网站审图号为GS2020(4619)的标准地图制作,底图无修改。

为进一步探究新的侵蚀聚集中心形成的驱动因素,对宁夏回族自治区大武口区、惠农区和内蒙古自治区乌拉特后旗三个区域的因子探测结果进行分析。由表3可知,2020年,新的侵蚀热点的形成主要受坡度、地貌类型和LUCC影响。分析LUCC、地貌与遥感影像数据可知,三个区域的LUCC类型均以草地为主,大武口区与惠农区地处贺兰山丘陵地区,乌拉特后旗地处阴山山脉区域,均为土壤侵蚀易发区,加之人类活动的影响,极易造成土壤侵蚀的发生。Li等[4]的研究发现,黄河上游净土壤侵蚀速率的空间格局主要由坡度控制,这也从一定程度上印证了所得结论。

表3 大武口区、惠农区和乌拉特后旗因子探测结果Tab.3 Factor detection results of Dawukou, Huinong and Urad Rear Banner

3 讨 论

由于不同学者采用的数据源不同,而且RUSLE模型的因子计算方法也不尽相同,因此,各文献对黄河流域或相关区域土壤侵蚀模数的估算也存在较大差异。此外,考虑到黄河流域土壤侵蚀情况较为复杂,从水利部《中国河流泥沙公报》仅能获知黄河部分水文站点的输沙模数,而在泥沙输移比未知时无法估算侵蚀模数,故对比了当前研究与其他研究来验证本文结果的准确性(表4)。当前研究结果与其他研究结果具有较好的线性拟合(R2=0.96),但模型精度仍有待进一步提高。

表4 其他研究的土壤侵蚀模数Tab.4 Soil erosion modulus obtained from previous studies

综上,黄河流域土壤侵蚀时空特征明显,水土流失严重的区域主要位于黄土高原地区,流域土壤侵蚀状况有所好转,但近年来阴山山脉与贺兰山区的水土流失较为严重。各因子对黄河流域土壤侵蚀的解释力均有限,植被覆盖度始终是黄河流域土壤侵蚀的重要控制因素,该结论与周璐红[44]、贾磊[45]等人的分析结果一致。黄河流域各县土壤侵蚀驱动因素的空间分异明显,但空间分布规律不显著。

不同季节的土壤侵蚀及其主要控制因素不同,Li等[4]研究发现,黄河上游6～8月土壤水蚀最为强烈,且NDVI在夏季对土壤侵蚀的解释力更强,降水在春秋季节、河源地区对土壤侵蚀的解释力更强。在后期研究中,还应探讨不同季节、不同地貌特征(梯田等)下的黄河流域土壤侵蚀及驱动因素的时空变化,以期为黄河流域水土保持精细化管理提供帮助。

4 结 论

1) 2000年、2010年和2020年黄河流域平均土壤侵蚀模数分别为1 994、1 860和1 384 t/(km2·a),侵蚀热点呈东北—西南走向,并呈带状集中于黄土高原地区。在坡度、地貌类型与LUCC的影响下,阴山西部与贺兰山区逐渐成为新的侵蚀热点区域。

2) 植被覆盖度为黄河流域土壤侵蚀最重要的控制因素。在县域尺度下,LUCC成为仅次于植被覆盖度的影响因素,人为因素在局部区域发挥了更大作用。各县土壤侵蚀驱动因素差异性特征明显,但除植被覆盖度在高植被覆盖区域解释力较强外,其他驱动因素无明显空间分布特征。

3) 2000—2020年,黄河流域土壤侵蚀状况好转,但侵蚀驱动因素的解释力呈现减弱趋势,土壤侵蚀驱动机理更加复杂。