(新疆昌源水務(wù)科學(xué)研究院有限公司，新疆 烏魯木齊 830000)

20世紀以來，隨著人口的增長和社會經(jīng)濟的發(fā)展，人類對水資源的需求不斷增加，許多地區(qū)不僅水資源日益短缺，而且還面臨一系列的環(huán)境問題。天然地表徑流作為水資源的重要組成部分及水循環(huán)過程中的關(guān)鍵環(huán)節(jié)，其變化蘊含著突變性、周期性、混沌特性等非線性變化特征[1]。詳細刻畫和揭示天然徑流的內(nèi)在規(guī)律是實現(xiàn)水資源可持續(xù)利用的前提，特別在干旱半干旱區(qū)內(nèi)流河流域，由于平原區(qū)降水稀少，水資源主要依靠源流山區(qū)由降水和冰川融水所形成的徑流補給，對氣候波動的響應(yīng)更為敏感[2]。IPCC根據(jù)氣候模式預(yù)測結(jié)果指出，到2100年，全球平均氣溫將上升1.1～6.4℃[3]，氣候變暖將加速水循環(huán)的過程，導(dǎo)致極端水文事件發(fā)生頻率增加和水資源的時空再分配[4]。因此，研究氣候變化背景下流域徑流的時空變化已成為水資源管理者和科學(xué)家們關(guān)注的重點。

塔里木河流域是新疆最重要的糧食、棉花和石油生產(chǎn)基地，具有自然資源豐富與生態(tài)環(huán)境脆弱并存的特點。塔里木河干流下游，受長期斷流和生態(tài)輸水的影響，已成為我國乃至世界范圍內(nèi)人類干擾受損極為嚴重的稀有區(qū)和人為促進生態(tài)修復(fù)極具成效的典型區(qū)，具有重要的科學(xué)研究和實踐借鑒價值。而保證塔里木河干流下游實現(xiàn)長期的生態(tài)輸水，則需要對源流天然地表徑流進行深入研究。

本文采用基于樣本熵的小波閾值去噪方法對集合經(jīng)驗?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)進行改進，分析了塔里木河“三源一干”(阿克蘇河、葉爾羌河、和田河與塔里木河干流)1960—2015年地表徑流的變化特征，研究成果可為流域水資源管理提供參考依據(jù)。

1 研究區(qū)概況

塔里木河流域總面積102萬km2，是中國最大的內(nèi)陸河流域，由于地處全封閉的內(nèi)陸區(qū)域，流域內(nèi)降水稀少，蒸發(fā)強烈。塔里木河流域盆地南部、西部和北部為阿爾金山、昆侖山和天山，山勢高峻，是河川徑流的產(chǎn)水區(qū)。高山之上有冰川和永久積雪，起到高山水庫作用；高原山坡攔截濕潤氣流，形成降水徑流，其中冰雪融水補給占河流補給量的45%～60%，降水占河流補給量的18%～33%。

塔里木河流域?qū)儆诘湫偷膬?nèi)陸耗散型河流，水資源形成于山區(qū)，消耗于平原區(qū)、荒漠區(qū)，消失于沙漠。自阿克蘇河、葉爾羌河、和田河三河匯合口的肖夾克以下河流稱之為塔里木河干流。塔里木河干流自身不產(chǎn)流，徑流全靠源流補給，干流全長1321km，尾閭位于塔里木盆地東部的臺特瑪湖。歷史上，塔里木河流域內(nèi)九大水系均有水匯入塔里木河干流。但受人類活動與氣候變化等的影響，目前僅有阿克蘇河、葉爾羌河、和田河補給塔里木河，稱為“三源流”。其中：阿克蘇河常年有水補給塔里木河，和田河在每年洪水期有水注入塔里木河，葉爾羌河在特大洪水期才有水進入塔里木河。本文研究區(qū)域為塔里木河三源流和干流(簡稱塔里木河“三源一干”)，見圖1。

圖1 塔里木河流域“三源一干”水系

2 研究方法

2.1 集合經(jīng)驗?zāi)B(tài)分解方法

集合經(jīng)驗?zāi)B(tài)分解(EEMD)方法利用高斯白噪聲具有頻率均勻分布的統(tǒng)計特性，自適應(yīng)地將數(shù)據(jù)分解為有限個本征模函數(shù)(Intrinsic Mode Function，IMF)，所分解出來的各IMF分量包含了原始數(shù)據(jù)的不同時間尺度的局部特征信號[5-6]。因此，本文采用EEMD方法對塔里木河“三源一干”徑流數(shù)據(jù)進行多尺度分解，分解結(jié)果可以由式(1)求得，重構(gòu)結(jié)果由式(2)求得。

(1)

(2)

式中Cij(t)——第i次加入白噪聲后分解所得的第j個IMF分量;

t——時間；

Cj(t)——分解得到的第j個IMF分量；

N——添加白噪聲序列的數(shù)目；

X(t)——原始信號；

rn——分解后的殘留余項；

n——分解次數(shù)。

2.2 改進的集合經(jīng)驗?zāi)B(tài)分解方法

天然水文過程和動態(tài)水文系統(tǒng)易受到與之相關(guān)的多種物理因素(如氣候、流域、地理特征等)影響，為分析這種非穩(wěn)定的水文時間序列變化規(guī)律，通常將其分為確定性成分和隨機成分，其中確定性成分包括周期性過程和暫態(tài)過程(趨勢、奇異點等)，而隨機成分就是通常所說的噪聲，是受到許多不確定和隨機的因素的影響而產(chǎn)生的。目前，在水文研究領(lǐng)域，應(yīng)用較為廣泛的去噪方法為基于樣本熵的小波閾值去噪方法，該方法計算流程詳見文獻[7]。

本文首先采用基于樣本熵的小波閾值方法對塔里木河“三源一干”徑流時間序列數(shù)據(jù)去噪，再對去噪后的序列數(shù)據(jù)做鏡像延伸處理，以便抑制直接采用EEMD方法時序列兩端產(chǎn)生的端點效應(yīng)，對延伸后的序列數(shù)據(jù)進行距平處理之后再利用EEMD方法得到各序列數(shù)據(jù)的固態(tài)模態(tài)函數(shù)以及趨勢項，并使用最大熵譜分析法對各固態(tài)模態(tài)函數(shù)進行周期分析。

3 結(jié)果與分析

3.1 基于樣本熵小波閾值方法的“三源一干”徑流數(shù)據(jù)去噪分析

分別對阿克蘇河、葉爾羌河、和田河與塔里木河干流1960—2015年徑流序列進行小波閾值去噪分析，結(jié)果如下：

a.阿克蘇河去噪分析。結(jié)果見圖2(a)，當閾值為3.1時，噪聲序列的樣本熵值達到最大，然后隨著閾值的再增加，樣本熵值開始減小，表明此刻的閾值為最適當去噪閾值。該閾值下去噪序列與原序列的變化見圖2(b)。

圖2 阿克蘇河年徑流序列去噪

對去噪結(jié)果進行評價(見表1)，阿克蘇河協(xié)合拉站多年平均徑流為49.52億m3，經(jīng)過小波閾值去噪后，去噪序列均值變?yōu)?9.54億m3，與原徑流序列均值相差-0.02億m3；原徑流序列的方差為57.88，經(jīng)過去噪后的徑流序列方差變?yōu)?9.80，小于原徑流序列，表明去除噪聲干擾后，原年徑流序列的復(fù)雜度降低；原年徑流序列和去噪后的年徑流序列的偏態(tài)系數(shù)分別為1.03和1.09，表明經(jīng)去噪后去噪序列保留了原年徑流序列的偏態(tài)特征；去噪序列的一階自相關(guān)系數(shù)相對于原年徑流序列增加了0.20，而分離出的噪聲序列的一階自相關(guān)系數(shù)為-0.06，接近于0，表明經(jīng)過去噪后，年徑流序列的自我相關(guān)性有所增加，而基于噪聲隨機性的特點(無自我相關(guān)性)，因此接近于0。經(jīng)過小波閾值去噪后，年徑流序列和分離出的噪聲的各項指標均符合要求，表明去噪結(jié)果合理，滿足要求。

表1 阿克蘇河徑流序列去噪結(jié)果評價

b.葉爾羌河去噪分析。結(jié)果見圖3(a)，當閾值為1.8時，噪聲序列樣本熵值達到最大，然后隨著閾值的增加，樣本熵值開始減小，表明此刻的閾值為最適去噪閾值。該閾值下去噪序列與原序列的變化見圖3(b)。

對去噪結(jié)果進行評價(見表2)，葉爾羌河卡群站多年平均徑流量為65.790億m3，經(jīng)小波閾值去噪后，去噪序列均值變?yōu)?5.850億m3，與原徑流序列相差0.033億m3，與分離出來的噪聲序列均值相等；原年徑流序列的方差為129.06，經(jīng)去噪后的年徑流序列均值為96.51，小于原年徑流序列，表明去除噪聲干擾后，原年徑流量序列的復(fù)雜度降低；原年徑流量序列和去噪后的年徑流量序列的偏態(tài)系數(shù)分別為0.165和0.126，表明經(jīng)去噪后去噪序列保留了原年徑流序列的偏態(tài)特征；去噪序列的一階自相關(guān)系數(shù)相對于原年徑流量序列的一階自相關(guān)系數(shù)分別為-0.157和-0.218，有所下降，而分離出的噪聲序列的一階自相關(guān)系數(shù)為-0.066，接近于0，表明經(jīng)去噪后，年徑流序列的自我相關(guān)性有所增加，而基于噪聲完全隨機的特點，無自我相關(guān)性，因此接近于0。經(jīng)過小波閾值去噪后，去噪后的年徑流序列和分離出的噪聲各項指標均符合要求，表明去噪結(jié)果合理，滿足要求。

圖3 葉爾羌河年徑流序列去噪

表2 葉爾羌河徑流序列去噪結(jié)果評價

c.和田河去噪分析。結(jié)果見圖4(a)，當閾值為30時，噪聲序列的樣本熵值達到最大，然后隨著閾值的增加，樣本熵值開始減小，表明此刻的閾值為最適去噪閾值，該閾值下去噪序列與原序列的變化見圖 4(b)。

圖4 和田河年徑流序列去噪

對去噪結(jié)果進行評價(見表3)，和田河多年平均徑流為22.26億m3，經(jīng)過小波閾值去噪后，去噪序列均值變?yōu)?2.29億m3，與原徑流序列的均值相差-0.03億m3，與分離出來的噪聲序列均值相等；原年徑流序列方差為26.36，經(jīng)去噪后的年徑流序列方差為13.78，小于原徑流序列，表明去除噪聲干擾后，原徑流序列的復(fù)雜度降低；原徑流序列和去噪后的徑流序列的偏態(tài)系數(shù)分別為0.61和0.73，經(jīng)去噪后偏態(tài)系數(shù)與原徑流偏態(tài)系數(shù)相近，表明去噪序列保留了原徑流序列的偏態(tài)特征；去噪序列的一階自相關(guān)系數(shù)與原徑流序列的一階自相關(guān)系數(shù)分別為-0.10和0.10，有所下降，而分離出的噪聲序列一階自相關(guān)系數(shù)為-0.30，接近于0，表明經(jīng)去噪后，年徑流序列的自我相關(guān)性有所增加，而基于噪聲完全隨機的特點(無自我相關(guān)性)，因此接近于0。經(jīng)過小波閾值去噪后，去噪后的年徑流序列和分離出的噪聲各項指標均符合要求，表明去噪結(jié)果合理，滿足要求。

d.塔里木河干流去噪分析。結(jié)果見圖5(a)，當閾值為36時，噪聲序列的樣本熵值達到最大，之后隨著閾值的再增加，樣本熵值開始減小，表明該閾值為最適去噪閾值。該閾值下去噪序列與原序列的變化見圖5(b)。

表3 和田河徑流序列去噪結(jié)果評價

對去噪結(jié)果進行評價(見表4)，塔里木河干流阿拉爾站多年平均徑流為45.40億m3，經(jīng)過小波閾值去噪后，去噪序列均值變?yōu)?5.20億m3，與原徑流序列均值相差0.20億m3，與分離出來的噪聲序列的均值相等；原徑流序列方差為143.90，而經(jīng)去噪后的年徑流序列方差為16.93，小于原徑流序列，表明去除噪聲干擾后，原徑流序列復(fù)雜度降低；原徑流序列和去噪后的徑流序列的偏態(tài)系數(shù)分別-0.10和0.10，經(jīng)去噪后偏態(tài)系數(shù)與原徑流量偏態(tài)系數(shù)相近，表明去噪序列保留了原徑流序列的偏態(tài)特征；去噪序列的一階自相關(guān)系數(shù)相對于原徑流序列的一階自相關(guān)系數(shù)分別為-0.09和0.10，有所上升，而分離出的噪聲序列的一階自相關(guān)系數(shù)為-0.15，接近于0，表明經(jīng)去噪后，徑流序列的自我相關(guān)性有所增加。經(jīng)過小波閾值去噪后，去噪后的徑流序列和分離出的噪聲各項指標均符合要求，表明去噪結(jié)果合理，滿足要求。

圖5 塔里木河干流徑流去噪

特征值年均值/億m3方差偏態(tài)系數(shù)一階自相關(guān)系數(shù)原序列45.40143.90-0.10-0.09去噪系列45.2016.930.100.10噪聲系列0.19122.610.40-0.15

3.2 基于EEMD的“三源一干”徑流周期性分析

對去噪后的“三源一干”年均徑流序列分別作鏡像延伸處理，以便抑制當使用經(jīng)驗?zāi)B(tài)分解時序列兩端產(chǎn)生的端點效應(yīng)，然后對序列進行距平處理，再對經(jīng)距平處理后的序列用集合經(jīng)驗?zāi)B(tài)分解(EEMD)處理，結(jié)果如下：

a.阿克蘇河徑流周期分析。當添加白噪聲為0.2，采樣頻率為1000時，結(jié)果見圖6，其中IMF1～IMF4分別為從高頻到低頻的固有模態(tài)函數(shù)，代表時間序列不同尺度的周期，Res為殘余項，代表序列的趨勢項(下同)。

利用最大熵譜分析法(Maximum Entropy Spectrum Analysis，MESA)[8]對IMF分量進行處理，得到各分量周期。其中：IMF1的周期為3.3年，方差貢獻率59%，表明其周期變化對序列的周期變化起主要作用；IMF2的周期為7.5年，方差貢獻率14%；IMF3的周期為15年，方差貢獻率為16%；IMF4的周期為30年，方差貢獻率為11%。由各IMF分量的年際變化可以看出，各分量在20世紀70年代波動劇烈，能量集中，徑流豐枯變化頻繁，進入80年代后振幅減小，波動減緩，進入90年代后，各分量振幅顯著增加，徑流豐枯變化加劇。由趨勢項可知，協(xié)合拉站年徑流呈增加趨勢，在2000年以后增速有所減緩。

圖6 阿克蘇河年徑流EEMD分解和方差貢獻率

b.葉爾羌河徑流周期分析。結(jié)果見圖7，利用MESA對各IMF分量進行處理，得到各分量周期。其中：IMF1的周期為4.3年，方差貢獻率為49%，表明其周期變化對序列的周期變化起主要作用；IMF2的周期為15年，方差貢獻率為18%；IMF3與IMF4的周期均為30年，其方差貢獻率分別為28%和5%。IMF1的周期分量在1960—1975年波動劇烈，能量集中，徑流豐枯變化頻繁，之后振幅開始減小，波動減緩，直到進入20世紀90年代后，各個分量振幅才顯著增加，徑流豐枯變化加劇，但在20世紀90年代后期又進入平穩(wěn)狀態(tài)，IMF2周期分量與其相似；由趨勢項可以看出，卡群站年徑流呈增加趨勢，2000年以后增速呈上升趨勢。

圖7 葉爾羌河年徑流EEMD分解和方差貢獻率

c.和田河徑流周期分析。利用MESA對各IMF分量進行處理，得到各分量周期(見圖8)。其中：IMF1的周期為3.3年，方差貢獻率為68%，表明其周期變化對序列的周期變化起主要作用；IMF2的周期為10年，方差貢獻率為28%；IMF3的周期與IMF4的周期均為30年，方差貢獻率分別為1.5%和2.2%。由各IMF的年際變化可知，IMF1在1960—1985年波動劇烈，能量集中，徑流豐枯變化頻繁，之后振幅開始減小，波動減緩，直到進入20世紀90年代后，振幅才顯著增加，徑流豐枯變化加劇，但在20世紀90年代后期又進入平穩(wěn)狀態(tài)。由趨勢項可以看出，和田河年徑流呈增加趨勢，在2000年以后增速顯著上升。

圖8 和田河年徑流EEMD分解和方差貢獻率

d.塔里木河干流徑流周期分析。采用MESA對各IMF分量進行處理，得到其周期(見圖9)。其中：IMF1周期為6年，方差貢獻率為16%；IMF2和IMF3的周期均為15年，方差貢獻率分別為32%和48%，表明這兩個周期分量起主要作用；IMF4為30年，其方差貢獻率僅為4%。由各IMF分量年際變化可以看出，IMF1在2005年以前，基本沒有波動，2005年后先下降后又在2010年上升；由趨勢項可知，阿拉爾站徑流呈遞減趨勢，在2009年以后有增加趨勢。IMF1和IMF2周期分量均在2009年后顯著上升，表明塔里木河干流進入豐水期。

圖9 塔里木河干流徑流EEMD分解和方差貢獻率

4 結(jié) 語

將EEMD應(yīng)用于河川徑流時間序列，可提取可靠真實的徑流變化信號，得到徑流變化的固有時間尺度。本文借助基于樣本熵的小波閾值去噪方法對徑流時間序列去噪，再對去噪后的序列作鏡像延伸處理，以便抑制使用EEMD時序列兩端產(chǎn)生的端點效應(yīng)，結(jié)合最大熵譜分析法對EEMD產(chǎn)生的固態(tài)模態(tài)函數(shù)進行周期分析，得到塔里木河“三源一干”徑流序列(1960—2015年)周期演變和變化趨勢。結(jié)果表明：

a.阿克蘇河年徑流序列存在3.3年、7.5年、15年與30年4個周期，主周期為3.3年；近56年徑流呈增加趨勢，在2000年以后增速有所減緩。

b.葉爾羌河年徑流序列存在4.3年、15年與30年3個周期，主周期為4.3年；近56年徑流呈增加趨勢，2000年以后上升趨勢明顯。

c.和田河年徑流序列存在3.3年、10年和30年3個周期，主周期為3.3年；近56年徑流呈增加趨勢，2000年以后增速顯著上升。

d.塔里木河干流徑流序列存在6年、15年和30年3個周期，主周期為15年；近56年徑流呈遞減趨勢，在2009年以后有增加趨勢，表明塔里木河干流進入豐水期。