水力侵蚀是引起土壤侵蚀的主要因素，水力侵蚀分为径流侵蚀和降雨侵蚀两种类型，其中降雨侵蚀力是指降雨引起土壤侵蚀的潜在能力, 降雨侵蚀力R值是评估土壤侵蚀潜在能力的一个动力指标^[1]。研究降雨侵蚀力R值的时空分布特征对区域的水土保持和防治具有重要的现实意义。目前，降雨侵蚀力的研究主要基于USLE模型, 以研究区的小区实测数据及该研究区现状对USLE模型进行改进，分析降雨侵蚀力R值的时空分布规律^[2]。在降雨侵蚀力R值的计算方法上，对其较早进行研究的是外国学者Wischmeier，他通过大量数据的分析后发现降雨的总动能(∑=E)及其30 min最大雨强(I₃₀)的乘积与土壤流失量的关系是最为紧密, 得出R=∑E·I₃₀的降雨侵蚀力表达式^[3]。随后世界上很多国家都利用该表达式计算R值。然而该公式是基于美国的观测数据求得的，在具体运用到其他国家和地区时，由于地质地貌、土壤类型等因素存在较大差异其降雨特征也存在很大的不同，该表达式很难满足实际需要。对此，不少的学者针对研究区的地形、地貌、土壤等自然特征，制定出适合该研究区降雨侵蚀力的最佳算法。例如：日本学者种田行男等^[4]通过对各地区进行试验分析研究，最终选择以E·I₁₀和E·I₆₀作为降雨侵蚀力R值的指标, 即R=∑E·I₁₀和R=∑E·I₆₀。中国学者针对西北黄土地区提出的降雨侵蚀力R值的计算方法R=∑E₆₀·I₁₀或R=∑E·I₃₀^[5], R=∑E·I₁₀^[6], R=P·I₃₀^[7]，这些方法运用到相应地区都有较好的精度。但是一般情况下很难获得某一地区的长时间的降雨强度、最大雨量、持续时间、降雨动能等详细资料, 且由于站点多，时间顺序长，使得整理过程十分繁琐, 基于对以上情况的考虑，在计算降雨侵蚀力时一般都会建立降雨侵蚀力的简易算法, 从而省去常规计算中难以获取的因子。例如：周伏健等^[8]提出的适宜福建地区的降雨侵蚀力R值计算的方法R=$\sum\limits_{i=1}^{12}$(-2.639 8+0.304 6P_i)；吴素业等^[9]，针对大别山地区提出的年降雨侵蚀力R值的估算公式R=$\sum\limits_{i=1}^{12}$(0.012 5P_i^{1.629 5})，章文波等^[10-11], 伍育鹏等^[12]提出的利用不同降雨资料对降雨侵蚀力进行计算的简易方法等。通过实际检验都有较高的精度，能够满足对土壤侵蚀计算的要求。

贵州省以喀斯特地貌为主，山高谷深，石灰岩发育使得多数地区土层普遍较薄，加上长时间大规模的破坏性拓荒，使得地表植被、土壤的覆盖遭到破坏，导致水土流失和石漠化成为贵州省生态环境最突出的问题^[13]。遥感调查数据显示，水土流失面积占贵州省总面积的41.62%^[14]。因此，对贵州省的降雨侵蚀力R值进行估算及时空特征进行分析显得尤为重要。本文通过贵州省33个气象站实测逐日降雨数据，反演计算近52 a贵州省降雨侵蚀力，并分析其时空变化特征，以期为贵州省评估降雨对土壤侵蚀的潜在作用、防治、制定水土保持措施及农业生产规划提供参考。

1 资料与方法 1.1 资料准备及选用数据来源

降雨资料主要来源于中国气象数据共享网(http://data.cma.cn/)提供的贵州省内33个气象监测站1966—2017年的日降雨数据，部分地区年降雨侵蚀力数据来源于贵州省水土保持监测站。各站点分布状况见封2附图 1，基本信息如表 1所示。

1.2 研究方法 1.2.1 降雨侵蚀力R值的计算方法

由于数据资料的局限性，本研究选用章文波等^[11]提出的利用日降雨量计算降雨侵蚀力R值的方法对贵州省1966—2017年的降雨侵蚀力R值进行计算，然后再对求得的数据进行时空特征分析，其中降雨侵蚀力R值的计算方法为：

$M=\alpha \sum\limits_{i=1}^{k}\left(P_{j}\right)^{\beta} $

(1)

式中：M为某个半月时段的降雨侵蚀力值MJ·mm/(hm²·h)；k为半月时段内侵蚀性降雨的天数(d); P_j为某个半月时段内第j天的日降雨量(mm)，要求其日降雨量必须大于12 mm(12 mm是降雨量是否为侵蚀性降雨的阈值)，否则以0计算，α与β为待定参数，这里α与β选取基于章文波等人提出的日雨量的最优参数进行计算，α=1.219 7，β=1.631 5^[10]。

1.2.2 EOF分析法

经验正交函数(empirical orthogonal function, EOF)分析法, 又称主成分分析法，是气候气象特征研究中较为常用的方法，它可以对某一区域内观测站点的值进行时空分析，能准确的反映相应的空间模态随时间的权重变化，即反映其空间分布特征。

1.2.3 突变检测及周期检验方法

突变检测和周期检验是分析长时间序列事物变化发展的重要内容之一，可以通过突变检测和周期变化清晰直观的了解事物的发展过程及其变化规律。当前检测突变的方法有许多，通过研究各种检验方法发现，H.B.Mann和M.G.Kendall提出的Mann-Kendall检验法有较好的理论基础和实际应用效果，目前该方法已经在气温，降雨等方面有了广泛的使用^[15-16]。其计算过程如下：

设气候序列为，x₁, x₂, …, x_n, S_k表示第i个样本，定义统计量：

$S_{k}=\sum\limits_{i=1}^{n} r_{i}, \quad r_{i}=\left\{\begin{array}{ll}1 & \left(x_{i}>x_{j}\right) \\ 0 & \left(x_{i}=x_{j}\right) \\ 0 & \left(x_{i}<x_{j}\right)\end{array}\right.\\ (其中 i<j, j=1, 2, 3, \cdots, n, k=1, 2, 3, \cdots, n $

(2)

在时间序列随机独立的假定下, S_k的均值和方差为：

$E\left[S_{k}\right]=\frac{k(k-1)}{4}, \\ \operatorname{var}\left[S_{k}\right]=\frac{k(k-1)(2 k+5)}{72} \quad(1 \leqslant k \leqslant n) $

(3)

将S_k标准化：$\mathrm{UF}_{k}=\frac{\left(S_{k}-E\left[S_{k}\right]\right)}{\sqrt{\operatorname{var}\left[S_{k}\right]}}$，UF₁=0反序列UB_k=-UF_k，k=n，n-1, …, 1，UB₁=0

当UF_k或者UB_k大于0时呈上升趋势，反之则呈下降趋势。给定显著性水平α，若|UF_k|＞U_α则存在明显的变化趋势，反之变化不明显。如果UF_k和UB_k出现相交情况，并且在两条显著性水平线之间，则交点即是存在突变的开始时间^[15-17]。

在周期变化检验上主要采用在水文分析方面有着显著成效的小波分析法，小波分析法是在傅里叶变化的基础上发展而来的，其对获取一个复杂时间序列的调整规律, 分辨时间序列在不同尺度上的演变特征非常有效.小波分析方法用于检验信号突变可以较准确的检验出信号突变时间点^[17]。本文选用水文研究中常用的复Morlet小波分析方法对贵州省的降雨侵蚀力R值进行研究，其表达为：

$\varphi(t)=e^{i \omega_{0} t} e^{-t^{2} / 2} $

(4)

式中：φ(t)为高斯包络下的复指数函数；i表示虚数；ω₀为常数^{[15, 18]}；t为时间。

1.3 降雨侵蚀力R值计算精度评估

由于本研究所获得的降雨数据资料没有详细的降雨细节记录。通过对比最终选择了基于日降雨量数据计算降雨侵蚀力R值的简易算法对本研究区进行计算，经过前人的多次验证之后发现，该方法在降雨较丰富的地区有较好的精度，而在降雨量较小的地区该方法计算的降雨侵蚀力相对误差较大。为了检验此方法在贵州地区的适用程度，抽取4个不同地区前人所做的结果与本次计算降雨侵蚀力的数据进行验证，验证结果如表 2所示, 通过对比可以看出毕节、贵阳、兴义、织金4个站均与前人通过经典算法或简易算法计算的值非常相近，相对误差均小于5.3%。因此，章文波等^[10]提出的利用日降雨量计算R值的算法和其他经典算法以及简易算法在计算贵州地区降雨侵蚀力时有较好的契合度。

由于贵州省属于喀斯特地区，地理环境与其他地区有所不同，因此需对本文所用的简易方法进行进一步验证，以证明该方法在贵州的适用性。为了验证该方法的适用性，本文收集了2012年至2017年期间六盘水(龙贵地)、安顺(三股水)、贵阳(修文)、黔南(贵定)、以及遵义(浒洋水)等各水土保持监测站点计算的年降雨侵蚀力值作为精度评估样本。各水土保持监测站计算的年降雨侵蚀力值是由次降雨过程计算出的值具有一定的精确性以及真实性。利用基于日降雨侵蚀力的方法估算的降雨侵蚀力R值与水土保持监测站点计算的年降雨侵蚀力R值进行拟合分析，从图 1可知，水土保持站计算值与本文估算值总体趋势保持一致，大部分时间点两者值相差不大。从图 2中也可以看出两者总体拟合情况较好，R²值达到了0.668 8。因此利用基于日降雨数据估算贵州省的降雨侵蚀力具有一定的科学性和合理性。

2 结果与分析 2.1 贵州省降雨侵蚀力R值分布特征

通过公式(1)对贵州省33个气象检测站点52 a的日降雨数据进行计算，得到全省各气象检测站的多年平均降雨侵蚀力R值，通过ArcGIS软件的克里金插值方法对贵州全省进行插值得到贵州省年平均降雨侵蚀力R值空间分布图、汛期年平均分布图以及枯水期年平均分布图。

从附图 2(见封2)可知贵州地区降雨侵蚀力R值分布的地域差异较明显的，贵州省各地区的多年平均降雨侵蚀力R值介于2 767.46~6 921.73 MJ·mm/(hm²·h)之间，平均R值为4 994.22 MJ·mm/(hm²·h)；贵州省降雨力R值地域分布总体上呈由北向南递增的趋势，在西南部、南部、东部边缘多年平均降雨侵蚀力R值较大，而西北部、北部多年平均降雨侵蚀力R值偏低。其中，在兴义、望谟、安顺、织金以及都匀地区一带形成了降雨侵蚀力R值的高值中心，R值均高于6 000 MJ·mm/(hm²·h)，相反在威宁以及毕节地区一带形成降雨侵蚀力R值的低值中心，其值均低于3 000 MJ·mm/(hm²·h)；形成这种分布特征的原因主要是由于贵州省的地形地貌造成的，从附图 1(见封2)中可以看出中西部、西北部地势较高，而兴义、安顺、都匀地区恰好处在夏季风迎风坡一带容易形成降雨，反之在毕节一带地势较高，最高海拔达2 800 m左右，当湿热水汽从南方北上时，大部分水汽已在迎风坡上升遇冷形成地形雨，因此到达西北高海拔地区的水汽较少，降雨也较少。

2.2 经验正交函数EOF分析 2.2.1 贵州省降雨侵蚀力空间分布特征

为了进一步分析贵州省降雨侵蚀力R值的空间分布特征，利用Matlab2017软件对贵州省内33个气象监测站点测得的降雨资料计算的1966—2017年的降雨侵蚀力R值进行EOF分析。经检验，KOM值约为0.665，Bartlett球形检验发现，近似卡方值为2 036.261，自由度36，p=0.000，说明适合做EOF分析。从表 3中可以看出，前5个特征向量贡献值的累计贡献率达69%，但是只有前两个特征根误差范围没有重叠通过了North显著性检验，累计贡献率达52%，基本上能解释贵州省降雨侵蚀力R值的两种空间分布特征。

由表 3可以看出，第一特征向量的方差贡献率为37%，远高于其他模态贡献率，是贵州省降雨侵蚀力R值的主要空间分布形式。由图 3可知，第一模态的特征值均为正值，这表明了贵州省1966—2017年间，贵州省的降雨侵蚀力R值变化趋势具有一致性，导致这一现象的主要原因是由于贵州省较大的降雨受西太平洋副热带高压的影响，其与该气压的面积指数与强度呈正相关，由此导致贵州省的降雨侵蚀力在大尺度范围上的变化有明显的一致性。高值区位于安顺市、都匀市等中部偏南地区，反映了该地区年降雨侵蚀力R值变化幅度较大。低值区主要在东北部，降雨侵蚀力R值变化幅度较小。现有研究表明在中部偏南地区多极端降雨，对枯丰都比较敏感，而东北部年降雨量较为稳定变化不大^[22]。

第二特征向量方差贡献率为15%，其表现出降雨侵蚀力R值东西呈现相反的分布模式，零值先出现在105°E左右，这体现了贵州省东部降雨侵蚀力R值较大时，西部就较小，而西部较大时，东部则较小。正值中心在西南部，负值中心出现在都匀附近。造成这种现象的主要原因主要是地形所致，以105°E分为东、西两部分，东部海拔较低，最低处仅有287.4 m，而西部的平均海拔高度在1 600 m左右，高度差异较大，所受得天气系统影响不同，造成降水差异，进而导致降雨侵蚀力差异。

2.2.2 降雨侵蚀力R值EOF时间系数分析

通过对EOF时间系数进行分析，贵州省的降雨侵蚀力R值主要有4种表现类型，第一模态决定了全年全省降雨侵蚀力R值较大或者全年全省降雨侵蚀力R值较小两种类型；第二模态决定了全年贵州省东部降雨侵蚀力R值偏大还是东部偏小，或者西部偏大东部偏小两种类型。EOF时间系数代表了对应特征向量空间分布模态的时间变化特征，系数的正负决定了模态的方向，正好表示与模态同向，负号表示与模态方向相反，并且系数的绝对值越大越能决定这一时刻的对应模态越典型。对贵州52 a的降雨侵蚀力R值的4种空间模态分布类型进行统计分析，取每年EOF时间系数绝对值最大值对应的特征向量作为对应年的降雨侵蚀力空间分布模态，结果如图 4所示，得出贵州省降雨侵蚀力R值偏大年为19 a，全省降雨侵蚀力R值偏小年为18 a，东部降雨侵蚀力R值偏大年为7 a，西部降雨侵蚀力偏大年为8 a。由此可看出贵州省52 a降雨侵蚀力R值，以第一向量的空间分布为主，达到了37 a，占总年数的71.15%，全省全年降雨侵蚀力偏大和偏小年相当。东部和西部差异型的降雨侵蚀力年数共15 a，占总年数的28.85%。这也与一、二模态贡献率反映的典型程度基本表现一致。

2.3 降雨侵蚀力R值年内变化趋势

通过对贵州省1966—2017年各月降雨量比例(见图 5)进行分析确定贵州省的汛期为5—8月，枯水期为8月至翌年5月。通过计算贵州地区汛期和枯水期降雨侵蚀力R值并对其空间插值后分析发现，贵州各地区的降雨侵蚀力R值在年内的变化十分显著。从附图 3(见封2)中可以看出贵州省的汛期降雨侵蚀力R值分布与多年平均降雨侵蚀力R值在空间分布上十分相似，最小值依然出现在西北部降雨侵蚀力R值为2 145 MJ·mm/(hm²·h)，最大值也同样出现在西南部的兴义附近其值为5 195 MJ·mm/(hm²·h)。从附图 4(见封2)可以看出，贵州省各地区多年平均枯水期降雨侵蚀力R值在582~1 833 MJ·mm/(hm²·h)之间，空间特征整体呈现由西向东递增的变化趋势，最小值出现在毕节一带，最大值出现在松桃、天柱、黎平一带。

从图 6可以看出大雨和中雨造的降雨侵蚀力占比达到了33.09%和44.36%，通过对贵州地区多年平均降雨进行分析发现，贵州省的大雨和暴雨主要集中在汛期，因此可以得出汛期降雨是导致降雨侵蚀力R值变化的主要原因。

贵州省属于亚热带季风气候，夏季高温多雨，冬季温和少雨。夏季贵州省的东部处于湿润的东南季风区内，西部则处于无明显的干湿季之分的东南季风向干湿明显的西南季风区的过渡地带，加之地形西高东低，整体呈中部向北、东、南3面倾斜，易在东南季风的迎风坡形成地形雨，由此导致了夏季西南部及东部降雨侵蚀力R值较大，西部和西北部则相对较低。而在冬半年由于受北方秦巴山系阻挡，南下冷空气多半绕道两湖盆地由偏东北方向入侵，在贵阳昆明一带形成准静止锋，西部六盘水、毕节一带经常处于锋前位置，故冬季多晴朗天气，在中部、东部正好处于锋后，因此多阴雨天气。

尽管冬季东部降雨天数较多，但是整体的降雨强度较低，达侵蚀性降雨的标准天数较少，因此造成了降雨侵蚀力R值冬春季节整体偏低。

2.4 降雨侵蚀力R值年际变化趋势

通过Mann-Kendall检验法对贵州省的多年年均降雨侵蚀力R值与汛期多年年均R值进行突变检测，检测结果如图 7a, 7b所示，分析发现贵州省降雨侵蚀力R值在0.05置信度线区间内UF与UB两条统计量曲线存在多个交点，分别是1971，1974，1975，1978，1981，1991，1992，2003，2014，2016年，这说明贵州省的降雨侵蚀力多突变，震荡频繁，尤其是在1971年到1981年之间突变频率最为频繁。

从UF统计量曲线来看，除了少数年份外贵州省的降雨侵蚀力年际R值多处于负值，但均没有通过0.05显著性水平线，这也进一步证明了其降雨侵蚀力变化趋势总体呈现下降趋势，但趋势不明显。从汛期年际降雨侵蚀力R值与年际降雨侵蚀力R值图的对比中可以发现两者的变化趋势基本保持一致，这体现了多年年均降雨侵蚀力R值主要受汛期降雨侵蚀力R值的影响。

从图 8贵州省年降水量小波实部时频分布特征可以清晰的看出贵州省降雨侵蚀力R值的变化过程中存在多时间尺度特点，在降雨侵蚀力的演变过程中存在3~8，8~21，21~32 a共3类时间尺度的周期变化规律。其中在3~8 a时间尺度上出现了升—降交替的准11.5次震荡；在8~21 a时间尺度经历了升—降交替的准6次震荡；在22~32 a时间尺度上存在准2.5次震荡。经图 9小波方差分析可以看出存在3个较为明显的波峰，依次对应的是6，12，28 a的时间尺度，其中28 a的峰值最大为其震荡主周期，其次是12 a和6 a的时间尺度，从主周期来看，贵州省的降雨侵蚀力R值正处于升高的范围内，并未达到峰值阶段，因此，未来几年内贵州省的降雨侵蚀力R值仍有升高趋势。

3 结论

本文基于贵州省内33个气象站监测站48 a的逐日降雨数据，运用降雨侵蚀力R值的相关计算方法，计算贵州省的多年平均降雨侵蚀力R值并对其时空特征进行研究分析，主要结论有：

(1) 贵州省的多年平均降雨侵蚀力R值在2 767.46~6 921.73 MJ·mm/(hm²·h)之间，地域分布较为明显，R值总体上呈由北向南递增的趋势，最低值出现在威宁以及毕节一带，导致其分布的主要原因是地势所造成的。

(2) EOF分析表明，贵州省多年年均降雨侵蚀力空间是主要分为两种类型，全局型以及东西相反型，这两种空间分布类型的累计方差贡献率达到了52%，基本上可以反映贵州省降雨侵蚀力的空间分布特征。贵州省降雨侵蚀力EOF时间系数分析表明，贵州省降雨侵蚀力偏大年份数量与降雨侵蚀力偏小年份量相当，分别是19 a和18 a。

(3) 通过分析得出贵州省的降雨侵蚀力R值主要受汛期降雨R值影响，侵蚀性降雨主要集中在汛期(5—9月)，全省各地汛期多年平均降雨R值均占多年平均降雨侵蚀力R值的60%以上，最大占比达到了80%左右。

(4) 运用Mann-Kendall检验法对贵州省的年降雨侵蚀力R值和汛期多年平均R值进行突变检测，发现年际降雨侵蚀力R值存在贵州省的降雨侵蚀力多突变，震荡频繁现象，变化趋势总体呈现下降趋势，但下降趋势不明显。

(5) 贵州省降雨侵蚀力R值变化的主周期为28 a左右，目前降雨侵蚀力R值正处于逐渐升高的范围内，并未达到峰值阶段，因此，未来几年内贵州省的降雨侵蚀力R值仍有升高趋势。