1960-2017年贵州省不同水系降雨侵蚀力时空变异特征

2021-04-16朱大运李少男

水土保持通报 2021年1期

朱大运, 杨倩, 陈海, 陈静, 李少男

(1.贵州师范大学喀斯特研究院, 贵州贵阳 550001; 2.国家喀斯特石漠化防治工程技术研究中心, 贵州贵阳 550001)

雨滴溅蚀和径流剥蚀是降雨作用于土壤地表侵蚀的主要途径[1]，可通过降雨侵蚀力指标进行表达，与降雨量、降雨强度、降雨历时、地形植被、人类活动等因子要素密切相关[2-3]。当前全球气候变化大背景下，极端气候事件发生的概率明显增多[4]，某种程度上增加了土壤侵蚀的风险。科学量化评估区域降雨侵蚀力，对水土流失定向阻控和水保措施制定具有重要指导意义。

以贵州高原为中心的西南喀斯特地区是中国水土流失最严重的区域之一[5]，其降雨侵蚀问题已引起广泛关注。基于EI30模型[6]、次降雨模型[7]、日降雨模型[8]、月降雨模型等[9]多种降雨侵蚀力计算模型，学者对贵州降雨侵蚀力时空特征开展了不同程度的研究。许月卿等[10]首先利用19个站点的气象数据对贵州降雨侵蚀力时空特征进行了初步分析，指出降雨侵蚀力在空间上由南到北递减，时间上呈增加趋势，且年内分配不匀主要集中在夏季。这与戴海伦等[11]基于贵州罗甸小区实测数据和多站点雨量资料得到的贵州降雨力时空特征相一致。阮欧等[12]基于正交函数、Mann-Kendall检验、小波分析等分析了贵州降雨侵蚀力的空间分布、突变性和周期性规律。Zhu等[13]研究了贵州省长时间尺度降雨侵蚀力与厄尔尼诺—拉尼娜事件的关系，结果表明ENSO持续时间与降雨侵蚀力呈正比。上述研究成果在一定程度上丰富了对贵州地区降雨侵蚀力时空分布特征的认知。但是，关于降雨侵蚀力水系分布特征却鲜有报道，贵州喀斯特生境特征脆弱，石漠化广泛分布且横跨长江、珠江2个流域，局地气候与土壤侵蚀特征较为复杂，因此揭示长时间尺度不同水系降雨侵蚀力的内在差异十分必要。本文以日降雨资料为基础，从水系的角度对贵州省1960—2017年降雨侵蚀力时空分布特征开展对比研究，以期为区内水土流失防治、生态保护等工作提供科学参考。

1 数据与方法

1.1 研究区概况

贵州省地处云贵高原东部，位于103°36′—109°35′E,24°37′—29°13′N之间，属亚热带湿润季风气候，受西南季风和东亚季风的双重影响。年平均气温15 ℃左右；多年平均降水量1 100～1 400 mm，其中47%的降水集中在夏季；植被类型丰富，以亚热带常绿阔叶林为主。水资源分布以苗岭为分水岭，南部属珠江流域，境内面积6.00×104km2，下辖南盘江水系、北盘江水系、红水河水系和柳江水系；北部属长江流域，境内面积1.16×105km2，下辖牛栏江水系、乌江水系、赤水河水系和沅江水系(图1)。贵州地貌以高原山地为主，山地和丘陵面积占92.5%。地势西高东低，起伏较大，岩溶发育典型，石漠化与水土流失十分严重。据水利部全国水土流失遥感调查结果和贵州省水土保持规划(2016—2030)，20世纪80年代贵州省水土流失面积占土地总面积的43.52%，经过多方持续不断的努力，水土流失得到有效控制，到2015年，水土流失面积占比下降到27.71%，防治形势依然严峻。

图1 贵州省境内水系与气象站点分布

1.2 数据源

贵州省1960—2017年33个国家地面标准观测站的逐日降水量数据来源于中国气象科学数据共享服务网(http:∥data.cma.cn/)。遵循欧洲气候评估数据审订标准(数据不得少于40 a，单个站点数据总缺失量不得超过10%，年数据缺失量不得超过20%或超过连续3个月)对原始数据进行了检验[14]，各站点气象数据均通过严格的质量控制和均一性检验。对个别缺失数据采用相邻站点线性回归方法进行插补，经过处理修正后的数据序列具有很好的连续性。

1.3 研究方法

1.3.1 降雨侵蚀力计算采用章文波等[15]修正的Richardson日降雨侵蚀力模型计算降雨侵蚀力，该算法已经在国内大范围尺度和区域尺度上得到广泛应用。具体公式如下：

(1)

式中：Ri表示第i个半月时段的侵蚀力值〔MJ·mm/(hm2·h)〕；k表示该半月时段内的侵蚀性降雨日数(d)；Pj表示半月时段内第j天≥12 mm的日雨量；α与β为模型参数,分别根据以下公式计算:

α=21.586β-7.198 1

(2)

(3)

式中：Pd12表示日雨量≥12 mm的日平均雨量(mm)；Py12表示日雨量≥12 mm的年平均雨量;不同区域的模型参数α与β存在差异,经计算研究区参数α的数值范围为0.219 1～3.879 1,平均值为1.146 0;参数β的数值范围为1.269 7～1.893 7,平均值为1.528 4。式(1)降雨侵蚀力是以逐年各半月为基本统计单元,经过累加计算得到年降雨侵蚀力值。

1.3.2 时空变化分析方法降雨侵蚀力年际、年代际变化趋势通过线性倾向估计来进行表达，使用最小二乘法进行倾向率估计，研究其变化率。线性倾向估计法是一种非常实用有效的表征序列变化趋势和倾向性的分析方法，在气候变化研究中应用广泛。采用ArcGIS地统计模块的反距离加权插值法，分析多年平均降雨侵蚀力空间分布特征。该插值法基于相近相似的原理，主要依赖反距离的幂值和插值点与样本点的距离，与任何实际的物理过程均不关联，是土壤侵蚀点状数据到连续面状数据转换的一种有效方法[16]。

Mann-Kendall突变分析是世界气象组织推荐的一种非参数的突变检验法，具有人为干扰性少、较验范围宽的优点，可以对时间序列变化趋势的显著性和突变性进行检验。它主要基于统计量序列的分析，正序列与逆序列曲线在显著性水平边界线内的交点即为统计序列发生突变的时间[17]。运用变异系数(CV)量化描述研究区各水系降雨侵蚀力空间分布的不均匀性，衡量序列离散程度，其公式如下：

(4)

降雨侵蚀力年内空间变化采用重心模型来表达，它基于牛顿力学，动态权衡各个地区作用力的大小并指示变量变化方向，地理位置和属性变化是决定重心的主要因素。降雨侵蚀力重心模型参考钟业喜等[18]，计算公式如下：

(5)

(6)

1.3.3 聚类与相关性分析方法运用聚类分析法对降雨侵蚀力空间区域进行划分。聚类分析是对样本或变量进行分析的一种统计方法，在多个领域的区域研究中均有应用。它根据事物本身的特性将相似的事物归类，被归为一类的事物具有较高的相似性，而不同类间的事物有着很大差异。本研究借助主成分分析法获取分类数目，通过isodata算法获取像元分组特征，从而实现对研究区降雨侵蚀力的聚类划分。

降雨侵蚀力与地形区位因子的相关性分析以像元点为基本单位。在ArcGIS 10.2中，将降雨侵蚀力空间插值数据转化为矢量像元点，并提取各个水系每一个像元点的降雨侵蚀力值。同时，基于SRTM-DEM和坡度、坡向数据赋予每一个像元点经纬度坐标、高程、坡度、坡向值，而后以水系为分类单元，采用Person相关性分析法计算降雨侵蚀力与有关因子的相关系数，涉及矢量像元点累计47 300多个。

2 结果与分析

2.1 降雨侵蚀力时间变化特征对比

2.1.1 不同流域降雨侵蚀力时间变化特征 1960—2017年贵州省境内长江流域降雨侵蚀力范围3 049.77～8 112.10 MJ·mm/(hm2·h)，多年平均值5 823.42 MJ·mm/(hm2·h)；珠江流域降雨侵蚀力范围5 933.61～7 747.38 MJ·mm/(hm2·h)，多年平均值6 893.22 MJ·mm/(hm2·h)，珠江流域整体降雨侵蚀能力高于长江流域。由图2可知，过去58 a贵州省境内长江流域年降雨侵蚀力呈现出线性上升趋势，上升速率为4.13 MJ·mm/(hm2·h·a)；在年代际变化上较为复杂，经历了升高—下降—升高—下降—升高的变化过程。珠江流域降雨侵蚀力线性趋势以下降为主，速率为-10.66 MJ·mm/(hm2·h·a)；年代际变化上表现出波动特征。然而，长江流域与珠江流域年降雨侵蚀力线性趋势均未达到显著性水平。距平分析表明，珠江流域降雨侵蚀力在时间序列变化上比长江流域更加剧烈。

图2 贵州省境内珠江流域及长江流域年均降雨侵蚀力变化与距平值

2.1.2 不同水系降雨侵蚀力时间变化特征利用线性回归法对比分析了8个水系降雨侵蚀力变化趋势，除沅江水系和红水河水系表现出统计意义上的上升趋势，其余水系均为下降趋势，反映出贵州大部分区域降雨侵蚀风险呈降低态势，其中沅江水系上升速率最快达18.44 MJ·mm/(hm2·h·a)，北盘江水系下降率最快为-24.44 MJ·mm/(hm2·h·a)(图3)。各水系降雨侵蚀力波动变化特征明显且差异显著，牛栏江水系、乌江水系、赤水河水系、沅江水系、南盘江水系、北盘江水系、红水河水系、柳江水系极差比分别为：6.56，3.32，4.49，3.11，9.01，3.60，3.95，5.26。通过5 a滑动平均线和极差比结果可知，珠江流域下辖水系降雨侵蚀力波动幅度明显高于长江流域，并且整体在趋势线起伏变化上更为一致，在特定时段尤其明显，如1965—1970,2005—2015年等。长江流域下辖的4个水系降雨侵蚀力波动幅度相对较小，但在趋势变化的同向性上差异较大。

图3 贵州省境内不同水系降雨侵蚀力5 a滑动变化趋势(k为趋势系数)

2.1.3 降雨侵蚀力时间序列突变分析利用Mann-Kendall法对不同水系降雨侵蚀力时间序列突变特征进行检测，并基于移动T检验法和Yamamoto指数法综合对比剔除虚假突变点。检测结果表明，过去58 a各水系降雨侵蚀力突变特征不明显，乌江、南盘江等水系虽然检测到突变节点，但是均未通过显著性检验，仅牛栏江水系在1976年达到p<0.05显著性水平；而沅江水系、红水河水系以及境内长江流域则未检测到可靠的突变点，说明其历年降雨侵蚀力不存在自然突变(表1)。

表1 贵州省境内不同水系降雨侵蚀力时间序列突变分析

2.2 降雨侵蚀力空间变化特征对比

2.2.1 不同水系降雨侵蚀力空间变化特征图4为贵州省境内不同水系年均降雨侵蚀力空间分布图。从图4可知，长江流域与珠江流域范围内降雨侵蚀力在空间变化上表现出截然不同的趋势，前者降雨侵蚀力从西北向东南递增，后者降雨侵蚀力则从西向东递减。长江流域以沅江水系降雨侵蚀最为严重，明显高于其他水系，3个降雨侵蚀高值中心分别出现在都匀、松桃、织金；珠江流域降雨侵蚀最严重的区域为南盘江水系，兴义、望谟、安顺是该流域降雨侵蚀高值中心。

图4 贵州省境内不同水系年均降雨侵蚀力空间分布

各水系年均降雨侵蚀力呈现出强烈的空间变异性，其中珠江流域下辖水系变异系数差异较小，流域平均值为0.22略高于长江流域。变异系数值最大的牛栏江水系，最小的乌江水系均属长江流域，二者相差近1.8倍，反映出该流域强烈的区域空间变异性(表2)。此外，对比发现年均降雨侵蚀力的大小与空间变异性的强弱不存在显著的对应关系，虽然在空间变异系数最高、年均降雨侵蚀力最低的牛栏江水系表现出反向对应关系，但在其他水系均未发现这种现象，这可能是受地形、土壤等多种因素的影响。

表2 贵州省境内不同水系降雨侵蚀力变异性分析

2.2.2 基于聚类分析分区的降雨侵蚀力对比水系是根据河道集水范围而形成的区域地理单元，主要受地形地貌的影响。为了减少下垫面的影响，进一步分析降雨侵蚀力在自然降雨条件下空间集聚变化特征，基于聚类分析法对贵州降雨侵蚀力时空变异特征进行区域划分和对比。范俊甫等[19]在ArcGIS 10中对贵州省58 a降雨侵蚀力进行空间插值，运用主成分分析方法实现对空间插值数据的降维，第一主成分至第六主成分量累积贡献率达90.16%，本文以此作为isodata聚类分析中类数量划分的依据。最后基于最大似然法将isodata聚类分析计算的特征值进行像元归类，划分为6个空间区域。

由图5可知，基于聚类分析分区与流域分区的多年平均降雨侵蚀力在空间界线上呈现截然不同的特征，聚类分区更加倾向于各中心点降雨侵蚀力高低的类别特征集聚。年均降雨侵蚀力最低值为Ⅰ区3 907.59 MJ·mm/(hm2·h)，最高值为Ⅳ区7 188.76 MJ·mm/(hm2·h)，极差值为3 281.87 MJ·mm/(hm2·h)(表3)。与水系分类统计结果相比，虽然聚类分区的年均降雨侵蚀力极差值小于水系分区，但是两种分类法极值倍差相同均为1.8倍，进一步表明贵州省降雨侵蚀力较强的空间差异性。变异系数除了一区值较低为0.19，其他几个区域则相对一致，介于0.24～0.26之间(表3)。总体上，基于聚类分区的各区域变异系数比基于水系划分的变异系数要小。

图5 基于聚类分区的贵州省年均降雨侵蚀力空间分布

表3 1960-2017贵州省境内不同聚类分区降雨侵蚀力变异性分析

进一步对6个区域降雨侵蚀力的时间变异性进行了趋势分析，在长时间尺度上第Ⅲ区域、第Ⅴ区域降雨侵蚀力呈现不显著增加趋势，增加速率分别为10.87 MJ·mm/(hm2·h·a)，7.64 MJ·mm/(hm2·h·a)；其余4个区域均为不显著减少趋势，第Ⅳ区域递减速率最高达25.68 MJ·mm/(hm2·h·a)(图6)。1960—2017年各区域年均降雨侵蚀力变化曲线走向上整体上较为一致，在急骤转换节点表现尤为明显，如1990，2015年前后的降雨侵蚀力快速增加趋势，但是波动幅度存在差别，反映出各区域对气候响应的差异性。

图6 贵州省基于聚类分区的降雨侵蚀力时间变化(k为趋势系数)

2.3 降雨侵蚀力重心时空迁移特征

基于重心模型计算了贵州境内各月降雨侵蚀力中心和大雨中心的位置，以及重心迁移方向(图7)。如图7所示，从地理空间上来看，降雨侵蚀力重心主要分布在贵阳辖区与周边的龙里、贵定县，降雨侵蚀力重心滞留在龙里县境内时间最长达5个月，其次为贵阳市辖区。1—12月份降雨侵蚀力重心位置在这3个区域呈不规则方向迁移，累计迁移距离191 km，平均每月移动15.9 km,9—10月份重心迁移距离最长，单月向东平移了31.5 km，7—8月份重心迁移缓慢仅2.8 km。降雨侵蚀最为强烈的6月至9月其重心坐标全部位于贵阳市境内，与大雨(≥25 mm)的月重心位置一致。大雨月重心位置分布相对集中，除了降水较少的12，1，2月，全部集聚在贵阳东北的乌当区境内。

图7 贵州降雨侵蚀力与大雨月重心迁移

2.4 降雨侵蚀力与地理因子相关性分析

通过降雨量推演的降雨侵蚀力仅是雨水对土壤侵蚀潜在能力的反映，它与地形地貌、人类活动等要素共同作用决定了区域水土流失状况。贵州是中国典型喀斯特区，地貌特征脆弱而复杂，相关性分析表明地形、区位等要素对局地降雨侵蚀力的影响作用明显(表4)。各水系年均降雨侵蚀力与经度、纬度相关性均达到显著水平(p<0.05)，与经度的相关性主要表现为正相关，其中柳江水系相关系数最高为0.89；与纬度的相关性主要表现为负相关，北盘江水系相关系数最高为-0.79。海拔与降雨侵蚀力相关性密切，大部分水系表现为显著负相关，其中区域海拔较高的牛栏江水系和南、北盘江水系与降雨侵蚀力的负相关系数最高，这可能是受下垫面的影响。李维杰等[17]研究发现降雨侵蚀力与海拔负相关主要受地形以及地势起伏的影响，海拔较高、坡面较陡、地势起伏较大的区域较为明显，这与牛栏江、南、北盘江水系下垫面地貌特征基本吻合。

由表4可知，坡向与降雨侵蚀力的相关性不显著，但是坡度与降雨力的相关性全部通过显著性水平检验，其中西北部牛栏江水系相关程度最高，存在较为严重的降雨侵蚀风险。

表4 不同水系降雨侵蚀力与影响因素相关性分析

3 讨论

降雨侵蚀力不仅直接影响土壤侵蚀的时空分布，而且作为重要因子参与侵蚀模型的计算，为水土流失预测提供关键支撑[3,15]。本文基于贵州省33个气象站点逐日降水资料，分析了贵州省境内8个水系近60 a降雨侵蚀力时空变化特征，结果表明各水系年降雨侵蚀力呈波动变化，在长时间尺度上以线性减少趋势为主。芦鑫[20]研究发现长江流域上游年均降雨侵蚀力呈下降趋势；赖成光等[21]对珠江流域降雨侵蚀力的分析发现，南、北盘江水系呈减少趋势，红水河水系呈增加趋势，与本研究结果趋于一致。

全球气候变化与各水系降雨侵蚀力时空演变关系密切，当气温上升，蒸发旺盛，大气环流与降水空间格局在全球尺度上发生重新调整，进而导致侵蚀外营力降雨量与降雨强度发生改变[22]。贵州位于中国第二阶梯向第三阶梯过渡地带，高海拔、低纬度的区位特征，使其气候变化复杂且受季风影响明显。赵志龙等[23]分析结果表明，1960—2016年贵州省降水量年际变化剧烈，呈不显著减少趋势，与本文降雨侵蚀力变化特征十分吻合，说明降雨量减少是引起侵蚀力呈递减变化的关键因素。然而，降水量的减少只是区域气候变化的一种直观反映。

从深层次原因来看，过去几十年中国近海海温升高和青藏高原热源作用减弱，季节性海陆温差变化异常，导致东亚季风环流减弱、夏季西太平洋副热带高压位置偏南，来自东部海洋的水汽无法深入长江流域西部地区[24]；同时期西南季风的减弱现象也被研究发现[25]。因此，东亚季风与西南季风减弱的叠加影响，导致贵州地区降水减少，可能是引起贵州降雨侵蚀力呈线性递减的主要原因。

此外，关于降雨侵蚀力在长时间尺度上的波动变化特征，厄尔尼诺—南方涛动(ENSO)的影响不容忽视。有学者研究发现贵州年降雨侵蚀力与ENSO指代因子海洋尼诺指数(ONI)、多变量ENSO指数(MEI)、南方涛动指数(SOI)关系密切，并且降雨侵蚀力大小与ENSO的持续时间呈正比[13]。然而，降雨侵蚀力变化是多种影响要素共同作用的结果，后续还需要从气候变化、人类活动、地理因子等综合影响的角度进行更多的分析。