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

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

0 引 言

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

1 研究區(qū)域概況

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

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

2 數(shù)據(jù)來(lái)源

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

3 研究方法

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

3.1 趨勢(shì)檢驗(yàn)法

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

3.1.1 線性趨勢(shì)回歸檢驗(yàn)法

(1)

(2)

3.1.2 Kendall秩次檢驗(yàn)法

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

(3)

(4)

(5)

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

3.1.3 Spear man秩次檢驗(yàn)法：

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

(6)

(7)

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

3.2 突變點(diǎn)檢驗(yàn)法

3.2.1 滑動(dòng)T檢驗(yàn)法

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

(8)

(9)

3.2.2 有序聚類

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

(10)

(11)

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

Sn(τ)=Vτ+Vn-τ

(12)

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

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

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

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

(13)

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

(14)

極差R定義為：

(15)

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

(16)

則有：

(17)

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

對(duì)上式取對(duì)數(shù)：

(18)

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

3.2.4 M-K檢驗(yàn)法

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

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

(19)

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

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

(20)

(21)

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

(22)

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

3.3 周期性分析法

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

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

4.1 趨勢(shì)性

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

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

4.2 突 變

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

圖3 突變檢驗(yàn)圖Fig. 3 Test chart of mutation

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

4.3 周 期

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

從年最大洪峰流量小波實(shí)部等值線圖(圖5)可以看出民和站歷年最大洪峰流量在時(shí)序演化過(guò)程中存在著多時(shí)間尺度特征?？偟膩?lái)說(shuō)，存在著25～32 a，15～25 a、7～15 a以及3～7 a的4類尺度的周期變化規(guī)律。縱觀整個(gè)時(shí)域，在25～32 a尺度上出現(xiàn)了準(zhǔn)3次震蕩；在15～25 a時(shí)間尺度上有4次震蕩。同時(shí)，還可以看出上述兩個(gè)尺度的周期變化在整個(gè)時(shí)域內(nèi)相對(duì)穩(wěn)定，尤其25～32 a尺度的震蕩周期具有全時(shí)域性； 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時(shí)間尺度的能量強(qiáng)、周期顯著。

從湟水流域民和站歷年最大洪峰流量小波方差圖中(圖6)也可以看出，周期中心時(shí)間尺度32 a對(duì)應(yīng)著第一主峰也是第一主周期，其他峰值分別對(duì)應(yīng)著20、10和4 a周期中心時(shí)間尺度，與上述的小波實(shí)部等值線和能量譜的表現(xiàn)相對(duì)應(yīng)。這4個(gè)時(shí)間尺度的周期控制著年最大洪峰流量在整個(gè)時(shí)域內(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í)間尺度小波實(shí)數(shù)過(guò)程圖Fig.6 Wavelet variance and 32 years time scale wavelet real process

4.4 變化原因分析

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

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

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

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

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

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

4.4.1 水利工程

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

圖10 湟水流域水庫(kù)庫(kù)容累積曲線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è)實(shí)施，植被覆蓋增加，通過(guò)landsat8[9]衛(wèi)星2016年5月與1995年6月兩幅影像圖制作的NDVI對(duì)比圖顯示，湟水河谷地靠近城區(qū)兩側(cè)近山上的植被指數(shù)增加范圍明顯(圖12橘黃色部分)，說(shuō)明裸地減少，植被覆蓋態(tài)勢(shì)向好。

4.4.3 人類活動(dòng)影響估算

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

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

5 結(jié) 論

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

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

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

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

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

□