莫崇勛，班華珍，謝燕平，何嘉奇，阮俞理，孫桂凱

(1. 廣西大學(xué)土木建筑工程學(xué)院，南寧 530004；2. 工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室，南寧 530004；3. 廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室，南寧 530004)

0 引 言

目前，關(guān)于洪峰流量的研究成果多聚焦于場(chǎng)次洪峰流量[1-5]，較少涉及洪峰流量序列演變規(guī)律方面?？v有研究，也多集中于非巖溶區(qū)，如任健等[6]運(yùn)用HHT分析法研究黃河花園口站年最大洪峰流量序列的多時(shí)間尺度特征；唐權(quán)輝等[7]基于北江干流4個(gè)水文站資料建立年最大洪峰流量序列，利用Mann-Kendall和EMD法對(duì)序列進(jìn)行特征分析；孫思瑞等[8]采用水文變異診斷系統(tǒng)對(duì)洞庭湖三口的年最大洪峰流量和年最高洪峰水位序列進(jìn)行變異診斷，結(jié)果表明各站點(diǎn)洪峰流量多呈現(xiàn)下降的趨勢(shì)變異。這些研究都是針對(duì)非巖溶區(qū)進(jìn)行的，而巖溶地貌作為世界上最脆弱的地貌類型之一，對(duì)環(huán)境變化的響應(yīng)更為敏感。隨著氣候變化和人類活動(dòng)影響的加劇，對(duì)巖溶區(qū)流域的洪峰流量序列進(jìn)行特征分析具有一定的研究意義。鑒于此，論文基于巖溶區(qū)壩首水文站1963-2014年實(shí)測(cè)逐日入庫(kù)流量，運(yùn)用水文變異診斷系統(tǒng)中的多種統(tǒng)計(jì)方法和Morlet小波分析法，研究年最大洪峰流量序列的演變規(guī)律，以期為巖溶區(qū)流域澄碧河水庫(kù)制定防洪減災(zāi)對(duì)策和水環(huán)境治理措施提供科學(xué)依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)來(lái)源

澄碧河隸屬于珠江流域西江水系右江干流的一級(jí)支流，發(fā)源于廣西百色市凌云縣青龍山北麓一帶，途經(jīng)弄林、浩坤、平塘、百色右江區(qū)等地，水流呈明伏流交替出流，至那東村為澄碧河水庫(kù)壩址。水庫(kù)壩址以上流域集雨面積2 000 km2，流域地貌以弄林為界分為兩個(gè)部分，浩坤～弄林以上為典型的巖溶峰林區(qū)，浩坤以下為非巖溶區(qū)。巖溶區(qū)流域，由于石灰?guī)r裸露，在長(zhǎng)期風(fēng)化作用下，形成大量的裂隙、溶溝、漏斗和落水洞，雨洪過(guò)程中，地面徑流通過(guò)這些通道大量涌入地下，形成地下徑流，短歷時(shí)內(nèi)對(duì)洪峰起到削峰作用。流域內(nèi)的浩坤溶洞伏流河形成天然水庫(kù)，庫(kù)容約3 億m3，對(duì)澄碧河水庫(kù)起調(diào)節(jié)作用。流域形狀近似矩形，地勢(shì)西北高而東南低，流域概況圖如圖1。

圖1 澄碧河流域水系圖Fig.1 Chengbi river Basin

澄碧河流域內(nèi)共有4個(gè)水文站，分別為壩首、平塘、浩坤和下甲水文站。從地理位置看，浩坤站和下甲站較能體現(xiàn)巖溶區(qū)的特性，但因這兩個(gè)水文站建站時(shí)間較晚，水文序列均不足30年，平塘站水文序列不夠完整。壩首水文站位于澄碧河流域出口斷面處，站點(diǎn)數(shù)據(jù)系列足夠長(zhǎng)，可靠性高，且一定程度上能表征巖溶區(qū)特性。故論文基于廣西巖溶區(qū)流域澄碧河水庫(kù)壩首水文站1963-2014年共52 a的日均入庫(kù)流量數(shù)據(jù)，采用年最大值法對(duì)洪峰流量進(jìn)行選樣，得到52 a的年最大洪峰流量序列。

2 研究方法

2.1 趨勢(shì)和突變檢驗(yàn)方法

目前，可用于檢驗(yàn)水文時(shí)間序列變異情況的方法很多，但因水文時(shí)間序列(特別是洪水序列)變異形式相對(duì)復(fù)雜，單一檢驗(yàn)方法往往可靠性較低。故本文采用謝平等[9]提出的水文變異診斷系統(tǒng)中多種統(tǒng)計(jì)方法，對(duì)年最大洪峰流量序列的趨勢(shì)成分和突變(跳躍)成分進(jìn)行檢驗(yàn)分析。檢驗(yàn)過(guò)程分為3個(gè)階段：①初步診斷，選取線性趨勢(shì)法、滑動(dòng)平均法和Hurst系數(shù)法對(duì)序列的趨勢(shì)和突變進(jìn)行檢測(cè)，定性分析序列的變異情況；②詳細(xì)診斷，選用線性回歸法、Spearman秩次相關(guān)法和Man-Kendall秩次相關(guān)法共3種方法分析序列的趨勢(shì)顯著性，結(jié)合Mann-Kendall檢驗(yàn)法、Lee-Heghinan檢驗(yàn)法、滑動(dòng) T 檢驗(yàn)法、滑動(dòng)秩和檢驗(yàn)法、有序聚類檢驗(yàn)法、累積距平檢驗(yàn)法、滑動(dòng)游程檢驗(yàn)法和滑動(dòng)F檢驗(yàn)法共8種變異診斷方法對(duì)序列的突變進(jìn)行分析；③綜合診斷，對(duì)詳細(xì)診斷結(jié)果進(jìn)行匯總，趨勢(shì)檢驗(yàn)結(jié)果一般較統(tǒng)一，對(duì)顯著性進(jìn)行求和即可，若結(jié)果大于等于1則顯著，反之，則不顯著；突變檢驗(yàn)結(jié)果相對(duì)復(fù)雜，除了進(jìn)行顯著性綜合外，還需進(jìn)行突變權(quán)重綜合，論文以所有可能突變點(diǎn)權(quán)重最大者作為最可能突變點(diǎn)。上述方法的具體過(guò)程見相關(guān)文獻(xiàn)[10-13, 15]。

2.2 Morlet小波分析法

小波分析是廣泛應(yīng)用于研究水文序列周期性規(guī)律的方法，Morlet小波分析的基本原理及計(jì)算步驟[14]如下。

小波函數(shù)ψ(t)∈L2(R)且滿足：

(1)

式中：ψ(t)為小波基函數(shù)，其作用是通過(guò)尺度伸縮和在時(shí)間軸上的平移構(gòu)成一簇函數(shù)系。

(2)

式中：ψa,b(t)為基小波；a為尺度因子；b為時(shí)間平移因子。

由于年最大洪峰流量序列是離散的，設(shè)函數(shù)f(kΔt)(k=1，2，…，N)，Δt為時(shí)間間隔。則離散序列小波變化表示為：

(3)

式中：Wf(a,b)為小波變化系數(shù)；其余同上。

對(duì)時(shí)間域上關(guān)于a的全部小波系數(shù)的平方進(jìn)行積分，即為小波方差：

(4)

3 結(jié)果與分析

3.1 趨勢(shì)性分析

采用一元線性法和5 a滑動(dòng)平均法對(duì)年最大洪峰流量序列的趨勢(shì)性進(jìn)行初步檢驗(yàn)，結(jié)果見圖2。由圖2可知，年最大洪峰流量序列整體呈微弱下降趨勢(shì)，年均下降率為6.04 m3/s；1979年之前，洪峰流量序列大多處于均值上方，其中1967年達(dá)到最大值，為3 399 m3/s，1977年次之，為3 075 m3/s；1979年之后，洪峰流量序列在均值附近平穩(wěn)波動(dòng)，其中2013年達(dá)到最小值，為301.77 m3/s。為進(jìn)一步檢驗(yàn)序列趨勢(shì)變化的顯著性，采用線性回歸法、Spearman秩次相關(guān)法及Man-Kendall秩次相關(guān)法進(jìn)行詳細(xì)診斷，顯著性α均取為0.05，結(jié)果見表1(圖略)，經(jīng)計(jì)算，線性回歸法和Spearman秩次相關(guān)法的統(tǒng)計(jì)量T值分別為-1.1和-0.15，絕對(duì)值均小于α=0.05顯著性水平臨界值1.64；Man-Kendall秩次相關(guān)法的統(tǒng)計(jì)量U為-0.209，絕對(duì)值小于α=0.05顯著性水平臨界值1.96；記不顯著性為“-1”，顯著性為“1”，則綜合診斷趨勢(shì)顯著性為“-3”，表明年最大洪峰流量序列的變化趨勢(shì)不顯著，即序列呈不顯著的下降趨勢(shì)。

圖2 年最大洪峰流量序列趨勢(shì)線Fig.2 The trend of annual maximum peak discharge series

診斷方法顯著值(|T|/|U|)檢驗(yàn)結(jié)果趨勢(shì)診斷線性回歸法T=-1.1不顯著(-1)Spearman秩次相關(guān)檢驗(yàn)法T=-0.15不顯著(-1)Kendall秩次相關(guān)檢驗(yàn)法U=-0.209不顯著(-1)綜合顯著性--不顯著(-3)

3.2 突變性分析

采用R/S分析原理計(jì)算出年最大洪峰流量序列的Hurst系數(shù)[15]h值為0.69，表明序列呈現(xiàn)中變異，初步判定該序列可能存在趨勢(shì)或突變(跳躍)變異。由3.1節(jié)知，趨勢(shì)變化不顯著，故主要對(duì)突變(跳躍)成分變異情況進(jìn)行詳細(xì)診斷。采用Man-Kendall檢驗(yàn)法、Lee-Heghinan法等8種變異檢測(cè)方法對(duì)序列進(jìn)行分析，顯著性α均取為0.05，Man-Kendall檢驗(yàn)法結(jié)果如圖3所示。由圖3可知，統(tǒng)計(jì)量UF和UB在α=0.05顯著性水平臨界線內(nèi)存在多個(gè)交點(diǎn)，分別為1966年、1971年、1980年、1994年、1998年和2011年，因1966、1971、1998和2011年過(guò)于接近序列兩端，可靠性相對(duì)較低，故選取1980年和1994年為可能突變點(diǎn)；其余方法的檢測(cè)結(jié)果見表2。由表2可知，各方法檢測(cè)結(jié)果不一致，突變點(diǎn)及其顯著性均有差異，故進(jìn)行綜合診斷分析，結(jié)果見表3。由表3可知，突變(跳躍)點(diǎn)為1979年，其綜合權(quán)重達(dá)到最大值為0.375，綜合顯著性為3，表明突變顯著。

圖3 Mann-Kendall變異檢測(cè)Fig.3 Mann-Kendall Mutation detection

綜上可知，52 a來(lái)，年最大洪峰流量序列整體呈不顯著的下降趨勢(shì)，并于1979年發(fā)生以跳躍(突變)形式為主的變異。

3.3 周期性分析

利用Morlet小波分析序列的周期，結(jié)合Matlab軟件繪制出如圖4和圖5所示。圖4(a)為小波系數(shù)實(shí)部等值線圖，小波系數(shù)正相位表示洪峰流量處于偏高值，負(fù)相位表示處于偏低值；由圖4可知，年最大洪峰流量主要存在5～10、11～20和25～30 a的時(shí)間尺度周期性變化規(guī)律。對(duì)應(yīng)于時(shí)間尺度的交替，序列出現(xiàn)明顯變化，其中11～20 a尺度震蕩最強(qiáng)烈，并以15 a為中心時(shí)間尺度，正負(fù)交替持續(xù)到1985年前后；5～10 a表現(xiàn)較明顯，其中心時(shí)間尺度為8 a左右，正負(fù)交替持續(xù)到1975年；25～30 a尺度出現(xiàn)得較微弱，其中心時(shí)間尺度為28 a左右，正負(fù)交替持續(xù)到1985年左右。1985年之后，洪峰流量未出現(xiàn)明顯的時(shí)間尺度周期性規(guī)律，但其相位幾乎全為正值，表明洪峰流量于1985年之后幾乎全部大于均值。圖4(b)為小波系數(shù)模平方圖，反映不同變化周期對(duì)應(yīng)的能量密度在時(shí)間域的分布強(qiáng)弱，對(duì)應(yīng)值越大，能量越強(qiáng)。由圖可知，序列存在兩個(gè)相對(duì)顯著的能量聚集中心：5～10和11～20 a的時(shí)間尺度，表明洪峰流量變化對(duì)應(yīng)于這兩個(gè)尺度周期性比較顯著。其中以11～20 a時(shí)間尺度震蕩能量變化最強(qiáng)，主要發(fā)生在1960-1985年，震蕩中心在1967年；5～10 a時(shí)間尺度十分突出，主要發(fā)生在1960-1975年，震蕩中心也為1967年。圖5表示波動(dòng)能量尺度分布的小波方差圖，由圖可知，年最大洪峰流量序列存在3個(gè)比較明顯的峰值。其中，最大峰值對(duì)應(yīng)15a的時(shí)間尺度，表明15 a的周期震蕩最強(qiáng)烈，即為年最大洪峰流量變化過(guò)程的主周期；8和28 a對(duì)應(yīng)著第二峰值，即為序列的次周期。

表2 年最大洪峰流量序列詳細(xì)診斷結(jié)果Tab.2 The detail diagnosis of annual maximum peak discharge series

表3 年最大洪峰流量序列綜合診斷結(jié)果Tab.3 The comprehensive diagnosis of annual maximum peak discharge series

圖4 年最大洪峰流量序列Morlet小波變化系數(shù)時(shí)頻分布Fig.4 The time-frequency distribution of the Morlet wavelet coefficients of annual maximum peak discharge

圖5 年最大洪峰流量序列小波方差圖Fig.5 The wavelet variance of annual maximum peak discharge

4 討論與小結(jié)

4.1 討 論

影響徑流量的因素眾多，其中氣候變化和人類活動(dòng)是兩個(gè)重要因素[16]，洪水是徑流的一個(gè)極端表現(xiàn)，洪峰是洪水過(guò)程的一個(gè)重要指標(biāo)[17]，降雨、徑流、氣溫等水文氣象要素變化一定程度上會(huì)引起洪峰流量的改變。澄碧河流域是典型的巖溶區(qū)與非巖溶區(qū)結(jié)合的流域，弄林以上巖溶區(qū)以峰從洼地為主，其特點(diǎn)為峰林峰從連座、洼地閉塞，土壤覆蓋層稀薄，巖溶裂隙、漏斗、溶溝和落水洞很發(fā)育，形成具有不規(guī)則性、多層性和連通性的地下河網(wǎng)。大部分洼地的消水溶洞與地下暗河、管道連接，形成暢排型的地下河網(wǎng)通道，能快速補(bǔ)給地下徑流；少部分洼地雖不存在消水洞，但土壤覆蓋層松散薄弱，垂向補(bǔ)給速率快，雨洪過(guò)程，洼地蓄水，短歷時(shí)調(diào)節(jié)中，起到削減洪峰作用[18, 19]。整個(gè)澄碧河流域近50 a年來(lái)氣候變化較明顯，澄碧河流域降雨量、氣溫和蒸發(fā)量分別呈下降、上升和下降趨勢(shì)，降雨量在20世紀(jì)60年代、80年代和 21世紀(jì)發(fā)生強(qiáng)變異[20, 21]， 這與壩首站年最大洪峰序列呈不顯著下降趨勢(shì)，并于1979年發(fā)生顯著突變基本一致。此外，20世紀(jì)70年代至80年代，澄碧河流域無(wú)巨大人類活動(dòng)影響(如建設(shè)大型水利工程、大面積植被砍伐等)。據(jù)此，可推斷影響年最大洪峰流量序列發(fā)生變異的主要原因可能是氣候變化。

相關(guān)研究指出，太陽(yáng)黑子活動(dòng)與全球氣候異常，與大氣環(huán)流和熱帶氣旋，與厄爾尼諾和拉尼娜現(xiàn)象，與各地區(qū)的降水異常、洪水災(zāi)害等均有良好的響應(yīng)關(guān)系[22]。太陽(yáng)黑子峰谷值出現(xiàn)時(shí)易發(fā)大洪水，從變化周期來(lái)看，巖溶區(qū)流域壩首站年最大洪峰流量存5～10、11～20和25～30 a的震蕩周期，與太陽(yáng)黑子的震蕩周期(11、20 a)[23]有一定的吻合度。

4.2 小 結(jié)

論文利用水文變異診斷系統(tǒng)中的多種統(tǒng)計(jì)方法和Morlet小波分析法對(duì)巖溶區(qū)澄碧河水庫(kù)52 a年最大洪峰流量進(jìn)行變化特征分析，得到了如下結(jié)論。

(1)1963-2014年，巖溶區(qū)流域年最大洪峰流量整體呈不顯著的下降趨勢(shì)，年均下降率為6.04 m3/s。

(2)水文序列于1979年發(fā)生突變，其綜合權(quán)重為0.375，綜合顯著性為3，說(shuō)明52 a來(lái)，巖溶區(qū)流域年最大洪峰流量存在跳躍形式為主的變異。

(3)1985年之前，年最大洪峰流量存在15 a的主周期，8和28 a的次周期；1985年之后，不存在明顯的周期變化規(guī)律。

(4)建議對(duì)年最大洪峰流量變化特征的影響因素分析時(shí)，結(jié)合氣候因素深入研究，探討1985年之后序列不存在周期性規(guī)律的原因。

(5)建議采用多個(gè)站點(diǎn)的數(shù)據(jù)對(duì)巖溶區(qū)流域進(jìn)行系統(tǒng)的研究分析。