河南省近十年来土壤侵蚀时空变化分析

2021-03-23黄硕文李健张欣佳邓联文张金萍

农业资源与环境学报 2021年2期

黄硕文，李健*，张欣佳，邓联文，张金萍

（1.郑州大学地球科学与技术学院，郑州 450001；2.郑州大学水利科学与工程学院，郑州 450001）

2019年9月，黄河流域生态保护和高质量发展被列为重大国家战略，水土流失相关研究是生态保护的重要组成部分。河南省处于黄河中下游，是农业大省和粮食主产区，但部分地区在气候、地形等因素的综合作用下极易发生土壤侵蚀，对耕地产生破坏，直接影响粮食作物的产量。为了更好地保护有限的耕地资源，实现区域高质量发展，有必要对河南省土壤侵蚀状况进行定量评价，分析其时空变化特征，探究引发土壤侵蚀的主导因素，并制定相关的保护措施。

对土壤侵蚀进行评价主要有实地观测和建立土壤侵蚀模型两种方法。在实际应用中，由于某些地区的客观条件限制难以进行实地观测，建立土壤侵蚀模型进行侵蚀量预报是常用的方法。根据建立的方法用途和所模拟的侵蚀过程可以将土壤侵蚀模型分为经验模型、物理模型和分布模型[1]。代表性的物理模型有水蚀预报模型WEPP[2]、欧洲土壤侵蚀预报模型EUROSEM[3]、非点源地区流域环境反应模型ANSWERS[4]等，分布模型以SHE 模型[5]最为典型。物理模型和分布模型由于面向过程建模，所需的资料和参数设置过多，难以在特定区域外进行推广。土壤侵蚀经验模型基于大量实测数据进行统计分析，结合土壤侵蚀过程的机理建立，国际上应用较广的主要有通用水土流失方程USLE 和其修正形式RUSLE。RUSLE 对USLE 中各因子含义及计算方法进行了必要的修正，增加了复杂坡度下土壤流失量的计算方法，具有结构简明、计算因子含义明显的特点，且容易与GIS 结合进行分析，在国内外土壤侵蚀定量评价中应用广泛[6-9]。土壤侵蚀过程中同时受到多种因素的影响，因此分析其内在驱动力对于水土保持政策的制定至关重要。地理探测器可用于分析不同自变量对因变量的解释能力，在土壤侵蚀领域的应用也得到了学者们的认可[10-11]。

本研究利用RUSLE模型与GIS技术相结合，对河南省全境土壤侵蚀强度进行定量评价，分析其时空变化情况，并利用地理探测器分析土壤侵蚀影响因素，为河南省水土流失防治重点保护区的划定和生态环境保护工作提供参考。

1 材料与方法

1.1 研究区概况

河南省（31°23′～36°22′N，110°21′～116°39′E）地处中原地区，是国家重要的粮食生产基地（图1），由于季风气候的影响以及跨海河、黄河、淮河、长江四大水系的地理位置，研究区内常常出现局部土壤强侵蚀，对生态环境造成破坏，但近几年河南加强了水土保持措施的制定与实施，水土流失情况得到了一定改善。

图1 河南省行政区划及数字高程模型Figure 1 Administrative division and digital elevation model of Henan Province

1.2 数据来源

本研究所选用的数据为河南省行政区划图、降雨数据（2008、2013、2018年）、ETM+卫星遥感数据、数字高程模型（DEM）数据、地貌类型数据和土壤数据。其中降雨数据为中国气象数据网（http：//data.cma.cn）的中国地面气候资料日值数据集；土壤数据（质地、有机质等）来源于黑河计划数据管理中心（http：//westdc.west⁃gis.ac.cn）的“基于世界土壤数据库的中国土壤数据集（V1.1）”，从中截取了河南省地区数据；DEM 90 m分辨率原始高程数据和ETM+卫星遥感影像来源于地理空间数据云平台（http：//www.gscloud.cn）；地貌形态类型数据来源于中国科学院资源环境数据中心（http：//www.resdc.cn）。将各个因子数据坐标系转化为高斯-克里格投影坐标系，栅格数据分辨率统一为90 m×90 m。

1.3 RUSLE模型

RUSLE模型充分考虑了影响土壤侵蚀的多个因子，对复杂、分布广泛的土壤侵蚀及模型参数变异性等问题均提出较好的解决方案，其公式为：

式中：A为模型计算得到的土壤侵蚀模数，t·hm-2·a-1，本研究乘以100，将单位换算为t·km-2·a-1；R为降雨侵蚀力因子；K为土壤可侵蚀性因子；L为坡长因子；S为坡度因子；C为作物覆盖和管理因子；P为水土保持措施因子[12]。RUSLE 模型通过计算各个因子的指标值对土壤侵蚀强度进行定量评价。

降雨侵蚀力因子（R）：降雨通过对土壤进行冲击引发水力侵蚀，是进行土壤侵蚀预报的重要因子之一。本研究采用基于月平均降雨量和年平均降雨量数据的Wischmeier公式[13]进行降雨侵蚀力的计算，公式为：

式中：R为年降雨侵蚀力，MJ·mm·hm-2·h-1；Pi、P分别为月平均、年平均降雨量，mm。

土壤可蚀性因子（K）：土壤可蚀性因子反映了土壤自身理化性质对于土壤侵蚀的抵抗能力。同等条件下，K值越大，土壤被冲蚀的可能性就越大。本研究根据已有参数采用Williams 等[14]在EPIC 模型中的算法，利用土壤机械组成和土壤有机质进行K值计算。具体计算公式为：

式中：Sn=1-Sd/100；Sd为砂粒含量，%；Si为粉粒含量，%；Cl为黏粒含量，%；C为有机质含量，%。

地形因子（LS）：地形因子直接影响坡面侵蚀的速率，是土壤侵蚀的动力因子。在RUSLE 模型中通常将坡长L因子与坡度S因子合并为LS 地形因子进行计算。国内学者基于不同地貌类型区和径流小区的实测数据对模型中地形因子计算参数进行了修正[15]，本研究采用第四次土壤侵蚀普查算法[16]计算LS因子，公式如下：

表1 不同土地利用类型的水土保持因子（P值）Table 1 Water and soil conservation factors（P values）of different land use types

式中：λ为坡长；θ为坡度；m为坡长系数。

植被覆盖和管理因子（C）：植被覆盖与管理因子（0≤C≤1）反映了植被覆盖对土壤侵蚀的抑制作用，为无量纲数，其值越大则植被发挥减蚀的作用越小。C因子的值通常通过植被覆盖度计算，本研究使用归一化植被指数NDVI 来计算地表植被覆盖度c，计算公式如下：

C因子的计算采用蔡崇法等[17]算法的修正形式，避免了当c值在0～0.1之间时异常值（C＞1）的出现，如式（6）所示：

水土保持因子（P）：水土保持措施因子P（0≤P≤1），其值越小代表水土保持措施越好，土壤侵蚀越不易发生。目前P的赋值主要由实地考察资料结合经验决定，本研究参照了美国农业部手册537 号[18]和覃杰香等[19]的相关研究，对不同土地利用类型区域和不同坡度的耕地分别按表1和表2进行赋值。

1.4 地理探测器

地理探测器由因子探测器、风险探测器、生态探测器和交互作用探测器四部分组成。它通过探测要素的空间分层异质性来揭示其背后驱动力[20]。

因子探测器可用于探测某个因素X对因变量Y的解释程度，其大小可用q值来度量（0≤q≤1），其值越大表示自变量X对因变量Y的解释力越强[20]。以此可以识别出影响土壤侵蚀的主导因素。

表2 不同坡度耕地的水土保持因子（P值）Table 2 Water and soil conservation factors（P values）of cultivated land on different slopes

交互探测器用于识别不同自变量因子之间的交互作用[20]，其交互方式见表3。

表3 自变量对因变量的交互作用方式Table 3 Types of interaction between two covariates

生态探测器可以对两个因子进行比较，判断它们的空间分布影响是否有显著差异[20]。风险探测器判断单个影响因子不同层间属性是否有显著差异，可用于识别土壤侵蚀高风险区。

本研究选取土壤侵蚀强度作为因变量，坡度、海拔、年平均降雨量、植被覆盖度、土地利用类型、地形地貌作为土壤侵蚀的影响因子，并对6 类因子进行离散化处理。多年平均降雨量等间距分成9 类，植被覆盖度分为8 类（＜0.3、0.3～0.4、0.4～0.5、0.5～0.6、0.6～0.7、0.7～0.8、0.8～0.9、0.9～1），坡度分为8 类（＜5°、5°～10°、10°～15°、15°～20°、20°～25°、25°～30°、30°～35°、＞35°），海拔分为4 类（＜500 m、500～1 000 m、1 000～1 500 m、＞1 500 m），土地利用类型分为林地、草地、水域、居民用地、建筑用地、未利用地6 个类别，地形地貌数据使用《中华人民共和国地貌图集（1∶100 万）》中的类别编号。然后以河南省为范围进行随机采样，共选取10 000个样本点，每个样本点赋予以上各类属性值作为地理探测器的输入数据。

图2 河南省2008—2018年降雨侵蚀力因子分布Figure 2 Distribution of rainfall erosivity factors in Henan Province from 2008 to 2018

2 结果与讨论

2.1 RUSLE模型土壤侵蚀分布计算结果

RUSLE 模型中降雨侵蚀力因子、土壤可侵蚀性因子、地形因子、植被覆盖与管理因子、水土保持因子计算结果如图2～图6所示。

将以上各因子结果相乘，得到河南省2008—2018年的土壤侵蚀分布情况（图7）。

图3 河南省土壤可侵蚀性因子分布Figure 3 Distribution of soil erosibility factors in Henan Province

图4 河南省地形因子分布Figure 4 Distribution of topographic factors in Henan Province

2.2 土壤侵蚀时空变化特征

根据《土壤侵蚀分类分级标准》（SL 190—2007），并结合河南省土壤侵蚀量实际分布情况对研究区土壤侵蚀强度进行分级：平均侵蚀模数0～500 t·km-2·a-1划分为微度侵蚀区；平均侵蚀模数500～2 500 t·km-2·a-1划分为轻度侵蚀区；平均侵蚀模数2 500 t·km-2·a-1以上区域统一划分为中度及以上侵蚀区。河南省各类侵蚀区的面积占比见表4。

图7 和表4 显示，河南省2008—2018 年土壤侵蚀分布广、程度低：整体以微度侵蚀为主，其面积占比超过95%，且2008 年以来逐年递增，分布在东部及中部冲积平原区；轻度侵蚀地区主要分布在西部、南部山区及中部风沙平原区，其面积占比呈逐年下降趋势，从图7 可以看出豫西部地区下降尤为显著；中度及以上侵蚀区域集中在西北及北部黄河流域，多年来总体变化不大。

图5 河南省2008—2018年植被覆盖与管理因子分布Figure 5 Distribution of vegetation cover and management factors in Henan Province from 2008 to 2018

图6 河南省2008—2018年水土保持因子分布Figure 6 Distribution of water and soil conservation factors in Henan Province from 2008 to 2018

图7 河南省2008—2018年土壤侵蚀模数分布Figure 7 Distribution of soil erosion modulus in Henan Province from 2008 to 2018

表4 2008—2018年河南省不同土壤侵蚀等级面积比例（%）Table 4 Proportion of area of different soil erosion levels in Henan Province from 2008 to 2018（%）

对各年份土壤侵蚀分级图进行叠加分析得到2008—2013、2013—2018 年两个阶段的土壤侵蚀强度转移矩阵（表5、表6）和转移图（图8）。

根据图表显示，2008—2013年河南省共有2 738.01 km2地区土壤侵蚀强度降低，243.88 km2侵蚀强度增加。侵蚀强度降低主要集中在轻度转化为微度，转化面积达2 324.77 km2。从前文分析可知侵蚀下降主要集中在西部地区，这与河南省近年来在西部山地生态区实施的一系列保护性措施有关。这些措施增加了植被面积，有效抑制了土壤侵蚀。侵蚀强度增加主要集中在微度侵蚀增强为轻度，主要是因为豫西北地区以山地为主，同时受河流侵蚀及耕地分布不合理等自然及人为因素的共同影响，治理难度较大。

2013—2018 年河南省共有946.57 km2地区土壤侵蚀强度减小，843.89 km2地区侵蚀强度加大。与2008—2013 年相比，侵蚀强度明显加强。这可能是由于2018年河南省降雨量明显增加，水力侵蚀增强。

表5 河南省2008—2013年土壤侵蚀强度转移矩阵（km2）Table 5 Transfer matrix of soil erosion intensity in Henan Province from 2008 to 2013（km2）

表6 河南省2013—2018年土壤侵蚀强度转移矩阵（km2）Table 6 Transfer matrix of soil erosion intensity in Henan Province from 2013 to 2018（km2）

2.3 土壤侵蚀影响因素分析

前述研究结果揭示了河南省土壤侵蚀的时空演变特征，借助地理探测器进一步探讨其驱动因子。

图8 土壤侵蚀动态转移分布图Figure 8 Dynamic transfer map of soil erosion

表7 土壤侵蚀影响因子的解释力（q值）Table 7 The q value of influencing factors of soil erosion

地理探测器中因子探测器结果如表7 所示，不同影响因子对土壤侵蚀强度的解释力有明显差异，排序为坡度＞土地利用类型＞地貌类型＞年平均降雨量＞海拔＞植被覆盖度。其中坡度对土壤侵蚀强度的解释力最强，高达39.63%，可以看作是影响河南省土壤侵蚀的主导因子；其次是土地利用类型，解释力为11.27%；植被覆盖度的q值最小，解释力不足1%，单因子作用下对土壤侵蚀的影响并不显著。

根据交互探测器结果（表8），不同影响因子通过交互作用对土壤侵蚀空间分布的解释力远大于单因子。坡度与其他因子的交互作用明显大于剩余因子间的交互作用，其中与土地利用类型的协同作用是显著控制因子，解释力高达57.14%，明显高于坡度的单因子作用。此外，植被覆盖度与其他因子的交互作用大大增强了单因子对土壤侵蚀的解释力，特别是与坡度因子的交互作用，q值增幅极大，说明了在地形起伏较大的地区有必要加强植被建设，通过因子间交互作用抑制土壤侵蚀的发生。

通过风险探测器对不同影响因子各分类的土壤侵蚀平均强度进行统计，识别出高风险区，结果（表9）表明坡度＞35°区域发生土壤侵蚀的风险最高。河南省土壤侵蚀强度随着坡度等级的增加而增大，坡度范围为0～25°时增长缓慢，大于25°时土壤侵蚀强度急剧增加，这部分区域主要分布在河南省西部的三门峡以及南阳市，土地利用类型以耕地为主，植被覆盖度较低，极易引发水土流失。降雨量和地貌类型因子不同分类间平均土壤侵蚀强度无显著差异，植被覆盖度的高风险区分布在0～0.3 区间，且随着植被覆盖度的升高，土壤侵蚀强度明显降低，说明植被对研究区土壤侵蚀有重要影响。当海拔大于1 500 m 时，侵蚀强度比前3 个分类有明显提升，通过与降雨量图层的叠置分析，发现这些区域也对应着强降雨区，两因子间协同作用造成了强烈的侵蚀。

表8 土壤侵蚀影响因子交互作用下的解释力（q值）Table 8 The q values of dominant interactions between soil erosion influencing factors

表9 各影响因子侵蚀高风险区域及平均侵蚀模数Table 9 High risk areas of soil erosion and its mean value of erosion modulus

生态探测器结果表明，降雨和坡度、降雨和土地利用类型、植被覆盖度和土地利用类型这三对因子对土壤侵蚀的空间分布影响具有显著差异，与前面分析相结合，可以推断出坡度、土地利用类型、植被覆盖度和降雨量因子对河南省土壤侵蚀具有主导作用，影响着土壤侵蚀空间分布格局。

总体来说，影响河南省土壤侵蚀的因素多种多样。河南地形复杂，虽大部分属于平原地区，但北部、西部、南部分别与太行山脉、秦岭余脉和大别山脉相接，这些区域多为丘陵和山地，坡度变化幅度大，降雨频繁且集中，两者间的协同作用加速了土壤侵蚀的过程。此外，河南作为农业大省，土地类型多为旱地，林地和草地占比较少，植被可以有效降低雨滴的动能，减少雨水的冲刷作用，避免对土壤的破坏，如果不注重耕作区的造林防护，势必会造成生态系统的破坏。为了推动黄河流域高质量发展，有必要转变土地利用方式：对于坡耕地，应加强退耕还林还草工程建设，提高植被覆盖率；对于暴雨易发地区，应设立监测站进行预警，避免由局部强降雨而引发的土壤侵蚀恶性循环；对于西部黄河流经的黄土区，可考虑建设大型的拦泥蓄水工程，控制泥沙下泄，同时又为耕作区提供水源保证。