范明元，李福林，仕玉治，陳華偉

(1.山東省水利科學(xué)研究院，山東 濟(jì)南 250013； 2.山東省水資源與水環(huán)境重點(diǎn)實(shí)驗(yàn)室，山東 濟(jì)南 250013)

黃河是中國(guó)第二大河流，也是世界上最著名的水少沙多河流。受到氣候變化和人類活動(dòng)雙重影響下，河道徑流、泥沙等水文要素儼然發(fā)生變化，尤其近幾十年入海水沙通量呈現(xiàn)銳減趨勢(shì)，甚至發(fā)生斷流，對(duì)黃河入?？谏鷳B(tài)環(huán)境產(chǎn)生了巨大影響，由此引起國(guó)內(nèi)許多學(xué)者的廣泛關(guān)注。對(duì)黃河入海口水沙通量要素的變化特征識(shí)別，可為認(rèn)知水沙動(dòng)態(tài)、變化驅(qū)動(dòng)力以及演變規(guī)律等方面提供數(shù)據(jù)支持，同時(shí)，對(duì)實(shí)現(xiàn)黃河入?？谏鷳B(tài)環(huán)境可持續(xù)發(fā)展具有重要的意義。

對(duì)于黃河入海口水沙通量的變化特征分析，已有許多成果，比較已有研究成果發(fā)現(xiàn)：①所應(yīng)用的研究時(shí)間序列長(zhǎng)度不同[1-3]；②特征分析的側(cè)重點(diǎn)不同且不全面，僅分析水文要素的變化趨勢(shì)、變點(diǎn)或者周期中的一個(gè)方面或者兩方面，如侯志軍等[2]根據(jù)1950—2003年利津站水沙月序列資料，運(yùn)用統(tǒng)計(jì)分析方法研究了黃河入?？趨^(qū)水沙通量的季節(jié)性變化、年際變化以及水沙通量之間的關(guān)系；樊輝等[3]基于1950—2007年利津站實(shí)測(cè)水沙年序列資料，采用Mann-Kendall方法分析黃河入海水沙通量的變化規(guī)律；③部分研究?jī)H對(duì)徑流或者輸沙量單一水文要素研究，且研究的尺度不同，如丁艷峰等[1]以黃河利津站1950—2004年的實(shí)測(cè)、天然月徑流序列為基礎(chǔ)資料，采用回歸分析、相關(guān)分析和隨機(jī)水文學(xué)的方法，僅對(duì)黃河入海徑流量的變化特征進(jìn)行了分析。此外，2002年以來(lái)，小浪底水庫(kù)先后實(shí)施了9次調(diào)水、調(diào)沙，這為徑流、輸沙量時(shí)間序列的變化規(guī)律識(shí)別，提供新的數(shù)據(jù)資料和人類活動(dòng)響應(yīng)。因此，筆者以1950—2009年水沙月時(shí)間序列為研究對(duì)象，采用應(yīng)用較為廣泛的小波分析方法及Mann-Kendall等多種方法，對(duì)時(shí)間序列的變化趨勢(shì)、變點(diǎn)及周期等多方面，綜合識(shí)別水沙通量的變化特征。

1 黃河入海水沙變化特征識(shí)別方法

國(guó)內(nèi)外對(duì)水文序列趨勢(shì)分析方面的研究很多，目前進(jìn)行水文序列趨勢(shì)分析的方法主要有Mann-Kendall秩次相關(guān)檢驗(yàn)法(簡(jiǎn)稱 MK)、Spearman秩次相關(guān)檢驗(yàn)法(簡(jiǎn)稱SP)、線性趨勢(shì)回歸檢驗(yàn)法等[4-7]。通過(guò)模型的比較分析，可以定性識(shí)別水文變量的趨勢(shì)變化。本文采用較為常用且經(jīng)檢驗(yàn)較為有效的MK方法、SP方法、有序聚類方法和小波分析方法，其中采用小波分解與重構(gòu)技術(shù)方法識(shí)別時(shí)間序列的變點(diǎn)，通過(guò)小波方差法識(shí)別時(shí)間序列周期。限于篇幅，僅給出MK方法和小波分析方法計(jì)算原理，其他方法可參考文獻(xiàn)[5]。MK方法是用來(lái)識(shí)別序列變化趨勢(shì)應(yīng)用廣泛的方法，而小波分析法是時(shí)、頻分析的有效方法，可通過(guò)小波分解與重構(gòu)進(jìn)行序列變點(diǎn)分析，也通過(guò)小波變換進(jìn)行序列周期變化特征識(shí)別。

1.1 MK方法

與參數(shù)統(tǒng)計(jì)檢驗(yàn)法相比，非參數(shù)MK檢驗(yàn)法更適用于非正態(tài)分布的資料。MK檢驗(yàn)是提取序列變化趨勢(shì)的有效工具，被廣泛應(yīng)用于氣候參數(shù)和水文序列的分析。

對(duì)于時(shí)間序列X={x1，x2，…，xn}，構(gòu)造非參數(shù)MK統(tǒng)計(jì)量S為

(1)

當(dāng)n≥10時(shí)，統(tǒng)計(jì)量S近似服從正態(tài)分布，其正態(tài)分布均值E(S)和方差Var(S)分別為

E(S)=0Var(S)=n(n+1)(2n+5)/18

(2)

正態(tài)分布檢驗(yàn)統(tǒng)計(jì)量Z為

(3)

假設(shè)原序列為無(wú)趨勢(shì)，當(dāng)給定顯著水平α后，若|Z|

1.2 小波分析法

1.2.1 小波方差法

小波分析(wavelet analysis，簡(jiǎn)稱WA)是一種時(shí)、頻多分辨率分析方法，即時(shí)間域和頻率域均可改變的分析方法，是傅立葉分析思想的延伸與拓展，其核心是小波變換[8]。

令L2(R)表示定義在實(shí)軸上、可測(cè)的平方可積函數(shù)空間，則對(duì)于時(shí)間序列f(t)∈L2(R)，其連續(xù)小波變換(continues wavelet transform，簡(jiǎn)稱CWT)公式為

(4)

通過(guò)連續(xù)小波變換，可獲得不同尺度下的小波變換系數(shù)Wf(a，b)。Wf(a，b)能夠同時(shí)反映時(shí)、域參數(shù)b和頻域參數(shù)a的特性。通過(guò)小波系數(shù)分布圖可以清楚地看出時(shí)間序列的周期變化情況，若是要準(zhǔn)確的判斷序列的周期，需利用小波方差進(jìn)行判斷。

對(duì)于長(zhǎng)度為n的離散時(shí)間序列，小波方差的計(jì)算公式為

(5)

式中：V(a)為尺度伸縮因子a、時(shí)間xj處的小波系數(shù)平方W2(a，xj)的均值，對(duì)于復(fù)小波函數(shù)W2(a，xj)則為系數(shù)模的平方，小波方差的各個(gè)峰值分別對(duì)應(yīng)顯著周期。當(dāng)小波系數(shù)達(dá)到最大值時(shí)，小波函數(shù)的尺度與序列周期吻合最好。

1.2.2 小波重構(gòu)與分解法

通常時(shí)間序列是趨勢(shì)項(xiàng)、周期項(xiàng)和隨機(jī)項(xiàng)三者的線性疊加，其中趨勢(shì)項(xiàng)對(duì)應(yīng)于小波分解后最大尺度的低頻重構(gòu)序列，而隨機(jī)項(xiàng)對(duì)應(yīng)高頻部分，除去趨勢(shì)項(xiàng)與隨機(jī)項(xiàng)的影響后所得到的就是時(shí)間序列的周期部分，而通過(guò)小波分解與重構(gòu)，隨機(jī)成分就會(huì)被分離出去，分解后的低頻系數(shù)重構(gòu)序列即可代表時(shí)間序列的變化趨勢(shì)，通過(guò)重構(gòu)序列折點(diǎn)可識(shí)別時(shí)間序列的變點(diǎn)。通常根據(jù)樣本容量，最多可以把序列分解成log2N個(gè)頻率級(jí)，一般情況下，N>2，對(duì)黃河入海水沙時(shí)間序列x，采用db3小波函數(shù)進(jìn)行分辨率為5的快速小波分解，可分別獲得近似系數(shù)CA1和詳細(xì)系數(shù)[CD1，CD2，…，CD5]，然后對(duì)其低頻系數(shù)進(jìn)行單支重構(gòu)。小波分解的結(jié)構(gòu)示意圖如圖1所示。

圖1 小波分解示意圖

本文采用Matlab中小波分析方法軟件包，應(yīng)用CWT進(jìn)行小波變換，用wavedec進(jìn)行小波分解，wrcoef進(jìn)行單支重構(gòu)，提取低頻信號(hào)系數(shù)，通過(guò)Matlab應(yīng)用函數(shù)的調(diào)用，可方便實(shí)現(xiàn)相應(yīng)功能運(yùn)算。

2 黃河下游水沙通量的趨勢(shì)變化識(shí)別

近些年來(lái)，受到氣候變化和人類活動(dòng)的影響，黃河入海口水沙通量發(fā)生了較大變化，據(jù)1950—2009年水沙通量序列數(shù)據(jù)全年、汛期、非汛期的不同年代統(tǒng)計(jì)分析(表1)，自20世紀(jì)60年代以來(lái)，不同年代均值呈現(xiàn)減少變化，尤其80年代以后，其水沙通量年代均值均小于多年平均值，說(shuō)明黃河入?？谒惩砍尸F(xiàn)減少變化趨勢(shì)。但2002年以來(lái)，隨著小浪底水庫(kù)調(diào)水調(diào)沙的運(yùn)行，黃河口徑流量有所回升，輸沙量變化較小。年均值、汛期及非汛期水沙通量時(shí)間序列過(guò)程如圖2、3所示。

表1 1950—2009年利津站水沙通量統(tǒng)計(jì)

為進(jìn)一步檢驗(yàn)黃河入海口水沙通量變化趨勢(shì)，采用MK和SP方法兩種方法對(duì)利津站1950—2009年全年、汛期、非汛期水沙通量序列趨勢(shì)進(jìn)行檢驗(yàn)，檢驗(yàn)結(jié)果見(jiàn)表2。由表2可知，全年、汛期、非汛期水沙通量均通過(guò)兩種方法99.5%置信檢驗(yàn)，呈現(xiàn)顯著下降變化趨勢(shì)。

圖2 利津站1950—2009年徑流變化趨勢(shì)

圖3 利津站1950—2009年輸沙量變化

表2 利津站1950—2009年MK和SP方法趨勢(shì)識(shí)別結(jié)果

3 黃河下游水沙通量的變點(diǎn)識(shí)別

采用小波分解與重構(gòu)方法和有序聚類方法對(duì)利津站1950—2009年不同尺度水沙通量進(jìn)行變點(diǎn)識(shí)別，結(jié)果見(jiàn)表3。由表3知：①對(duì)于不同尺度的徑流量，小波分解與重構(gòu)方法識(shí)別結(jié)果顯示不同尺度下較為明顯的變點(diǎn)為1985年和2002年，而有序聚類法識(shí)別結(jié)果顯示全年和汛期一致，最有可能的變點(diǎn)為1985年，非汛期有所不同，為1971年；②對(duì)于不同尺度的輸沙量，小波分解與重構(gòu)方法識(shí)別結(jié)果顯示年均值和汛期的變點(diǎn)為1969、1985、2002年，非汛期為1985，2002年，而有序聚類方法識(shí)別結(jié)果顯示年均值和汛期最有可能變點(diǎn)為1985年。

表3 利津站1950—2009年水沙通量不同尺度變點(diǎn)識(shí)別結(jié)果

圖4為兩種方法識(shí)別年徑流量的變點(diǎn)過(guò)程圖，圖5為兩種方法識(shí)別年均輸沙通量的變點(diǎn)過(guò)程圖。

圖4 兩種方法識(shí)別年徑流量的變點(diǎn)

圖5 兩種方法識(shí)別年均輸沙通量變點(diǎn)

對(duì)比兩種方法識(shí)別結(jié)果，以20世紀(jì)80年代為分界線，前后水沙通量變化較為明顯，1985年為水沙通量的主要變化轉(zhuǎn)折點(diǎn)，其次為2002年。原因是龍羊峽水庫(kù)1986—1989年連續(xù)4年蓄水導(dǎo)致黃河入海水沙通量在1985年出現(xiàn)明顯轉(zhuǎn)折，2002年為小浪底水庫(kù)初運(yùn)行起始點(diǎn)，2002年以來(lái)利用小浪底水庫(kù)調(diào)蓄功能實(shí)行9次調(diào)水調(diào)沙后，使全年、非汛期入海徑流量有所抬升，呈現(xiàn)較為明顯徑流變點(diǎn)，此外，水庫(kù)汛期調(diào)節(jié)洪峰時(shí)可致使河流輸沙能力降低，大量泥沙滯留導(dǎo)致入海泥沙通量減少，出現(xiàn)變點(diǎn)。

4 黃河下游水沙通量的周期識(shí)別

采用小波分析法對(duì)利津站1950—2009年不同尺度的水沙通量進(jìn)行周期識(shí)別，結(jié)果見(jiàn)表4。對(duì)兩種不同水文要素，其年均、汛期第一、第二主周期基本相同；對(duì)于不同時(shí)間尺度的，其非汛期水沙通量的第一、第二主周期不同于全年和汛期的周期變化特征。圖6為年均徑流量的小波周期分析圖和小波方差過(guò)程圖；圖7為年均輸沙量的小波周期分析圖和小波方差過(guò)程圖。

表4 利津站1950—2009年水沙通量不同尺度周期識(shí)別結(jié)果

圖6 年均徑流量的周期分析和小波方差過(guò)程

圖7 年均輸沙量的小波周期分析和小波方差過(guò)程

5 結(jié) 論

通過(guò)多種方法綜合分析1950—2009年黃河入海水沙通量的不同尺度的變化趨勢(shì)、變點(diǎn)及周期等三方面的變化特征，并對(duì)黃河入海通量水沙通量變化特征形成原因進(jìn)行了分析，不同方法識(shí)別結(jié)果顯示：①采用MK，SP兩種方法識(shí)別的黃河入海水沙通量的不同時(shí)間尺度的變化趨勢(shì)均呈現(xiàn)顯著下降；②采用小波重構(gòu)與分解法和有序聚類方法識(shí)別水沙通量的主要變點(diǎn)為1985年，次變點(diǎn)為2002年；不同時(shí)間尺度變點(diǎn)識(shí)別結(jié)果有所不同，但差別不大，兩種方識(shí)別結(jié)果有所不同，同樣差別不大；③采用小波分析法(小波方差法)識(shí)別水沙通量的不同時(shí)間尺度的變化周期，年均、汛期徑流量的第一主周期為8a，第二主周期為22a；非汛期徑流周期有所不同，第一主周期為6a，第二主周期為17a；年均、汛期輸沙量的第一主周期為3a或者9a，第二主周期為24a或者25a，非汛期輸沙量第一主周期為8a，第二主周期為18a。

[1] 丁艷峰，潘少明.近50年黃河入海徑流變化特征及影響因素分析[J].第四紀(jì)研究，2007，27(5)：709-716.

[2] 侯志軍，由寶宏，茹玉英.黃河入海口水沙通量變化規(guī)律[J].泥沙研究，2007(5)：60-67.

[3] 樊輝，劉艷霞，黃海軍.1950-2007年黃河入海水沙通量變化趨勢(shì)及突變特征[J].泥沙研究，2009(5)：9-16.

[4] MANN H B. Nonparametric tests against trend[J]. Econometrica，1945，13：245-259.

[5] 仕玉治.氣候變化及人類活動(dòng)對(duì)流域水資源的影響及實(shí)例研究[D].大連：大連理工大學(xué)，2011.

[6] 熊立華，周芬，肖義，等.水文時(shí)間序列變點(diǎn)分析的貝葉斯方法[J]. 水電能源科學(xué)，2003，21(4)：39-41.

[7] PETTITT A N. A non-parametric approach to the change point problem[J]. Applied Statistics，1979，28(2)：126-135.

[8] 桑燕芳，王棟.水文序列小波分析中小波函數(shù)選擇方法[J].水利學(xué)報(bào)，2008，39(3)：295-306.