湟水流域水沙情势演变及驱动力分析

余 欣，李其江，刘希胜，闵 敏，杨 阳，时 璐，马 蕊

(1.黄河水利委员会 黄河水利科学研究院，河南 郑州 450003;2.青海省水利厅，青海 西宁 810001;3.青海省水文水资源测报中心，青海 西宁 810001）

湟水流域位于我国西北地区的生态过渡带,干流西宁站以下的黄土丘陵区水土流失严重[1］。近年来,受人类活动和气候变化的影响,湟水流域水沙关系发生了很大变化。

本文以湟水流域民和站以上区域为研究对象,基于水沙长时间序列资料,利用趋势分析、突变检验等方法,从侵蚀性降水量(日降水量大于12 mm)、NDVI指数和水土保持措施等三个方面,分析了水沙情势演变特征及其驱动因素,以期为水资源开发利用、水旱灾害防治、水土流失治理等提供参考。

1 数据与方法

1.1 数据来源

收集了湟水干流湟源、西宁和乐都3个青海省属水文站实测径流、泥沙数据以及黄河水利委员会水文局所属湟水干流民和水文站实测径流、泥沙数据。数据序列为1956—2020年,资料质量较好、序列完整。

1.2 研究方法

1.2.1 趋势检验

趋势检验采用Sen斜率估计法和Mann-Kendall[1］检验法(M-K法)。Sen斜率估计法是一种稳健的非参数统计趋势计算方法,该方法计算效率高,对测量误差和离群数据不敏感,常被用于长时间序列数据的趋势分析中。Sen斜率β公式为

式中:median为中值函数；Xi、Xj为序列中数据值。

β值为正表示上升趋势,β值为负表示下降趋势。

M-K[3-5］法是一种非参数统计检验方法,变量可以不具有正态分布特征,适用于水文变量的趋势检验。假设时间序列(X1,X2,…,Xn),统计量S为

方差Var(S):

检验指标Z:

在双边趋势检验中,对于给定的置信水平α,若则在置信水平α上时间序列数据存在明显的上升或下降趋势。Z为正值表示上升趋势,Z为负值表示下降趋势。表示其通过了置信度为95%的显著性检验。

以上两种方法结合可以增强研究方法的抗噪性,并在一定程度上提高检验成果的准确性。

1.2.2 突变检验

突变检验采用Pettitt检验法[2,6］。Pettitt检验法是非参数检验法,该方法在检验连续序列时是高效的,能够较好地识别水文时间序列的突变点,在突变点检验中应用较多,且物理意义清晰。该检验基于Mann-Whitney的统计函数认为两个序列(X1,X2,…,Xt)和(Xt＋1,Xt＋2,…,XT)均来自于同一序列,对于连续的序列,Ut,T和Vt,T由下式计算:

式中:Ut,T、Vt,T为不同时段检验值统计量。

当｜Ut,T｜最大时,对应的Xt是可能的突变点。当突变点Ut,T＞0时,该序列具有向下突变趋势,反之则具有向上突变的趋势。可能突变点的显著性水平由下式计算:

当POA(t)≤0.5时,该点为有效突变点。

1.2.3 周期检验

Morlet小波[7-9］作为一种分解工具和检测长时间序列中多尺度突变的方法,具有时域和频域的良好局部化特征,而且随着信号不同频率成分在时域取样的疏密,调节其放大倍数,可以清楚地分析信号在各个层次上的变化,进而对源信号在不同时间尺度上的变化趋势进行判断和甄别。

2 结果与分析

2.1 水沙变化分析

2.1.1 趋势性

采用湟水干流湟源、西宁、乐都、民和4个水文站的径流、泥沙数据,分析1956—2020年湟水流域从上游到中上游、中游、中下游主要水文站的天然年水沙时序,见图1、图2。

图1 湟水干流各水文站天然年径流量过程线

图2 湟水干流各水文站天然年输沙量过程线

通过Sen斜率估计法和M-K趋势检验法,分析各站径流量、输沙量变化趋势,结果见表1。1956—2020年各水文站径流量呈增加趋势,增加速率在0.2亿～0.5亿m3/10 a之间,M-K检验指标Z在2.04～3.40之间,呈显著上升趋势。同期各站输沙量呈减小趋势,减少速率在3.3万～386.2万t/10 a之间,M-K检验指标Z在-6.8～-2.2之间,呈显著下降趋势。

表1 湟水干流水沙趋势分析结果

通过Pettitt突变检验可得,各站径流量突变发生在1982年和2002年左右,见表2。1982—2001年与1956—1981年相比多年平均径流量增加3.9%～11.6%,2002—2020年与1982—2001年相比多年平均径流量增加7.5%～19.7%。各站输沙量在2000年左右发生突变,2001—2020年与1956—2000年相比多年平均输沙量减少45.6%～79.3%。

表2 湟水干流水沙突变分析结果

2.1.2 周期性

用小波分析法识别湟水流域1956—2020年径流量、输沙量序列周期成分。湟水径流量上下游呈现相对一致的周期变化,在1956—2020年主要存在中心时间尺度为9、22、34、47 a的周期分布,在34 a的周期尺度上对应小波方差最大,说明周期震荡最强,为主要周期,见图3。

图3 湟水干流各水文站径流量小波周期分析

从湟水干流各站输沙量小波周期分析可以看出,上游湟源站与下游3站周期变化不同(见图4)。上游的湟源站以30 a尺度为主周期；中下游西宁、乐都、民和站的小波变换周期类似,存在11～25 a尺度的震荡,尤其在11、14、21 a尺度的小波震荡能量较强。

图4 湟水干流各水文站输沙量小波周期分析

2.2 水沙关系分析

2.2.1 水沙异源

按照水文站区间控制单元,将流域划分为湟源以上、湟源—西宁、西宁—乐都、乐都—民和4个区间。经过分析,流域产水、产沙存在区间差异性,水沙存在明显异源现象。

流域水量主要来源于湟源至西宁区间,该区间贡献水量和沙量分别占民和站实测的41.4%和19.1%；而沙量主要来源于乐都至民和区间,该区间贡献水量和沙量分别占民和站实测的16.1%和42.7%,两个区间水量和沙量贡献率大小相反,见表3。

表3 湟水流域各区间水沙量对比

考虑各水文站控制面积不同,分别计算各区间的产水模数与产沙模数,分析可知,4个区间的产水模数接近,在9.48万～11.31万m3/km2之间。因为产水模数差异不大,所以区间面积的大小是产水量多少的主导因素,由于湟源—西宁区间面积最大,因此其成为主要产水区间。

从上游至下游,4个区间的产沙模数成倍增加,由湟源以上的134.16 t/km2增加到乐都—民和区间的2289.86 t/km2,在产水模数接近的前提下,4个区间的水土流失程度呈现自上游至下游成倍递增趋势。

2.2.2 上下游水沙关联性

湟水干流上下游相邻两站年径流量相关性较好,由上游至下游两两相关关系增强,说明年径流自上游至下游具有良好的协同关系,见图5。

图5 湟水干流各站实测年径流量上下游相邻两站相关关系(单位:亿m3)

在相邻两站实测年输沙量相关关系中,除湟源与西宁站的相关关系较差外,中下游各站相关性较好,说明湟水流域中下游主要产沙地区各站的协同关系好于上游,见图6。

图6 湟水干流各站实测年输沙量上下游相邻两站相关关系(单位:万t)

2.2.3 水沙时序变化协同作用

湟水干流各站径流量与输沙量的关系散乱,相关性较差,呈现不同变化趋势。时序上,各站年径流量整体呈现上升态势,输沙量则呈现下降态势。

综合径流输沙突变时间、历年过程以及下垫面变化等情况,将长程时序分为1956—1980年、1981—2002年以及2003—2020年3个时段进行纵深分析。各站实测年径流量总体呈现上升态势,在2002年以后,各站均值、中位数、四分位值(上、下四分位值)、最大最小值均有所上升。3个阶段,各站实测年输沙量呈现下降态势,尤其是西宁、乐都、民和站输沙量持续下降,其均值、中位数、四分位值(上、下四分位值)、最大最小值均有所下降,见图7。

图7 湟水干流各站实测年径流量、输沙量不同时段对比

2.3 驱动力分析

2.3.1 降水指标

降水是湟水流域水资源的主要补给来源,但侵蚀性降水同时会导致土壤水蚀。本文选用日降水量大于12 mm的年侵蚀性降水量(PRCPTOT12)、侵蚀性降水天数(R12)、侵蚀性降水强度(SDII12)三个降水指标,分析湟水1981—2020年PRCPTOT12、R12、SDII12的长期变化趋势(见图8)。湟水流域PRCPTOT12、R12趋势表现基本一致,年代际呈现“先增后减”态势。2001—2010年年平均PRCPTOT12为研究时段最大,达到211.0 mm,PRCPTOT12以9.00 mm/10 a速率呈不显著上升趋势,最高值为2009年的264.2 mm,最低值为1982年的109.8 mm。湟水流域研究时段内多年平均R12以0.43 d/10 a速率呈不显著上升趋势,总体波动较大,变化范围在6.0～14.2 d之间,1991—2000年时段均值上升,至2000年后转为偏多趋势,在2001—2010年为研究时段最大,R12达到11.5 d,2001—2010年、2011—2020年两时期分别较均值偏多0.8、0.3 d。湟水流域研究时段内多年平均SDII12为18.3 mm/d,变化范围在16.2～20.1 mm/d之间,总体波动不大,年代际呈现“增—减—增”特征。黄土高原地区降水量显著减少,减少幅度自东南向西北逐渐变小,侵蚀性降水量与年降水量变化相似,但湟水流域降水量和侵蚀性降水量呈增加趋势,不同于黄土高原变化趋势[10］,这与本文研究结果类似。

图8 侵蚀性降水指标变化趋势

采用相关分析法进一步分析湟水干流民和站径流量、输沙量与降水指标的响应关系(见表4)。从相关性分析结果看,湟水流域PRCPTOT12、R12与民和站径流量呈极显著正相关关系,SDII12与输沙量呈显著正相关关系,其余相关关系不显著。这表明PRCPTOT12、R12的增加是湟水流域径流量增加的主要因素,SDII12是引起输沙量变化的重要因素。这也意味着尽管增加趋势不显著,但随着极端天气事件的增多,由水蚀导致的土壤侵蚀总体上还将继续增多,这必然对湟水流域高质量发展产生深远影响。

表4 民和站径流量、输沙量与降水指标的相关系数

2.3.2NDVI指数

植被变化受人类活动、气候条件等因素共同作用,不同下垫面条件导致水热条件空间分布不均,进而造成NDVI指数空间异质性[11］。植被覆盖也是衡量水沙变化的关键因素之一[12］。湟水流域NDVI年均值与5—9月均值空间分布见图9,NDVI年均值介于0.13～0.68之间,区域性差异明显,表现为由西北向东南方向递减,从上游向中下游递减,表明西北部植被覆盖率高于东南部,这与王兮之等[13-14］研究结论相似。流域西北部多为海拔较高、气候寒冷区域,植被以林地和草地为主,区域内气温低,热量少,蒸发量小,城乡居民建设用地较少,受人类活动干扰小,植被覆盖度相对较高,降水的侵蚀产沙能力较弱。被浅山区包围的中下游河谷地带,地势相对较低,水热条件较好,人类活动频繁,沿河两岸遍布城乡居民建设用地、工矿企业和耕地,NDVI较上游偏低,下垫面截留固沙能力较弱。

图9 湟水流域NDVI指数空间分布

湟水流域降水主要集中在5—9月,在此期间NDVI均值较年均值有所提高,尤其中上游区域变化最为明显。夏季降水较多,上游海拔较高,气温相对下游偏低,降水蒸发较慢,雨水经下渗进入土壤促进植物生长,因此在一定程度上可以减少雨水侵蚀、增强地表渗透性与抗蚀性[15］,但局部强降水易冲刷植被覆盖度较低的湟水谷地坡面和沟道,造成输沙量增大,进一步印证中下游是输沙的主要来源。

2.3.3 水土保持措施

湟水流域人口占青海省总人口的60%左右,人类活动频繁,下垫面条件复杂。自2000年以来,青海省水利部门通过新建基本农田、营造水保林和经济林、封育治理等措施有效改善湟水流域水土流失状况,重点实施了小流域综合治理、淤地坝坝系建设、坡耕地水土流失综合防治、巩固退耕还林区基本口粮田工程等一系列生态修复和预防保护工程,对于减少水土流失、增加土壤入渗、拦蓄坡沟泥沙和地表径流在一定程度上发挥了重要作用。根据全国第一次水利普查数据,湟水水利工程数量增长不明显,但随着黑泉水库建成累计总库容陡增,这也是西宁—民和区间输沙量减少的原因之一。累计建设淤地坝661座,小型蓄水工程62864个。湟水中下游黄土丘陵区水土流失得到抑制,2000年之后输沙量显著下降,治理效果、生态效益显著提高。

3 结论

(1)湟水干流4站年径流量和输沙量分别呈现不同变化趋势。各站径流量均呈增加趋势,速率为0.2亿～0.5亿m3/10 a,突变发生在1982年和2002年左右。各站输沙量则呈下降趋势,下游减小幅度较大,在2000年左右发生突变。

(2)湟水流域径流量上下游呈现相对一致的周期变化,主要存在中心时间尺度为9、22、34、47 a的周期分布；输沙量湟源站存在30 a尺度的主周期,而下游西宁、乐都、民和站较为相似,在11、14、21 a尺度的小波方差较为突出。

(3)湟水流域径流纵向空间具有连贯协同性,中下游泥沙协同关系好于上游。水沙趋势各异,产水、产沙存在明显异源现象,湟源—西宁区间面积最大致使其成为主要来水区间,乐都—民和区间为主要产沙区间。

(4)PRCPTOT12、R12对径流量变化具有促进作用,SDII12变化引起输沙量的变化。得益于湟水流域生态修复和预防保护工程的实施,流域内植被覆盖度提高,水土流失程度降低、土壤入渗增加、拦蓄坡沟泥沙和地表径流成效显著。