1901-2016年黄土高原降雨侵蚀力时空变化

2022-07-03陈剑南刘益麟李朋飞胡晋飞高健健党恬敏

水土保持研究 2022年4期

陈剑南, 刘益麟, 李朋飞, 胡晋飞, 高健健, 党恬敏

(1.西安科技大学测绘科学与技术学院, 西安 710054; 2.黄河水利委员会绥德水土保持科学试验站, 陕西榆林 719000; 3.黄河流域水土保持生态环境监测中心, 西安 710021)

黄土高原是全球土壤侵蚀最严重的地区之一[1]。水力侵蚀是黄土高原最重要的侵蚀类型，而降水及其产生的径流是水力侵蚀的主要驱动力。降雨侵蚀力用于表征降雨引起土壤侵蚀的潜在能力[2]，目前已广泛用于气候变化监测、土壤侵蚀模拟等研究中[3-4]。因此，厘清降雨侵蚀力时空特征是开展黄土高原土壤侵蚀研究与防治的基础。

20世纪50年代以来，中外学者针对降雨侵蚀力时空变化开展了大量研究。例如，Panagos等[5]估算了希腊1965—1994年月平均降雨侵蚀力；Liu等[6]分析了1980—2017年全球降雨侵蚀力的时空变化；刘斌涛等[7]估算了1960—2009年中国590个气象站的降雨侵蚀力，并通过插值得到中国降雨侵蚀力的时空分布；殷水清等[8]绘制了1961—2016年全国1 km降雨侵蚀力等值线图，分析了中国降雨侵蚀力的时空分布及重现期。黄土高原是我国降雨侵蚀力研究的重点区域。例如，穆兴民等[9]绘制了1956—2002年陕北黄土高原降雨侵蚀力的等值线并分析了其时空变化特征；孙从建等[10]研究了黄土高原塬面保护区内降雨侵蚀力的时空变化趋势及影响因素；李静等[11]研究了1971—2000年黄土高原不同地貌分区的降雨侵蚀力的时空特征，并得到降雨侵蚀力存在2.7 a的波动周期；KEO等[12]利用Kringing内插法得到黄土高原空间连续分布的多年平均降雨侵蚀力分布并讨论其在1960—2011年的变化趋势。然而，已有研究多关注20世纪50年代以后黄土高原降雨侵蚀力变化，对于更长时段(如横跨20世纪百年尺度)的降雨侵蚀力研究依然缺乏，限制了对黄土高原土壤侵蚀长时段变化及其与气候变化关系的理解。

综上所述，本文基于高分辨率地表气候格网数据集的逐月降水数据，评估1901—2016年黄土高原降雨侵蚀力时空变化特征。首先采用逐月实测降水量验证CHELSAcrust数据集的精度，并基于该数据集估算黄土高原1901—2016年降雨侵蚀力；其次，利用Mann-Kendall非参数检验法、经验模态分解法(Empirical Mode Decomposition， EMD)、累积距平、逐像元趋势检验等方法分析黄土高原1901—2016年降雨侵蚀力的时空变化特征，为黄土高原土壤侵蚀侵蚀防治与预报提供参考。

1 材料与方法

1.1 研究区概况

黄土高原地处中国中部偏北，东起太行山，西至乌鞘岭，南接秦岭，北抵长城(100.8°—114.6°E，33.7°—41.3°N)，可分为沙地沙漠区、农灌区、丘陵沟壑区、高原沟壑区、河谷平原区、土石山区(图1)[13]。黄土高原沟壑纵横，地形复杂，总面积约64万km2。该区域属于半湿润、干旱和半干旱区，年均气温4～12℃，年均降水量400～600 mm，降水以短历时暴雨为主，且集中于7—9月[14]。区域内土壤结构松散，抗蚀能力较差，植被覆盖率较低，加之汛期暴雨多发，造成了严重的土壤侵蚀。20世纪70年代以来，黄土高原实施了大规模水土保持措施，尤其自1999年以来，实施了退耕还林(草)工程，有效缓解了区域内的土壤侵蚀。然而，近年来极端降水事件多发[15-16]，为已有水土保持措施带来了挑战。

图1 研究区位置

1.2 数据来源

1901—2016年降水数据来自高分辨率地表气候格网数据集(Climatologies at high resolution for the earth′s land surface areas，CHELSAcrust)。该数据集基于英国东英吉利大学气候研究中心生产的CRU TS 4.01数据集，采用增量更改算法生产而成[17]，提供20世纪以来的气候数据资料，空间分辨率为1 km，数据集获取网址为https:∥chelsa-climate.org/chelsacruts/。此外，由国家气象信息中心编制的《中国地面气候资料日值数据集》，获取了黄土高原内国家气象站1980—2016年的实测降水量数据，用于验证CHELSAcrust数据集的精度，获取网址为：http:∥data.cma.cn/。

1.3 研究方法

1.3.1 降雨侵蚀力估算方法降雨侵蚀力计算方法较多，初期研究多基于次降水强度及降雨动能估算降雨侵蚀力[18-19]。然而，次降水信息较难获取，限制了该类方法的应用。2003年章文波等[20-21]建立了基于不同类型降水资料(如日降水量、月降水量、年降水量)的降雨侵蚀力估算方法，为利用气象站整编降水资料评估降雨侵蚀力提供了方法基础，已在黄土高原上得到广泛应用。本文采用逐月降雨侵蚀力模型计算区域的降雨侵蚀力，计算方法如下:

(1)

(2)

式中:Pi,j为第i年j月的降水量(mm)；N为年数;R为多年平均降雨侵蚀力[MJ·mm/(hm2·h)]；α4，β4为模型参数。模型参数选用逐月雨量的模型参数,即α4=0.1833,β4=1.9957。

1.3.2 降水数据集精度验证方法为验证CHELSA crust数据集精度，选取黄土高原内均匀分布的13个气象站点(图1)，利用1980—2016年逐月降水量观测数据，以Nash效率系数(Nash-Sutcliffe efficiency coefficient，NSE)、决定系数(Coefficient of determination，R2)的结果来评定逐月降水数据集精度。R2取值范围0～1，Nash取值范围为-∞～1，其值越接近1表明模拟精度越高，降雨数据集越可靠，一般认为，R2>0.6，Nash>0.5时，结果可信[22]。

1.3.3 降雨时空分布特征分析方法

(1) 经验模态分解法(EMD)。利用经验模态分解法分析1901—2016年降雨侵蚀力年际周期性。经验模态分解法(Empirical mode decomposition，EMD)是Huang等[23]提出的一种自适应信号时频处理方法。该方法依据数据系列的瞬时特征，将原始时间序列信号按照频率特征，从高到低分解为若干个本征模函数(Intrinsic model function，IMF)，所分解出来的IMF分量包含原信号不同时间尺度的局部特征信号，以反映原始时间序列的周期与振幅[24]。

(2) 非参数Mann-Kendall突变检验。Mann-Kendall突变检验法是检验气候、水文等要素时间序列发生突变的非参数检验方法。该方法通过绘制UF，UB曲线，当UF>0，则表明序列呈现上升趋势，反之同理。给定显著性水平α，若|UF|>Uα，则表明序列存在显著的趋势变化。若UF和UB相交，且交点在临界线Uα之内，则交点对应的时刻可能为突变开始的时间[25]。本文使用该法评估1901—2016年降雨侵蚀力的突变特征，置信水平设置为95%(Uα=±1.96)。

(3) 累积距平法。使用累积距平法探明1901—2016年黄土高原和各个分区降雨侵蚀力的年际变化的趋势和阶段。其计算方法如下：

(3)

(4) 逐像元趋势性检验。通过逐像元一元线性回归来判断1901—2016年黄土高原地区降雨侵蚀力时空变化特征及显著性。使用ArcGIS 10.3.1栅格计算器，计算一元线性回归中各项参数，得到逐栅格降雨侵蚀力线性回归方程，并进行显著性F检验，基于给定的显著性α=0.05，通过查表判断线性变化的显著性，通过ArcGIS 10.3.1重分类工具实现降雨侵蚀力空间变化分布和显著性分布。

2 结果与分析

2.1 降水数据集精度验证结果

基于1980—2016年各气象站点的降水量实测值，验证CHELSAcrust数据集逐月降水量的精度(表1)。结果表明，各气象站点Nash系数最小为0.65，最大为0.91，均值为0.79，R2最大为0.89，最小为0.65，均值为0.82，说明该数据集模拟精度较高，能够满足长时间序列分析需求。

表1 CHELSAcrust数据集逐月降水量精度验证结果

2.2 降雨侵蚀力空间分布特征分析

2.2.1 年均降雨侵蚀力空间分布 1901—2016年间黄土高原年均降雨侵蚀力差距较大，从东南向西北递减(图2)，东部大部分地区降雨侵蚀力超过3 000 MJ·mm/(hm2·h)，而西北地区降雨侵蚀力较小，均在1 000 MJ·mm/(hm2·h)以下。降雨侵蚀力大于4 000 MJ·mm/(hm2·h)的地区主要集中在土石山区和丘陵沟壑区东部。各分区的降雨侵蚀力差异明显，土石山区降雨侵蚀力较大，平均降雨侵蚀力为2 414 MJ·mm/(hm2·h)。农灌区的降雨侵蚀力均值较小，为699.30 MJ·mm/(hm2·h)。不同分区降雨侵蚀力由大到小依次为：土石山区>河谷平原区>丘陵沟壑区>高原沟壑区>沙地沙漠区>农灌区。

在20世纪20年代之前，黄土高原平均降雨侵蚀力为1 233.35 MJ·mm/(hm2·h)，其中土石山区的降雨侵蚀力最大，为5 253.75 MJ·mm/(hm2·h)。降雨侵蚀力<2 000 MJ·mm/(hm2·h)的面积占整个黄土高原面积的80.92%～88.44%，而降雨侵蚀力>4 000 MJ·mm/(hm2·h)的面积仅占总面积的0.34%～1.43%；表明该时期黄土高原的降雨侵蚀力总体维持在较低水平；20世纪30年代开始，黄土高原的降雨侵蚀力激增，尤其是1930—1939年，土石山区以及丘陵沟壑区东部的极强降雨侵蚀力[>5 000 MJ·mm/(hm2·h)]显著增加，面积占比达到3.10%(图2)。1930—1979年，降雨侵蚀力在1 000～3 000 MJ·mm/(hm2·h)的面积占比达到60%以上，降雨侵蚀力>4 000 MJ·mm/(hm2·h)的面积占比为1.99%～6.53%，比1901—1929年增加4.5倍以上，各分区在这一时期的降雨侵蚀力达到峰值，且增大趋势向西北地区的高原沟壑区和沙地沙漠区延伸。20世纪80年代—21世纪初，降雨侵蚀力逐渐下降，降雨侵蚀力<2 000 MJ·mm/(hm2·h)的面积占到整个黄土高原面积的79.04%～88.59%，而降雨侵蚀力>4 000 MJ·mm/(hm2·h)的面积仅占0.57%～0.16%，且降雨侵蚀力较大的地区逐渐退至黄土高原东部的土石山区局部地区。2010—2016年，降雨侵蚀力>4 000 MJ·mm/(hm2·h)面积占比增至3.06%，表明黄土高原降雨侵蚀力高值区域再次增加。

图2 1901-2016年黄土高原降雨侵蚀力时空分布

如图3所示，1901—2016年，黄土高原降雨侵蚀力增加的区域主要集中在黄土高原中部及西部部分区域，减少的区域主要分布在黄土高原西部和东北部。降雨侵蚀力的显著性变化在空间分布上呈现出边缘地区多为非显著变化而中部地区多为显著变化的规律。其中，丘陵沟壑区中部、沙地沙漠区中部、高原沟壑区西部和东部等区域的降雨侵蚀力显著增长；沙地沙漠区东部和高原沟壑区部分地区的降水侵蚀力显著减少。

图3 1901-2016年黄土高原降雨侵蚀力变化趋势

2.3 降雨侵蚀力时间演变特征

2.3.1 降雨侵蚀力年际变化 1901—2016年，黄土高原平均降雨侵蚀力呈非显著下降趋势(图4)。年均降雨侵蚀力为1 462.17 MJ·mm/(hm2·h)，其中，1925年降雨侵蚀力为研究时间序列最大值，为3 497.3 MJ·mm/(hm2·h)，较平均值高出约139.2%。1991年降雨侵蚀力出现最小值，为545.4 MJ·mm/(hm2·h)，相较于平均值减少了约62.7%。对各分区而言，土石山区降雨侵蚀力明显高于其他五大分区，沙地沙漠区及农灌区的降雨侵蚀力相比之下较小。沙地沙漠区和河谷平原区的降雨侵蚀力非显著上升，其他4个分区降雨侵蚀力非显著下降。各分区降雨侵蚀力随时间的变化相似，在1901—1916年期间，各区降雨侵蚀力维持在较低水平[<4 000 MJ·mm/(hm2·h)]，1917年降雨侵蚀力激增，丘陵沟壑区达到峰值5 242.12 MJ·mm/(hm2·h)；1917—1960年，降雨侵蚀力明显增大且各分区较大降雨侵蚀力集中出现，1939年土石山区达到峰值9 001.02 MJ·mm/(hm2·h)，1940年沙地沙漠区达到峰值2 037.20 MJ·mm/(hm2·h)，1949年河谷平原区和高原沟壑区均达到峰值，分别为5 072.44，2 497.41 MJ·mm/(hm2·h)，1959年农灌区达到峰值2 030.51 MJ·mm/(hm2·h)。1961—2000年间，降雨侵蚀力波动下降，除河谷平原区其他各区均达到谷值，沙地沙漠区、农灌区、丘陵沟壑区谷值在1965年，高原沟壑区出现在1972年，土石山区出现在1997年。2010年后降雨侵蚀力出现反弹，2013年出现峰值。

黄土高原和各分区的年降雨侵蚀力累积距平结果变化趋势类似(图5)，可将降雨侵蚀力变化划分为1901—1930年、1931—1980年、1980—2016年3个阶段，1930年前累积距平值为负，说明降雨侵蚀力低于平均值。自1930年起，降雨侵蚀力与多年均值差值为正，累积距平值呈现持续上升趋势，到1957年，累积距平值达到0，之后持续上升，直到1979年达到峰值。1980年开始，累积距平值逐步下降，表明降雨侵蚀力低于多年平均水平。各分区的降雨侵蚀力距平与黄土高原降雨侵蚀力累积距平的变化趋势基本一致。值得注意的是，黄土高原及各分区的降雨侵蚀力累积距平变化于2010年出现较大反复，说明2010年之后的降雨侵蚀力有增强迹象。

利用Mann-Kendall法对黄土高原及6个分区进行突变诊断(图6)，整个黄土高原降雨侵蚀力UF和UB序列均存在多个交叉点，6个分区的降雨侵蚀力与黄土高原的UF和UB曲线波动基本一致，交点大多分布在20世纪20年代以前和90年代以后。结合年降雨侵蚀力变化趋势(图4)及年降雨侵蚀力累积距平(图5)，说明黄土高原的降雨侵蚀力不存在明显突变。黄土高原及6个分区降雨侵蚀力的正序列在1930年之前均表现为先升后降，UF值在0刻度线上下浮动，无持续性趋势；1930年之后，黄土高原、高原沟壑区及河谷平原区的正序列曲线表现为先升后降，但UF值始终>0，表明降雨侵蚀力在增加，在1980年之后增加幅度减小，且均在60—80年代左右突破0.05显著性水平线，说明在1960—1980年降雨侵蚀力增加较显著。土石山区与以上3个区域相似，但在2007年之后UF值在0值上下波动，说明2007年之后该区降雨侵蚀力反复变化但不显著。农灌区、沙地沙漠区、丘陵沟壑区在2000年之后UF<0且未超过置信区间，表明该时段内降雨侵蚀力非显著减小。

2.3.2 降雨侵蚀力周期性分析利用EMD法对黄土高原地区1901—2016年的116 a降雨侵蚀力进行分析，得到了5个本征模分量(IMF)及最低频分量(RES)(图7)。表2计算了各个IMF分量、方差以及方差贡献率。可以看出5个本征模分量(IMF)累计方差贡献达到96.03%，而其中IMF1贡献率最大，为56.53%，其次为IMF2分量，为32.77%。IMF5分量的贡献率最小，为1.48%。由此得出，黄土高原地区降雨侵蚀力变化存在周期性规律，2.75 a变化周期最显著，其次为5.33 a，再其次为10.29，24.5，44.5 a，这3个周期的方差贡献率差异不明显，均为弱周期。降雨侵蚀力的周期变化受制于该区域气候要素的波动周期。太阳黑子和大气环流的异常活动可以影响水热平衡，从而改变降水时序和空间分布。研究表明，太阳黑子存在11 a的活动周期[26]，这与本研究得出的本征模分量IMF3(10.29 a)的周期基本吻合。大气环流异常因子厄尔尼诺现象也存在2～7 a的变化周期[27]，本研究所得出的本征模分量IMF1(2.75 a)，IMF2(5.33 a)的周期与其基本一致。

注：△为降雨侵蚀力最大值；▽为降雨侵蚀力量小值。

图5 1901-2016年黄土高原及各分区年均

3 结论

(1) CHELSAcrust数据集降水量数据精度较高，可用于黄土高原长时段研究需求。

(2) 黄土高原降雨侵蚀力总体为东部大于西部。黄土高原地区各个分区降雨侵蚀力差距较大，不同分区降雨侵蚀力由大到小依次为：土石山区>河谷平原区>丘陵沟壑区>高原沟壑区>沙地沙漠区>农灌区。1901—2016年，降雨侵蚀力增加的区域主要集中在黄土高原中部及西部部分区域，减少的区域主要分布于黄土高原西部和东北部。降雨侵蚀力的显著性变化在空间分布上呈现出边缘多为非显著变化而中部多为显著变化的规律。

(3) 1901—2016年黄土高原地区年平均降雨侵蚀力约为1 462.17 MJ·mm/(hm2·h)。时间序列上6个分区与黄土高原的变化趋势基本一致。研究时段内，黄土高原降雨侵蚀力变化不显著且无明显突变点，降雨侵蚀力变化总体可分为1901—1930年、1930—1980年和1980—2016年3个阶段。