丁志宏，張金萍，趙焱

（1.天津市中水科技咨詢有限責任公司，天津300170；2.鄭州大學水利與環(huán)境學院，河南鄭州450001；3.黃河水利科學研究院水資源研究所，河南鄭州450003）

基于CEEMDAN的黃河源區(qū)年徑流量多時間尺度變化特征研究

丁志宏1，張金萍2，趙焱3

（1.天津市中水科技咨詢有限責任公司，天津300170；2.鄭州大學水利與環(huán)境學院，河南鄭州450001；3.黃河水利科學研究院水資源研究所，河南鄭州450003）

基于EMD、EEMD和互補EEMD方法的不足之處，將CEEMDAN方法應用于黃河源區(qū)唐乃亥水文站1956—2015年的年徑流量序列多時間尺度分析，采用納什效率系數(shù)定量評價了各階模態(tài)的重構(gòu)信號對年徑流量序列模擬精度的貢獻程度，指出了提高年徑流量預測精度的工作方向，揭示了黃河源區(qū)水文水資源系統(tǒng)變化的復雜年際特征，解釋了20世紀90年代至21世紀初黃河連續(xù)枯水時段的成因。

黃河源區(qū)；唐乃亥；CEEMDAN；年徑流量；多時間尺度；納什效率系數(shù)；枯水時段

DOI

黃河是中華民族的母親河，以占全國河川徑流總量2.2%的天然徑流量支撐著全國總?cè)丝诘?2%、耕地總面積的15%，除了承擔著本流域的供水任務外，還承擔著向流域外的河北省和天津市引黃供水的任務，黃河水資源量的變化對于相關(guān)區(qū)域社會經(jīng)濟發(fā)展具有重大影響。

黃河源區(qū)是指黃河從河源至唐乃亥水文站之間的高寒草甸草原區(qū)，唐乃亥水文站控制流域面積12.19萬km2。黃河源區(qū)以占黃河流域面積13%的匯水面積貢獻了黃河年徑流量的33%，是黃河流域最重要的產(chǎn)流區(qū)，素有“黃河水塔”之稱，該區(qū)域徑流量的變化對于整個黃河流域水資源的變化具有至關(guān)重要的影響和控制性作用。

徑流量的年際變化過程具有多時間尺度性。所謂多時間尺度性，是指徑流量的變化在某一時間段內(nèi)不是只以一種固定的頻率（周期、時間尺度）在運動，而是同時包含著各種頻率（周期、時間尺度）的變化和局部波動，是包括氣象、水文、土壤、植被、社會等各子系統(tǒng)在內(nèi)的多種動力學機制同時發(fā)揮作用的結(jié)果，是徑流量在時域中呈現(xiàn)復雜變化的根本原因。

綜上所述，為深入分析黃河源區(qū)年徑流量的多時間尺度波動特征，本文運用CEEMDAN方法對黃河源區(qū)的控制性水文站——唐乃亥水文站1956—2015年的年徑流量序列進行分解，探討其在各個時間尺度層次上的波動特征，以期為黃河水資源的合理開發(fā)、高效利用、有效保護、科學調(diào)度等工作提供技術(shù)參考和科學依據(jù)。

1 CEEMDAN基本理論

為了克服小波變換的諸多不足之處，Huang等人于1998年提出了經(jīng)驗模態(tài)分解（Empirical Mode Decomposition，簡稱EMD）[1]。EMD是一種可以用于分析非線性系統(tǒng)產(chǎn)生的非平穩(wěn)序列的數(shù)據(jù)驅(qū)動型的適應性方法，它將一個序列分解為局部的、完全數(shù)據(jù)驅(qū)動的、具有快速和慢速波動周期的一系列分量，這些幅度和頻率經(jīng)過調(diào)制的分量被稱為本征模態(tài)函數(shù)（Intrinsic Mode Function，簡稱IMF）。EMD已在水文水資源領(lǐng)域得到了廣泛應用[2-4]。

但是，EMD算法的局部特性可能會產(chǎn)生一種被稱為模態(tài)混淆（混頻）的現(xiàn)象：在一個模態(tài)中存在具有完全不同的尺度的波動或者在不同的模態(tài)中產(chǎn)生具有相似尺度的波動，而理想的情況是每一個模態(tài)中的尺度是相似的。為了減輕EMD的模態(tài)混淆現(xiàn)象，Huang等人于2009年提出了集合經(jīng)驗模態(tài)分解[5]（EnsembleEmpiricalModeDecomposition，簡稱EEMD），該方法是對帶有噪聲的原始序列的集合進行EMD分解。所謂集合，即是添加了白噪聲的原始序列的若干副本，通過求集合的平均值來得到最終分解結(jié)果。通過添加白噪聲來減少模態(tài)混淆是利用了EMD的二值濾波器組特性以及遍布整個時間—頻率空間的噪聲[6]，以此來求得更多的在整個時間跨度內(nèi)具有相似尺度的更規(guī)則的模態(tài)。盡管EEMD被證明可以大幅減少混頻現(xiàn)象并在水文水資源等領(lǐng)域得到了廣泛應用[7-9]，然而，該方法在解決舊問題的同時也產(chǎn)生了新問題：在EEMD的重構(gòu)序列（即所有模態(tài)之和）中存在著殘留噪聲，另外，EEMD的每一次EMD分解所產(chǎn)生的模態(tài)的個數(shù)可能是不同的，這使得最終在求集合平均值時變得困難。為了解決EEMD存在的問題，Huang等人于2010年提出了互補EEMD[10]（Complementary Ensemble Empirical Mode Decomposition，簡稱互補EEMD），即通過使用互補噪聲對（即增加和減去相反的白噪聲）減輕了重構(gòu)序列存在的噪聲殘留問題。然而，互補EEMD的數(shù)學完整性不能被證明，而且最終在求集合平均值時存在的問題依然沒有得到解決，因為互補EEMD依舊會產(chǎn)生每一次EMD分解所產(chǎn)生的模態(tài)的個數(shù)可能會不同這一問題。

2011年，Torres等人提出了具適應性噪聲的完全集合經(jīng)驗模態(tài)分解[11]（Complete Ensemble Empiri?cal Mode Decomposition with Adaptive Noise，簡稱CEEMDAN），對EEMD進行了重要改進，解決了EEMD存在的上述兩個問題并在多個領(lǐng)域得到了應用[12-14]。2014年，Torres等人又對CEEMDAN進行了改進[15]，完美解決了CEEMDAN初始算法[11]所存在的個別模態(tài)包含殘留噪聲以及分解早期可能存在虛假模態(tài)等兩個問題。

綜上所述，CEEMDAN是在EMD和EEMD的基礎(chǔ)上發(fā)展而來的，本文按照技術(shù)發(fā)展歷程來對CEEMDAN[15]的基本理論介紹如下：

EMD是把1個序列分解為若干數(shù)目的IMF，而IMF必須滿足2個條件：①極值點（極大值和極小值）的個數(shù)和跨零點的個數(shù)必須相等或者至多相差1個；②局部平均值，即上包絡線和下包絡線的平均值，必須為0。EMD算法具體參見文獻[1]。

EEMD把相應的IMF的平均值定義為“真實”模態(tài)，這些IMF是通過向原始序列中添加白噪聲后再進行EMD分解而得到的。設x() t為待分解序列，EEMD算法可描述如下：

在進行EMD分解時，均需進行不同次數(shù)的迭代，迭代終止與否的判斷指標采用限制標準差SD，SD定義為[1]：式中：h1l(t)為EMD分解時第l次篩選所得的數(shù)據(jù)；h1(l-1)(t)為EMD分解時第l-1次篩選所得的數(shù)據(jù)；T為序列長度。

SD的值一般取0.2～0.3，即滿足0.2＜SD＜0.3時分解過程即可結(jié)束。采用此標準的物理考慮是：既要使得d(i)k足夠接近IMF的要求，又要控制分解的次數(shù)，從而使得所得IMF分量保留原始序列中的幅值調(diào)制信息。

值得指出的是，在EEMD中，每個x(t)(i)都是獨立地被分解，對每一個x(t)(i)來說，每一次實現(xiàn)中的每一個分解階段得到的都是獨立的，其間沒有關(guān)聯(lián)。

使用噪聲輔助技術(shù)改進EMD的主要思想是往序列中添加一些可控噪聲，以創(chuàng)造新的極值點。使用這種方式，局部平均值被“強迫”吸引在原始序列中的新極值點被創(chuàng)造出來的那些部分，而同時原始序列中的沒有新的極值點被創(chuàng)造出來的那些部分沒有被改變，即該算法被強迫聚焦到尺度—能量空間的一些特別點上。取平均值就是為了更好地估計局部均值，這些局部均值在原始序列添加噪聲后的各個實現(xiàn)中是略有不同的。

然而，EEMD通過取平均值來估計的是模態(tài)而不是局部均值。這是因為EEMD是獨立地分解每一個具噪聲的原始序列，所以在每一次分解的第1個階段有1個局部均值和1個模態(tài)，則真實模態(tài)就是具噪聲原始序列的EMD分解所求得的模態(tài)的平均，這其中就包含著一些殘余噪聲。這就造成EEMD存在以下問題：①分解是不完全的，即存在重構(gòu)誤差；②每一次所得的模態(tài)個數(shù)可能會不同，造成最后求集合的平均值時存在困難。

在互補EEMD中，噪聲是成對地被添加到原始序列上（一個是正的，一個是負的），由此產(chǎn)生2個集合：式中：y(i)為具噪聲的原始序列的具有互補性的2個副本；x(t)和w(i)同前。

盡管這一方法顯著減小了重構(gòu)序列中的殘余噪聲，但是仍然不能保證和會產(chǎn)生相同數(shù)目的模態(tài)，使得最后的求平均值還是存在困難，同時模態(tài)中仍然存在噪聲殘余。

CEEMDAN的具體算法為：令Ek(?)為通過EMD產(chǎn)生第k階模態(tài)的算子，令M(?)為產(chǎn)生將要被進行分解的序列的局部均值的算子，令w(i)為均值為零、方差為1的白噪聲，是在實現(xiàn)中求取平均值的算子，可以看出E1(x)=x-M(x)，則：

（1）使用EMD計算x(i)=x+β0E1(w(i))（即x的第i次實現(xiàn)）的局部均值，以求得第1個殘差：

（4）對于k=3,…,K，計算第k個殘差：

（5）計算第k階模態(tài)：

（6）返回第4步計算下一個k。

重復進行第4步至第6步直到所求得的殘差滿足以下條件之一為止：①不能被EMD進一步分解為止；②滿足IMF條件；③局部極值點的個數(shù)小于3個。

綜上，經(jīng)CEEMDAN重構(gòu)，最終殘差滿足：

K是模態(tài)的總階數(shù)。因此，原始序列x可以表達為：

上述分解過程確保了CEEMDAN的完整性并因此保證了原始序列得以準確重構(gòu)。模態(tài)的最終階數(shù)只取決于原始序列數(shù)據(jù)和停止準則。系數(shù)βk=εkstdrk允許在每一個階段進行信噪比（SNR）選擇。在EEMD中，噪聲和殘差之間的信噪比SNR隨著階數(shù)k的增加而增加，這是因為當k＞1時，第k階殘差中的能量只是在計算開始時所添加的噪聲能量的若干分之一。為了模擬這種現(xiàn)象，CEEMDAN將ε0設置為初始噪聲和原始序列的理想信噪比SNR的倒數(shù)：若將SNR表達為標準差的商數(shù)，則有β0=ε0std(x)/std[E1(w()i)]，為了獲得后續(xù)分解階段中的具有較小波動幅度的噪聲實現(xiàn)，在剩余模態(tài)中，直接使用其前一步通過EMD進行分解時得到的噪聲，不用其標準差來進行歸一化處理，即βk=ε0std(rk)，k≥1。

所有EMD類方法的由細到粗進行篩分的理論本質(zhì)決定了其分解所得到的第1階模態(tài)所揭示的是原始序列中變化最快（頻率最高、周期最短）的序列分量。

2 實際應用

2.1 基本資料選用

本文分析選用唐乃亥水文站1956—2015年的年徑流量序列，如圖1所示。

圖1 唐乃亥站年徑流量序列

2.2 年徑流量的CEEMDAN分解

運用CEEMDAN方法[15]對圖1所示的唐乃亥站年徑流量序列進行多時間尺度分解，實現(xiàn)次數(shù)取200，限制標準差SD的取值為0.25，分解結(jié)果如圖2—5所示。

由圖2—5可知：

（1）CEEMDAN將唐乃亥站1956—2015年的年徑流量序列分解為4階模態(tài)，其中包括3個IMF分量（圖2—4）和1個趨勢項Res分量（圖5），反映了黃河源區(qū)產(chǎn)匯流系統(tǒng)變量演化的復雜多時間尺度性。

（2）第1階模態(tài)IMF1是振幅最大、周期最短、頻率最高的一個波動，依次下去的其他各階模態(tài)的振幅逐漸減小、周期逐漸變長、頻率逐漸降低。

（3）第1階模態(tài)IMF1具有準2～7 a波動周期，以準2～4 a波動周期為主，在60 a的觀測時段內(nèi)其平均振幅34.35億m3，最大振幅75.19億m3，最小振幅1.90億m3。

圖2 唐乃亥站年徑流量序列的IMF1分量

圖3 唐乃亥站年徑流量序列的IMF2分量

圖4 唐乃亥站年徑流量序列的IMF3分量

圖5 唐乃亥站年徑流量序列的Res分量

（4）第2階模態(tài)IMF2具有準5～10 a波動周期，以準7 a波動周期為主，在60 a的觀測時段內(nèi)其平均振幅36.79億m3，最大振幅57.27億m3，最小振幅2.64億m3；自1967—1982年為高幅振蕩時段，平均振幅51.86億m3，最大振幅57.27億m3，最小振幅44.37億m3；自1982—1998年為低頻振蕩時段，最大波動周期為準10 a。

（5）第3階模態(tài)IMF3具有準18 a和準28 a波動周期，在60 a的觀測時段內(nèi)其振幅呈增加趨勢，平均振幅26.40億m3，最大振幅33.59億m3，最小振幅18.74億m3。

（6）第4階模態(tài)Res分量顯示的是唐乃亥站年徑流量的整體變化趨勢，1977和2008年是Res曲線的2個拐點，1951—1977年唐乃亥站年徑流量整體呈增加趨勢，增幅為14.36%；1977—2008年唐乃亥站年徑流量整體呈減小趨勢，減幅為16.54%；自2008年開始唐乃亥站年徑流量整體又呈增加趨勢，進入新一輪的增長周期。按照目前所展現(xiàn)出來的半個波動周期（1977—2008年）推論，Res分量具有準62 a波動周期。

（7）20世紀90年代至21世紀初，黃河源區(qū)出現(xiàn)的、進而導致黃河流域出現(xiàn)的連續(xù)枯水時段主要是由于第3階模態(tài)IMF3的異常波動造成的，這是因為若按照90年代之前的波動周期為準18 a的波動規(guī)律，該模態(tài)應該在1993年出現(xiàn)谷值、2002年出現(xiàn)峰值，但是卻在1993年之后繼續(xù)下行直至在2002年出現(xiàn)谷值，形成一個完整的“逆峰”，又與Res分量的同期下降趨勢相疊加，導致黃河出現(xiàn)長達10 a左右的枯水期。

2.3 模態(tài)重構(gòu)精度評價

根據(jù)CEEMDAN的分解理論及實際計算，圖2—5所示的4階模態(tài)可以完全重構(gòu)圖1所示的實測徑流量序列。為了定量評價各階模態(tài)在重構(gòu)中的作用，采用納什效率系數(shù)（Nash-Sutcliffe efficiency co?efficient，簡稱NSE）來表征不同模態(tài)組合和實測徑流量序列之間的模擬精度，納什效率系數(shù)的定義為：

NSE的取值范圍為負無窮至1，E越接近1，表示模擬質(zhì)量越好，模型可信度越高；E等于1，表示模型模擬值與實測值完全一致，誤差為0；E越接近0，表示模擬結(jié)果越接近觀測值的平均值水平，即模型總體結(jié)果可信，但過程模擬誤差大；E遠遠小于0，則表明模型的模擬結(jié)果是不可信的。

以圖1所示徑流量序列作為實測值、圖2—5所示的4階模態(tài)的依次組合作為模擬值，分別計算其NSE值，結(jié)果如圖6—8所示。

圖6 唐乃亥站年徑流量實測值與2、3、4階模態(tài)重構(gòu)值的比較

圖7 唐乃亥站年徑流量實測值與3、4階模態(tài)重構(gòu)值的比較

圖8 唐乃亥站年徑流量實測值與模態(tài)4的比較

從圖6—8可知：

（1）第2、3、4階模態(tài)重構(gòu)所得的序列與實測徑流量序列之間的NSE為0.56，則第1階模態(tài)對重構(gòu)實測值的模擬精度貢獻率為0.44。

（2）第3、4階模態(tài)重構(gòu)所得的序列與實測徑流量序列之間的NSE為0.27，則第2階模態(tài)對重構(gòu)實測值的模擬精度貢獻率為0.29。

（3）第4階模態(tài)與實測徑流量序列之間的NSE為0.06，則第4階模態(tài)對重構(gòu)實測值的模擬精度貢獻率為0.06，同時也說明第4階模態(tài)作為實測徑流量序列的平均值是有統(tǒng)計意義的，符合NSE的定義，證明CEEMDAN分解理論和計算結(jié)果的正確性。

（4）周期越短、振幅越大、頻率越高的模態(tài)在重構(gòu)中的作用越是突出，對模擬精度的提高貢獻越大；但是，必須指出的是，這種貢獻的作用是相對的，是建立在作為平均值的第4階模態(tài)得以準確模擬的基礎(chǔ)之上的，若無第4階模態(tài)以及第3階模態(tài)作為基礎(chǔ)，第1階和第2階模態(tài)等高頻模態(tài)對原始序列的模擬精度將會出現(xiàn)數(shù)量級的誤差。由此可知，高頻模態(tài)是刻畫序列變化的細節(jié)和局部，低頻模態(tài)則掌控著序列變化的全局和趨勢。

3 結(jié)語

在分析EMD、EEMD和互補EEMD方法不足之處的基礎(chǔ)上，將CEEMDAN方法應用于黃河源區(qū)唐乃亥水文站1956—2015年的年徑流量時序的多時間尺度分析并采用納什效率系數(shù)定量評價了各階模態(tài)在模擬實測徑流信號的重構(gòu)中的精度貢獻率，結(jié)果表明唐乃亥站年徑流量分別具有準2～7 a、準5～10 a、準18 a、準28 a、準62 a的波動周期；黃河源區(qū)在20世紀90年代至21世紀初的連續(xù)枯水時段是由于第3階模態(tài)的異常波動造成的；在掌握序列趨勢變化的基礎(chǔ)上提高第1階模態(tài)的預測精度是提高唐乃亥站年徑流量預測精度的工作方向；按照各階模態(tài)的周期和振幅所顯示的變化趨勢，預計自21世紀10年代至中葉的這一時期內(nèi)，唐乃亥站年徑流量將呈現(xiàn)在波動中增加的趨勢，黃河源區(qū)將進入豐水期，這對于黃河流域水資源開發(fā)、利用與保護工作是有利的。

具適應性噪聲的完全集合經(jīng)驗模態(tài)分解（Com?plete Ensemble Empirical Mode Decomposition with Adaptive Noise，簡稱CEEMDAN）作為EMD和EEMD的最新改進方法，是具有分解完整性、模態(tài)精準性等優(yōu)點的一種全新的非線性、非平穩(wěn)序列多時間分辨率分析方法，在水文水資源領(lǐng)域的分析、建模和預測中具有廣闊的應用前景。

[1]Huang N E，Shen Z，Long S R，et al.The empirical mode de?composition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sci?ences，1998，454：903-995.

[2]馮平，丁志宏，韓瑞光，等.基于EMD的降雨徑流神經(jīng)網(wǎng)絡預測模型[J].系統(tǒng)工程理論與實踐，2009，29（1）：152-158.

[3]Zhang Jinping，Ding Zhihong，Yuan Wenlin，et al，Research on the relationship between rainfall and reference crop evapo?transpiration with multi-time scales[J].Paddy and Water En?vironment，2013，11（1-4）：473-482.

[4]趙文剛，邢旭光，馬孝義.基于EMD方法的土壤入滲空間異質(zhì)性及其影響因素研究[J].灌溉排水學報，2016，35（3）：61-67.

[5]Wu Z，Huang N E.Ensemble empirical mode decomposition: a noise-assisted data analysis method[J].Advances in Adap?tive Data Analysis，2009，1（1）:1-41.

[6]Flandrin P，Rilling G，Goncalvès P.Empirical mode decompo?sition as a filter bank[J].IEEE Signal Process，LETT，2004，11（2）：112-114.

[7]王兵，李曉東.基于EEMD分解的歐洲溫度序列的多尺度分析[J].北京大學學報（自然科學版），2011，47（4）：627-635.

[8]安學利，潘羅平，張飛.基于EEMD和近似熵的水電機組擺度去噪方法[J].水力發(fā)電學報，2015，34（4）：163-169.

[9]Ouyang Qi，Lu Wenxi，Xin Xin，et al.Monthly Rainfall Fore?castingUsingEEMD-SVRBasedonPhase-SpaceReconstruc?tion[J].WaterResourcesManagement，2016，30：2311-2325.

[10]Yeh J R，Shieh J S，Huang N E.Complementary ensemble empirical mode decomposition：a novel noise enhanced data analysis method[J].Advances in Adaptive Data Analysis，2010，2（2）：135-156.

[11]Torres M E，Colominas M A，Schlotthauer G，et al.A com?plete ensemble empirical mode decomposition with adaptive noise[A].Prague，Czech：IEEE International Conference on Speech and Signal Processing（ICASSP），2011：4144-4147.

[12]Antico Andrés，Torres María E，Diaz Henry F.Contributions of different time scales to extreme Paraná floods[J].Clim?mate Dynamics，2016，46：3785–3792.

[13]Han J，Bann M.Van der.empirical mode decomposition for seismic time-frequency analysis[J].Geophysics，2013，78（2）：9-19.

[14]韓慶陽，孫強，王曉東，等.CEEMDAN去噪在拉曼光譜中的應用研究[J].激光與光電子學進展，2015，52（11）：274-280.

[15]Colominas M A，Schlotthauer G.，Torres M.E.Improved com?plete ensemble EMD：A suitable for biomedical signal pro?cessing[J].Biomedical Signal Processing and Control，2014，14（11）：19-29.

TV121

A

1004-7328（2016）06-0001-06

2016—07—12

國家自然科學基金項目（51309202）；水利部公益性行業(yè)科研專項經(jīng)費項目（201101016）

丁志宏（1979—），男，博士，高級工程師，主要從事水文水資源研究與咨詢工作。