基于小波变换的近70年燕山水库径流演变特征分析

2023-05-24王雪彦范思聪王富强石家豪

华北水利水电大学学报(自然科学版) 2023年2期

王雪彦, 范思聪, 王富强,3, 石家豪,3

(1.华北水利水电大学,河南郑州 450046; 2.河南省燕山水库管理局,河南平顶山 467224; 3.河南省黄河流域水资源节约集约利用重点实验室,河南郑州 450046)

受气候、地貌、植被等自然条件以及人类活动耦合作用的影响,流域水文循环有所加速。径流是水文循环不可或缺的一部分,是地表水资源的主要形式,对人类生存和发展起至关重要的作用[1-2]。变化环境下水库径流及地表水资源量的变化将对区域水资源管理带来严峻挑战。因此,开展水库径流演变特征研究对区域水资源合理分配及管控具有重要意义。

目前,关于径流演变特征的研究主要集中于径流的年代际变化趋势、突变特征和周期性分析等方面,相应的研究方法颇多,其中:趋势性主要利用基尼系数[3]、重标极差(R/S)[4]、M-K趋势分析[5]和Spearman秩次相关检验[6]等方法进行分析;突变特征主要通过滑动T检验法和累计距平分析方法[7]等进行分析;周期性主要利用最大熵谱分析法[8]、人工神经网络[9]和奇异谱分析[10]等进行分析。但上述分析方法自身存在一定缺陷:①只在频域上具备局部化性质;②在诊断突变点过程中缺乏数学上的严谨性[11-12]。20世纪80年代初,以Fourier分析方法为基础发展起来的小波分析技术,选用一种窗口大小可改变、位置可动的变窗进行谱分析,从而同时满足信号时域和频域局部化的要求,很好地克服了上述传统谱分析方法的不足,故此研究方法在水文时间序列周期性分析中被广泛应用[13-14]。受天气系统、流域下垫面和其他因素的影响,径流、暴雨、洪水等水文要素时间序列几乎都是非平稳的,内部变化极为复杂,表现出多时间尺度的变化特征,利用小波变换分析技术可揭示水文要素时间序列在时域和频域上的局域变化特征,为径流分析和径流预报提供理论依据。

国内学者对沙颍河流域径流演变特征进行了大量的研究,如:甘容等[15]以沙颍河流域为研究对象,选取改进的SWAT模型模拟了径流过程,评价了模拟效果,分析了沙颍河降雨及径流变化过程;戴韵秋等[16]运用线性趋势法等分析了沙颍河上游降雨及径流序列的演变特征,量化了气候、人类活动这两个因素对沙颍河流域径流变化的影响。以往的研究多集中于沙颍河流域径流量演变特征方面,而有关沙颍河主要支流上的燕山水库径流演变过程的研究较为单一。从目前已有研究成果来看,有关燕山水库径流的研究主要是分析降水、人类活动等影响因子对径流的影响以及其多年趋势性变化特征[17],关于径流的季节性变化和周期性研究相对薄弱。鉴于此,本文结合基尼系数及Mann-Kendall检验法分析了燕山水库入库径流量的年内分配、年际和各季节变化趋势,并基于Morlet小波分析法分析了多时间尺度下的变化特征和演变规律,以期为水库水资源合理开发与利用提供依据。

1 数据来源与研究方法

1.1 数据来源

燕山水库位于沙颍河主要支流澧河上游的干江河上,坝址位于官寨水文站和官寨(二)水文站之间,控制流域面积约为1 170 km2,官寨水文站控制流域面积约为1 130 km2,官寨(二)水文站控制流域面积约为1 200 km2。本文基于官寨水文站1955—1993年共计39年的实测径流资料、官寨(二)水文站1994—2008年共计14年的实测径流资料、燕山水库水文站2009—2018年共计10年的实测径流资料进行分析。为使资料一致,对燕山水库实测径流数据进行还原,采用还原后的1955—2018年的逐月径流数据,分析水库的径流变化特征。

1.2 研究方法

1.2.1 基尼系数

1912年意大利统计与社会学者Corrado Gini提出了一种用来量化收入分配平等程度的指标——基尼系数,其被广泛应用于经济领域分布不均匀度的量化评价中。基尼系数GI值的取值范围为[0,1],当GI<0.2时认为收入分配程度绝对均匀,当GI>0.5时认为收入分配程度不均匀,其值越大,表示不均匀程度越高[18-19]。大量学者将其引入到水文序列的年内分配均匀度变异分析与径流时空分布均匀度的研究中。本文利用基尼系数GI值判断径流年内分配程度。

1.2.2 Mann-Kendall非参数检验法

Mann-Kendall检验法(简称M-K法)是世界气象组织推荐并广泛用于研究水文和气候长时间序列的一种非参数检验法,检测范围大,定量化程度高,受人为干扰影响较小,易于反映数据的突变性,是水文、气象要素时间序列趋势检验中应用较多且具有理论意义的方法[20-21]。本文选用M-K法诊断径流年际变化趋势的显著性和检验径流年际变化的突变点。

1)趋势检验法。假设径流时间序列x1、x2、…、xn服从独立同分布,构造统计量Z:

(1)

其中:

(2)

(3)

sgn为符号函数,若统计量Z>0认为时间序列有增加趋势;若统计量Z<0认为时间序列有减少趋势。给定显著性水平α,当|Z|≥Zα/2时,表示水文时间序列有显著的上升或下降趋势。

2)突变检验法。设径流时间序列为x1、x2、…、xn,构造一秩序列Sk:

(4)

其中:

(5)

式中:1≤i≤n;1≤j≤i,j=1、2、…、n。

在时间序列服从独立同分布的假设下,Sk的均值、方差分别为:

(6)

将Sk标准化为:

(7)

式中:UF1=0;给定显著性水平α,若|UF|>Uα,表明序列存在明显的变化趋势。反向序列的顺序变化为xn、xn-1、…、x1,反向序列标准化为UBk,这里设定显著性水平值α=0.05,临界值U0.05=±1.96。UFk和UBk组成2条曲线,UFk和UBk的交点称为该径流序列的突变点。

1.2.3 小波分析

小波函数是指具有震荡特性、能够快速衰减到零的一类函数。于20世纪80年代发展起来的小波分析技术同时在时域和频域上具有良好的局部化功能特性,可满足对水文时间序列的局部化分析要求,解析其内部精细结构,在分析序列的多时间尺度特性过程中,可有效识别主要频率成分和读取局部化信息,有利于径流序列的多种频率成分分析,可将其广泛应用于分析水文时间序列周期变化特性和水文预测等方面[22-23]。设时间序列f(t)∈L2(R)(L2(R)表示能量有限)的连续小波变换为:

(8)

(9)

小波方差反映的是径流时间序列的能量波动随时间尺度a的变化情况,以此来精准判断信号中不同时间尺度的主周期[24]。小波方差计算公式为:

(10)

因水文水资源系统具有时变特性,所以水文时间序列常表现出非平稳特性。为详细刻画水文时间序列的非平稳性,需要一种可同时在时域和频域上进行分析的函数。Morlet提出的小波分析满足这一要求,故本研究选用Morlet小波函数对径流周期性特征进行分析。

2 结果分析与讨论

2.1 径流年内变化特征

基于燕山水库1955—2018年入库径流序列资料,首先运用四元一次方程拟合出洛伦兹曲线(图1),得出曲线的函数表达式,用定积分求出洛伦兹曲线与横坐标轴围成的面积之后,计算得出基尼系数为0.62,表明燕山水库的年内径流量分配极度不均匀。

图1 径流量年内分配洛伦茨曲线

不同时段径流量的年内变化如图2所示。由图2可知:①1—6月和10—12月的径流量较小,变化缓慢,汛期(7—9月)的径流量较大且变化幅度明显;燕山水库径流量存在季节性分配不均的特点,其中春、夏、秋、冬季径流量占全年径流量的比例分别为11%、63%、21%、5%,夏季径流量占比最大,冬季径流量占比最小。②不同时段的月径流量过程线均呈单峰型,20世纪80年代后月径流量峰值的发生时间出现前移。1955—1979年的径流量峰值发生在8月,而1980—2000年和2001—2018年的径流量峰值均出现在7月。③不同时段的年径流量总体呈减少趋势,1980—2000年和2001—2018年的径流量较1955—1979年的分别减少了9%、27%,且汛期各月径流量的变化最为明显。

图2 不同时段径流量的年内变化

2.2 径流年际变化特征

2.2.1 径流变化趋势分析

图3为不同时间尺度下的燕山水库径流量变化过程。从图3(a)可以看出,燕山水库径流量以0.026亿m3/年的线性倾向率下降。从图3(b)—3(e)可以看出,春、夏、秋、冬季径流量呈现不同程度的下降趋势,其倾向率分别为0.004亿、0.017亿、0.005亿、0.000 4亿m3/年。1975年8月淮河上游的洪汝河和沙颍河流域发生了历史少见的特大暴雨,使燕山水库入库径流量达到峰值。因此,夏季径流量的峰值和年径流量的峰值出现时间点吻合。从5年滑动平均曲线可以看出:研究区全年和季节径流量均在趋势线附近上下波动,呈现“增—减”的波动变化;5年滑动平均曲线均有1个明显的峰值,20世纪50年代至20世纪末期的波动较大,径流量变化剧烈;2000年后的波动较小,径流量变化缓慢。

图3 1955—2018年燕山水库径流量变化过程

采用M-K检验法进一步诊断燕山水库全年和各季节径流量变化趋势的显著性,得到Z值依次为-1.84、-1.67、-1.25、-0.89和-0.64。由此可以看出,Z值均小于0,且均未超过0.05置信水平阈值(±1.96)。因此,全年和各季节径流量均呈非显著性下降趋势。

2.2.2 径流突变分析

采用M-K突变检验法对燕山水库径流序列进行突变性分析,其结果如图4所示。由图4可得,UF和UB大多未超出显著性水平α=0.05临界值。

图4 1955—2018年燕山水库年径流突变点检验结果

由图4(a)可知:UF值在1955—2018年基本小于零,只有在1984年大于零;UF和UB交点较少,明显的突变点只出现在1963年和1993年,且1993年后波动变化并超出临界线,说明燕山水库径流在1963年和1993年发生突变。由图4(b)可知:1955—1974年和1981—2018年的径流量呈下降趋势,1975—1980年的径流量呈上升趋势,但上升趋势不明显;1957—1963年超出临界值,UF和UB的交点出现在1963年和1981年,说明春季径流量在1963年和1981年出现突变点。由图4(c)可知:夏季径流量在2003—2007年呈上升趋势,其余时间段均呈下降趋势;UF和UB的交点出现在1964年和1986年,说明夏季径流量在1964年和1986发生突变。由图4(d)可知:秋季径流量在1963—1977年和1981—1991年呈上升趋势,在1955—1962和1982—2018年呈下降趋势;UF和UB的交点出现在1963年和1986年,说明秋季径流量在1963年和1986年发生突变。由图4(e)可知:冬季径流量在分界点上下波动变化,UF和UB的交点出现在1963年,说明冬季径流量在1963年发生突变。

从以上分析可知:燕山水库1955—2018年的年径流量和季节径流量均存在多个突变点,年径流量突变点大致发生在1963年和1993年,各季径流量突变点发生时间不同,大致发生在1963年和1986年。

2.3 径流周期性特征

本文应用小波分析法计算得到燕山水库不同时间尺度下的径流序列的小波变换等值线图和小波方差图,分别如图5—9所示。小波变换实部分布图(图5—9的(a)图)的下半部等值线相对密集,为高频曲线,其对应小尺度周期的震荡;而上半部等值线相对稀疏,为低频曲线,其对应大尺度的周期震荡。

图5 1955—2018年全年径流量小波分析结果

由图5(a)可知,年径流量存在明显的周期变化,在18～24年、12～15年和3～8年时间尺度上分别有5次、6次和12次的“丰—枯”交替变化。由图6(a)可知,春季径流量在26～28年尺度上经历4次“枯—丰”交替变化,在12～18年以及3～11年尺度上分别有6次和10次“丰—枯”交替变化。由图7(a)可知,夏季径流量在26～30年尺度上经历4次“枯—丰”变化,在18～24年尺度上经历6次“丰—枯”交替变化,对于更小的时间尺度其“丰—枯”交替变化更频繁。由图8(a)可知,秋季径流量在24～30年尺度上有4次“枯—丰”变化,8～16年尺度上有6次“丰—枯”交替变化,对于较小的尺度交替变化更频繁。由图9(a)可知,冬季径流量在2～8年时间尺度上经历7次“丰—枯”交替变化。

图6 1955—2018年春季径流量小波分析结果

图7 1955—2018年夏季径流量小波分析结果

图8 1955—2018年秋季径流量小波分析结果

图9 1955—2018年冬季径流量小波分析结果

由图5—9的(b)图可以看出,燕山水库径流量存在显著的多时间尺度变化特点,年径流量存在15年的主周期,春、夏、秋、冬季径流量分别存在14年、12年、11年和3年的主周期。

综上可知,燕山水库四季的径流量变化并非以一种固定周期的形式发生,而是以各种时间尺度相互嵌套的结构出现,1955—2018年燕山水库径流量变化周期总体上呈现出18～24年的大尺度变化特点和3～6年、8～16年的小尺度变化特点;年径流量存在15年的主周期,春、夏、秋、冬季径流量分别存在14年、12年、11年和3年的主周期。