唐驕宇 盧 洋

(1.長安大學,陜西 西安 710000;2.黃河水利委員會三門峽庫區(qū)水文水資源局,河南 三門峽 472000)

渭河流域為干旱半干旱區(qū)域,在氣候變化和人類活動的影響下,渭河上游降水、徑流和泥沙均有不同程度的變化。渭河上游流域下墊面為黃土溝壑,水土流失嚴重,水資源供需矛盾突出,降水、徑流和泥沙的年際和季節(jié)變化規(guī)律對于研究地區(qū)水資源開發(fā)利用有重要影響。因此,研究區(qū)域內降水、徑流、泥沙變化規(guī)律,為渭河上游防汛抗旱、生態(tài)環(huán)境保護、水資源規(guī)劃和開發(fā)利用提供技術支撐,對于分析水循環(huán)規(guī)律、預測極端水文事件等具有重要意義。

本文以渭河上游天水水文勘測區(qū)(以下簡稱“天水測區(qū)”)為研究區(qū)域。天水測區(qū)上起渭河武山,下至天水市麥積區(qū)東岔鎮(zhèn),包括區(qū)間支流散渡河、葫蘆河、藉河、牛頭河流域,測區(qū)內有黃土高原、秦嶺、六盤山,年降水量地域分布極不均勻。采用的水文資料系列為1990—2020年,年降水量采用測區(qū)各雨量站實測數據,據此計算出流域平均面雨量;年徑流量、年輸沙量采用北道水文站實測值,利用線性趨勢法、肯德爾秩次檢驗法、啟發(fā)式分割、有序聚類法、季-海哈林法、雙累積曲線法等方法研究了降水、徑流和泥沙的變化特點、趨勢性、突變性,進而分析變化的影響原因,分析成果可為流域水資源合理開發(fā)利用與配置提供參考。

1 分析計算方法

1.1 趨勢性分析

分別采用線性回歸法、非參數檢驗法、肯德爾秩次檢驗法等進行趨勢性分析。

1.1.1 線性回歸法

線性回歸法構建水文變量x與其時間t的一元線性回歸方程:

x(t)=a+bt

(1)

回歸系數b的顯著性可通過t檢驗進行判斷,構建的T統(tǒng)計量:

(2)

服從自由度為(n-2)的t分布,給定顯著性水平α,可查得tα/2,如果T>Tα/2,則認為回歸效果顯著,存在線性趨勢;否則,線性趨勢不顯著。

1.1.2 非參數檢驗法

非參數檢驗法是對原始數據的秩進行分析,不需要樣本服從一定的分布,根據樣本信息對總體分布情況進行推斷和處理,并針對其不同條件下的參數進行優(yōu)化和整理。

1.1.3 肯德爾秩次檢驗法

利用肯德爾秩次檢驗法(Mann-Kendall法,簡稱M-K)[1]對水文序列的趨勢性采取顯著性檢驗方式,是一種非參數方法(亦稱無分布檢驗),其優(yōu)點是樣本不需要遵從一定的分布,也不受少數特異值的干擾,計算過程相對簡單,被廣泛應用于水文變量、要素等非正態(tài)分布的趨勢分析中,計算出相應的檢驗值,對水文序列的趨勢性和突變性進行有效的檢驗,其原理見表1。

表1 M-K趨勢判斷表(顯著性水平α=0.05)

時間序列X以(X1,X2,X3,…,Xn)表示,建立標準正態(tài)分布統(tǒng)計量Z:

(3)

其中

(4)

(5)

(6)

式中:當n>10時,s為近似服從正態(tài)分布;Var(s)為方差;m為數據相同的組數;tk為與第k組的數據相同的個數。

給定顯著性水平α,依據表1可以判斷序列的上升或下降趨勢顯著情況。

以Kendall傾斜度β量化序列,β可表示為

(7)

式中:l0時,呈現上升趨勢,反之呈下降趨勢。

1.2 跳躍成分分析

跳躍成分分析分別采用有序聚類法、季-海哈林法、啟發(fā)式分割法對降水、徑流進行突變點診斷,年輸沙量除采用前述三種方法外,還采用了雙累積曲線法進行突變點診斷分析。

1.2.1 有序聚類法[2]

有序聚類法利用最優(yōu)分割點來推求突變點,基本原則是使同要素之間的離差平方和較小。設有水文序列x1,x2,…,xn,假設可能的突變點為τ(2≤τ≤n-1),則突變點前后的離差平方和分別為

(8)

(9)

有序聚類法常用的目標函數為

(10)

當式(10)中S取極小值時,對應的τ為最優(yōu)分割點,可推斷為突變點。

1.2.2 季-海哈林檢驗法[2]

對于水文時間序列x(t),假定總體為正態(tài)分布,且可能變異點τ的先驗分布為均勻分布的情況下,推得τ的后驗分布為

(1≤τ≤n-1)

(11)

(12)

式中:k為比例常數;n為樣本容量。當統(tǒng)計值f(τ)取最大值時,滿足后驗分布條件,此時對應的τ值為最可能分割點,由此確定最可能變異年份。

1.2.3 啟發(fā)式分割法

相關研究[3-4]表明啟發(fā)式分割能在不同時間尺度高效地把非平穩(wěn)時間序列劃分成多個不同均值的子序列,并且相比多種常用突變檢測方法,能有效地排除虛假的變異點。對于時間序列X[X1,X2,…,Xn],其各點的合并偏差可采用下式計算:

(13)

式中:SD為合并偏差;N1、N2分別為該點左邊部分和右邊部分的點數;S1、S2分別為該點左邊部分和右邊部分的標準差。各點的統(tǒng)計量T計算公式為

(14)

式中:μ1、μ2分別為該點左邊部分和右邊部分的均值。

統(tǒng)計量T最大值Tmax的置信度P(Tmax)的計算式為

P(Tmax)≈{1-I[v/(v+Tmax2](δv,δ)}η

(15)

當P(Tmax)達到指定的置信度時,則該點為顯著的均值突變點,對原序列進行分割,之后對分割后的子序列繼續(xù)檢測。若P(Tmax)未能達到指定的置信度區(qū)間,則不再分割。

2 計算結果分析

2.1 降水

2.1.1 降水趨勢變化

天水測區(qū)平均年降水量年際變化見圖1,其多年平均年降水量為534.9mm,最大為755.4mm(2003年),最小為361.2mm(1997年)。年降水量趨勢線表明31年來降水總體呈平穩(wěn)微弱上升趨勢。線性回歸法、非參數檢驗法、肯德爾秩次檢驗法均一致診斷為趨勢不顯著,但計算檢驗值與標準值較接近(見表2),說明上升與非上升趨勢處于臨界狀態(tài)。

圖1 天水測區(qū)年降水量年際變化

表2 天水測區(qū)年降水量序列趨勢診斷結果

2.1.2 降水突變點診斷

對測區(qū)平均年降水量進行M-K突變檢驗,并設置0.05顯著性水平,即得到置信區(qū)間,臨界值為-1.96和1.96,分析結果見圖2。根據分析,年降水量UF、UB曲線在2010年后出現多個交點,相交于置信期間內,且兩曲線趨勢一致,通過UB、UF交叉點可知,突變點在2017年前后。

圖2 天水測區(qū)年降水量M-K秩次檢驗

采用有序聚類法和季-海哈林法進行同步驗證,得出相同檢驗值,|T|=2.42,T(0.05/2)=1.96,|T|>T(0.05/2),根據檢驗值趨勢圖,有序聚類法分析出其2017年資料系列達到最小值,季-海哈林法分析出其2017年資料系列達到最大值,根據突變點判別,在2017年前后均值發(fā)生顯著跳躍,見圖3、圖4。

圖3 天水測區(qū)年降水量有序聚類法檢驗

圖4 天水測區(qū)年降水量季-海哈林法檢驗

2.2 徑流

2.2.1 徑流趨勢變化

對北道水文站1990—2020年年徑流數據進行統(tǒng)計分析,得到徑流變化序列,見圖5。可知渭河上游多年平均年徑流量為7.32億m3,最大年徑流量為18.97億m3(2020年),最小年徑流量為1.29億m3(1997年)。由圖5可知,北道水文站年徑流量總體呈現升高趨勢,線性變化率為0.0952億m3/a,5年滑動平均顯示

圖5 天水測區(qū)年徑流量年際變化

近31年來年徑流量呈現緩慢上升趨勢。線性回歸法、非參數檢驗法、肯德爾秩次檢驗法均一致診斷為“趨勢不顯著”,但計算檢驗值與標準值較接近(見表3),說明上升與非上升趨勢處于臨界狀態(tài)。

表3 天水測區(qū)年徑流量序列趨勢診斷結果

2.2.2 徑流突變點診斷

北道水文站年徑流量M-K突變檢驗結果見圖6,設置0.05顯著性水平,即得到置信區(qū)間臨界值為-1.96和1.96,根據分析,得出2013年年徑流量UF<0,說明年徑流量處于下降趨勢,2013—2020年處于上升趨勢。20世紀90年代UF、UB值超過置信區(qū)間臨界值,表明在此期間年徑流量變化顯著。整體來看,采用有序聚類法和季-海哈林法同步驗證,得出檢驗值|T|=2.42>U(0.05/2)=1.96,根據檢驗值趨勢圖可分析出,在2017年資料系列達到最值,根據方法的突變點判別,均在2017年前后均值發(fā)生顯著跳躍,見圖7、圖8,與M-K檢驗結果一致。

圖6 天水測區(qū)年徑流量M-K秩次檢驗

圖7 天水測區(qū)年徑流量有序聚類法檢驗

圖8 天水測區(qū)年徑流量季-海哈林法檢驗

根據以上方法,設定置信度為0.95,最小分割尺度為25,根據年徑流變異分析結果,天水測區(qū)年徑流量啟發(fā)式分割檢驗見圖9。由圖9可知,2018年對應的T值最大,并且相對應的P(Tmax)為0.61,小于臨界值0.95,說明2018年為第一突變點[5-6]。

圖9 天水測區(qū)年徑流量啟發(fā)式分割檢驗

2.2.3 降水徑流EMD結果分析

運用EMD方法對天水測區(qū)多年來降水和徑流數據進行多尺度分解,結果見圖10、圖11,其中包含有多個具有物理意義的平穩(wěn)固有模態(tài)函數(Intrinsic Mode Functions,IMF)和具有單一性的趨勢項(Re-sidual,Res),從中可得出如下重要結論:

圖10 年降水量序列EMD分解

圖11 年徑流量序列EMD分解

a.根據EMD原理,渭河上游天水測區(qū)降水量,原始分量除外,可分解成3個不同波動周期的振蕩分量,徑流量可分解成4個不同波動周期的振蕩分量,反映了該測區(qū)降水和徑流在不同變化程度上的多時間尺度性。

b.20年代初及2015年前后降水和徑流IMF1分量變動幅度大,降水是徑流變化的主要影響因素,周期較為一致。IMF1是第一個本征模態(tài)函數,波動振幅最大,頻率最高,波長最短,周期2～5年,隨后的分解過程中,本征模態(tài)函數振幅逐漸減小,頻率逐漸降低,波長逐漸變大,其對原始時間序列的影響程度逐漸降低。

c.Res分量顯示的是天水測區(qū)年降水量和年徑流量的整體變化趨勢,近31年來整體呈急劇上升趨勢。EMD通過信號逐步去噪和趨勢剔除[7],實現了序列的平穩(wěn)化處理。

2.3 泥沙變化

2.3.1 輸沙量趨勢變化

對北道水文站1990—2020年年輸沙量數據進行統(tǒng)計分析,得到變化序列,見圖12。渭河上游多年平均年輸沙量為3490萬t,最大年輸沙量為15600萬t(1992年),最小年輸沙量為298萬t(2014年)。北道水文站年輸沙量總體呈現下降趨勢,線性變化率為198萬t/a,5年滑動平均顯示近31年來輸沙量呈現持續(xù)下降趨勢。采用線性回歸法、非參數檢驗、肯德爾秩次檢驗方法分析結果一致,均為下降趨勢顯著,見表4。

圖12 天水測區(qū)年輸沙量年際變化

表4 天水測區(qū)年輸沙量序列趨勢診斷結果

2.3.2 泥沙突變點診斷

對北道水文站年輸沙量進行M-K突變檢驗,見圖13。設置0.05顯著性水平,即得到置信區(qū)間臨界值為-1.96和1.96。根據分析得出,在1990—2020年,UF<0,說明年輸沙量呈顯著下降趨勢。20世紀90年代以來,UF和UB曲線存在較多交點,整體來看,相關檢驗值[|U|=3.651]<[U(0.05/2)=1.96],趨勢顯著。采用有序聚類法和季-海哈林法同步驗證,得出檢驗值[|T|=4.67]>[U(0.05/2)=1.96],根據檢驗值趨勢圖分析出,在1992年資料系列達到最值,根據方法的突變點判別,資料系列在1992年前后均值發(fā)生顯著跳躍,見圖14和圖15,與M-K檢驗結果一致。

圖13 天水測區(qū)年輸沙量M-K秩次檢驗

圖14 天水測區(qū)年輸沙量有序聚類法檢驗

圖15 天水測區(qū)年輸沙量季-海哈林法檢驗

2.3.3 徑流-泥沙關系分析

通過繪制徑流-泥沙雙累積曲線(見圖16),可以看出累積曲線顯著分為3個階段,1990—1992年為一個變化時段,該時段相關系數為0.9994;1993—2002年為一個變化時段,該時段相關系數為0.9935;2003—2020年為一個變化時段,該時段相關系數為0.9563;每個階段相關性均較好。由圖16可知,1992年、2002年前后均為輸沙量的突變點,累積輸沙量的趨勢線存在變緩趨勢,說明在徑流增加的情況下,輸沙量逐漸減少,這與流域內的生態(tài)景觀工程和近些年的水沙治理成效都有著密不可分的關系。

圖16 徑流-泥沙雙累積曲線

3 結 語

經過分析,天水測區(qū)屬干旱半干旱區(qū),徑流主要由降水補給,南北兩岸支流有不同水文特性,南岸面積小,地形陡峻,地表多為森林覆蓋或為剝蝕山地,降水較多,植被好,支流距渭河干流距離短,水多沙少;北岸支流水少沙多,水量較貧,主要為暴雨洪水,陡漲陡落。流域降水、徑流分別采用線性回歸、非參數檢驗、肯德爾秩次檢驗三種方法分析,結果一致,均為趨勢性不明顯,但統(tǒng)計值與判斷值較接近,說明處在臨界狀態(tài),而采用滑動平均曲線均有微弱上升趨勢。分別采用M-K秩次檢驗、啟發(fā)式分割、有序聚類法、季-海哈林突變診斷方法進行同步驗證,得出降水和徑流均在2017年前后產生突變一致性結論。對于年輸沙量分別采用線性回歸、非參數檢驗、肯德爾秩次檢驗三種方法分析,結果一致,均為下降趨勢明顯。采用啟發(fā)式分割、有序聚類、季-海哈林突變診斷方法得出,在1992年前后產生突變的一致結論,采用徑流-泥沙雙累積曲線法分析,明顯發(fā)現年輸沙量存在1992年和2002年兩個突變點。

天水測區(qū)在降水、徑流無明顯減少的情況下,年輸沙量減少趨勢明顯,說明流域內土地利用、生態(tài)保護、植被恢復等人類活動劇烈,特別是水土保持工作成效顯著。

采用徑流-泥沙雙累積曲線法發(fā)現年輸沙量存在兩個突變點,但基于數理統(tǒng)計分析判斷的啟發(fā)式分割、有序聚類、季-海哈林突變診斷方法年輸沙量診斷為一個,可見數理統(tǒng)計分析判斷方法對于多個突變點的診斷有待進一步完善。