1 引言

土壤侵蚀已成为全球性的重大生态环境问题之一^[1]。定量地研究评价土壤侵蚀的时间和空间分异规律,明确其发生强度,区域类型,对有效地实施水土保持工作和改善土地利用方式具有重要科学价值。为此,国内外众多学者做了许多杰出的研究。基于景观生态学的“尺度—格局—过程”原理,建立了多尺度土壤侵蚀^[2]与生态安全评价指数^[3]。研究了不同土地利用方式^[4,5]、不同土壤属性^[6]、及不同植被类型景观^[7]条件下的水土流失效应。利用RULSE模型及GIS技术研究了土壤侵蚀的时空变化^[8,9,10]。利用WEPP模型研究了坡度^[11]、坡长^[12]及耕作制度^[13]对土壤侵蚀的影响,并对WEPP模型在不同地区的适用性进行评价^{[14,15,16,17]}。

但是,喀斯特关键带是个特殊的生态系统^[18],传统经典方法在喀斯特关键带并不适用。喀斯特关键带的特殊性体现在以下几各方面,首先不同岩性下垫面产生侵蚀性径流的降水临界值存在差异。在非喀斯特区地表径流系数一般在20%以上,甚至达到30%以上^[19]。喀斯特区的径流系数介于0.01%~14.72%^{[20,21,22,23]}。喀斯特区的地表径流系数远低于非喀斯特区,将非喀斯特地区的降雨侵蚀力算法直接运用于喀斯特地区是不合适的,需要对其加以改进。其次是碳酸盐岩成土物质的先天不足使得喀斯特地区的成土速度十分缓慢,喀斯特山区的土层通常仅有几厘米到几十厘米厚^[24],其实际的土壤存有量远低于理论侵蚀量,土壤侵蚀被高估。第三是喀斯特洼地因有一定量缓坡耕地的存在,通常被当作侵蚀区,但从区域整体上来看,喀斯特洼地是区域物质与水流的汇聚区,在土壤侵蚀过程中取到汇流与物质累积的作用。因此,从局部上来看,喀斯特洼地坡面是土壤侵蚀的源,但从整个洼地系统来看,其实质却是土壤侵蚀的汇。

有鉴于此,本文以世界上喀斯特最为集中连片的中国西南为研究对象,通过对土壤侵蚀的多次改进,更加精准高效地评估喀斯特关键带土壤侵蚀,为长江与珠江的上游生态屏障建设提供科学参考,为中国西南喀斯特石漠化的治理提供理论依据。

2 研究区与数据来源

2.1 研究区概况

本文所涉及的喀斯特关键带在行政区划上包含云南省、贵州省、四川省、重庆市、湖北省、湖南省、广西区、广东省8个省区（20°13'N~34°19'N、97°21'E~117°19'E）,总面积约193.82万km²（图1）,约占中国国土面积的1/5。地跨中国地貌的三级阶梯,地势西高东低,包含川西高原、云贵高原、横断山脉、四川盆地、洞庭湖平原、两广丘陵多个一级地貌单元,地貌条件复杂,土壤水蚀作用强烈^[25]。同时是世界三大喀斯特集中连片分布区之一。气候兼具低纬气候、季风气候、山原气候的特点,年降水量500~2500 mm之间,降水主要集中于4-9月,受地形、气候因素的制约,降水由东南向西北递减,空间差异较大。地带性植被为亚热带、热带常绿阔叶林。土壤类型丰富,地带性土壤以红壤和黄壤为主,非地带性土壤以石灰土为主。区域生态环境脆弱,是中国水土流失风险最严重的区域之一。

显示原图| 下载原图ZIP| 生成PPT

图1   研究区位置及高程

Fig. 1   Location and elevation of the study area

2.2 数据来源

研究采用的基础数据包括：研究区数字高程模型（DEM）,国家级气象站点日降水数据,土地利用数据,土壤类型及土壤质地数据,土壤有机质数据,NDVI月尺度MODIS中国合成产品数据（表1）。

3 研究方法

3.1 降雨侵蚀力优化和改进

降雨侵蚀力的计算方法采取降雨动能E和最大30 min降雨强度的乘积来表示时,即EI₃₀法,称为R值经典算法。Richardson等在研究中发现EI₃₀与降雨量P之间存在幂指函数关系^[26],形如EI₃₀ = aP^b,其中P为日降水量,a、b为待定参数;Elsenbeer等的研究也证实了此关系的存在,并对该方法进行了验证^[27]。CREAMS模型中降雨侵蚀力的计算便采取了幂指函数的基本形式,系数a = 1.03,b = 1.51。胡续礼等通过对CREAMS模型、Richardson模型、郭新波修正模型的验证对比,表明CREAMS模型有着更好的稳定性^[28],能达到更高的精度要求。在为喀斯特^{[23, 29]}与非喀斯特^[30]区域限定有效降水侵蚀阈值的基础上,本文选取CREAMS日降水模型进行降雨侵蚀力计算：

Rnk=1.03Pi1.51(Pi≥d12)（1）

Rk=1.03Pi1.51(Pi≥d30)（2）

式中：R_nk为非喀斯特地区的降雨侵蚀力（MJ mm ha^-1 h^-1 yr^-1）;R_k为喀斯特地区的降雨侵蚀力（MJ mm ha^-1 h^-1 yr ^-1）;P_i为日降水量（mm）;d₁₂为非喀斯特区域侵蚀性降雨标准,即日降水量≥12 mm;d₃₀为喀斯特区域侵蚀性降雨标准,即日降水量≥ 30 mm。

3.2 岩石风化速率

喀斯特关键带的土壤来源于母岩风化^[31],其土壤允许流失量的大小取决于成土速率的高低^[32,33]。参照白晓永等^[34]提出的喀斯特地区不同岩石类型及其组合方案的成土速率计算公式：

W=v×Q×ρ×M+N×(1-M)（3）

式中：W为岩石风化成土速率（t ha^-1 yr ^-1）;v为碳酸盐岩溶蚀速率（mm a^-1）,换算为m³ ha^-1 yr^-1;Q为酸不溶物含量（%）;ρ为碳酸盐岩容重（tm^-3）;M为碳酸盐岩含量（%）;N为非碳酸盐岩成土速率（t ha^-1yr ^-1）。

以碳酸盐岩的泥质含量作为岩石组合类型的划分指标^[32],将碳酸盐岩地区的岩石组合类型分为三大类（表2）。曹建华等^[35,36]的监测记录表明,碳酸盐岩的溶蚀速率介于30~89.7 mm/ka之间,白云岩平均溶蚀速度为55.27 mm/ka,石灰岩的平均溶蚀速度为96.72 mm/ka。根据陈晓平^[37]的研究结果,石灰岩方解石密度为2.75 t/m³,白云石密度2.86 t/m³,碎屑岩的成土速率以5 t ha^-1 yr^-1计算^[38]。

3.3 喀斯特洼地空间分布信息的提取与计算

由于DEM的误差以及一些真实地形的存在,使得DEM表面存在着一些凹陷区域。首先对原始DEM数据进行填洼,得到无洼地的DEM。通过计算中心网格与邻域网格的最大距离权落差来确定水流的流向,在ArcGIS中通过Hydrology工具集中的Flow Direction工具进行提取。其次是利用Hydrology工具集中的Sink工具,通过输入水流方向数据可计算得到洼地空间分布,Watershed工具则用于计算洼地的贡献区域。Zonal工具集中的Zonal Statistic工具用于计算每个洼地所形成的贡献区域的最低高程,Zonal Fill工具则用于计算每个洼地贡献区域出水口高程^[39]。依据洼地贡献区域的最低高程和出水口高程,运用Raster Calculator工具计算出区域洼地深度。

3.4 RUSLE模型

本文采用目前应用最为广泛的修正土壤流失模型（RULSE）进行研究区土壤侵蚀的量化估算,方程式如下：

A=R×K×L×S×C×P（4）

式中：A为土壤侵蚀模数（t ha^-1 yr^-1）;R为降雨侵蚀力因子（MJ mm ha^-1 h^-1 yr^-1）;K为土壤可蚀性因子（t ha h ha^-1MJ^-1mm^-1）;L为坡长因子;S为坡度因子;C为地表植被覆盖与管理因子;P为水土保持措施因子。

降雨侵蚀力（R）的改进计算详见3.1节。

土壤可蚀性因子（K）本文选择Williams等^[40]提出的EPIC模型进行计算,公式为：

K=0.1317{0.2+0.3exp[-0.0256SAN1-SIL100]×SILCLA+SIL0.3×[1-0.25CC+exp(3.72-2.95C)]×[1-0.7SNSN+exp(22.9SN-5.51)]}（5）

式中：SAN为砂粒含量百分比;SIL为粉砂含量百分比;CLA粘粒含量百分比;C为有机质含量百分比;SN=1-SAN/100。

坡度因子（S）参照Liu^[41]提出的分段修正公式：

S=10.8sinθ+0.03(θ<5°)16.8sinθ-0.5(5°≤θ≤14°)21.9sinθ-0.96(14°≤θ)（6）

坡长因子（L）参照Wischmeier等^[42]和Mocool等^[43]提出的算法公式：

L=(λ/22.13)α（7）

α=β/(β+1)β=(sinθ/0.0896)/(3(sinθ)0.8+0.56)（8）

式中：θ为坡度;λ为坡长（m）;α为坡长因子指数。

植被覆盖与管理因子（C）参照Van der Knijff等^[44]提出的计算公式：

C=exp-αNDVIβ-NDVI（9）

式中：α和β是决定NDVI-C关系曲线的参数,α = 2,β = 1。

参照陈思旭等^[45]、许月卿等^[46]的研究成果,确定研究区水土保持措施因子（P）（表3）。K、LS、C、P各因子图层如图2所示。

显示原图| 下载原图ZIP| 生成PPT

图2   RULSE模型因子图层

Fig. 2   RULSE model factor layers

4 结果分析

4.1 喀斯特关键带不同岩性降雨侵蚀力

未考虑岩性差异,将全年降雨都计算为降雨侵蚀力时,非喀斯特区域降雨侵蚀力的平均值为7787.23 MJ mm ha^-1 h^-1 yr^-1,喀斯特区域降雨侵蚀力的平均值为7321.58 MJ mm ha^-1 h^-1 yr^-1。非喀斯特区域与喀斯特区域的平均值仅有465.65 MJ mm ha^-1 h^-1 yr^-1的差值,不同岩性背景下的降雨侵蚀力的差异未被区分开（图3a、3b）。

显示原图| 下载原图ZIP| 生成PPT

图3   不同算法降雨侵蚀力及其差值与改进幅度

Fig. 3   Rainfall erosivity of different algorithms and its minus and improvement rate

考虑到不同岩性下垫面产汇流差异,为不同的岩性区域限定不同的有效降雨侵蚀阈值。将日降水量≥ 12 mm作为非喀斯特区域的有效侵蚀降雨阈值时,非喀斯特区域平均降雨侵蚀力为6434.31 MJ mm ha^-1 h^-1 yr^-1。将日降水量≥ 30 mm作为喀斯特区域的有效侵蚀降雨阈值时,喀斯特区域平均降雨侵蚀力为3854.9 MJ mm ha^-1 h^-1yr^-1。非喀斯特区域与喀斯特区域的平均降雨侵蚀力差值为2579.41 MJ mm ha^-1 h^-1 yr^-1,喀斯特区域的平均降雨侵蚀力仅为非喀斯特区域的59.91%,喀斯特与非喀斯特不同岩性对降水产生侵蚀地表径流的差异得到充分地体现（图3c、3d）,与彭涛等^{[19,20,21,22,23]}的实验观测相符。与有效阈值算法相比（图3e、3f、3g、3h）,无阈值算法将非喀斯特区域的平均降雨侵蚀力高估了17.37%,将喀斯特区域的平均降雨侵蚀力高估了47.35%。

4.2 喀斯特关键带不同岩性背景区土壤允许流失量的计算与空间格局

根据碳酸盐岩不同组合中酸不溶物的平均含量和碳酸盐岩的平均风化溶蚀速率计算碳酸盐岩地区的成土速率并将其作为土壤允许流失量。结果表明,连续性碳酸盐岩组合分布区为0.21 t ha^-1 yr^-1,碳酸盐岩夹碎屑岩组合分布区为1.2 t ha^-1 yr^-1,碳酸盐岩与碎屑岩互层组合为2.89 t ha^-1 yr^-1。

连续性碳酸盐岩组合亚类中,连续性白云岩组合区域为0.19 t ha^-1 yr^-1,连续性石灰岩组合区域为0.24 t ha^-1 yr^-1,白云岩与石灰岩混合组合区域为0.21 t ha^-1 yr^-1;碳酸盐岩夹碎屑岩组合亚类中,白云岩夹碎屑岩分布区域为1.12 t ha^-1 yr^-1,石灰岩夹碎屑岩分布区域为1.29 t ha^-1 yr^-1;碳酸盐岩与碎屑岩互层组合亚类中,白云岩与碎屑岩互层分布区域为2.76 t ha^-1 yr^-1,石灰岩与碎屑岩互层分布区域为3.03 t ha^-1 yr^-1。

显示原图| 下载原图ZIP| 生成PPT

图4   不同岩性土壤允许流失量及其空间分布

Fig. 4   Soil loss tolerance and its spatial distribution of different lithology

由图4可知,喀斯特地区存在多个土壤允许流失量,具有多样性的特点。但是,同样是碳酸盐岩,在连续性碳酸盐岩、碳酸盐岩夹碎屑岩、碳酸盐岩与碎屑岩互层中,其土壤允许流失量却差别较大,这是因为在喀斯特地区,影响其土壤允许流失量大小的关键因素除了岩石中酸不溶物的含量外,还与区域地层中碎屑岩的含量有关。

4.3 喀斯特洼地空间格局与碳酸盐岩分布的匹配性分析

洼地大小与其汇集降雨面积与原始地形有关,其规模和深度随汇集面积而变化。从空间上看,洼地主要集中分布于广西北部、云南的东部、贵州除黔东南州外的大部分区域、重庆市全境、四川的东部、湖北西部及湖南西北部地区。研究区内的洼地总面积为70.82×10⁵ ha,占区域总面积的3.65%;其中喀斯特洼地面积为25.80×10⁵ ha,占总洼地面积的36.43%;非喀斯特洼地面积为45.01×10⁵ ha,占总洼地面积的63.57%。

显示原图| 下载原图ZIP| 生成PPT

图5   中国西南洼地空间分布及局部放大示意图

Fig. 5   Spatial distribution and partial enlargement of depressions

洼地在空间上连续性与碳酸盐岩的分布区基本重合（图5）,属喀斯特类型的洼地。这类型的洼地在成因上是由碳酸盐岩的可溶蚀性,加之长期的水流作用下溶蚀形成。此外,还有部分非喀斯特洼地的存在,主要集中分布于四川盆地中部以及零星分布于其他区域。这部分洼地的成因有别于喀斯特洼地,主要是由早期以地壳差异性抬升作用为主的构造运动形成^[47]。

4.4 传统算法与优化算法对土壤侵蚀量的测算与对比分析

对多次校正前后的土壤侵蚀进行统计（表4,图6）,结果显示,将全年降雨都计算为降雨侵蚀力时,研究区的土壤侵蚀面积为14.07×10⁷ ha,占总面积的73.19%,土壤侵蚀总量为90.9×10⁸ t yr^-1,平均侵蚀模数为47.3 t ha^-1 yr^-1。

显示原图| 下载原图ZIP| 生成PPT

图6   不同算法中国西南土壤侵蚀模数

Fig. 6   Soil erosion rating of different algorithms

将日降水量≥ 12 mm和≥ 30 mm分别做为非喀斯特与喀斯特区域的有效侵蚀降雨阈值,对降雨侵蚀力进行优化改进后,研究区的土壤侵蚀面积为13.39×10⁷ ha,占总面积的69.66%,土壤侵蚀总量为56.71×10⁸ t yr^-1,平均侵蚀模数为29.51 t ha^-1 yr^-1。与传统算法相比,第一次改进后的侵蚀面积、侵蚀量和平均侵蚀模数分别下降4.83%、37.61%、37.61%。

将碳酸盐岩地区的岩石风化速率作为最大土壤允许流失量阈值,对喀斯特地区的土壤侵蚀量进行校正后,研究区的土壤侵蚀面积为10.25×10⁷ ha,占总面积的53.32%,土壤侵蚀总量为47.63×10⁸ t yr^-1,平均侵蚀模数为24.79 t ha^-1 yr^-1。与传统算法相比,第二次改进后的侵蚀面积、侵蚀量和平均侵蚀模数分别下降27.14%、47.60%、47.58%。

洼地是区域的负地貌,是区域范围内物质汇聚的沉积区。将研究区范围内的洼地看作是一个独立的闭合系统,不会产生洼地范围以外的土壤侵蚀,因此洼地的真实侵蚀量应计算为0 t yr^-1。研究区的土壤侵蚀面积为10.16×10⁷ ha,占总面积的53.32%,土壤侵蚀总量为47.52×10⁸ t yr^-1,平均侵蚀模数为24.73 t ha^-1 yr^-1。与传统算法相比,第三次改进后的侵蚀面积、侵蚀量和平均侵蚀模数分别下降27.79%、47.72%、47.72%。

不同算法下的土壤侵蚀有着如下共同特点,就侵蚀面积而言,土壤侵蚀等级越高则占总面积的比例越低,其中强度、极强度及剧烈侵蚀面积均在10%以下。就侵蚀量而言,侵蚀等级越高则侵蚀量越大,对总侵蚀量的贡献率越高,其中以剧烈侵蚀最大,各算法均大于50%;极强度侵蚀次之,均大于13%。

5 讨论

高效精确地估算喀斯特关键带土壤侵蚀,明确其发生的类型及强度,对喀斯特关键带生态系统管理和社会经济良性发展具有重要的理论和运用价值。本文基于喀斯特关键带不同岩性差异性及正负地形源汇关系,对喀斯特关键带的土壤侵蚀进行多次改进（表5）。

喀斯特地区并不是每场降雨都产生径流,主要是大雨（25 mm ≤ x < 50 mm,x为雨强）尤其是暴雨（x ≥ 50 mm）才生产径流^[23]。设定合理的有效侵蚀性降水阈值,能够区分发生和不发生侵蚀的降雨事件,有效地剔除未产生侵蚀的降雨量^[30],使降雨侵蚀力的计算结果更加精确,从而解决因降雨侵蚀高估而造成土壤侵蚀高估的问题。本文将喀斯特与非喀斯特相区分,在未改进降雨侵蚀力时,喀斯特区域的平均降雨侵蚀力相当非喀斯特区域的94.02%;改进降雨侵蚀力后,喀斯特区域的平均降雨侵蚀力仅为非喀斯特区域的59.91%。

喀斯特关键带因其成土物质的先天不足,成土速度慢、土层薄,实际成土量远小于其理论侵蚀量。按水利部的传统侵蚀标准,喀斯特地区的土壤侵蚀强度以微度和轻度侵蚀为主,但喀斯特地区较低的侵蚀强度是由于土壤的限制^[32]。虽然喀斯特地区土壤侵蚀的绝对量小,但其土壤侵蚀的风险性往往非常严重。因此,急需建立一套适用于喀斯特地区的土壤侵蚀分类分级标准及风险评价方法。此外,在土壤侵蚀过程中喀斯特洼地系统应该被看作是完整的系统。虽然洼地在坡面尺度上存在坡面汇入径流沟蚀、超渗径流面蚀等水土流失过程,从其最终的侵蚀结果来看,被侵蚀的土壤仍处在整个洼地系统之中,只是其位置发现了变化。被侵蚀搬运的土壤在洼地底部富集形成可耕作层,在一定程度上增加了喀斯特区域的耕地面积。

通过表6的对比可知,通过对土壤侵蚀算法的多次改进,能够有效的剔除因降雨侵蚀力及实际土壤保有量高估造成的计算误差,这种考虑喀斯特土壤侵蚀可承受阈值的方法,实现了对喀斯特土壤侵蚀量的精确量化。

由于喀斯特特殊的二维三元地质结构,使得其土壤侵蚀存在地表流失与地下漏失两种形式。需要指出的是,本文主要是针对喀斯特地区的地表土壤侵蚀进行的算法改进,地下漏失部分将会作为进一步的研究,还有待完善。

6 结论

（1）将全年降雨都计算为降雨侵蚀力时,非喀斯特与喀斯特区域的降雨侵蚀力没有显著差异。设定不同的有效降雨侵蚀阈值后,非喀斯特与喀斯特平均降雨侵蚀力的差值为2579.41 MJ mm ha^-1 h^-1 yr^-1,喀斯特区域的平均降雨侵蚀力仅相当于非喀斯特区域的59.91%,不同岩性对降水产生侵蚀地表径流的差异得到充分地体现。

（2）在碳酸盐岩地区不同岩层组和类型的成土速率差别较大,在连续性碳酸盐岩、碳酸盐岩夹碎屑岩、碳酸盐岩与碎屑岩互层中其土壤允许流失量分别为0.21 t ha^-1 yr^-1、1.2 t ha^-1 yr^-1、2.89 t ha^-1 yr^-1。碳酸盐岩地区存在多个土壤允许流失量,具有多样性的特点。

（3）研究区内的洼地总面积为70.82×10⁵ ha,占区域总面积的3.65%,洼地在空间上的连续性与碳酸盐岩的分布区基本重合。其中喀斯特类型的洼地数量较多,空间分布相对集中,但是规模面积较小;而非喀斯特洼地的数量相对较少,空间分布较为零散单一,其面积规模通常较大,这是由于其成因的差异性导致。

（4）相比于全年降雨累计算法,经过多次校正后的土壤侵蚀面积降低了3.91×10⁷ ha,修正幅度达27.79%;土壤侵蚀总量减少了43.38×10⁸ t yr^-1,修正幅度达47.72%;平均侵蚀模数减低了22.57 t ha^-1 yr^-1,修正幅度达47.72%。

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子