張金萍，肖宏林，張 鑫

(鄭州大學(xué)水利科學(xué)與工程學(xué)院，鄭州 450001)

0 引 言

大型水庫的運行對河流水沙過程具有很強的調(diào)控作用，往往引起河流生態(tài)環(huán)境演化[1]、河床演變[2]以及水文情勢變化[3]，其影響甚至是深遠的。因此，研究水庫運行對水沙的調(diào)節(jié)作用，對于認(rèn)識水沙變化成因和演化機理，優(yōu)化水庫運行方式具有重要的意義。目前，有關(guān)水庫運行對徑流和泥沙的調(diào)控作用已有不少研究，例如，韓閃閃等采用水沙指標(biāo)變化幅度公式定量分析了小浪底水庫運行對黃河下游水沙變化的影響[4]，班璇等采用變化范圍法定量評價了三峽水庫蓄水后長江中游水沙時空變化[5]。實際上，與其他隨機變量時間序列一樣，徑流和泥沙時間序列也具有多時間尺度特征[6]。與常規(guī)日、月、季時間尺度不同，該多時間尺度性是指系統(tǒng)變量的變化在某一時間段內(nèi)不是只以一種固定的頻率(周期、時間尺度)在運動，而是同時包含著各種頻率(周期、時間尺度)的變化和局部波動，是復(fù)雜非線性系統(tǒng)的重要演化特性之一。因此，將多時間尺度概念引入水庫運行對水沙變化的研究中，分析水庫運行對水沙條件變化局部特征的影響，可以掌握水沙時間序列微觀和宏觀層面的本質(zhì)特征，使徑流和泥沙研究機制更加科學(xué)，更接近實際。

為了更清晰了解水庫運行對水沙的調(diào)控作用，本文以黃河上游龍羊峽水庫為例，利用互補集合經(jīng)驗?zāi)B(tài)分解方法分析水沙各分解序列的細部演化特征，采用多時間尺度熵和納什效率系數(shù)分析建庫前后水沙各分解序列的變化特征和各分量對原始序列的貢獻程度，以期辨識水庫對河流水沙微觀和宏觀不同層面的變化影響。

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

1.1 研究區(qū)概況

龍羊峽水庫是黃河上游具有多年調(diào)節(jié)性能的綜合利用型水庫，水庫設(shè)計蓄水位為2 600 m，總庫容247 億m3，調(diào)節(jié)庫容194 億m3，屬于大一型水庫，同時該水利工程兼顧防洪、發(fā)電、灌溉、防凌等功能，具有很大的綜合效益[7]。貴德水文站位于龍羊峽水庫大壩下游54.8 km處，是水庫下游控制性水文站。龍羊峽水庫和貴德水文站地理位置如圖1所示。本文獲取了貴德水文站1959-2013年共55 a實測年徑流量、年輸沙量數(shù)據(jù)，根據(jù)龍羊峽水庫修建時間對研究時段進行劃分，其中1959-1986年為建庫前階段，1987-2013年為建庫后階段。

圖1 龍羊峽水庫和貴德水文站位置圖

1.2 研究方法

1.2.1 CEEMD方法

經(jīng)驗?zāi)B(tài)分解方法(EMD方法)是黃鍔等人于1998年提出的一種新型自適應(yīng)信號時頻分析方法[8]，它是通過一個“篩選”過程從復(fù)雜的多分量原始信號中分離本征模態(tài)函數(shù)(IMF)，但在計算過程中可能會出現(xiàn)模態(tài)混淆現(xiàn)象，端點效應(yīng)等問題。WU 等人提出了一種集總經(jīng)驗?zāi)B(tài)分解(EEMD)方法，它是在EMD方法的基礎(chǔ)上向原始信號加入白噪聲，可以有效解決EMD方法的模態(tài)混淆問題，但會在信號重構(gòu)序列中存在殘余噪聲[9]； Yeh等對EEMD方法作了進一步的改進，加入的輔助噪聲采用正、負成對的形式，不但能夠消除重構(gòu)信號中的殘余噪聲，而且加入的噪聲集合次數(shù)可以很低，計算效率較高，這種方法被稱為互補集合經(jīng)驗?zāi)B(tài)分解(CEEMD)方法[10]。

CEEMD方法是在EMD方法基礎(chǔ)上進行的，具體步驟如下：

首先定義原始信號為X(t)；

向原始信號中加入n組正、負成對出現(xiàn)的輔助白噪聲，這樣得到的集合信號個數(shù)為2n，即產(chǎn)生兩套集合IMF分量。公式如下：

(1)

式中：B為輔助白噪聲；A1、A2是加入成對正、負輔助白噪聲后形成的信號。

對加入輔助白噪聲后的2n個信號進行EMD分解，每個信號通過分解均可得到一組IMF分量，則最終產(chǎn)生2n組IMF分量。

通過對2n組對應(yīng)IMF分量求取平均值，即可得到相應(yīng)的IMF分量。

(2)

式中:Cij為第i個信號分解的第j個IMF分量；Cj為第j個IMF分量。

通過分解最終可將原始信號分為m個IMF分量和一個趨勢項(RES分量)Rm，即：

(3)

1.2.2 多時間尺度熵

熵來源于熱力學(xué)，香農(nóng)于1948年將其引入信息論中，他把排除信息冗余后的平均信息量稱為信息熵，并給出了求解信息熵的公式[11]。若系統(tǒng)變量有N種取值：U1…Ui…UN，對應(yīng)概率為：P1…Pi…PN，且各種取值的出現(xiàn)彼此獨立。此時信源的平均不確定性，即每個符號不確定性-lnPi統(tǒng)計得到的數(shù)學(xué)期望稱為信息熵[12]，信息熵表示系統(tǒng)本身所擁有信息量多少、系統(tǒng)的復(fù)雜程度[13]。其計算公式為：

(4)

對于水文長時間序列而言，為了更好地評價建庫前后樣本容量N的取值對信息熵的影響，需要對信息熵公式進行改進。采用信息熵的變形公式，即考慮了建庫前后年份N值的不一致對熵值的影響，其公式如下：

(5)

在本研究中，首先采用CEEMD方法將徑流、泥沙原始序列進行分解，得到對應(yīng)的多個IMF分量，然后將采用信息熵變形公式計算得到的各IMF分量的熵值稱為多時間尺度熵。具體計算過程如下：設(shè)原始信號X(t)為貴德水文站實測徑流量或輸沙量序列，經(jīng)CEEMD分解可得到M個不同的IMF分量，每個分量包含N年的實測分解值, 則得有判定矩陣：

C=(Cjt)M×N(其中j=1,2,…,Mt=1,2,…,N)

(6)

式中：Cjt表示第j個IMF的第t年的實測分解值。

為了排除各IMF分量中負值的干擾，需要對M個IMF分量的所有數(shù)據(jù)作歸一化處理，得到歸一化矩陣B，歸一化處理公式為：

(7)

式中：Cmax,Cmin表示同一IMF分量在N年中的數(shù)據(jù)最大值和最小值。通過(7)式可將各IMF分量數(shù)據(jù)歸一化到[0,1]區(qū)間。

通過信息熵變形公式，則第j個IMF分量的多時間尺度熵為：

(8)

1.2.3 納什效率系數(shù)

納什效率系數(shù)NSE可用來表征經(jīng)CEEMD分解的徑流量、輸沙量各模態(tài)分量及不同模態(tài)分量組合對原始序列的模擬精度，即它能定量評價建庫前后徑流量、輸沙量各模態(tài)分量對原始序列的貢獻程度。求解納什效率系數(shù)的公式如下：

(9)

NSE的取值范圍為[-∞,1]，對于不同模態(tài)分量組合值而言，NSE越接近1，表示模擬效果越好，模型組合可信度越高，對于單獨各模態(tài)分量而言，表示對原始序列的貢獻度越大；NSE越接近0，表示模擬結(jié)果越接近實測值的平均值水平，但過程模擬誤差大；當(dāng)NSE遠遠小于0，則表示模型的模擬結(jié)果是不可信的[14]。

2 結(jié)果分析

2.1 CEEMD分解

利用CEEMD方法分別對建庫前后兩個時段的徑流量、輸沙量原始序列進行分解，徑流量、輸沙量原始序列及各分解序列數(shù)據(jù)如圖2和圖3所示。

由圖可以看：

(1)無論是建庫前還是建庫后，貴德水文站的徑流量和輸沙量均包含著復(fù)雜的多時間尺度周期性變化和宏觀趨勢走向。建庫前徑流量和輸沙量具有相同的多時間尺度準(zhǔn)周期，3個IMF分量的準(zhǔn)周期分別約為3、7和17 a，除振幅波動有差異外，二者基本上滿足同頻同步波動變化。建庫前RES分量表明二者宏觀趨勢走向基本一致，呈現(xiàn)增加的趨勢。

(2)建庫后徑流量各IMF分量對應(yīng)的準(zhǔn)周期分別約為4、9和20 a，輸沙量則分別約為4、8和10 a。可以看出，龍羊峽水庫運行以后，受水庫攔蓄和調(diào)度的影響，出庫的徑流量和輸沙量更趨于平穩(wěn)。與建庫前相比，徑流量和輸沙量的各時間尺度振幅均呈現(xiàn)減小趨勢，中高頻分量對應(yīng)的準(zhǔn)周期增大，而低頻分量對應(yīng)的徑流量準(zhǔn)周期增大，輸沙量準(zhǔn)周期減小。這表明水庫運行對水沙的影響并不一致，由此導(dǎo)致徑流量和輸沙量的多時間尺度變化不僅振幅差異較大，而且波動周期也差異較大，甚至在個別時間尺度上波動變化呈現(xiàn)異步狀態(tài)，最終導(dǎo)致宏觀趨勢走向呈現(xiàn)完全相反的態(tài)勢。

(3)RES分量反映的是徑流量、輸沙量整體的變化趨勢特征，即低頻長周期模態(tài)分量控制著序列變化的全局和趨勢。建庫前徑流量和輸沙量的RES分量呈現(xiàn)隨時間逐漸上升的趨勢，建庫后徑流量RES分量仍逐漸上升，但輸沙量卻趨于下降。這表明龍羊峽水庫運行對下游輸沙量的影響作用比徑流量大。水庫蓄水?dāng)r沙，合理調(diào)度，年徑流量雖然減少，但上升趨勢仍很明顯，而年輸沙量不僅在數(shù)量發(fā)生了絕對的減少，其趨勢性也發(fā)生了根本性改變，呈現(xiàn)下降態(tài)勢。

圖2 建庫前徑流量、輸沙量各序列變化圖

圖3 建庫后徑流量、輸沙量各序列變化圖

2.2 多時間尺度熵計算

利用多時間尺度熵計算方法，計算建庫前后徑流量、輸沙量分解序列的多時間尺度熵，結(jié)果如表1所示。根據(jù)其分解后的波動周期大小，將各IMF分量的準(zhǔn)周期分別稱為短周期、中周期和長周期，因RES分量缺乏周期波動，故不與IMF分量作對比分析。

表1 建庫前后徑流量、輸沙量各分量的熵值

通過上表1可知：

建庫前，徑流量、輸沙量在多時間尺度下的熵值變化規(guī)律并不明顯，但在建庫后，徑流量、輸沙量的熵值隨波動周期增大而不斷減小。此外，受水庫運行的影響，熵值在IMF分量和RES分量的變化明顯不同。與建庫前相比，建庫后徑流量、輸沙量各IMF分量的熵值變大，但RES分量卻變小。表明水庫運行后不同時間尺度的徑流量和輸沙量隨機性增強，可預(yù)測性降低，但就宏觀發(fā)展趨勢而言，可預(yù)測性升高。水庫運行、人類活動等因素雖然加劇了徑流量、輸沙量在不同時間尺度下的復(fù)雜程度，但卻促使其變化越來越有序，從而使整個水沙系統(tǒng)由不穩(wěn)定向穩(wěn)定逐漸過渡。

建庫前，徑流量IMF2分量熵值最大，輸沙量IMF1分量熵值最大，表明建庫前徑流量中頻分量和輸沙量高頻分量所攜帶的信息量最多。因此，建庫前對龍羊峽水庫下游徑流量、輸沙量的研究可集中于約7 a和3 a的時間周期上。徑流量、輸沙量對應(yīng)的RES分量熵值也較大，其攜帶信息量也較多，因此建庫前對徑流量和輸沙量的研究也應(yīng)注重對其整體發(fā)展態(tài)勢的研究和分析上。

建庫后，徑流量、輸沙量IMF1分量對應(yīng)的熵值最大，表明二者短周期分量所攜帶的信息量最多，因此，建庫后對水庫下游徑流量、輸沙量的研究可集中于約4 a的時間周期上。同時，建庫后徑流量、輸沙量的RES分量熵值大幅減少，特別是輸沙量的RES分量熵值減幅最大，說明水庫運行對河流輸沙量的影響程度要比徑流量的影響大。

2.3 納什效率系數(shù)計算

通過計算各分量的納什效率系數(shù)，用以分析建庫前后徑流量、輸沙量各分量對原始序列的貢獻程度，結(jié)果如表2所示。

由表2可知：對于徑流量IMF分量而言，建庫前IMF2分量對應(yīng)的NSE最大，為0.416，即IMF2分量對徑流原始序列的模擬貢獻程度最大；建庫后IMF1分量對應(yīng)的NES最大，為0.391，即IMF1分量對徑流原始序列的模擬貢獻程度最大。對于輸沙量IMF分量而言，建庫前后的IMF1分量對應(yīng)的NES均最大，其對輸沙量原始序列的模擬貢獻程度最大。由此可見，建庫前對徑流量的研究主要基于中周期的演變特征，輸沙量則主要集中于短周期上，而建庫后對徑流量和輸沙量的研究則集中于短周期，這與多時間尺度熵分析結(jié)果相對應(yīng)。

另一方面，建庫前后徑流量NSE的方差分別為0.025和0.023，輸沙量NSE的方差分別為0.073和0.039。可以看出，建庫后徑流量、輸沙量NSE的方差比建庫前小，且輸沙量NSE的方差減小程度較大。這表明建庫后徑流量、輸沙量各分量的NSE更接近平均水平，且由于水庫的調(diào)控作用，能使各時間尺度的徑流量、輸沙量分配更均勻，達到削減洪峰流量、調(diào)水調(diào)沙的目的。

3 結(jié) 論

本文將CEEMD方法與熵理論、納什效率系數(shù)結(jié)合，分析水庫運行對水沙調(diào)控的影響，結(jié)果表明：

(1)龍羊峽水庫運行對下游水沙條件變化均有顯著影響，但對輸沙量的影響要大于徑流量。同時，水庫的運行會加劇徑流量、輸沙量系統(tǒng)的復(fù)雜性，但是隨著多時間尺度波動周期的增大，水沙各分解序列的復(fù)雜性降低，可預(yù)測性升高。

(2)建庫前，徑流量與輸沙量有著很強的相關(guān)性，這種相關(guān)性不僅體現(xiàn)在宏觀趨勢上，而且還體現(xiàn)在微觀局部特征上，具體表現(xiàn)為各分量的水沙變化基本滿足同頻同步。水庫運行后水沙各分解序列同頻同步變化被打破，多時間尺度下水沙變化和波動周期均存在較大差異，甚至在個別時間尺度上波動變化呈現(xiàn)異步狀態(tài)。

(3)建庫前徑流量中頻IMF分量、輸沙量高頻IMF分量所攜帶的信息量較多，建庫后二者的信息量都集中于高頻IMF分量，并且高頻IMF分量對應(yīng)的周期增大，同時建庫前后對應(yīng)分量的NSE也較大。所以綜合考慮建庫前后波動周期的變化，建議今后黃河上游對徑流量和輸沙量的相關(guān)監(jiān)測和研究周期分別為4～7 a和3～4 a。

□