渭河上游水沙變化特征分析

唐驕宇 盧 洋

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

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

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

1 分析計算方法

1.1 趨勢性分析

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

1.1.1 線性回歸法

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

x(t)=a+bt

(1)

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

(2)

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

1.1.2 非參數檢驗法

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

1.1.3 肯德爾秩次檢驗法

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

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

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

(3)

其中

(4)

(5)

(6)

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

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

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

(7)

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

1.2 跳躍成分分析

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

1.2.1 有序聚類法[2]

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

(8)

(9)

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

(10)

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

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

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

(1≤τ≤n-1)

(11)

(12)

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

1.2.3 啟發式分割法

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

(13)

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

(14)

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

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

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

(15)

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

2 計算結果分析

2.1 降水

2.1.1 降水趨勢變化

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

圖1 天水測區年降水量年際變化

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

2.1.2 降水突變點診斷

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

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

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

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

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

2.2 徑流

2.2.1 徑流趨勢變化

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

圖5 天水測區年徑流量年際變化

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

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

2.2.2 徑流突變點診斷

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

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

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

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

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

圖9 天水測區年徑流量啟發式分割檢驗

2.2.3 降水徑流EMD結果分析

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

圖10 年降水量序列EMD分解

圖11 年徑流量序列EMD分解

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

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

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

2.3 泥沙變化

2.3.1 輸沙量趨勢變化

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

圖12 天水測區年輸沙量年際變化

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

2.3.2 泥沙突變點診斷

對北道水文站年輸沙量進行M-K突變檢驗,見圖13。設置0.05顯著性水平,即得到置信區間臨界值為-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年前后均值發生顯著跳躍,見圖14和圖15,與M-K檢驗結果一致。

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

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

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

2.3.3 徑流-泥沙關系分析

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

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

3 結 語

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

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

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