唐驕宇 盧 洋

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

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

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

1 分析計(jì)算方法

1.1 趨勢性分析

分別采用線性回歸法、非參數(shù)檢驗(yàn)法、肯德爾秩次檢驗(yàn)法等進(jìn)行趨勢性分析。

1.1.1 線性回歸法

線性回歸法構(gòu)建水文變量x與其時(shí)間t的一元線性回歸方程:

x(t)=a+bt

(1)

回歸系數(shù)b的顯著性可通過t檢驗(yàn)進(jìn)行判斷,構(gòu)建的T統(tǒng)計(jì)量:

(2)

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

1.1.2 非參數(shù)檢驗(yàn)法

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

1.1.3 肯德爾秩次檢驗(yàn)法

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

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

時(shí)間序列X以(X1,X2,X3,…,Xn)表示,建立標(biāo)準(zhǔn)正態(tài)分布統(tǒng)計(jì)量Z:

(3)

其中

(4)

(5)

(6)

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

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

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

(7)

式中:l0時(shí),呈現(xiàn)上升趨勢,反之呈下降趨勢。

1.2 跳躍成分分析

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

1.2.1 有序聚類法[2]

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

(8)

(9)

有序聚類法常用的目標(biāo)函數(shù)為

(10)

當(dāng)式(10)中S取極小值時(shí),對(duì)應(yīng)的τ為最優(yōu)分割點(diǎn),可推斷為突變點(diǎn)。

1.2.2 季-海哈林檢驗(yàn)法[2]

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

(1≤τ≤n-1)

(11)

(12)

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

1.2.3 啟發(fā)式分割法

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

(13)

式中:SD為合并偏差;N1、N2分別為該點(diǎn)左邊部分和右邊部分的點(diǎn)數(shù);S1、S2分別為該點(diǎn)左邊部分和右邊部分的標(biāo)準(zhǔn)差。各點(diǎn)的統(tǒng)計(jì)量T計(jì)算公式為

(14)

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

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

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

(15)

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

2 計(jì)算結(jié)果分析

2.1 降水

2.1.1 降水趨勢變化

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

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

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

2.1.2 降水突變點(diǎn)診斷

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

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

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

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

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

2.2 徑流

2.2.1 徑流趨勢變化

對(duì)北道水文站1990—2020年年徑流數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析,得到徑流變化序列,見圖5?？芍己由嫌味嗄昶骄陱搅髁繛?.32億m3,最大年徑流量為18.97億m3(2020年),最小年徑流量為1.29億m3(1997年)。由圖5可知,北道水文站年徑流量總體呈現(xiàn)升高趨勢,線性變化率為0.0952億m3/a,5年滑動(dòng)平均顯示

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

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

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

2.2.2 徑流突變點(diǎn)診斷

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

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

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

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

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

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

2.2.3 降水徑流EMD結(jié)果分析

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

圖10 年降水量序列EMD分解

圖11 年徑流量序列EMD分解

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

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

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

2.3 泥沙變化

2.3.1 輸沙量趨勢變化

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

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

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

2.3.2 泥沙突變點(diǎn)診斷

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

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

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

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

2.3.3 徑流-泥沙關(guān)系分析

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

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

3 結(jié) 語

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

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

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