基于GIS的白洋淀流域降雨侵蚀力时空分布特征分析

王佳欢, 张晨星, 杨新兵, 杨文姬, 米秋菊

(1.河北农业大学 林学院, 河北 保定 071000; 2.山合林(北京)水土保持技术有限公司, 水利部沙棘开发管理中心(水利部 水土保持植物开发管理中心), 北京 100038; 3.河北金洋环境科技有限公司, 河北 保定 071000)

土壤侵蚀作为当前我国面临的众多生态问题中最严重的一项，其引发因素包括降雨、地形、土壤、植被等，其中降雨引起的土壤水蚀是诱发土壤侵蚀的主要驱动力[1]。降雨侵蚀在于地表产生的点式击溅以及径流冲刷产生的综合效应，除降雨的时间长短、强度大小、次数多少外，坡面径流的流速、流量以及流态也会对降雨侵蚀的综合效应产生影响[2]。关于土壤侵蚀的研究，国内外学者致力于土壤侵蚀模型的构建，1965年国外学者在通用土壤流失方程(USLE)应用中建立以小时为单位降雨强度与动能的数学关系，其结果可表征降雨的侵蚀强度，因广泛应用于侵蚀计算而被定义为经典算法[3]。由于经典算法所需数据获取与记录难度较大，为简化算法，学者们开始致力于新模型的研究[4]，就国内而言，在空间尺度上，王万中[5-6]、章文波[7]、刘赋涛[8]等先后研究了中国降雨侵蚀力的计算与分布特征；在时间尺度上，赖成光[9]、Huang[10]等分别对珠江流域和长江流域的降雨侵蚀力时间特征进行了研究。

2017年国家成立雄安新区，提出“绿色生态，宜居新城”的生态定位，在这种情况下，防洪建设、水土保持、水资源调节等一系列生态工程成为新区规划建设的重点，由于新区地处白洋淀流域腹地，上游河流均汇聚于此，生态工程建设无一不受制于大环境，因此，白洋淀流域生态建设成为新区发展的必要前提。高彦春[11]、温静[12]、张晨星[13]等先后对白洋淀流域环境因子变化规律、景观特征及其动态变化进行了研究，为流域治理提供决策性支持。目前，涉及白洋淀流域的研究多集中于生态涵养及景观变化方面，关于降雨侵蚀的研究鲜有报道，针对白洋淀流域生态战略地位以及对于打造雄安新区所发挥的重要生态价值，本文应用日雨量模型开展白洋淀流域降雨侵蚀力研究，通过分析降雨侵蚀力空间、时间分布特征，阐明其分布原因，从而揭示其变化规律，旨在为流域内开展土壤水蚀研究、地质灾害预警、水土保持建设提供科学参考，为雄安新区生态建设提供有力保障。

1 材料与方法

1.1 研究区概况

白洋淀流域位于中纬度地区(北纬38°01′—40°04′、东经113°40′—116°30′)，地跨京冀晋三地，为海河流域的二级流域，总面积33 877.64 km2。流域自西北向东南倾斜，西部为山地丘陵，区间山峦起伏，沟深谷远，东部为平原，地势平坦，高程跨度0～2 781 m。该流域地处温带季风气候区，大陆性较强，多年平均气温为13.3℃，1月份气温较低，平均气温达-4.5℃，7月份气温较高，平均气温达26.7℃，年均降雨量约500 mm，多集中于6—8月份，全年无霜期138～192 d。白洋淀素有“九河下梢”之称，源于大清河水系的9条河流均经淀区汇入渤海。土壤类型多样，以潮土、褐土为主，主要植被类型为温带落叶阔叶林、针叶林以及针阔混交林等，优势树种包括侧柏(Platycladusorientalis(L.) Franco)、油松(PinustabuliformisCarr.)、白桦(BetulaplatyphyllaSuk.)、山杨(PopulusdavidianaDode)、槲树(QuercusdentataThunb.)等。流域内交通发达、道路完备，其67%地区属保定辖区，人口多分布于京广线以东地区，2020年保定市生产总值3.4×103亿元。

1.2 研究方法

降雨数据源于中国气象数据网(www.data.cma.cn)，总共涉及流域内及周边气象站点64个，选取2003—2018年的逐日降雨数据。DEM高程数据源自地理空间数据云平台(www.gscloud.cn)。各站点位置见图1。

1.2.1 降雨侵蚀力计算 降雨侵蚀力计算方法源于日雨量模型[7]，具体计算方法如下：

(1)

式中:R半月为侵蚀力半月值;m为半月;Pk为半月内第k天12.0 mm以上的日雨量;m为半月内日雨量12.0 mm以上的天数;a与β均为日雨量模型参数。

α= 21.586β-7.1891

(2)

β= 0.836 3 + 18.144/Pd12+ 24.455+Py12

(3)

式中:Pd12为12 mm以上的日雨量平均值;Py12为日雨量12 mm以上的年降雨平均值,根据公式(1—3)可计算得到年降雨侵蚀力值。

1.2.2 降雨侵蚀力空间分布 基于DEM数据和《河北省水利志》，采用ArcGIS 10.2软件中Spatial Analyst模块的填洼、流向、流量、流域出水口、分水岭、提取边界线等水文分析工具，提取白洋淀流域边界。通过对64个站点降雨侵蚀力插值分析得到空间上连续分布的栅格数据，根据生成的白洋淀流域边界进行掩膜提取(栅格输出分辨率30 m×30 m)，生成降雨侵蚀力空间分布图。计算得到的气象站点降雨侵蚀力采用SPSS 24.0进行线性逐步回归分析，输入因素包括纬度、经度、海拔。

1.2.3 降雨侵蚀力时间分布 时间分布包括年内分布和年际分布，年内分为春季(3—5月)、夏季(6—8月)、秋季(9—11月)和冬季(1月、2月和12月)，有研究表示，白洋淀流域存在周期为2～6年的降水变化规律[14]，结合前人研究成果，年际变化以4年为一周期，即2003—2006年、2007—2010年、2011—2014年、2015—2018年，通过倾向率[15]表示16年的年内变化和年际变化(具体计算方法如下)，采用ArcGIS克里金插值功能将分析结果反映到图层上。

(4)

式中:yi为降雨侵蚀力值；xi为时间；a为降雨侵蚀力倾向率；b为计算所得回归系数。

1.2.4 Mann-Kendall突变检验 Mann-Kendall法常用于检验时间序列下降雨突变特征显著情况，具体计算方法如下：

x1，x2，…，xn为气候序列，Sk为样本i的xi>xj累计数(i≥j≥1)，定义统计量：

(5)

Sk均值与方差(假定随机时间序列独立)：

E(Sk)=k(k-1)/4

(6)

var(Sk)=k(k-1)(2k+5)/72， 1≤k≤n

(7)

Sk标准化：

(8)

式中:UF1=0,反序列中UBk=-UFk(k=xn,xn-1,…，x1),UB1=0。

根据UF和UB绘制趋势图，气候序列上升时UF或UB大于0，气候序列下降时UF或UB小于0，UF与UB交点表示气候序列发生突变的时间，交点值越过1.96水平线则上升趋势显著，交点前后所在区间表示突变区域，突变区域越过1.96水平线则突变显著。

2 结果与分析

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

由图2可知，2003—2018年白洋淀流域年降雨侵蚀力范围为909.62～3 358.38 (MJ·mm)/(hm2·h)，平均2 284.54 (MJ·mm)/(hm2·h)，空间分布上呈东南和西北偏低、西南和东北偏高趋势，从东南向西北延伸形成一条廊道分割两个降雨极点(西南极点为阜平县；东北极点为霞云岭)，廊道上降雨侵蚀力呈先增后降趋势，降雨侵蚀力极点位于望都县、唐县和定州市范围内，其中，西南地区降雨侵蚀力大于2 800 (MJ·mm)/(hm2·h)，东北地区降雨侵蚀力大于3 000 (MJ·mm)/(hm2·h)，而东南—西北方向上降雨侵蚀力范围为909.62～2 800 (MJ·mm)/(hm2·h)。结合图3分析，白洋淀流域高程范围为0～2 781 m，东南地区地势平坦，高程跨度0～100 m，西北地区峰峦叠嶂，高程跨度400～2 781 m，东南平原与西北山地交界的海拔起伏带多为低山丘陵，高程跨度100～400 m，降雨极点多分布于低山丘陵地带及其紧邻的平原区，东南地区平原地带高于西北地区山岳地带，降雨侵蚀力在东南—西北高程变化梯度上先增后降。

图2 研究区降雨侵蚀力空间分布

图3 研究区DEM分布

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

由表1可知，单因子输入中经度和海拔的降雨侵蚀力回归方程均存在极显著相关关系(p<0.001)，而纬度与降雨侵蚀力没有明显显著关系(p>0.05)，双因子输入中纬度、经度、海拔任意双因子的降雨侵蚀力回归方程都达极显著相关关系(p<0.001)，三因子输入中因子与降雨侵蚀力同样具有极显著相关关系(p<0.001)。根据回归方程R2值判断，经度和海拔与回归方程偏离较小，纬度偏离较大，随着输入因子的增多，R2值越大，说明多因子可更好地与降雨侵蚀力产生关联。

表1 降雨侵蚀力线性逐步回归分析

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

2.3.1 年内分布特征 由图4可知，白洋淀流域春季(3—5月)降雨侵蚀力范围为83.91～262.03 (MJ·mm)/(hm2·h)，平均190.94 (MJ·mm)/(hm2·h)，东南边缘地区降雨侵蚀力较大，其他地区较小。夏季(6—8月)降雨侵蚀力范围为650.77～2 719.00 (MJ·mm)/(hm2·h)，平均1 826.75 (MJ·mm)/(hm2·h)，西北边缘地区降雨侵蚀力较大，其他地区较小。秋季(9—11月)降雨侵蚀力范围为184.43～445.29 (MJ·mm)/(hm2·h)，平均325.75 (MJ·mm)/(hm2·h)，东南地区、西北地区和中部地区降雨侵蚀力较大，西南和东北地区较小。冬季(1月、2月和12月)降雨侵蚀力范围为0.16～3.91 (MJ·mm)/(hm2·h)，平均1.77 (MJ·mm)/(hm2·h)，由西南向东北逐渐增大。总的来说，从春季到秋季，降雨侵蚀力峰值线由东南向西北再到东南移动，而冬季降雨侵蚀力呈不均匀分布。

2.3.2 年际分布特征 由图5可知，白洋淀流域2003—2006年降雨侵蚀力范围为862.91～2 726.92 (MJ·mm)/(hm2·h)，平均2 053.79 (MJ·mm)/(hm2·h)，西南地区定州市、无极县、行唐县和东北地区易县范围内降雨侵蚀力较大，其他地区较小。2007—2010年降雨侵蚀力范围为1 027.75～3 031.42 (MJ·mm)/(hm2·h)，平均2 266.32 (MJ·mm)/(hm2·h)，西南地区曲阳县、行唐县、新乐市、灵寿县、正定区和东北地区高碑店市、涿州市范围内降雨侵蚀力较大，其他地区较小。2011—2014年降雨侵蚀力范围为870.40～3 185.66 (MJ·mm)/(hm2·h)，平均2 584.82 (MJ·mm)/(hm2·h)，东北地区高碑店市、涿州市范围内降雨侵蚀力较大，其他地区较小。2015—2018年降雨侵蚀力范围为1 060.13～3 104.70 (MJ·mm)/(hm2·h)，平均2 475.96 (MJ·mm)/(hm2·h)，东北地区霞云岭、房山区范围内降雨侵蚀力较大，其他地区较小。总的来说，2003—2010年降雨侵蚀力峰值地区偏向西南和东北地区，2011—2018年偏向东北地区。

图4 不同季节降雨侵蚀力分布

图5 不同年份降雨侵蚀力分布

2.4 降雨侵蚀力时间变化特征

2.4.1 年内变化特征 降雨侵蚀力的变化趋势常以线性趋势表示，即倾向率，通过对白洋淀流域2003—2018年各季节降雨侵蚀力倾向率进行插值得到流域内季节侵蚀力的变化分布图。由图6可知，白洋淀流域春季倾向率整体为正，范围为34.49～69.55 (MJ·mm)/(hm2·h)，平均49.81 (MJ·mm)/(hm2·h)，无极县、深泽县、安国市、蠡县和肃宁县范围内倾向率较大，其他地区相对较小。夏季倾向率整体为正，范围为34.60～293.88 (MJ·mm)/(hm2·h)，平均162.84 (MJ·mm)/(hm2·h)，流域北方的灵丘县、涞源县和霞云岭地区倾向率较大，保定市、高阳县和蠡县范围内出现倾向率高于周边地区的峰值结构，其他地区相对较小。秋季倾向率整体为负，范围为-133.49～-60.30 (MJ·mm)/(hm2·h)，平均-89.47 (MJ·mm)/(hm2·h)，阜平县范围内倾向率最大，易县范围内倾向率次之，其他地区相对较小。冬季倾向率整体为负，范围为-0.13～-0.01 (MJ·mm)/(hm2·h)，平均-0.07 (MJ·mm)/(hm2·h)，冬季倾向率变化甚微。整体上看，春夏两季降雨侵蚀力逐年增加，而秋季逐年下降，冬季变化不明显。

图6 降雨侵蚀力倾向率分布

2.4.2 年际变化特征 由图7可知，2003—2018年白洋淀流域降雨侵蚀力倾向率范围为-49.48～129.22 (MJ·mm)/(hm2·h)，平均24.54 MJ·mm/(hm2·h)，说明16年来流域的降雨侵蚀力逐年递增。流域降雨侵蚀力倾向率呈东北高、西南低的分布特征，其中，东北方向为61.22～129.22 (MJ·mm)/(hm2·h)，而西南方向为-49.48～0.56 (MJ·mm)/(hm2·h)，且多分布于石家庄辖区内，流域内其他区域倾向率基本为正值，并在唐县、望都县、定州市和阜平县附近出现倾向率峰值结构。

2.4.3 Mann-Kendall突变检验 对白洋淀流域各年降雨侵蚀力进行分析(图8)，得到其16年的变化趋势，2003—2018年间UF统计量均≥0，说明降雨侵蚀力整体上为上升势态，但UF统计量低于U0.05显著性水平，说明流域的上升势态不显著。UF统计量和UB统计量在2003—2007年、2013—2016年两个时间段内出现交点，交点绝对值均小于1.96，低于U0.05显著性水平，因此，16年间白洋淀流域降雨侵蚀力不存在显著性突变点以及突变区间。

图7 研究区降雨侵蚀力倾向率分布

图8 研究区降雨侵蚀力Mann-Kendall突变检验

3 讨 论

3.1 降雨侵蚀力空间分布差异

白洋淀流域地形复杂多样，气候多变，各区域降雨侵蚀力存在不同程度变化，同时，白洋淀流域作为环京津水土流失重点防治区，深入了解流域降雨侵蚀力特征，可为及时落实水土保持措施提供依据。白洋淀流域地势西北高、东南低，在保定阜平至北京霞云岭一线上形成海拔起伏带，东南暖气流随地势升高而被抬升，伴随气团垂直向上运动的同时温度降低，水汽凝结形成降雨，在风速和坡度增加的共同作用下，降雨效应更加明显[16-17]，西北地区，东南暖气流沿河谷向西北腹地深入，但西北方向海拔升高和太行山区背风坡增加的地形变化[18]抑制了气流的前进，东部地区为广袤平原，东南暖气流受地形因素的抬升作用降低而向西北快速推进，降雨多汇聚于海拔起伏带上，因此，本研究中西南、东北以及东南—西北方向上的中部地区降雨侵蚀力较大，西北山区和东南平原相对较小。降雨侵蚀力影响因素众多，因而不同地区侵蚀力值大小也存在差异，赖成光等[9]对珠江流域展开研究，结果显示，降雨侵蚀力随经度增加、纬度降低而逐渐增大。刘斌涛等[19]对西南地区降雨侵蚀力与海拔关系进行研究，结果显示，降雨侵蚀力变化趋势在高海拔地区明显增加，在低海拔地区明显降低。本研究中，降雨侵蚀力虽与不同输入因子相关关系存在差异，但与纬度、经度、海拔等三因子的相关关系最紧密，且与三者分别呈负、正、负的相关关系，出现该结果的原因在于白洋淀流域轮廓东西、南北距离差异小，其纬度跨度2°3′、经度跨度2°50′，且地形起伏地带位于中部的东北—西南分界线上，致使降雨线呈东南—西北走向，随纬度变大、经度变小、海拔升高，东南季风逐渐降低，侵蚀性降雨锐减。基于上述的线性逐步回归分析，可为白洋淀流域降雨侵蚀力测算提供依据。

3.2 降雨侵蚀力时间分布差异

白洋淀流域气候特点为雨热同期，年内降水分配不均匀，年际降水同样存在不确定性，由于该流域内存在洪灾史，因此，研究降雨侵蚀力的时间分布特征可预防降水带来的灾害。年内分布中，流域降雨侵蚀力春季东南地区大、增长趋势高，夏季西北地区大、增长趋势高，秋季东南—西北方向大、降低趋势低，冬季东北地区大、降低趋势低，原因在于春季为雨季前季，东南暖气流率先到达流域东南区域，雨季的到来使西北半干旱区相比其他地区雨量骤增程度大，而秋季受西伯利亚寒流影响降雨线向东南移动，冬季主要以降雪形成产生。在霞云岭—阜平沿线上出现的秋季降雨侵蚀力以及倾向率变化可能是由局部地区地形地貌引起的水汽差异造成的，钟科元等[20]对松花江流域降雨侵蚀力的研究结果显示，山地区域的水汽差异可引起局部异于整体的现象。研究结果显示，冬季降雨侵蚀力以及春季、冬季降雨侵蚀力倾向率均呈不均匀分布，原因在于流域内降雨侵蚀力及其倾向率的局部差异在大范围内极易引起地域上明显的变化，这种明显变化造成分布上的不均匀。年际变化中，西南地区降雨侵蚀力整体上呈降低趋势，东北地区整体上呈增加趋势，造成西南地区降雨侵蚀力倾向率为负，东北地区倾向率为正，原因可能是东北方向受季风气候、海平面高压、太平洋气流以及温室气体的影响大于西南地区[11]，同时，东北地区一侧为低山丘陵半环形环绕，一侧为山前平原，山间河谷呈喇叭口向东南开放，利于季风气流长驱直入，地形因素也是导致东北地区多年降雨侵蚀力增大的重要原因。整体上流域的降雨侵蚀力呈逐年增加趋势，但Mann-Kendall突变检验结果显示该趋势不显著，说明流域降雨带来的侵蚀变化平稳、波动小，干旱、洪涝等灾害性年份少见。

4 结 论

2003—2018年白洋淀流域年降雨侵蚀力范围为909.62～3 358.38 (MJ·mm)/(hm2·h)，东南—西北方向上变化趋势为先增后降，而西南和东北地区降雨侵蚀力变化剧烈程度高于其他地区，通过分析站点降雨侵蚀力在空间上的分布规律发现，其值大小与纬度、经度、海拔的相关关系分别为负、正、负。降雨侵蚀力峰值线季节变化基本呈东南—西北方向移动，年际变化基本呈西南—东北方向移动；倾向率峰值线春季和夏季在东南—西北方向上移动，秋季、冬季和年际变化在西南—东北方向上变化。Mann-Kendall检验结果显示，16年间的降雨侵蚀力变化不明显。根据本研究结果，流域内海拔起伏带的降雨侵蚀力以及西北山区夏季倾向率的变化均较大，该地区应提高水土流失防范意识，汛期积极落实相关措施，减少降雨侵蚀。