黃河下游利津站入海水沙通量時序變化特征識別

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

(1.山東省水利科學研究院，山東 濟南 250013； 2.山東省水資源與水環(huán)境重點實驗室，山東 濟南 250013)

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

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

1 黃河入海水沙變化特征識別方法

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

1.1 MK方法

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

對于時間序列X={x1，x2，…，xn}，構造非參數MK統(tǒng)計量S為

(1)

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

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

(2)

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

(3)

假設原序列為無趨勢，當給定顯著水平α后，若|Z|

1.2 小波分析法

1.2.1 小波方差法

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

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

(4)

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

對于長度為n的離散時間序列，小波方差的計算公式為

(5)

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

1.2.2 小波重構與分解法

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

圖1 小波分解示意圖

本文采用Matlab中小波分析方法軟件包，應用CWT進行小波變換，用wavedec進行小波分解，wrcoef進行單支重構，提取低頻信號系數，通過Matlab應用函數的調用，可方便實現(xiàn)相應功能運算。

2 黃河下游水沙通量的趨勢變化識別

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

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

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

圖2 利津站1950—2009年徑流變化趨勢

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

表2 利津站1950—2009年MK和SP方法趨勢識別結果

3 黃河下游水沙通量的變點識別

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

表3 利津站1950—2009年水沙通量不同尺度變點識別結果

圖4為兩種方法識別年徑流量的變點過程圖，圖5為兩種方法識別年均輸沙通量的變點過程圖。

圖4 兩種方法識別年徑流量的變點

圖5 兩種方法識別年均輸沙通量變點

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

4 黃河下游水沙通量的周期識別

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

表4 利津站1950—2009年水沙通量不同尺度周期識別結果

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

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

5 結 論

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

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

[2] 侯志軍，由寶宏，茹玉英.黃河入?？谒惩孔兓?guī)律[J].泥沙研究，2007(5)：60-67.

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

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

[5] 仕玉治.氣候變化及人類活動對流域水資源的影響及實例研究[D].大連：大連理工大學，2011.

[6] 熊立華，周芬，肖義，等.水文時間序列變點分析的貝葉斯方法[J]. 水電能源科學，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] 桑燕芳，王棟.水文序列小波分析中小波函數選擇方法[J].水利學報，2008，39(3)：295-306.