基于小波變換的近70年燕山水庫徑流演變特征分析

王雪彥, 范思聰, 王富強(qiáng),3, 石家豪,3

(1.華北水利水電大學(xué),河南 鄭州 450046; 2.河南省燕山水庫管理局,河南 平頂山 467224; 3.河南省黃河流域水資源節(jié)約集約利用重點(diǎn)實(shí)驗(yàn)室,河南 鄭州 450046)

受氣候、地貌、植被等自然條件以及人類活動耦合作用的影響,流域水文循環(huán)有所加速。徑流是水文循環(huán)不可或缺的一部分,是地表水資源的主要形式,對人類生存和發(fā)展起至關(guān)重要的作用[1-2]。變化環(huán)境下水庫徑流及地表水資源量的變化將對區(qū)域水資源管理帶來嚴(yán)峻挑戰(zhàn)。因此,開展水庫徑流演變特征研究對區(qū)域水資源合理分配及管控具有重要意義。

目前,關(guān)于徑流演變特征的研究主要集中于徑流的年代際變化趨勢、突變特征和周期性分析等方面,相應(yīng)的研究方法頗多,其中:趨勢性主要利用基尼系數(shù)[3]、重標(biāo)極差(R/S)[4]、M-K趨勢分析[5]和Spearman秩次相關(guān)檢驗(yàn)[6]等方法進(jìn)行分析;突變特征主要通過滑動T檢驗(yàn)法和累計(jì)距平分析方法[7]等進(jìn)行分析;周期性主要利用最大熵譜分析法[8]、人工神經(jīng)網(wǎng)絡(luò)[9]和奇異譜分析[10]等進(jìn)行分析。但上述分析方法自身存在一定缺陷:①只在頻域上具備局部化性質(zhì);②在診斷突變點(diǎn)過程中缺乏數(shù)學(xué)上的嚴(yán)謹(jǐn)性[11-12]。20世紀(jì)80年代初,以Fourier分析方法為基礎(chǔ)發(fā)展起來的小波分析技術(shù),選用一種窗口大小可改變、位置可動的變窗進(jìn)行譜分析,從而同時(shí)滿足信號時(shí)域和頻域局部化的要求,很好地克服了上述傳統(tǒng)譜分析方法的不足,故此研究方法在水文時(shí)間序列周期性分析中被廣泛應(yīng)用[13-14]。受天氣系統(tǒng)、流域下墊面和其他因素的影響,徑流、暴雨、洪水等水文要素時(shí)間序列幾乎都是非平穩(wěn)的,內(nèi)部變化極為復(fù)雜,表現(xiàn)出多時(shí)間尺度的變化特征,利用小波變換分析技術(shù)可揭示水文要素時(shí)間序列在時(shí)域和頻域上的局域變化特征,為徑流分析和徑流預(yù)報(bào)提供理論依據(jù)。

國內(nèi)學(xué)者對沙潁河流域徑流演變特征進(jìn)行了大量的研究,如:甘容等[15]以沙潁河流域?yàn)檠芯繉ο?選取改進(jìn)的SWAT模型模擬了徑流過程,評價(jià)了模擬效果,分析了沙潁河降雨及徑流變化過程;戴韻秋等[16]運(yùn)用線性趨勢法等分析了沙潁河上游降雨及徑流序列的演變特征,量化了氣候、人類活動這兩個(gè)因素對沙潁河流域徑流變化的影響。以往的研究多集中于沙潁河流域徑流量演變特征方面,而有關(guān)沙潁河主要支流上的燕山水庫徑流演變過程的研究較為單一。從目前已有研究成果來看,有關(guān)燕山水庫徑流的研究主要是分析降水、人類活動等影響因子對徑流的影響以及其多年趨勢性變化特征[17],關(guān)于徑流的季節(jié)性變化和周期性研究相對薄弱。鑒于此,本文結(jié)合基尼系數(shù)及Mann-Kendall檢驗(yàn)法分析了燕山水庫入庫徑流量的年內(nèi)分配、年際和各季節(jié)變化趨勢,并基于Morlet小波分析法分析了多時(shí)間尺度下的變化特征和演變規(guī)律,以期為水庫水資源合理開發(fā)與利用提供依據(jù)。

1 數(shù)據(jù)來源與研究方法

1.1 數(shù)據(jù)來源

燕山水庫位于沙潁河主要支流澧河上游的干江河上,壩址位于官寨水文站和官寨(二)水文站之間,控制流域面積約為1 170 km2,官寨水文站控制流域面積約為1 130 km2,官寨(二)水文站控制流域面積約為1 200 km2。本文基于官寨水文站1955—1993年共計(jì)39年的實(shí)測徑流資料、官寨(二)水文站1994—2008年共計(jì)14年的實(shí)測徑流資料、燕山水庫水文站2009—2018年共計(jì)10年的實(shí)測徑流資料進(jìn)行分析。為使資料一致,對燕山水庫實(shí)測徑流數(shù)據(jù)進(jìn)行還原,采用還原后的1955—2018年的逐月徑流數(shù)據(jù),分析水庫的徑流變化特征。

1.2 研究方法

1.2.1 基尼系數(shù)

1912年意大利統(tǒng)計(jì)與社會學(xué)者Corrado Gini提出了一種用來量化收入分配平等程度的指標(biāo)——基尼系數(shù),其被廣泛應(yīng)用于經(jīng)濟(jì)領(lǐng)域分布不均勻度的量化評價(jià)中。基尼系數(shù)GI值的取值范圍為[0,1],當(dāng)GI<0.2時(shí)認(rèn)為收入分配程度絕對均勻,當(dāng)GI>0.5時(shí)認(rèn)為收入分配程度不均勻,其值越大,表示不均勻程度越高[18-19]。大量學(xué)者將其引入到水文序列的年內(nèi)分配均勻度變異分析與徑流時(shí)空分布均勻度的研究中。本文利用基尼系數(shù)GI值判斷徑流年內(nèi)分配程度。

1.2.2 Mann-Kendall非參數(shù)檢驗(yàn)法

Mann-Kendall檢驗(yàn)法(簡稱M-K法)是世界氣象組織推薦并廣泛用于研究水文和氣候長時(shí)間序列的一種非參數(shù)檢驗(yàn)法,檢測范圍大,定量化程度高,受人為干擾影響較小,易于反映數(shù)據(jù)的突變性,是水文、氣象要素時(shí)間序列趨勢檢驗(yàn)中應(yīng)用較多且具有理論意義的方法[20-21]。本文選用M-K法診斷徑流年際變化趨勢的顯著性和檢驗(yàn)徑流年際變化的突變點(diǎn)。

1)趨勢檢驗(yàn)法。假設(shè)徑流時(shí)間序列x1、x2、…、xn服從獨(dú)立同分布,構(gòu)造統(tǒng)計(jì)量Z:

(1)

其中:

(2)

(3)

sgn為符號函數(shù),若統(tǒng)計(jì)量Z>0認(rèn)為時(shí)間序列有增加趨勢;若統(tǒng)計(jì)量Z<0認(rèn)為時(shí)間序列有減少趨勢。給定顯著性水平α,當(dāng)|Z|≥Zα/2時(shí),表示水文時(shí)間序列有顯著的上升或下降趨勢。

2)突變檢驗(yàn)法。設(shè)徑流時(shí)間序列為x1、x2、…、xn,構(gòu)造一秩序列Sk:

(4)

其中:

(5)

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

在時(shí)間序列服從獨(dú)立同分布的假設(shè)下,Sk的均值、方差分別為:

(6)

將Sk標(biāo)準(zhǔn)化為:

(7)

式中:UF1=0;給定顯著性水平α,若|UF|>Uα,表明序列存在明顯的變化趨勢。反向序列的順序變化為xn、xn-1、…、x1,反向序列標(biāo)準(zhǔn)化為UBk,這里設(shè)定顯著性水平值α=0.05,臨界值U0.05=±1.96。UFk和UBk組成2條曲線,UFk和UBk的交點(diǎn)稱為該徑流序列的突變點(diǎn)。

1.2.3 小波分析

小波函數(shù)是指具有震蕩特性、能夠快速衰減到零的一類函數(shù)。于20世紀(jì)80年代發(fā)展起來的小波分析技術(shù)同時(shí)在時(shí)域和頻域上具有良好的局部化功能特性,可滿足對水文時(shí)間序列的局部化分析要求,解析其內(nèi)部精細(xì)結(jié)構(gòu),在分析序列的多時(shí)間尺度特性過程中,可有效識別主要頻率成分和讀取局部化信息,有利于徑流序列的多種頻率成分分析,可將其廣泛應(yīng)用于分析水文時(shí)間序列周期變化特性和水文預(yù)測等方面[22-23]。設(shè)時(shí)間序列f(t)∈L2(R)(L2(R)表示能量有限)的連續(xù)小波變換為:

(8)

(9)

小波方差反映的是徑流時(shí)間序列的能量波動隨時(shí)間尺度a的變化情況,以此來精準(zhǔn)判斷信號中不同時(shí)間尺度的主周期[24]。小波方差計(jì)算公式為:

(10)

因水文水資源系統(tǒng)具有時(shí)變特性,所以水文時(shí)間序列常表現(xiàn)出非平穩(wěn)特性。為詳細(xì)刻畫水文時(shí)間序列的非平穩(wěn)性,需要一種可同時(shí)在時(shí)域和頻域上進(jìn)行分析的函數(shù)。Morlet提出的小波分析滿足這一要求,故本研究選用Morlet小波函數(shù)對徑流周期性特征進(jìn)行分析。

2 結(jié)果分析與討論

2.1 徑流年內(nèi)變化特征

基于燕山水庫1955—2018年入庫徑流序列資料,首先運(yùn)用四元一次方程擬合出洛倫茲曲線(圖1),得出曲線的函數(shù)表達(dá)式,用定積分求出洛倫茲曲線與橫坐標(biāo)軸圍成的面積之后,計(jì)算得出基尼系數(shù)為0.62,表明燕山水庫的年內(nèi)徑流量分配極度不均勻。

圖1 徑流量年內(nèi)分配洛倫茨曲線

不同時(shí)段徑流量的年內(nèi)變化如圖2所示。由圖2可知:①1—6月和10—12月的徑流量較小,變化緩慢,汛期(7—9月)的徑流量較大且變化幅度明顯;燕山水庫徑流量存在季節(jié)性分配不均的特點(diǎn),其中春、夏、秋、冬季徑流量占全年徑流量的比例分別為11%、63%、21%、5%,夏季徑流量占比最大,冬季徑流量占比最小。②不同時(shí)段的月徑流量過程線均呈單峰型,20世紀(jì)80年代后月徑流量峰值的發(fā)生時(shí)間出現(xiàn)前移。1955—1979年的徑流量峰值發(fā)生在8月,而1980—2000年和2001—2018年的徑流量峰值均出現(xiàn)在7月。③不同時(shí)段的年徑流量總體呈減少趨勢,1980—2000年和2001—2018年的徑流量較1955—1979年的分別減少了9%、27%,且汛期各月徑流量的變化最為明顯。

圖2 不同時(shí)段徑流量的年內(nèi)變化

2.2 徑流年際變化特征

2.2.1 徑流變化趨勢分析

圖3為不同時(shí)間尺度下的燕山水庫徑流量變化過程。從圖3(a)可以看出,燕山水庫徑流量以0.026億m3/年的線性傾向率下降。從圖3(b)—3(e)可以看出,春、夏、秋、冬季徑流量呈現(xiàn)不同程度的下降趨勢,其傾向率分別為0.004億、0.017億、0.005億、0.000 4億m3/年。1975年8月淮河上游的洪汝河和沙潁河流域發(fā)生了歷史少見的特大暴雨,使燕山水庫入庫徑流量達(dá)到峰值。因此,夏季徑流量的峰值和年徑流量的峰值出現(xiàn)時(shí)間點(diǎn)吻合。從5年滑動平均曲線可以看出:研究區(qū)全年和季節(jié)徑流量均在趨勢線附近上下波動,呈現(xiàn)“增—減”的波動變化;5年滑動平均曲線均有1個(gè)明顯的峰值,20世紀(jì)50年代至20世紀(jì)末期的波動較大,徑流量變化劇烈;2000年后的波動較小,徑流量變化緩慢。

圖3 1955—2018年燕山水庫徑流量變化過程

采用M-K檢驗(yàn)法進(jìn)一步診斷燕山水庫全年和各季節(jié)徑流量變化趨勢的顯著性,得到Z值依次為-1.84、-1.67、-1.25、-0.89和-0.64。由此可以看出,Z值均小于0,且均未超過0.05置信水平閾值(±1.96)。因此,全年和各季節(jié)徑流量均呈非顯著性下降趨勢。

2.2.2 徑流突變分析

采用M-K突變檢驗(yàn)法對燕山水庫徑流序列進(jìn)行突變性分析,其結(jié)果如圖4所示。由圖4可得,UF和UB大多未超出顯著性水平α=0.05臨界值。

圖4 1955—2018年燕山水庫年徑流突變點(diǎn)檢驗(yàn)結(jié)果

由圖4(a)可知:UF值在1955—2018年基本小于零,只有在1984年大于零;UF和UB交點(diǎn)較少,明顯的突變點(diǎn)只出現(xiàn)在1963年和1993年,且1993年后波動變化并超出臨界線,說明燕山水庫徑流在1963年和1993年發(fā)生突變。由圖4(b)可知:1955—1974年和1981—2018年的徑流量呈下降趨勢,1975—1980年的徑流量呈上升趨勢,但上升趨勢不明顯;1957—1963年超出臨界值,UF和UB的交點(diǎn)出現(xiàn)在1963年和1981年,說明春季徑流量在1963年和1981年出現(xiàn)突變點(diǎn)。由圖4(c)可知:夏季徑流量在2003—2007年呈上升趨勢,其余時(shí)間段均呈下降趨勢;UF和UB的交點(diǎn)出現(xiàn)在1964年和1986年,說明夏季徑流量在1964年和1986發(fā)生突變。由圖4(d)可知:秋季徑流量在1963—1977年和1981—1991年呈上升趨勢,在1955—1962和1982—2018年呈下降趨勢;UF和UB的交點(diǎn)出現(xiàn)在1963年和1986年,說明秋季徑流量在1963年和1986年發(fā)生突變。由圖4(e)可知:冬季徑流量在分界點(diǎn)上下波動變化,UF和UB的交點(diǎn)出現(xiàn)在1963年,說明冬季徑流量在1963年發(fā)生突變。

從以上分析可知:燕山水庫1955—2018年的年徑流量和季節(jié)徑流量均存在多個(gè)突變點(diǎn),年徑流量突變點(diǎn)大致發(fā)生在1963年和1993年,各季徑流量突變點(diǎn)發(fā)生時(shí)間不同,大致發(fā)生在1963年和1986年。

2.3 徑流周期性特征

本文應(yīng)用小波分析法計(jì)算得到燕山水庫不同時(shí)間尺度下的徑流序列的小波變換等值線圖和小波方差圖,分別如圖5—9所示。小波變換實(shí)部分布圖(圖5—9的(a)圖)的下半部等值線相對密集,為高頻曲線,其對應(yīng)小尺度周期的震蕩;而上半部等值線相對稀疏,為低頻曲線,其對應(yīng)大尺度的周期震蕩。

圖5 1955—2018年全年徑流量小波分析結(jié)果

由圖5(a)可知,年徑流量存在明顯的周期變化,在18～24年、12～15年和3～8年時(shí)間尺度上分別有5次、6次和12次的“豐—枯”交替變化。由圖6(a)可知,春季徑流量在26～28年尺度上經(jīng)歷4次“枯—豐”交替變化,在12～18年以及3～11年尺度上分別有6次和10次“豐—枯”交替變化。由圖7(a)可知,夏季徑流量在26～30年尺度上經(jīng)歷4次“枯—豐”變化,在18～24年尺度上經(jīng)歷6次“豐—枯”交替變化,對于更小的時(shí)間尺度其“豐—枯”交替變化更頻繁。由圖8(a)可知,秋季徑流量在24～30年尺度上有4次“枯—豐”變化,8～16年尺度上有6次“豐—枯”交替變化,對于較小的尺度交替變化更頻繁。由圖9(a)可知,冬季徑流量在2～8年時(shí)間尺度上經(jīng)歷7次“豐—枯”交替變化。

圖6 1955—2018年春季徑流量小波分析結(jié)果

圖7 1955—2018年夏季徑流量小波分析結(jié)果

圖8 1955—2018年秋季徑流量小波分析結(jié)果

圖9 1955—2018年冬季徑流量小波分析結(jié)果

由圖5—9的(b)圖可以看出,燕山水庫徑流量存在顯著的多時(shí)間尺度變化特點(diǎn),年徑流量存在15年的主周期,春、夏、秋、冬季徑流量分別存在14年、12年、11年和3年的主周期。

綜上可知,燕山水庫四季的徑流量變化并非以一種固定周期的形式發(fā)生,而是以各種時(shí)間尺度相互嵌套的結(jié)構(gòu)出現(xiàn),1955—2018年燕山水庫徑流量變化周期總體上呈現(xiàn)出18～24年的大尺度變化特點(diǎn)和3～6年、8～16年的小尺度變化特點(diǎn);年徑流量存在15年的主周期,春、夏、秋、冬季徑流量分別存在14年、12年、11年和3年的主周期。

3 結(jié)論

本文選取燕山水庫1955—2018年還原后的徑流資料,運(yùn)用基尼系數(shù)、M-K檢驗(yàn)法和小波分析法對徑流趨勢性和周期性變化進(jìn)行了深入研究,主要結(jié)論如下:

1)入庫徑流量的年內(nèi)分配和季節(jié)分配具有分配不均的特點(diǎn),其中春、夏、秋、冬季徑流量占全年的比例分別為11%、63%、21%、5%。1980年后,月徑流量峰值的發(fā)生時(shí)間由1979年之前的8月前移至7月。1980—2000年和2001—2018年的徑流量較1955—1979年的分別減少9%和27%。

2)1955—2018年燕山水庫全年和各季徑流均以一定的速率下降,且均未通過α=0.05的顯著性檢驗(yàn),呈非顯著性下降趨勢;年徑流量在1963年、1993年發(fā)生突變,季節(jié)徑流突變點(diǎn)大致發(fā)生在1963年和1986年。

3)燕山水庫徑流量年際變化具有多重時(shí)間尺度上的變化特征,總體上呈現(xiàn)出18～24年的大尺度和3～6年、8～16年的小尺度變化特征;年徑流存在15年的主周期,春、夏、秋、冬季徑流量分別存在14年、12年、11年和3年的主周期。