湟水流域洪水極值時間演變特征及趨勢歸因分析

楊 陽，時 璐，王 崗，段 水強(qiáng)

(青海省水文水資源勘測局，西寧 810001)

0 引 言

在氣候轉(zhuǎn)變的大趨勢下，人類社會活動受到自然環(huán)境的制約，同時自然環(huán)境也受到人類社會活動的改造和考驗，自然環(huán)境的變化致使其內(nèi)的自然現(xiàn)象也隨之發(fā)生變化，其中水文現(xiàn)象的變化是敏感的指針之一，通過對水文要素在時間和空間上變化的研究并對其引起變化的原因進(jìn)行分析，提升對有利與不利因素的甄別，使其變化正向循環(huán)發(fā)展就顯示出重要意義。由此，人類社會活動、氣候變化、水資源演變等問題都受到了社會各界的廣泛關(guān)注。青海省境內(nèi)除去大通河的湟水流域是青海省人口分布最密集、經(jīng)濟(jì)發(fā)展程度最高的地方，同時也是受人類社會活動影響最為頻繁的地區(qū)，研究分析該流域水文系統(tǒng)的變化就顯的意義重大，本文力圖從湟水流域控制站民和站洪峰流量時序的趨勢、突變以及周期變化的角度揭示湟水流域極值內(nèi)在的演變特征(或規(guī)律)。

1 研究區(qū)域概況

湟水河(青海境內(nèi)除去大通河，下同)青海人民的母親河，黃河上游重要一級支流，位于青藏高原與黃土高原銜接地帶，發(fā)源于海晏縣包呼圖山對望青海湖，北至大通山、達(dá)坂山與南界拉脊山比肩相對，兩山之間形成了湟水流域，其峽谷將盆地連接串成葫蘆狀一字排列，流域成條狀，向東延伸至甘肅最終匯入黃河，境內(nèi)全長374 km，其流域面積16 091 km2，流域平均海拔為2 500 m，大地構(gòu)造屬祁連山褶皺帶，地質(zhì)條件復(fù)雜，山區(qū)占比較大；流域偏北、西風(fēng)帶過境頻繁，盛夏受西南季風(fēng)影響，河西走廊冷空氣沿河谷倒灌湟水谷地，普遍存在山谷風(fēng)環(huán)流，風(fēng)向隨晝夜交替變換，白天從山谷吹向山坡為谷風(fēng)，夜間山坡吹向谷地為山風(fēng)，當(dāng)大尺度大氣環(huán)流較弱時地方性的山谷風(fēng)現(xiàn)象明顯[1]，流域多年平均降水深527 mm，多年平均徑流深為134 mm。流域控制站為民和水文站，控制著15 342 km2的流域面積，面積占比達(dá)95%。湟水流域水系見(圖1)。

圖1 湟水流域水系分布圖Fig.1 Water system distribution map of Huangshui Basin

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

本文數(shù)據(jù)來源于青海省氣象局、水文局、水土保持局和民和站1950-2016年的歷史最大洪峰流量實測數(shù)據(jù)，其數(shù)據(jù)已在往年的水文年鑒中刊布，資料符合客觀事實，精度滿足分析要求。

3 研究方法

本文應(yīng)用了多種方法對民和站洪峰流量時序的趨勢、突變和周期性進(jìn)行了對比分析，同時也是互相驗證的手段，如下文逐一介紹。

3.1 趨勢檢驗法

本文分別采用了線性趨勢回歸檢驗法[2]、Kendall秩次檢驗法[2]和spearman秩次檢驗法[2]對民和站歷史洪峰流量的趨勢性進(jìn)行檢驗分析，通過構(gòu)建和計算統(tǒng)計量U、T與相應(yīng)的臨界值比較后判斷民和站歷史洪峰流量的變化趨勢是否顯著。

3.1.1 線性趨勢回歸檢驗法

(1)

(2)

3.1.2 Kendall秩次檢驗法

在序列x1,x2,…,xn中，計算統(tǒng)計量U：

(3)

(4)

(5)

式中：p為序列某一元素xi大于前面xi-1,xi-2,…,x2,x1累積數(shù)的總和；n為系列長度。

3.1.3 Spear man秩次檢驗法：

分析序列x1,x2,…,xi與時序i的相關(guān)關(guān)系，秩次相關(guān)系數(shù)為：

(6)

(7)

式中：n為序列長度；Ri表示序列xi的秩次，di=Ri-i。統(tǒng)計量T服從自由度為(n-2) 的t分布。

3.2 突變點檢驗法

3.2.1 滑動T檢驗法

滑動T檢驗法[3,4]是在滑動點前后，分別抽取容量相同的n1、n2兩個子序列，觀察兩個子序列均值的差異是否顯著來檢驗突變。定義統(tǒng)計量：

(8)

(9)

3.2.2 有序聚類

有序聚類分析法[3,4]就是尋找突變前后系列離差平方和的總和最小的點即最優(yōu)分割點。 突變前后兩個系列離差平方和分別表示為：

(10)

(11)

式中：xτ為突變點τ前的水文序列均值；xn-τ為突變點τ后的水文序列均值。

Sn(τ)=Vτ+Vn-τ

(12)

式中：Sn(τ)為總離差平方和。Sn(τ)取最小值時對應(yīng)的τ即為最優(yōu)分割點。

3.2.3 R/S重標(biāo)極差分析法

R/S重標(biāo)極差分析法又稱R/S檢驗法是赫斯特在大量實證研究的基礎(chǔ)上提出的一種時間序列統(tǒng)計方法[3], 多用于突變點識別和未來趨勢預(yù)測。本文利用R/S定量識別突變年份，基本原理和方法如下：

考慮一個時間序列{ξ(t)},t=1,2,…,n對于任意正整數(shù)τ≥1,定義均值序列：

(13)

用X(t)表示累積離差：

(14)

極差R定義為：

(15)

標(biāo)準(zhǔn)差S定義為：

(16)

則有：

(17)

式中：c為某常數(shù)；H為赫斯特指數(shù)。

對上式取對數(shù)：

(18)

對于不同區(qū)間長度τ，可得到不同的R(τ)/S(τ)值，以ln(τ)為自變量、lnR(τ)/S(τ)為因變量作散點圖，采用最小二乘估計擬合直線，直線斜率即為時間序列的赫斯特指數(shù)。為能定量地識別變異點，在ln(τ)與lnR(τ)/S(τ)散點圖中，對前后兩點作斜率以此滑動類推，取斜率差異絕對值最大的點作為最可能變異點。

3.2.4 M-K檢驗法

Mann-Kendall檢驗法是世界氣象組織推薦的非參數(shù)檢驗方法，并已廣泛地用來分析水文要素時間序列的變化情況，由于該檢驗方法不要求樣本遵從一定的分布，且檢驗結(jié)果基本不受少數(shù)異常值的干擾，適合水文、氣象等非正態(tài)分布的數(shù)據(jù)[3-8]。

它的檢驗原理是假設(shè)有一時間序列如下：x1,x2,x3,…,xn，構(gòu)造一秩序列mi，mi表示xi>xj,(1≤j≤i)的樣本累積數(shù)。構(gòu)造統(tǒng)計量：

(19)

式中：mi為某一元素xi大于前面xi-1，xi-2，…,x2,x1的累積數(shù)。

在原序列隨機(jī)獨立的情況下，dk的均值、方差分別為：

(20)

(21)

然后將dk標(biāo)準(zhǔn)化，構(gòu)造統(tǒng)計量：

(22)

再計算反序列，上面計算過程，得到UBk，k=n,n-1,…,1。同樣，UB1=0。

3.3 周期性分析法

小波分析[5]是通過典型小波經(jīng)過頻率窗口的伸縮和時間窗口位置平移后，將源信號分解成由低頻到高頻在時間域組成的信號組，從而能夠在時間域內(nèi)對源信號進(jìn)行精細(xì)化的周期分析和頻率的時間定位，揭示出隱藏在時間序列中的多種變化周期，進(jìn)而能夠?qū)υ葱盘栐诓煌瑫r間尺度上的未來變化趨勢進(jìn)行判斷和甄別。本文利用Morlet小波對湟水流域民和站歷年洪峰流量數(shù)據(jù)進(jìn)行內(nèi)在的變化周期分析。

4 結(jié)果與原因分析

4.1 趨勢性

對湟水流域民和站1950-2016洪峰流量分析表明，其變化范圍在164～1 300 m3/s，最大值出現(xiàn)在1952年，最小值出現(xiàn)在2011年(圖2)，利用線性趨勢回歸檢驗、Kendall秩次檢驗和Spearman秩次檢驗三種方法對湟水流域民和站點洪峰流量分析，其統(tǒng)計量都通過了顯著水平α=0.05的信度檢驗，說明其1950-2016年洪峰流量有明顯下降趨勢，變化速率為-5.4 m3/(s·a)(表1)。

圖2 民和站1950-2016最大洪峰流量過程線Fig.2 Maximum flood discharge process of Minhe station from 1950 to 2016

4.2 突 變

采用滑動T(n1=n2=5)、有序聚類、R/S檢驗法和Mann-Kendall檢驗法[11-13]，對湟水流域民和站洪峰流量系列進(jìn)行突變分析。從圖(3)可知采用滑動T和有序聚類檢測到突變年份為1999年，R/S檢驗法檢測到突變年份為1997年和1999年，Mann-Kendall檢驗法檢測到突變年份為2000年；從圖4洪峰流量跳躍圖中可以得出突變年份在2000年左右，由湟水流域民和站歷年最大洪峰流量過程線綜合判斷，認(rèn)為湟水流域民和站洪峰流量由1997年增加到減少再到1999年的增加呈現(xiàn)周期性，而真正在2000年后呈現(xiàn)穩(wěn)定狀態(tài)。

圖3 突變檢驗圖Fig. 3 Test chart of mutation

圖4 湟水流域民和站最大洪峰流量序列跳躍變異圖Fig.4 Jump variation of maximum flood discharge sequence of Minhe station in Huangshui basin

4.3 周 期

通過Morlet復(fù)小波尺度變換和位移因子后獲得小波系數(shù)，利用該小波系數(shù)的實部、小波模、小波模方及方差等數(shù)據(jù)[8]，對湟水流域民和站最大洪峰流量系列進(jìn)行內(nèi)在周期性分析，進(jìn)而揭示最大洪峰流量在時間序列上的周期變化[7,14]。

從年最大洪峰流量小波實部等值線圖(圖5)可以看出民和站歷年最大洪峰流量在時序演化過程中存在著多時間尺度特征?？偟膩碚f，存在著25～32 a，15～25 a、7～15 a以及3～7 a的4類尺度的周期變化規(guī)律。縱觀整個時域，在25～32 a尺度上出現(xiàn)了準(zhǔn)3次震蕩；在15～25 a時間尺度上有4次震蕩。同時，還可以看出上述兩個尺度的周期變化在整個時域內(nèi)相對穩(wěn)定，尤其25～32 a尺度的震蕩周期具有全時域性； 7～15 a尺度的震蕩周期不穩(wěn)定，1950-1970年震蕩強(qiáng)烈歷經(jīng)2次周期震蕩，之后的周期不顯著。

小波系數(shù)的模和模方相當(dāng)于小波能量譜，它可以分析出不同周期的震蕩能量[7]。從年最大洪峰流量小波模和模方等值線圖(圖5)可以看出25～32 a、18～25 a、7～15 a時間尺度的能量強(qiáng)、周期顯著。

從湟水流域民和站歷年最大洪峰流量小波方差圖中(圖6)也可以看出，周期中心時間尺度32 a對應(yīng)著第一主峰也是第一主周期，其他峰值分別對應(yīng)著20、10和4 a周期中心時間尺度，與上述的小波實部等值線和能量譜的表現(xiàn)相對應(yīng)。這4個時間尺度的周期控制著年最大洪峰流量在整個時域內(nèi)的變化特征。

圖5 湟水流域民和站最大洪峰流量小波系數(shù)、模、模方等值線圖Fig.5 The isoline of wavelet coefficient, modulus, square of the modulus of Maximum flood discharge of Minhe station

圖6 小波方差和32年時間尺度小波實數(shù)過程圖Fig.6 Wavelet variance and 32 years time scale wavelet real process

4.4 變化原因分析

根據(jù)郭生練[17]研究成果，氣候變化將對洪峰流量產(chǎn)生重大影響，氣溫的升高伴隨著降水量也隨之增加，將顯著的增大洪水的頻次和洪峰流量。我們收集湟水流域9個氣象站點氣溫資料和85個具有30年以上實測資料的雨量站點，經(jīng)預(yù)處理后插補(bǔ)延長至1956-2016系列，通過泰森多邊形法和降水等值線的繪制(圖7)，分別獲取了湟水流域面氣溫和面降水。通過線性趨勢分析(表1)發(fā)現(xiàn)，民和站洪峰流量在隨著氣溫顯著升高、降水不顯著增加趨于平穩(wěn)的狀態(tài)下反而呈現(xiàn)減小的趨勢(圖2和圖8)，由于洪峰流量是地表徑流縮短時間尺度的產(chǎn)物，其產(chǎn)生的物理機(jī)制一樣，都是降水通過下墊面調(diào)蓄后產(chǎn)生的，長時序或者大時間尺度下二者變化趨勢一致(表1)，所以我們通過徑流作為洪峰流量的媒介進(jìn)一步分析，經(jīng)過該流域最大的水庫黑泉水庫蓄水變量還原后，發(fā)現(xiàn)還原后的年徑流很敏感得開始呈現(xiàn)不顯著增加趨勢(表1和圖9)，與該區(qū)域氣溫、降水的變化趨勢一致，通過該還原分析得知影響湟水流域徑流、洪峰流量的主導(dǎo)因素是下墊面的變化或者是受人類活動影響為主要誘因。

圖7 湟水流域降水等值線Fig.7 Rainfall isoline of Huangshui Basin

圖8 湟水流域氣溫和降水過程線Fig.8 Temperature and Precipitation process of Huangshui Basin

圖9 民和站實測年徑流和經(jīng)還原后民和站年徑流量過程線圖Fig.9 Measured annual runoff process and restoration runoff process of Minhe Station

表1 湟水流域各水文要素參數(shù)統(tǒng)計表Tab.1 Statistic table of Hydrological parameters in Huangshui basin

湟水流域是青海省經(jīng)濟(jì)最發(fā)達(dá)、人口分布最密集的地區(qū)，所以該地區(qū)人類活動頻繁，工程措施較多，下墊面變化影響大[16]。

4.4.1 水利工程

來自水利普查水庫工程的數(shù)量規(guī)模統(tǒng)計，自2000年以來雖然工程數(shù)量增長不明顯但是規(guī)模增長顯著，從(圖10)可知，2000年左右累積總庫容驟然上升，比之前的庫容增加了進(jìn)2.8倍，主要是北川河上游黑泉水庫在這一時期的建成運行，其也是造成民和站洪峰流量突變的主要原因，北川河上游是湟水流域的降水高值區(qū)和主要洪水來源區(qū)，黑泉水庫有著削峰滯洪的作用。

圖10 湟水流域水庫庫容累積曲線Fig.10 Accumulation curve of reservoir capacity in Huangshui Basin

圖11 湟水流域降水洪峰雙累積曲線Fig.11 Double cumulative curve of precipitation flood peak in Huangshui Basin

圖12 2016年5月與1995年6月NDVI疊加圖Fig.12 Overlay map of NDVI between May 2016 and June 1995

4.4.2 水土保持工程

根據(jù)青海省水土保持局提供資料顯示，截止2016年青海省湟水流域水土保持治理累積總面積(包括農(nóng)牧、林業(yè))約為4 650 km2，占湟水流域面積的近1/3，隨著每年水土保持工程建設(shè)的落地，其效益也逐漸呈現(xiàn)，下墊面覆蓋條轉(zhuǎn)好，使其削減洪峰、調(diào)節(jié)徑流的功能增強(qiáng)。如我省湟水谷地水土保持工程的建設(shè)實施，植被覆蓋增加，通過landsat8[9]衛(wèi)星2016年5月與1995年6月兩幅影像圖制作的NDVI對比圖顯示，湟水河谷地靠近城區(qū)兩側(cè)近山上的植被指數(shù)增加范圍明顯(圖12橘黃色部分)，說明裸地減少，植被覆蓋態(tài)勢向好。

4.4.3 人類活動影響估算

人類活動實施后(下墊面變化)湟水流域降水與民和站最大洪峰流量雙累積曲線發(fā)生變化(圖11)，在2000年左右以后直線斜率發(fā)生明顯變化，將2000年之前的回歸線作為基準(zhǔn)線，把2000年之后回歸線的斜率(也就是受人類活動影響較大的階段)模擬到基準(zhǔn)線的斜率上，獲取1956-1999,2000-2016年兩階段的洪峰流量模擬均值，對比計算出人類活動因素影響率(表2)。

表2 人類活動影響率估算表Tab.2 Estimate table of the impact of human activities

5 結(jié) 論

通過對湟水流域控制站民和站1950-2016年最大洪峰流量進(jìn)行了趨勢、突變、周期及原因分析，主要得出以下結(jié)論。

(1)民和站歷年最大洪峰流量總體呈現(xiàn)顯著下降趨勢，其下降速率為-5.4 m3/(s·a)。

(2)湟水流域民和站歷年最大洪峰流量突變點為2000年左右,1999年后最大洪峰流量發(fā)生由大到小的突變,2000年后呈現(xiàn)穩(wěn)定狀態(tài)。

(3)民和站歷年最大洪峰流量時頻域中存在4個周期，控制著歷年最大洪峰流量在整個時間域內(nèi)的變化特征，有4 a左右時間尺度的短振蕩周期，10、20 a左右時間尺度的中振蕩周期和超過32 a左右時間尺度的長振蕩周期，并且以長振蕩為主要周期分布在整個時域，重現(xiàn)期為20 a左右。

(4)在大尺度氣候變化下，人類活動、下墊面的變化是影響流域徑流、洪峰流量的主要因素。由于下墊面的變化和水利、水土保持工程措施的實施在很大程度上直接或間接的改變了流域的水資源狀況，導(dǎo)致湟水流域民和站實測洪峰流量和實測徑流呈現(xiàn)減少趨勢，大型水庫黑泉水庫的建成運行成為民和站實測洪峰流量突變的主要原因，估算人類活動因素對實測洪峰流量影響率為49.7%。

□