基于GIS和RUSLE的淮河流域土壤侵蚀研究——以信阳市商城县为例

李亚平,卢小平,张 航，路泽忠，王舜瑶

(河南理工大学测绘与国土信息工程学院，焦作 454003)

0 引言

土壤侵蚀是世界三大环境问题之一，会造成土壤有机质流失、泥沙淤积和土地退化等，严重威胁区域的可持续发展[1]。随着地理信息系统(geographic information system,GIS)和遥感(remote sensing,RS)技术的日益成熟，结合RS/GIS技术和修正通用土壤流失方程(revised university soil loss equation,RUSLE)的区域土壤侵蚀定量评估，显著提高了监测及预测效率和空间可视化能力[2,3]，在世界中小流域或区域的土壤流失研究中得到了广泛应用[4]。

坡长L因子是RUSLE模型中的重要影响因子，而实际研究中大尺度的坡长很难通过实测获取，且利用DEM提取坡长的方法存在许多复杂性和不确定性，使L因子成为模型应用受到限制的主要因素[5,6]。Mahala[7]和Yue Li等学者[8]的众多研究都以规则格网的单元格代表坡面的一个坡段，忽略了单元格间的汇流过程，使其形成水文孤岛[9]。Desmet等[10]基于汇水面积的L因子提取方法得到了广泛应用。雍斌等[11]基于多流向算法对累积分配中单位等高线汇流面积进行了改进。然而这些研究都忽略了上坡地表覆盖/植被对汇流的影响，植被和地形的耦合作用必然会导致汇流变化，进而对侵蚀产生影响[12]。秦伟等[13]考虑了上坡土地利用/覆盖，但其实验仅针对黄土高原地区，且缺乏对比实验与验证，对多流向上坡累计汇流算法的介绍也不够详尽。针对现有研究中L因子提取算法不科学和对上坡汇流面积计算不准确等问题，本文改进了RUSLE模型中地形因子的提取算法，提出了考虑地表土地利用/植被对汇流影响的基于多流向的上坡汇流累积的L因子提取算法，并以淮河流域的商城县为例，选择适合该区域的其他因子算法获取土壤侵蚀量及其分布特征，通过对比证明了本文方法的有效性。

1 地形因子提取方法

1.1 通用算法

通用的RUSLE模型中提取地形因子的计算公式为

(1)

L=(λ/22.1)m，

(2)

(3)

式中：L和S分别为坡长和坡度因子;θ为利用数字高程模型(digital elevation model，DEM)数据提取的坡度值，°；m为坡长指数；λ为坡长，m。

通用的RUSLE模型提取的L因子忽略了上坡汇流作用，导致L因子小于真实值，影响土壤侵蚀评定的准确性。本文采用基于上坡汇流的L因子算法对L因子的提取进行了改进。

1.2 本文提出的改进算法

基于上坡汇流的L因子算法首先要计算上坡累积汇流面积。ArcGIS软件中汇流累积计算采用的是单流向算法，而真实的水流方向并不唯一，多流向算法更加符合真实三维地形中的汇流状况[14]。本研究采用Quinn等[15]提出的多流向算法，使用3×3的栅格窗口计算中心像元累积的上坡汇流面积，通过下坡方向加权后的有效等高线长度分配到相邻八像元中下坡方向像元的流量，权重由下坡方向决定，如图1。

图1 多流向径流分配方法Fig.1 Flow partition using the multiple flow direction algorithm

各下坡方向的累积汇流面积公式为

(4)

式中：Ai为第i个相邻下坡单元的累积汇流面积；A为总上坡面积；tanβi为第i个下坡方向的坡角，βi为坡度；Wi为第i个与流向垂直的有效等高线长，为栅格尺寸与权重的积(主方向上权重为0.5，对角方向为0.354)；k为下坡单元个数。

Griffin等[10]在考虑汇流累积及发散后，提出汇流计算L因子的算法，即

(5)

坐标(i，j)处L因子可以表示为

(6)

式中ASi,j-out和ASi,j-in分别为坐标(i，j)出水口和入水口单位汇水面积。

多流向的汇流累积算法只反映了地形对汇流的影响，忽略了地表植被对产流的影响，而事实上植被对降雨再分配和阻滞地表径流等影响显著。虽然RUSLE中的植被覆盖与管理因子C反映了植被覆盖对侵蚀的影响，但C因子仅代表所在单元格的影响，并无法体现上坡土地利用和植被对汇流的作用[13]。因此本文使用不同土地利用类型的汇流面积贡献率tn更加真实地反映上坡植被和地形的作用。首先，根据淮河流域有关研究成果得到不同土地利用类型的平均产流系数[16]，考虑到RUSLE模型是以休闲耕地为基准建立的，tn以裸地的产流系数为标准换算得到(表1)，得到单元上坡实际汇水面积，即

表1 不同土地利用汇流面积贡献率Tab.1 Contribution rate of confluence area for different land uses

(7)

Di,j=R·(sinθi,j+cosθi,j)，

(8)

式中：Ai,j-out和Ai,j-in分别为(i，j)处出水口和入水口的汇水面积；R为栅格宽度；Di,j为所在栅格的有效等值线长度；θi,j为坡向。将(7)代入(6)中得到

(9)

算法实现首先要对DEM进行填洼处理；其次要注意坡长的截断位置，本文选择水体作为坡长L因子的计算终点；最后使用MATLAB编程完成算法。技术路线如图2。

图2 L因子提取流程图Fig.2 Flow chart for extraction of L factor

2 研究区概况及研究方法

2.1 研究区概况

淮河流域大别山区是集山区、贫困区、水土流失严重区于一体的特殊地区。水土流失致使区域生态环境恶化，严重阻碍了经济社会发展[17]。本文选择淮河流域大别山区内的信阳市商城县作为研究区。该区域位于大别山北麓，地理位置E115°06′～115°37′，N31°23′～32°05′，土地总面积为2 109.66 km2，流域面积占全县总面积的60%，海拔高程44.5～1 584 m，地势由南向北倾斜，形成中低山、低山丘陵和丘陵垄岗3大自然区。

2.2 数据收集

研究使用的数据包括：①Landsat卫星2010年空间分辨率为30 m的TM遥感影像；②2010年空间分辨率为2.5 m的SPOT5数据；③空间分辨率为10 m的DEM数据；④淮河流域内及周边共58个气象站点1981—2010年间的月降雨数据及年降雨数据，来源于中国气象局气象中心；⑤中国土壤数据集来源于寒区旱区科学数据中心。

2.3 RUSLE模型

RUSLE模型充分考虑了影响土壤侵蚀的自然要素，具有较强的实用性。其表达式为

A=R·K·LS·C·P，

(10)

式中：A为年平均土壤侵蚀量，t·hm-2·a-1；R为降雨侵蚀力因子，MJ·mm·hm-2·h-1·a-1；K为土壤可蚀性因子，t·hm2·h·hm-2·MJ-1·mm-1;LS为地形因子，L和S分别为坡长和坡度因子；C为植被覆盖与管理因子；P为水土保持措施因子。

2.3.1 降雨侵蚀力因子R

降雨侵蚀力是指由降雨导致土壤侵蚀的潜在能力，是侵蚀的主要驱动力[18]。本文R因子的计算公式为

R=4.1×F-152，

(11)

(12)

式中:F为修正参数；Pj为第j月的降雨量，mm；P为年平均降雨量，mm。将计算所得的各气象站点的年降雨侵蚀力在ArcGIS中进行克吕格空间插值，得到研究区年降雨侵蚀力图。

2.3.2 土壤可蚀性因子K

土壤可蚀性是反映土壤抵抗雨滴击溅或径流冲刷作用的能力。计算公式为

(13)

Sn=1-Sa/100，

(14)

式中：Sa为砂粒的含量，%；Si为粉粒的含量，%；Cl为粘粒的含量，%；C0为有机碳的含量，%。

2.3.3 植被覆盖与管理因子C

植被覆盖与管理因子反映了植被覆盖或作物管理措施对土壤侵蚀量的影响[19]。计算公式为

(15)

(16)

(17)

式中：NDVI为归一化植被指数；ρNIR和ρR分别为近红外波段和红波段反射率；f为植被覆盖度，%；NDVImax和NDVImin分别为NDVI的最大和最小值。

2.3.4 水土保持措施因子P

运用研究区遥感影像在ENVI软件平台支持下建立解译标志，采用人机交互解译，将研究区分为耕地、林地、草地、建设用地、水体和湿地6类，对分类结果进行外业验证与内业修正，直到满足精度要求，得到研究区土地利用分类数据。结合土地利用数据和由DEM生成的坡度信息，参照方广玲等[20]和黄金良等[21]的研究成果对P因子赋值，见表2。

表2 不同土地利用类型P值Tab.2 P values of different land uses

3 结果与分析

3.1 对比验证

依照上述方法得到研究区各因子空间分布图，分别与通用算法和本文改进算法提取的LS因子叠加运算，得到两种算法的土壤侵蚀模数图，并依据水利部颁布的《土壤侵蚀分类分级标准》进行分级，得到研究区土壤侵蚀强度空间分布图(图3)。结果显示，通用算法研究区年平均侵蚀模数为12.2 t·hm-2·a-1，属轻度侵蚀；改进算法的年平均侵蚀模数为28.16 t·hm-2·a-1，属中度侵蚀，与信阳市水土保持监测站所提供数据结果吻合[22]。将2种算法的结果进行叠加分析，可知2种算法的侵蚀等级相同的面积居多，但整体上改进算法得到的结果侵蚀更严重。从2种结果异同区域均匀选取122个特征点结合SPOT5影像进行人工目视解译和野外核查验证，可知本文方法有118处与解译结果一致，验证结果混淆表如表3所示，说明本文提出的算法计算结果更为真实有效。如图4所示区域为采矿用地、坡地，无植被，斑块为亮白色，与周围地物特征差异明显，判定为剧烈侵蚀。本文方法结果为剧烈侵蚀，通用方法结果却为强烈侵蚀。其下坡方向的区域，因其带来的雨水冲刷，侵蚀程度也应较强，本文方法结果为极强烈或剧烈侵蚀，通用方法结果却为中度或强烈侵蚀，不及本文方法准确可靠。

(a)通用算法 (b)本文算法

图3 土壤侵蚀强度空间分布

Fig.3Soilerosiondistributionmap

表3 本文方法试验的混淆矩阵Tab.3 Confusion matrix of the method test in this paper

图4 两种结果对比验证

Fig.4Comparisonandverificationoftworesults

3.2 土壤侵蚀坡度分布特征

统计分析本文算法的土壤侵蚀分布结果，可知研究区总侵蚀面积达905.95 km2，占区域总面积的42.94%，以轻度侵蚀为主。将侵蚀等级与坡度等级叠加统计分析(见表4)可知随着坡度的增加，平均侵蚀模数显著增加，且包含的主要的侵蚀强度等级越高；极强烈和剧烈主要分布在[15°，35°)坡度带，该区域是受人类活动影响较大的易遭受侵蚀的坡耕地和林地的主要分布区域，所以侵蚀明显，应作为治理的重点区域，对于坡耕地最好坡改梯或等高耕种，坡度≥25°的要坚决退耕还林；侵蚀面积并没有随着坡度的增加持续增加，在≥35°区域反而减少了，这与胡世雄等[23]、王星等[24]等的研究结果一致。

表4 不同坡度的土壤侵蚀分布Tab.4 Soil erosion distribution of different slope degrees

3.3 土壤侵蚀土地利用分布特征

将土壤侵蚀结果与土地利用分布图叠加分析(见表5)，可知林地和草地年平均侵蚀模数较高；极强和剧烈侵蚀主要分布在林地。研究区森林覆盖率占49.46%，使得植被覆盖整体较高而侵蚀强度相对较低。但研究区林地多为实施退耕还林后的人工林、顺坡种植且无任何水土保持措施的经济林和疏林地，致使土壤侵蚀严重。草地包含封山育林形成的低覆盖度荒草地和人工种植草地，面积比例不高但侵蚀明显。

表5 不同土地利用的土壤侵蚀分布Tab.5 Soil erosion distribution of different land uses

4 结论

1)本文基于多流向汇流累积算法提取L因子，并考虑了地表土地利用/覆盖的影响加以改进。综合利用RS/GIS技术，基于RUSLE模型将通用算法和改进后算法的土壤侵蚀成果进行对比验证，实验证明了本文算法的有效性。

2)研究区年平均土壤侵蚀模数为28.16 t·hm-2·a-1，属于中度侵蚀，总侵蚀面积达905.95 km2，以轻度侵蚀为主。受地形影响，侵蚀主要分布在中南部地势较陡区域。

3)随着坡度升高，平均侵蚀模数和包含的主要土壤侵蚀等级也升高，极强烈和剧烈侵蚀主要分布在[15°，35°)坡度区域，应重点治理。加强经济林和疏林地改造，保护高密度草地和湿地，防止演化为水土易流失的其他土地利用/覆盖类型是研究区长期的任务。