王雪彥, 范思聰, 王富強,3, 石家豪,3
(1.華北水利水電大學,河南 鄭州 450046; 2.河南省燕山水庫管理局,河南 平頂山 467224; 3.河南省黃河流域水資源節(jié)約集約利用重點實驗室,河南 鄭州 450046)
受氣候、地貌、植被等自然條件以及人類活動耦合作用的影響,流域水文循環(huán)有所加速。徑流是水文循環(huán)不可或缺的一部分,是地表水資源的主要形式,對人類生存和發(fā)展起至關重要的作用[1-2]。變化環(huán)境下水庫徑流及地表水資源量的變化將對區(qū)域水資源管理帶來嚴峻挑戰(zhàn)。因此,開展水庫徑流演變特征研究對區(qū)域水資源合理分配及管控具有重要意義。
目前,關于徑流演變特征的研究主要集中于徑流的年代際變化趨勢、突變特征和周期性分析等方面,相應的研究方法頗多,其中:趨勢性主要利用基尼系數[3]、重標極差(R/S)[4]、M-K趨勢分析[5]和Spearman秩次相關檢驗[6]等方法進行分析;突變特征主要通過滑動T檢驗法和累計距平分析方法[7]等進行分析;周期性主要利用最大熵譜分析法[8]、人工神經網絡[9]和奇異譜分析[10]等進行分析。但上述分析方法自身存在一定缺陷:①只在頻域上具備局部化性質;②在診斷突變點過程中缺乏數學上的嚴謹性[11-12]。20世紀80年代初,以Fourier分析方法為基礎發(fā)展起來的小波分析技術,選用一種窗口大小可改變、位置可動的變窗進行譜分析,從而同時滿足信號時域和頻域局部化的要求,很好地克服了上述傳統(tǒng)譜分析方法的不足,故此研究方法在水文時間序列周期性分析中被廣泛應用[13-14]。受天氣系統(tǒng)、流域下墊面和其他因素的影響,徑流、暴雨、洪水等水文要素時間序列幾乎都是非平穩(wěn)的,內部變化極為復雜,表現(xiàn)出多時間尺度的變化特征,利用小波變換分析技術可揭示水文要素時間序列在時域和頻域上的局域變化特征,為徑流分析和徑流預報提供理論依據。
國內學者對沙潁河流域徑流演變特征進行了大量的研究,如:甘容等[15]以沙潁河流域為研究對象,選取改進的SWAT模型模擬了徑流過程,評價了模擬效果,分析了沙潁河降雨及徑流變化過程;戴韻秋等[16]運用線性趨勢法等分析了沙潁河上游降雨及徑流序列的演變特征,量化了氣候、人類活動這兩個因素對沙潁河流域徑流變化的影響。以往的研究多集中于沙潁河流域徑流量演變特征方面,而有關沙潁河主要支流上的燕山水庫徑流演變過程的研究較為單一。從目前已有研究成果來看,有關燕山水庫徑流的研究主要是分析降水、人類活動等影響因子對徑流的影響以及其多年趨勢性變化特征[17],關于徑流的季節(jié)性變化和周期性研究相對薄弱。鑒于此,本文結合基尼系數及Mann-Kendall檢驗法分析了燕山水庫入庫徑流量的年內分配、年際和各季節(jié)變化趨勢,并基于Morlet小波分析法分析了多時間尺度下的變化特征和演變規(guī)律,以期為水庫水資源合理開發(fā)與利用提供依據。
燕山水庫位于沙潁河主要支流澧河上游的干江河上,壩址位于官寨水文站和官寨(二)水文站之間,控制流域面積約為1 170 km2,官寨水文站控制流域面積約為1 130 km2,官寨(二)水文站控制流域面積約為1 200 km2。本文基于官寨水文站1955—1993年共計39年的實測徑流資料、官寨(二)水文站1994—2008年共計14年的實測徑流資料、燕山水庫水文站2009—2018年共計10年的實測徑流資料進行分析。為使資料一致,對燕山水庫實測徑流數據進行還原,采用還原后的1955—2018年的逐月徑流數據,分析水庫的徑流變化特征。
1.2.1 基尼系數
1912年意大利統(tǒng)計與社會學者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服從獨立同分布,構造統(tǒng)計量Z:
(1)
其中:
(2)
(3)
sgn為符號函數,若統(tǒng)計量Z>0認為時間序列有增加趨勢;若統(tǒng)計量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年代發(fā)展起來的小波分析技術同時在時域和頻域上具有良好的局部化功能特性,可滿足對水文時間序列的局部化分析要求,解析其內部精細結構,在分析序列的多時間尺度特性過程中,可有效識別主要頻率成分和讀取局部化信息,有利于徑流序列的多種頻率成分分析,可將其廣泛應用于分析水文時間序列周期變化特性和水文預測等方面[22-23]。設時間序列f(t)∈L2(R)(L2(R)表示能量有限)的連續(xù)小波變換為:
(8)
(9)
小波方差反映的是徑流時間序列的能量波動隨時間尺度a的變化情況,以此來精準判斷信號中不同時間尺度的主周期[24]。小波方差計算公式為:
(10)
因水文水資源系統(tǒng)具有時變特性,所以水文時間序列常表現(xiàn)出非平穩(wěn)特性。為詳細刻畫水文時間序列的非平穩(wěn)性,需要一種可同時在時域和頻域上進行分析的函數。Morlet提出的小波分析滿足這一要求,故本研究選用Morlet小波函數對徑流周期性特征進行分析。
基于燕山水庫1955—2018年入庫徑流序列資料,首先運用四元一次方程擬合出洛倫茲曲線(圖1),得出曲線的函數表達式,用定積分求出洛倫茲曲線與橫坐標軸圍成的面積之后,計算得出基尼系數為0.62,表明燕山水庫的年內徑流量分配極度不均勻。
圖1 徑流量年內分配洛倫茨曲線
不同時段徑流量的年內變化如圖2所示。由圖2可知:①1—6月和10—12月的徑流量較小,變化緩慢,汛期(7—9月)的徑流量較大且變化幅度明顯;燕山水庫徑流量存在季節(jié)性分配不均的特點,其中春、夏、秋、冬季徑流量占全年徑流量的比例分別為11%、63%、21%、5%,夏季徑流量占比最大,冬季徑流量占比最小。②不同時段的月徑流量過程線均呈單峰型,20世紀80年代后月徑流量峰值的發(fā)生時間出現(xiàn)前移。1955—1979年的徑流量峰值發(fā)生在8月,而1980—2000年和2001—2018年的徑流量峰值均出現(xiàn)在7月。③不同時段的年徑流量總體呈減少趨勢,1980—2000年和2001—2018年的徑流量較1955—1979年的分別減少了9%、27%,且汛期各月徑流量的變化最為明顯。
圖2 不同時段徑流量的年內變化
2.2.1 徑流變化趨勢分析
圖3為不同時間尺度下的燕山水庫徑流量變化過程。從圖3(a)可以看出,燕山水庫徑流量以0.026億m3/年的線性傾向率下降。從圖3(b)—3(e)可以看出,春、夏、秋、冬季徑流量呈現(xiàn)不同程度的下降趨勢,其傾向率分別為0.004億、0.017億、0.005億、0.000 4億m3/年。1975年8月淮河上游的洪汝河和沙潁河流域發(fā)生了歷史少見的特大暴雨,使燕山水庫入庫徑流量達到峰值。因此,夏季徑流量的峰值和年徑流量的峰值出現(xiàn)時間點吻合。從5年滑動平均曲線可以看出:研究區(qū)全年和季節(jié)徑流量均在趨勢線附近上下波動,呈現(xiàn)“增—減”的波動變化;5年滑動平均曲線均有1個明顯的峰值,20世紀50年代至20世紀末期的波動較大,徑流量變化劇烈;2000年后的波動較小,徑流量變化緩慢。
圖3 1955—2018年燕山水庫徑流量變化過程
采用M-K檢驗法進一步診斷燕山水庫全年和各季節(jié)徑流量變化趨勢的顯著性,得到Z值依次為-1.84、-1.67、-1.25、-0.89和-0.64。由此可以看出,Z值均小于0,且均未超過0.05置信水平閾值(±1.96)。因此,全年和各季節(jié)徑流量均呈非顯著性下降趨勢。
2.2.2 徑流突變分析
采用M-K突變檢驗法對燕山水庫徑流序列進行突變性分析,其結果如圖4所示。由圖4可得,UF和UB大多未超出顯著性水平α=0.05臨界值。
圖4 1955—2018年燕山水庫年徑流突變點檢驗結果
由圖4(a)可知:UF值在1955—2018年基本小于零,只有在1984年大于零;UF和UB交點較少,明顯的突變點只出現(xiàn)在1963年和1993年,且1993年后波動變化并超出臨界線,說明燕山水庫徑流在1963年和1993年發(fā)生突變。由圖4(b)可知:1955—1974年和1981—2018年的徑流量呈下降趨勢,1975—1980年的徑流量呈上升趨勢,但上升趨勢不明顯;1957—1963年超出臨界值,UF和UB的交點出現(xiàn)在1963年和1981年,說明春季徑流量在1963年和1981年出現(xiàn)突變點。由圖4(c)可知:夏季徑流量在2003—2007年呈上升趨勢,其余時間段均呈下降趨勢;UF和UB的交點出現(xiàn)在1964年和1986年,說明夏季徑流量在1964年和1986發(fā)生突變。由圖4(d)可知:秋季徑流量在1963—1977年和1981—1991年呈上升趨勢,在1955—1962和1982—2018年呈下降趨勢;UF和UB的交點出現(xiàn)在1963年和1986年,說明秋季徑流量在1963年和1986年發(fā)生突變。由圖4(e)可知:冬季徑流量在分界點上下波動變化,UF和UB的交點出現(xiàn)在1963年,說明冬季徑流量在1963年發(fā)生突變。
從以上分析可知:燕山水庫1955—2018年的年徑流量和季節(jié)徑流量均存在多個突變點,年徑流量突變點大致發(fā)生在1963年和1993年,各季徑流量突變點發(fā)生時間不同,大致發(fā)生在1963年和1986年。
本文應用小波分析法計算得到燕山水庫不同時間尺度下的徑流序列的小波變換等值線圖和小波方差圖,分別如圖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年的主周期。
綜上可知,燕山水庫四季的徑流量變化并非以一種固定周期的形式發(fā)生,而是以各種時間尺度相互嵌套的結構出現(xiàn),1955—2018年燕山水庫徑流量變化周期總體上呈現(xiàn)出18~24年的大尺度變化特點和3~6年、8~16年的小尺度變化特點;年徑流量存在15年的主周期,春、夏、秋、冬季徑流量分別存在14年、12年、11年和3年的主周期。
本文選取燕山水庫1955—2018年還原后的徑流資料,運用基尼系數、M-K檢驗法和小波分析法對徑流趨勢性和周期性變化進行了深入研究,主要結論如下:
1)入庫徑流量的年內分配和季節(jié)分配具有分配不均的特點,其中春、夏、秋、冬季徑流量占全年的比例分別為11%、63%、21%、5%。1980年后,月徑流量峰值的發(fā)生時間由1979年之前的8月前移至7月。1980—2000年和2001—2018年的徑流量較1955—1979年的分別減少9%和27%。
2)1955—2018年燕山水庫全年和各季徑流均以一定的速率下降,且均未通過α=0.05的顯著性檢驗,呈非顯著性下降趨勢;年徑流量在1963年、1993年發(fā)生突變,季節(jié)徑流突變點大致發(fā)生在1963年和1986年。
3)燕山水庫徑流量年際變化具有多重時間尺度上的變化特征,總體上呈現(xiàn)出18~24年的大尺度和3~6年、8~16年的小尺度變化特征;年徑流存在15年的主周期,春、夏、秋、冬季徑流量分別存在14年、12年、11年和3年的主周期。