南渡江流域降雨侵蚀力时空分布与变化趋势研究

卢旭东,齐 实,陈家东,郭建超,张 林,周 飘

(1.北京林业大学水土保持学院,北京 100083;2.海南省水利水电勘测设计研究院有限公司,海南 海口 571100)

水土流失会造成土壤养分、有机质、持水能力和渗透速率降低,威胁着生态系统平衡发展,已经成为全球最为严重的生态环境问题之一[1-2]。水土流失的影响因素主要包括降雨、土壤、植被、地形和人类活动,其中,降雨是诱发水土流失最根本的驱动力[3]。降雨加剧水土流失的方式主要为雨滴击溅和径流剥蚀,因此,降雨时间、强度和雨型特征可直接影响水土流失发生的程度和风险。

降雨侵蚀力作为通用土壤流失方程(USLE)及其修正方程(RUSLE)中估算土壤侵蚀的重要参数,是表征降雨侵蚀能力的直接指标,可用来定量反映降雨引发土壤侵蚀和水土流失的潜在风险[4-5]。降雨侵蚀力是基于侵蚀性降雨量由模型计算得到的,最初的计算模型由WISCHMEIER等[6]提出,模型计算公式为E·I30,其中,E为降雨动能,MJ·hm-2;I30为最大30 min降雨强度,mm·h-1。该模型计算结果精度较高,但长序列降雨数据需求通常难以得到满足,同时数据处理步骤较为繁琐、耗时[7-8]。因此,在该模型的基础上发展出不同时间尺度的简易降雨侵蚀力模型[9-10],可利用气象站常规观测数据实现降雨侵蚀力的计算。我国降雨侵蚀力研究从20世纪80年代开始,孙保平等[11]建立以年降雨量为基础数据的多年平均降雨侵蚀力模型并在宁夏西吉县成功应用;黄炎和等[12]提出月雨量模型;章文波等[4]提出日雨量模型,该模型数据获取的可行性得到大幅度提升,并简化了降雨侵蚀力计算过程,同时相较于年、月雨量模型,又可以保证多年平均降雨侵蚀力及其季节分布的结果精度。

海南岛位于中国最南端,地处热带,具有全年雨量丰富、雨季持续时间长和降雨强度大等气候特点,这些气候因素加大了该地区水土流失潜在风险[13]。南渡江流域是海南岛第一大流域,南渡江为海南省北部及东北部地区的生产和生活提供了主要水资源,同时上游径流汇入的松涛水库是海南西部重要的饮用水和农业灌溉用水水源[14],该流域水土流失动态变化直接关系着海南省北部、西部及东北部大部分地区水资源和水环境。鉴于此,利用南渡江流域及其邻近13个气象站1971—2020年逐日降雨数据,采用Mann-Kendall非参数趋势/突变检验、小波周期分析和反距离加权插值等方法,分析南渡江流域降雨侵蚀力时空分布特征及其变化趋势,为南渡江流域及其典型土壤侵蚀区风险性预测、评估及预警提供重要科学依据。

1 材料与方法

1.1 研究区概况

南渡江流域(19°09′～19°55′ N、109°36′～110°34′ E)位于海南省中北部(图1),流域总面积为5 333 km2,海拔范围为2～1 379 m,平均海拔高度为116 m,地势西南高东北低,这种地势造就了南渡江河流呈现自西向东的走向。

图1 南渡江流域地形及气象站位置分布

流域地处热带北部边缘,属于热带季风气候区,干湿两季明显,台风频繁,降雨丰富,流域1971—2020年年均降水量为2 137 mm,雨季(5—10月)和旱季(11—翌年4月)降水量分别占81%和19%,南部地区降雨量多于北部地区[15]。流域干流南渡江丰水少沙,此外集水面积大于100 km2的支流有20条,包括15条一级支流和5条二级支流。流域土壤类型有13种,以麻黄砂土为主,其后依次为淡麻砖土、肉泥田和湖光岩焦灰土。流域内植被覆盖度高(覆盖度≥75%),一旦被大规模开荒等人为扰动,热带雨林遭到破坏,容易造成水土流失加剧。

1.2 数据来源

气象数据来源于海南省水文水资源勘测局提供的南渡江流域及其周边13个气象站点逐日降雨量数据。各气象站点基本信息见表1。由于不同站点数据起始记录时间不同,部分站点存在数据缺测情况,所有数据均遵循欧洲气候评估数据集制作标准[16]进行筛选:(1)数据时序不少于40 a;(2)单个站点缺失数据占比不超过10%;(3)每个站点每年缺失数据不超过20%或连续3个月。对于单个站点局部缺失数据,采用临近2个及以上站点观测数据平均值进行补值,以得到完整序列降雨数据。

表1 南渡江流域气象站基本信息

1.3 研究方法

1.3.1降雨侵蚀力的计算

由于研究区缺少30 min高分辨率降雨数据,因此选择日降雨量模型[4]计算降雨侵蚀力,该模型已在中国的国家和地区尺度上得到广泛应用[17-18]。模型计算公式为

(1)

式(1)中,R半月,i为第i个半月内降雨侵蚀力,MJ·mm·hm2·h-1·a-1;k为某半月内侵蚀性降雨日数,d;Pk为侵蚀性降雨量,即半月内第k天≥12 mm的日雨量,mm;α、β为模型参数,能够反映不同区域降雨特征,利用日降雨资料进行估算。

α=21.586β-7.189 1,

(2)

(3)

式(2)～(3)中,Pd,12为日雨量≥12 mm的日平均雨量,mm;Py,12为日雨量≥12 mm的年平均雨量,mm。通过统计各半月降雨侵蚀力可得到流域月降雨侵蚀力、年降雨侵蚀力和多年平均降雨侵蚀力。

1.3.2时间变化分析

采用Mann-Kendall法和小波分析2种方法对时间序列的趋势变化、突变年份和周期变化进行分析。

(1)Mann-Kendall法。该方法在分析时间序列的变化趋势时,可以确定时间序列(X1,X2,…,Xn)及变化趋势是否显著。对于任何要检查的序列Xt(t=1,2,…,n),n为要检查的序列长度。统计量S定义为

(4)

式(4)中,Xj和Xk为时间序列对应的年份数据;n为时间序列长度;sgn(Xj-Xk)为符号函数。当Xj>Xk时,sgn值为1;当Xj

E(S)=0,

(5)

var(S)=n(n-1)(2n+5)/18。

(6)

式(5)～(6)中,E(S)为均值函数;var(S)为方差函数。

标准化检验统计量Z计算公式为

(7)

在双侧检验中,对于给定显著性水平α,如果|Z|≥Z1-α/2,则原假设不可接受,满足显著性水平,序列呈显著上升或下降趋势。该研究α为0.05,即Z1-α/2=1.96。当Z>0时,呈现上升趋势;当Z<0时,呈下降趋势。Mann-Kendall检验也可以通过正序列UF和逆序列UB统计量判断突变开始时间,UF和UB曲线的交点即为突变时间点[19]。

(2)小波分析。小波分析是一种时域和频域的分析方法,可以在时间变化上分析降雨侵蚀力的周期性变化规律。连续小波基本公式为

(8)

(9)

式(9)中,var(a)为小波方差。

小波方差随a变化的过程被称为小波方差图,小波方差图可以用来确定波动的能量强度和时间尺度分布。因此,小波方差图可以用来确定水文序列中的主时间尺度,即主周期。小波分析流程采用MATLAB计算小波系数和小波系数的实部,并绘制小波等值线图和方差图,由图中丰枯的交替来判断周期变化规律[20]。

1.3.3空间变化分析

空间变化分析主要采用反距离加权插值法进行。反距离加权插值法的原理是彼此相近的事物相较于远处的事物更具有相似性。因此,在计算时,距离观测者较近的观测物被赋予较大权重,较远的观测物则相反。反距离加权插值法基本公式为

(10)

(11)

(12)

式(10)～(12)中,Ri为第i点降雨侵蚀力,MJ·mm·hm-2·h-1;Rj为第j个气象站降雨侵蚀力计算值,MJ·mm·hm-2·h-1;n为气象站点数;λj为第j个气象站的权重;p为指数值,该研究中其值为2;dji为点j与气象站点i之间的距离,km。

采用ArcGIS 10.8对南渡江流域年侵蚀性降雨量和年降雨侵蚀力进行空间插值,对得到的空间插值分布图进行分析。

2 结果与分析

2.1 降雨侵蚀力时间分布特征

2.1.1降雨侵蚀力年际分布

南渡江流域1971—2020年多年平均降雨侵蚀力为16 497.67 MJ·mm·hm-2·h-1。如图2所示,50 a间,流域降雨侵蚀力最大值和最小值分别为23 692.14 和11 841.33 MJ·mm·hm-2·h-1,分别出现在2000和1983年;多年平均侵蚀性降雨量为1 692.41 mm,最大值为2 217.36 mm,最小值为1 392.67 mm,极值出现年份与降雨侵蚀力一致,这验证了侵蚀性降雨量对降雨侵蚀力的关键作用。

图2 南渡江流域侵蚀性降雨量与降雨侵蚀力年际变化

降雨侵蚀力最大值为最小值的2.01倍,侵蚀性降雨量最大值为最小值的1.59倍,南渡江流域年降雨侵蚀力和年侵蚀性降雨量的变异系数分别为0.29和0.21,均属于中等变异,但年降雨侵蚀力变异程度要大于年侵蚀性降雨量。这表明与侵蚀性降雨量相比,降雨侵蚀力在研究期间的年际变化更为剧烈。

图2显示,流域年降雨侵蚀力的变化趋势与年侵蚀性降雨量大体一致。1971—2020年流域降雨侵蚀力表现为波动上升变化趋势,总体分为4个阶段:1971—1983年为第1阶段,降雨侵蚀力表现为波动下降趋势;1984—2001年为第2阶段,降雨侵蚀力呈波动上升趋势;2002—2013年为第3阶段,降雨侵蚀力在2005年发生陡降后,转为上升趋势,整体则呈波动上升趋势;2014—2020年为第4阶段,降雨侵蚀力呈下降趋势。

2.1.2降雨侵蚀力年内分布

图3显示,南渡江流域降雨侵蚀力主要集中在7—9月,占全年降雨侵蚀力的51.26%;降雨侵蚀力最大值出现在8月,为2 789.75 MJ·mm·hm-2·h-1,占全年的17.20%;最小值出现在1月,为221.01MJ·mm·hm-2·h-1,仅占全年的1.36%,月季差异显著。

图3 南渡江流域月降雨侵蚀力分布

降雨侵蚀力季节分布变化见图4。除春季(3—5月)降雨侵蚀力呈现下降变化趋势外,夏季(6—8月)、秋季(9—11月)和冬季(12—翌年2月)均呈上升变化趋势。降雨侵蚀力主要集中在夏季和秋季,夏季降雨侵蚀力最大,为6 877.23 MJ·mm·hm-2·h-1,占全年降雨侵蚀力的42.28%;冬季降雨侵蚀力最小,为406.76 MJ·mm·hm-2·h-1,占全年的2.50%。

图4 南渡江流域降雨侵蚀力季节变化

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

2.2.1年降雨侵蚀力空间分布

各气象站点年均降雨侵蚀力变化范围为12 340.72～22 890.31 MJ·mm·hm-2·h-1,最大值出现在加报站,最小值出现在加来站,最大值为最小值的1.85倍,各站点降雨侵蚀力存在较大差异(图5)。

图5 南渡江流域各气象站点多年平均降雨侵蚀力

利用南渡江流域1971—2020年平均侵蚀性降雨量和降雨侵蚀力数据进行反距离加权插值得到流域年侵蚀性降雨量和降雨侵蚀力空间分布情况(图6)。流域年降雨侵蚀力空间分布特征表现为从东南部和西南部的高值区,向北部的低值区降低,整体呈从南向北逐渐降低趋势,这与年侵蚀性降雨量分布规律一致。根据降雨侵蚀力大小可将南渡江流域大体分为3个区域:南部及东南部包括加报、大丰、南坤、大陆坡和塘尾站,降雨侵蚀力基本达到17 000 MJ·mm·hm-2·h-1以上;中部地区包括南丰、昆仑、定安和新德站,降雨侵蚀力大致分布在15 000～17 000 MJ·mm·hm-2·h-1;北部及西北部地区包括加来、加潭、美亭和龙塘站,降雨侵蚀力大致分布在15 000 MJ·mm·hm-2·h-1以下。

图6 南渡江流域侵蚀性降雨量和降雨侵蚀力空间分布

2.2.2季节降雨侵蚀力空间分布

降雨侵蚀力在不同季节存在较大差异,对南渡江流域各季节降雨侵蚀力进行空间特征分析。如图7所示,各季节降雨侵蚀力空间分布整体与年降雨侵蚀力一致,呈现由南部向北部逐渐递减的趋势。其中,夏季降雨侵蚀力占全年比例最大,不同地区降雨侵蚀力变化范围为5 049.13～10 706.87 MJ·mm·hm-2·h-1,降雨侵蚀力最高值出现在南部和东南部地区,最低值为西北部地区;秋季降雨侵蚀力占全年比例仅次于夏季,不同地区降雨侵蚀力范围为3 765.97～9 058.60 MJ·mm·hm-2·h-1,极值区域分布与夏季一致;降雨侵蚀力全年占比最小的季节为冬季,降雨侵蚀力变化范围为154.43～634.65 MJ·mm·hm-2·h-1,最高值在东南部地区,最低值在西北部地区。

图7 南渡江流域季节降雨侵蚀力空间分布

2.3 降雨侵蚀力变化特征

2.3.1降雨侵蚀力时间变化分析

由南渡江流域降雨侵蚀力小波分析结果(图8)可知,研究区多年平均降雨侵蚀力在整个时间序列演化过程中,显著存在4～7和28～33 a 2个变化周期,且在1971—2020年稳定分布,具有全域性。小波方差图峰值为降雨侵蚀力序列演变的周期值。如图8所示,小波方差图峰值分别为6和30 a,且30 a峰值最高,说明年均降雨侵蚀力序列在30 a左右周期震荡最强,因此,30 a为流域降雨侵蚀力主变化周期,且在主周期内存在6 a的次变化周期。

图8 基于小波分析的年降雨侵蚀力周期变化

由Mann-Kendall突变检验结果(图9)可知,降雨侵蚀力正序列UF和逆序列UB曲线在置信区间内存在4个交点,最早相交于1971—1972年,交点处UF>0,但UF和UB曲线交点均位于±1.96内,即|Z1-α/2|<1.96,未满足0.05显著性水平,表明降雨侵蚀力在研究期限内未发生显著突变。

图9 年降雨侵蚀力Mann-Kendall突变分析

2.3.2降雨侵蚀力空间变化分析

通过对降雨侵蚀力变异系数CV值和Mann-Kendall统计量Z值进行反距离加权插值,得到南渡江流域降雨侵蚀力年际变化特征参数的空间分布情况(图10)。各气象站变异系数范围为0.24～0.43,地区差异比较明显,北部地区变异系数相对较大,南部地区变异系数相对较小,整体呈现为由北向南逐步递减趋势。流域南部地区降雨量丰富,年侵蚀性降雨量稳定,降雨侵蚀力年际变化相对较小,变异系数均小于0.28。由降雨侵蚀力Mann-Kendall统计量Z值分布可知,研究区除大丰、南丰和龙塘站降雨侵蚀力呈现降低趋势外,其余站点均呈现上升趋势,在呈上升趋势的地区需要加强水土保持措施建设,防止水土流失加剧。

图10 南渡江流域降雨侵蚀力变异系数及Mann-Kendall Z值空间分布

2.4 不同地形区降雨侵蚀力变化

为进一步分析地形对降雨侵蚀力的影响,根据南渡江流域地形特征将其划分为东、北部平原区和西、南部山地区2大类地形区(以400 m海拔为划分标准),并统计东、北部平原区和西、南部山地区降雨侵蚀力分布及其与侵蚀性降雨、经纬度和海拔的相关关系(表2)。

表2 南渡江流域不同地形区降雨侵蚀力及其影响因素

如表2所示,西、南部山地区侵蚀性降雨量、降雨侵蚀力均高于东、北部平原区;变异系数由高到低为东、北部平原区>全流域>西、南部山地区,西、南部山地区降雨侵蚀力较东、北部平原区更加集中稳定。东、北部平原区Z值为负数,该区域降雨侵蚀力呈下降趋势,西、南部山地区Z值为正数,降雨侵蚀力呈上升趋势,且西、南部山地区|Z|值比东、北部平原区大,说明其变化趋势更为明显。东、北部平原区和西、南部山地区降雨侵蚀力与侵蚀性降雨量和海拔均呈显著正相关,分别通过0.01和0.05显著性检验;东、北部平原区和西、南部山地区降雨侵蚀力与经度相关不明显,与纬度存在一定相关性,相关系数分别为0.76和0.42,分别通过0.01和0.05显著性检验。

3 讨论

降雨侵蚀力作为土壤侵蚀模型的关键参数,量化降雨侵蚀力的时空分布及其变化对模型应用以及准确评估水土保持措施的效果和预警水土流失发生风险至关重要[21]。笔者研究了1971—2020年南渡江流域降雨侵蚀力时空分布及其变化趋势,发现南渡江流域年均降雨侵蚀力为16 497.67 MJ·mm·hm-2·h-1,参考文献[22-23]分类标准,属于高降雨侵蚀力水平区域。南渡江流域降雨侵蚀力高于珠江流域[24]、赣江流域[25]和松花江流域[16]等诸多区域,主要原因在于各区域所在纬度和季风带不同。南渡江流域年内降雨侵蚀力峰值出现在7—9月,这是由于流域属于热带北部边缘地区,受热带季风气候影响,夏、秋两季台风、暴雨频发,造成侵蚀性降雨最为集中[13]。

1971—2020年流域降雨侵蚀力表现为波动上升变化趋势,这一结果与中国大陆降雨及降雨侵蚀力侵蚀力研究[21,26]中南部地区降雨侵蚀力变化趋势保持一致。流域侵蚀性降雨量与降雨侵蚀力空间分布相似,这与泾河流域[27]和大汶河流域[28]的研究结果一致,进一步说明侵蚀性降雨是影响降雨侵蚀力的关键因素。长期降雨侵蚀力数据是一项稳定的环境指标,可以反映降雨侵蚀力未来情景,对水土保持和区域环境治理具有重要参考价值。

由于地理位置和降雨特征的差异,降雨侵蚀力在不同地形地貌区呈现出不同变化趋势。许功伟等[28]对大汶河流域降雨侵蚀力时空变化的研究发现,大汶河流域降雨侵蚀力与高程存在较强相关性,东、北部山地地区降雨侵蚀力显著高于西、南部平原地区;此外,王佳欢等[29]通过建立白洋淀流域降雨侵蚀力与地形地貌指标的回归方程,发现区域降雨侵蚀力与纬度、海拔呈负相关关系,与经度呈正相关关系。笔者研究中南渡江流域降雨侵蚀力与海拔相关性最高,其次为纬度,与其他地区存在差异,这可能是大气环流和地理分布差异造成的影响。南渡江流域地貌差异较大,整体呈现西南高东北低的趋势,热带季风从东部由海洋吹向陆地,造成流域北部降雨侵蚀力由东向西逐渐减小,但流域南部地区由于地势抬高,使得暖湿气流上升,导致降雨量和降雨侵蚀力增加,同时这一地区地形起伏大,水土流失风险更大,应加强西、南山地区水土流失预防和治理。

降雨侵蚀力反映降雨对土壤侵蚀的潜在能力,但其值大小并不直接表示该地区水土流失强度大小,因为人为因素、土壤质地和植被因素均影响着水土流失的发生。南渡江流域处于热带岛屿,台风和暴雨的频繁发生也是导致土壤流失的主要因素[13],因此,需要关注气候变化条件下暴雨事件增多,进而导致降雨侵蚀力高值频现[30],这会进一步增加区域水土流失发生风险。此外,受研究区降雨数据限制,笔者选择日降雨侵蚀力模型进行计算,在今后研究中,可以考虑加强对数据的收集及整理,提高数据精度,利用降雨强度、降雨动能等数据对研究区降雨侵蚀力进行更深入对比研究。

4 结论

(1)南渡江流域1971—2020年年均降雨侵蚀力范围为11 841.33～23 692.14 MJ·mm·hm-2·h-1,均值为16 497.67 MJ·mm·hm-2·h-1,降雨侵蚀力存在30 a的主变化周期,在主周期内存在6 a的次变化周期,研究期内降雨侵蚀力未发生显著突变,与侵蚀性降雨量相比,降雨侵蚀力年际变化更为剧烈。

(2)南渡江流域降雨侵蚀力在7、8和9月最为突出,占比分别为14.59%、17.20%和15.22%,这3个月降雨侵蚀力分布集中且稳定;4个季节中,除春季外,夏季、秋季和冬季降雨侵蚀力均呈上升变化趋势。

(3)南渡江流域降雨侵蚀力空间分布特征与侵蚀性降雨量分布特征基本一致,呈现自南部区域向北部区域逐步递减的趋势;降雨侵蚀力高值和低值中心均分别与侵蚀性降雨量高值和低值中心基本吻合,各站点降雨侵蚀力变化趋势存在明显变异,空间上整体呈现从由南向北逐步递减的趋势。

(4)南渡江流域降雨侵蚀力与侵蚀性降雨量、海拔呈显著正相关,且分别通过0.01和0.05显著性检验;降雨侵蚀力与经度相关性不明显,与纬度存在一定相关性,相关系数分别为0.76和0.42,且分别通过0.01和0.05显著性检验。