巖溶區流域年最大洪峰流量變化特征分析

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

(1. 廣西大學土木建筑工程學院，南寧 530004；2. 工程防災與結構安全教育部重點實驗室，南寧 530004；3. 廣西防災減災與工程安全重點實驗室，南寧 530004)

0 引 言

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

1 研究區概況與數據來源

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

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

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

2 研究方法

2.1 趨勢和突變檢驗方法

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

2.2 Morlet小波分析法

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

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

(1)

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

(2)

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

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

(3)

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

對時間域上關于a的全部小波系數的平方進行積分，即為小波方差：

(4)

3 結果與分析

3.1 趨勢性分析

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

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

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

3.2 突變性分析

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

圖3 Mann-Kendall變異檢測Fig.3 Mann-Kendall Mutation detection

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

3.3 周期性分析

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

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

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

圖4 年最大洪峰流量序列Morlet小波變化系數時頻分布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 討論與小結

4.1 討 論

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

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

4.2 小 結

論文利用水文變異診斷系統中的多種統計方法和Morlet小波分析法對巖溶區澄碧河水庫52 a年最大洪峰流量進行變化特征分析，得到了如下結論。

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

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

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

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

(5)建議采用多個站點的數據對巖溶區流域進行系統的研究分析。