張金萍，肖宏林，張 鑫

(鄭州大學(xué)水利與環(huán)境學(xué)院，河南 鄭州 450001)

黃河以泥沙多而聞名，因泥沙淤積造成的黃河水患給國家和人民造成了極大危害。因此，研究黃河徑流-泥沙關(guān)系及其演化特征具有重要的意義。胥維坤[1]運(yùn)用統(tǒng)計(jì)學(xué)方法分析了小浪底水庫運(yùn)行后對(duì)黃河下游水沙變化的影響。方怒放[2]把三峽水利樞紐中的王家橋小流域作為研究區(qū)，利用多元相關(guān)分析和逐步回歸分析研究了降雨、徑流、泥沙三者之間對(duì)應(yīng)的關(guān)系性。吳敬東等[3]建立了小流域坡面不同治理措施下的徑流量與降雨量、降雨歷時(shí)之間的數(shù)學(xué)模型,確定了徑流量和輸沙量之間的關(guān)系。河川水沙常常受多種水文因素的影響，在不同時(shí)間尺度下呈現(xiàn)復(fù)雜的波動(dòng)特征和關(guān)系特性，因此對(duì)水沙研究應(yīng)包含時(shí)序多時(shí)間尺度細(xì)部演化特征分析[4]。本文針對(duì)流域產(chǎn)匯流系統(tǒng)在多時(shí)間尺度下的徑流-泥沙關(guān)系，通過經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)方法分析其原始序列的細(xì)部演變規(guī)律和特征，同時(shí)結(jié)合信息熵理論，辨識(shí)河流水沙關(guān)系及其演化特征。

1 數(shù)據(jù)來源及研究方法

1.1 數(shù)據(jù)來源

潼關(guān)水文站位于黃河中游，是黃河上第二大水文站，潼關(guān)斷面在黃河、北洛河、渭河3條河流匯流下游不遠(yuǎn)處，洪水常常在這里相遇，洪水此消彼長，水沙條件復(fù)雜。本研究獲取了潼關(guān)水文站1919—2015年實(shí)測(cè)徑流量和輸沙量資料作為研究數(shù)據(jù)。潼關(guān)水文站地理位置和研究時(shí)段內(nèi)的年徑流量、輸沙量原始序列分別見圖1、圖2。

圖1 潼關(guān)水文站地理位置

圖2 徑流、泥沙原始序列

1.2 研究方法

1.2.1EMD方法

EMD方法是Huang等[5]于1998年提出的一種新型自適應(yīng)信號(hào)時(shí)頻分析方法，近年來在水文水資源領(lǐng)域廣泛應(yīng)用[6]。該方法是通過一個(gè)“篩選”過程從被分析原始信號(hào)中分離本征模函數(shù)(intrinsic mode function, IMF)，并且提取的IMF分量需要滿足以下條件：①在整個(gè)時(shí)段內(nèi)包含所有數(shù)據(jù)信息的集合中，其數(shù)據(jù)極大值點(diǎn)與極小值點(diǎn)的數(shù)量必須等于經(jīng)過0點(diǎn)的數(shù)量或者數(shù)量差值為一個(gè)；②在任何時(shí)刻，其數(shù)據(jù)的極大值點(diǎn)連接形成的包絡(luò)線與極小值點(diǎn)連接形成的包絡(luò)線的總體均值最終為0[7]。對(duì)于復(fù)雜信號(hào)若不滿足IMF所要求的條件，可以利用EMD方法按照分解的步驟對(duì)復(fù)雜信號(hào)進(jìn)行分解形成一系列滿足條件的IMF分量[8]。EMD方法的步驟為：

步驟1：首先定義一個(gè)原始信號(hào)為X(t)。

步驟2：分析原始信號(hào)X(t)找出它全部的極大值點(diǎn)與極小值點(diǎn)及其對(duì)應(yīng)位置，并利用插值方法處理形成上下包絡(luò)線，求出極大值點(diǎn)和極小值點(diǎn)形成的包絡(luò)線的均值M1和原始信號(hào)X(t)與均值M1之間的差值H1，即H1=X(t)-M1。若H1滿足IMF的兩個(gè)條件，則進(jìn)行步驟4，否則進(jìn)行步驟3。

步驟3：把H1看作原始信號(hào)重復(fù)步驟2求得上下包絡(luò)線形成的均值M2，再通過上述方法計(jì)算出差值H2，并判斷其是否滿足IMF的條件，若還不滿足繼續(xù)重復(fù)直至找到滿足IMF條件的分量。

步驟4：將第一個(gè)IMF分量從X(t)中分離，從而得到R1，即R1=X(t)-H1。

步驟5：重新把分離后的差值R1看作原始信號(hào)并重復(fù)步驟2～4得到n個(gè)滿足條件的IMF分量和1個(gè)趨勢(shì)項(xiàng)(residual, RES)分量Rn，從而完成對(duì)原始信號(hào)的分解過程。

通過以上步驟可將潼關(guān)水文站徑流量、輸沙量進(jìn)行EMD分解，提取多個(gè)IMF分量。

1.2.2信息熵理論

Shannon[9]參考熱力學(xué)上的知識(shí)，把信息熵定義為排除了信息冗余后的平均信息量。若系統(tǒng)變量有U1，U2，…，Un共n種取值，對(duì)應(yīng)概率分別為：P1，P2，…，Pn且各種取值的出現(xiàn)彼此獨(dú)立，此時(shí)信源的平均不確定性應(yīng)為每個(gè)符號(hào)不確定性-lnPi統(tǒng)計(jì)得到的數(shù)學(xué)期望稱為信息熵[10-11]。信息熵H(U)的計(jì)算公式為

(1)

信息論中用信息熵來表示系統(tǒng)本身所擁有信息量的多少。一個(gè)系統(tǒng)表現(xiàn)得越穩(wěn)定、越有規(guī)律，則該系統(tǒng)對(duì)應(yīng)的信息熵值就越小，所以信息熵也用來表示系統(tǒng)的穩(wěn)定性程度、有序化程度[12-13]。

互信息在信息論中表示兩變量之間的相關(guān)性，它反映隨機(jī)變量之間共有的信息量或相互間的依賴關(guān)系[14]。若兩隨機(jī)變量間的聯(lián)系越緊密，則它們的互信息熵值就越大，同理在相反情況下，則兩者之間的互信息熵值就越小[15]。

分別計(jì)算徑流、泥沙原序列及分解后的IMF分量的信息熵，并通過計(jì)算徑流、泥沙的互信息從而確定黃河流域徑流-泥沙的關(guān)系。以潼關(guān)水文站實(shí)測(cè)徑流量或輸沙量為例，設(shè)原始信號(hào)X(t)為潼關(guān)水文站1919—2015年徑流量或輸沙量數(shù)據(jù)，經(jīng)EMD分解可得到n個(gè)不同的IMF分量，每個(gè)分量包含m年的實(shí)測(cè)分解值。則得有判定矩陣：

C=(Cst)n×m(s=1,2,…,n；t=1,2,…,m)

(2)

式中Cst為第s個(gè)IMF分量的第t年的實(shí)測(cè)分解值。

為了排除干擾，在計(jì)算各IMF分量的信息熵時(shí)，首先需要對(duì)n個(gè)IMF分量的所有數(shù)據(jù)作歸一化處理。得到歸一化矩陣B,歸一化處理公式為

(3)

式中：Bst為B的元素；Cmax、Cmin分別為同一IMF分量系列數(shù)據(jù)的最大值和最小值。

將區(qū)間[0,1]平均分成L組，每組長度為1/L，對(duì)應(yīng)的區(qū)間為[0，1/L)，[1/L，2/L)，…，[L-1/L，1]。對(duì)于某個(gè)IMF分量m年歸一化數(shù)據(jù)中第t年在相應(yīng)的某個(gè)分組區(qū)間內(nèi)的概率Pl為t/m[16]。由信息熵的定義和計(jì)算公式知，為保證lnPl始終有意義，需要對(duì)Pl進(jìn)行修正，即當(dāng)Pl為0時(shí)假設(shè)lnPl為0。則第s個(gè)IMF分量的信息熵為

(4)

根據(jù)式(4)可依次求出經(jīng)EMD分解的IMF分量信息熵值，同理求得趨勢(shì)項(xiàng)的信息熵值。在此基礎(chǔ)上，利用式(5)求出它們之間的互信息熵，得出相互間的關(guān)聯(lián)度。

MI(a,b)=H(a)+H(b)-H(a,b)

(5)

式中：MI(a,b)為互信息熵值；H(a)、H(b)分別為徑流量、輸沙量各分量的信息熵；H(a,b)為二者對(duì)應(yīng)各分量之間的聯(lián)合熵。

2 結(jié)果分析

2.1 EMD方法分析

應(yīng)用前文EMD分解步驟對(duì)潼關(guān)水文站1919—2015年徑流量、輸沙量原始序列進(jìn)行分解，分別得到徑流量和輸沙量的4個(gè)IMF分量與1個(gè)RES分量。徑流量、輸沙量各分量序列見圖3。

由圖3可知，徑流量和輸沙量均包含著復(fù)雜得多時(shí)間尺度周期性變化和宏觀趨勢(shì)走向。徑流量和輸沙量IMF1分量對(duì)應(yīng)的準(zhǔn)周期約為5 a和4 a，IMF2分量的準(zhǔn)周期約為8 a和6 a，IMF3分量的準(zhǔn)周期約

(a)IMF1分量

(b)IMF2分量

(c)IMF3分量

(d)IMF4分量

(e)RES分量

為16 a和15 a，IMF4分量的準(zhǔn)周期約為35 a和28 a，可見徑流和泥沙各IMF分量周期性關(guān)系在不同時(shí)間尺度上變化規(guī)律基本保持一致。這表明水沙相關(guān)性密切，水沙長時(shí)間序列在模糊特征下具有確定性特征，即存在波動(dòng)周期的確定性，可為不同時(shí)期水沙聯(lián)合預(yù)測(cè)、水沙變化特征研究提供參考。徑流量、輸沙量原始序列影響因素較多，經(jīng)過EMD分解后隨著分解層數(shù)的增多分解尺度的加大，徑流、泥沙的干擾信息被去除，影響因素減少，有效信息增多，在不同的分層下表現(xiàn)為越來越有規(guī)律性。在趨勢(shì)變化方面，表現(xiàn)為分解序列更有序，系統(tǒng)由無序到有序過渡的趨勢(shì)性。特別是徑流量、輸沙量的趨勢(shì)項(xiàng)呈現(xiàn)出較強(qiáng)的規(guī)律性，都呈現(xiàn)先上升后逐漸下降的趨勢(shì)，徑流上升速度比泥沙上升速度略快，二者趨勢(shì)走向基本同步。

潼關(guān)水文站徑流量、輸沙量的趨勢(shì)項(xiàng)在1950年左右分別出現(xiàn)拐點(diǎn)，在之后很長一段時(shí)間內(nèi)不斷減少，這是受自然和人類活動(dòng)等多重因素影響的結(jié)果。1949年以來，水利工程建設(shè)、河道外取水、水土流失綜合治理措施和修建淤地壩等生態(tài)工程建設(shè)的增加及其他人類活動(dòng)干預(yù)對(duì)水文系統(tǒng)產(chǎn)生重大影響，使潼關(guān)水文站徑流量，輸沙量出現(xiàn)明顯的減少[17]。整體而言，年徑流量與年輸沙量隨時(shí)間表現(xiàn)為顯著下降的趨勢(shì)與文獻(xiàn)[18]的結(jié)論相符。這表明原序列經(jīng)EMD分解的RES分量可以反映徑流量、輸沙量整體的變化趨勢(shì)，即低頻模態(tài)分量控制著序列變化的全局和趨勢(shì)。

2.2 信息熵分析

利用前文方法計(jì)算徑流量、輸沙量分解序列的信息熵結(jié)果見表1。徑流量、輸沙量IMF1分量對(duì)應(yīng)的波動(dòng)周期為短周期，IMF2分量對(duì)應(yīng)的波動(dòng)周期為中周期，IMF3分量對(duì)應(yīng)的波動(dòng)周期為中長周期，IMF4分量對(duì)應(yīng)的波動(dòng)周期為長周期，因RES分量缺乏周期波動(dòng)，故不與IMF分量作對(duì)比。

表1 徑流量與輸沙量分量信息熵

由表1可見，徑流量和輸沙量原序列熵值較大，徑流量熵值拐點(diǎn)首次出現(xiàn)在IMF1處，輸沙量熵值拐點(diǎn)出現(xiàn)在IMF2處。表示隨著分解層數(shù)的增多、分解尺度的加大，徑流、泥沙分解序列表現(xiàn)為越來越有序，影響因素越來越少，系統(tǒng)呈現(xiàn)由不穩(wěn)定向穩(wěn)定過渡、由無序到有序的過渡的趨勢(shì)性，因此信息熵值減小。但當(dāng)分層過多，分解尺度過大，使系統(tǒng)的其他影響因素體現(xiàn)出來，從而使得信息熵值增大。

原序列對(duì)應(yīng)的徑流熵值比泥沙熵值稍大，可以判定徑流系統(tǒng)較泥沙復(fù)雜，須密切關(guān)注徑流量的變化。徑流、泥沙各序列信息熵值變化基本同步，則表明徑流量、輸沙量之間關(guān)系應(yīng)表現(xiàn)為同步變化或者異步變化特征，由徑流量、輸沙量原序列經(jīng)EMD方法分析得知二者滿足同步變化關(guān)系。

2.3 互信息分析

為了深層次分析潼關(guān)水文站徑流-泥沙的關(guān)系，計(jì)算各序列與原序列之間的互信息熵值以及徑流量與輸沙量原序列之間的互信息熵值，結(jié)果見表2。

表2 徑流量和輸沙量原序列與其分解序列之間的互信息熵值

由表2可知，徑流量、輸沙量的原序列與RES分量之間的互信息最大，可見如果徑流、泥沙系統(tǒng)比較復(fù)雜，為了便于研究徑流、泥沙變化可采用RES分量數(shù)據(jù)信息闡述原始數(shù)據(jù)信息，并利用趨勢(shì)項(xiàng)變化特征觀察其整體發(fā)展態(tài)勢(shì)。徑流原序列與各IMF分量的互信息呈現(xiàn)遞增的趨勢(shì)性，即它們之間的共有信息量越來越多，原序列與IMF4的互信息最大，說明IMF4對(duì)徑流原序列的貢獻(xiàn)程度最大。對(duì)于徑流IMF分量而言，低頻分量與原序列的互信息較大，所以低頻分量所攜帶的信息量比較多，因此對(duì)徑流量宜采用長周期研究比較可靠。泥沙原序列與各IMF分量的互信息先減小后增大,泥沙原序列與IMF1互信息最大，說明IMF1對(duì)泥沙原序列的貢獻(xiàn)率最大。對(duì)于泥沙IMF分量而言，高頻分量與原序列的互信息較大，高頻分量所攜帶的信息量比較多，因此對(duì)輸沙量宜采用短周期研究比較可靠。

為了更進(jìn)一步研究徑流量、輸沙量之間的關(guān)系，計(jì)算了徑流與泥沙之間原序列及各分解序列的互信息熵值，原序列之間為0.689，IMF1分量之間為0.575，IMF2分量之間為0.338，IMF3分量之間為0.534，IMF4分量之間為0.675，RES分量之間為1.410。可以看出，趨勢(shì)項(xiàng)之間的互信息最大，這是由于經(jīng)EMD方法分解后排除了干擾，使徑流、泥沙的有效信息充分體現(xiàn)出來，從而使兩者之間的趨勢(shì)項(xiàng)互信息比原序列之間的互信息大。原序列之間互信息也較大，表示雖然徑流、泥沙原序列的影響因素比分解序列要復(fù)雜，系統(tǒng)比較無序，但是原序列信息量比較多。對(duì)于各分解序列來說，互信息基本呈現(xiàn)遞增的趨勢(shì)性，表示經(jīng)EMD分解排除干擾后各系列之間的關(guān)聯(lián)性增大，它們之間的共有信息量逐漸增大。對(duì)于徑流量與輸沙量經(jīng)EMD分解的各IMF分量而言，徑流量與輸沙量對(duì)應(yīng)的長周期IMF4分量互信息最大，因此徑流量與輸沙量關(guān)系及其變化特征宜采用長周期的研究。

3 結(jié) 論

潼關(guān)水文站年徑流量與年輸沙量隨時(shí)間表現(xiàn)為顯著下降的趨勢(shì)，且徑流量和輸沙量有著很強(qiáng)的相關(guān)性，這種相關(guān)性不僅體現(xiàn)在宏觀趨勢(shì)上，而且還體現(xiàn)在微觀局部特征上，具體表現(xiàn)為在多時(shí)間尺度下波動(dòng)特征基本同步。EMD方法與信息熵結(jié)合分析表明，徑流量與其長周期分量互信息最大，因此對(duì)徑流量變化特征宜采用長周期研究比較可靠；輸沙量與其短周期分量互信息最大，因此對(duì)輸沙量變化特征宜采用短周期研究比較可靠；而徑流量與輸沙量的長周期分量及趨勢(shì)項(xiàng)互信息較大，所以黃河中游潼關(guān)附近水沙關(guān)系可采用長周期研究，趨勢(shì)項(xiàng)分量數(shù)據(jù)信息闡述原始數(shù)據(jù)信息比較可靠。