灤河上游年降水量多時間尺度變化的EEMD分析

張小鵬

灤河上游年降水量多時間尺度變化的EEMD分析

張小鵬

（山西水利職業(yè)技術(shù)學(xué)院，山西運(yùn)城044004）

以多倫氣象站1956—2014年的年降水量序列為代表，運(yùn)用EEMD方法分析了灤河上游年降水量的多時間尺度變化特性。結(jié)果表明，灤河上游年降水量具有準(zhǔn)2～4年、準(zhǔn)4～6年、準(zhǔn)10年、準(zhǔn)15～20年波動周期，整體變化呈衰減趨勢。

灤河上游；多倫；年降水量；多時間尺度；EEMD；衰減；最嚴(yán)格水資源管理制度

灤河是海河流域的一條重要水系，發(fā)源于河北省豐寧縣大灘鎮(zhèn)，經(jīng)沽源縣向北流入內(nèi)蒙古多倫縣境，至外溝門子又進(jìn)入河北省境內(nèi)，蜿蜒于峽谷之間，到潘家口越長城，經(jīng)灤縣進(jìn)入平原，于樂亭縣境內(nèi)注入渤海，全長888 km，流域面積4.47萬km2。自河源至張百灣為上游，自張百灣至灤縣為中游，灤縣至入海口為下游。

氣象系統(tǒng)是復(fù)雜的非線性動力學(xué)系統(tǒng)，作為其重要輸出變量，降水量的年際變化存在著多時間尺度性，即年降水量序列在某一時間段內(nèi)不是只以一種固定的頻率（時間尺度、周期）在運(yùn)動，而是同時包含著各種頻率（時間尺度、周期）的變化和局部波動［1］，是多種動力機(jī)制同時在發(fā)揮作用，使得氣象系統(tǒng)變化在時域中存在著多層次的時間尺度和局部化特征。作為水文水資源系統(tǒng)中的最活躍的基本輸入變量，降水量的變化對于流域水資源量的形成和分布具有重要影響。

筆者采用EEMD方法分析灤河上游的多倫氣象站近60年來年降水量的多時間尺度變化特性，以期為評價灤河流域水資源情勢、制定水量分配方案以及貫徹落實(shí)最嚴(yán)格水資源管理制度等工作提供科學(xué)參考。

1　數(shù)據(jù)和方法

1.1數(shù)據(jù)

本分析采用位于灤河上游的多倫縣氣象站1956—2014年的年降水量序列，如圖1所示。

圖1　多倫站年降水量序列

1.2方法

就其實(shí)質(zhì)而言，小波變換是一種窗口可調(diào)的傅立葉變換，它要求在小波窗內(nèi)所分析的信號必須是平穩(wěn)的，小波變換會造成很多虛假的諧波，基函數(shù)的選擇對小波分解結(jié)果有顯著影響［2］。為此，Huang等人提出了經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）方法［3］，該方法是基于信號的局部特征時間尺度從原信號中提取一系列的幅度和頻率都經(jīng)過調(diào)制的函數(shù)，這些函數(shù)稱之為本征模態(tài)函數(shù)（Intrinsic Mode Function，IMF），原始信號可以用本征模態(tài)函數(shù)之和來進(jìn)行還原和表達(dá)，各IMF分量包含了原始信號的不同時間尺度局部特征信息，具有明顯的物理背景，而其中最低頻率的IMF分量通常代表原始信號的趨勢或均值。但是，EMD方法在分解過程中容易產(chǎn)生混頻即模態(tài)混疊或稱尺度混合現(xiàn)象，為此Huang等人又對EMD方法進(jìn)行了改進(jìn)，提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解（Ensemble Empirical Mode Decom-position，EEMD）方法［4］，具體是通過向原信號中多次添加不同的白噪聲之后再分別進(jìn)行EMD分解，然后再對經(jīng)過多次EMD分解所得的各個IMF分量分別求平均值而得到最終的實(shí)際分量。EEMD方法通過多次集合平均來抵消白噪聲的影響，可以有效地改善EMD方法所存在的模態(tài)混疊現(xiàn)象。

EEMD方法的具體計(jì)算步驟為：首先將給定振幅的白噪聲序列疊加在待分解的數(shù)據(jù)序列上，形成混合序列，然后再對此混合序列進(jìn)行EMD分解。如此反復(fù)，每次加入振幅相同的新的白噪聲序列都會得到不同的IMF分量，最后，將各次分解得到的IMF分量進(jìn)行集合平均以作為相應(yīng)IMF分量的最終分解結(jié)果。

EMD方法中的本征模態(tài)函數(shù)（IMF）要滿足以下2個條件：①在整個數(shù)據(jù)范圍內(nèi)，過零點(diǎn)和極值點(diǎn)的數(shù)量必須保持相等或至多相差1；②在任何點(diǎn)處，所有極大值點(diǎn)形成的上包絡(luò)線和所有極小值點(diǎn)形成的下包絡(luò)線的平均值始終保持為0。某一信號可以進(jìn)行EMD分解的前提為：①被分解的信號至少存在2個極值點(diǎn)：1個極大值點(diǎn)和1個極小值點(diǎn)；②局部特征時間尺度可定義為信號中兩臨近極大值點(diǎn)或極小值點(diǎn)的時間間隔；③若信號中不存在極值點(diǎn)，但包含若干拐點(diǎn)，可以先對信號進(jìn)行若干次微分，使得極值點(diǎn)顯露出來后，再對分解得到的分量進(jìn)行積分來求得最后結(jié)果。

EMD分解的基本思想是：若加入待分解數(shù)據(jù)序列的極小值或極大值數(shù)目比下跨零點(diǎn)（或上跨零點(diǎn)）的數(shù)目多2個或2個以上，則該數(shù)據(jù)序列需要進(jìn)行平穩(wěn)化處理。平穩(wěn)化處理時，首先，利用三次樣條函數(shù)把序列x（t）的局部極小值點(diǎn)和局部極大值點(diǎn)分別擬合成x（t）的下包絡(luò)線和上包絡(luò)線，然后再計(jì)算2條包絡(luò)線的平均值m1。再從原始數(shù)據(jù)序列x（t）中減去m1，即可得到一個移除了低頻信號的新的數(shù)據(jù)序列：

通常，h1并不是IMF分量，為此尚需對h1重復(fù)以上處理過程以進(jìn)行k次篩選直至所得到的包絡(luò)平均值趨于零為止，此時所得數(shù)據(jù)為：

式中：h1k為第k次篩選所得的數(shù)據(jù)；h1（）k-1為第k-1次篩選所得的數(shù)據(jù)。可使用限制標(biāo)準(zhǔn)差SD的值來判斷每次所得的篩選結(jié)果是否為IMF分量，SD定義為：

式中：T為數(shù)據(jù)序列長度；其他變量含義同前。

一般地，SD值取0.2～0.3，即滿足0.2〈SD〈0.3時，EMD分解過程即可結(jié)束，這樣既使得hk（）t足夠接近IMF的要求，又可以控制分解的次數(shù)，從而使所得IMF分量保留原始信號中幅值和頻率的調(diào)制信息［5］。

當(dāng)h1k滿足SD的要求時，令c1=h1k，即可得到信號x（t）的第一個IMF分量，它代表了原始信號序列中的頻率最高的組成成分。從原始數(shù)據(jù)序列x（t）中減去第一個IMF分量c1，就得到一個移除了高頻組分的差值數(shù)據(jù)序列：r1=x（t）-c1。若r1中仍包含x（t）的較長周期的局部特征時間尺度信息，可將r1作為待分解信號，再重復(fù)式（1）—（3）的過程，直到所剩信號r1中的信息對所研究目的而言意義已很小或者已是單調(diào)函數(shù)時即可停止分解運(yùn)算，此時的rn就代表著原始數(shù)據(jù)序列的趨勢或均值。至此，便得到了信號x（t）的一系列IMF分量：c1，c2，…，cn，且r1-c2=r2，r2-c3=r3，…，rn-1-cn=rn。原始數(shù)據(jù)序列即可由這些IMF分量以及1個均值或趨勢項(xiàng)表示：

EEMD方法是在每次開始初始分解前，在原始信號序列中添加具有一定幅值的白噪聲，然后再應(yīng)用EMD方法進(jìn)行分解，如此重復(fù)多次，直至分解完成，再取各IMF分量的各次分解值進(jìn)行求和以得到最終分解結(jié)果。EEMD方法在繼承EMD方法的自適應(yīng)分解特征的同時，通過引入白噪聲再進(jìn)行集合平均，使得最終分解得到的IMF分量保持了物理意義上的唯一性。

2　結(jié)果

運(yùn)用EEMD方法對圖1所示的多倫氣象站1956—2014年的年降水量序列進(jìn)行多時間尺度分解，擾動白噪聲與原始序列的信噪比取0.2，集合的樣本數(shù)取100，分解結(jié)果如圖2—6所示。

從中，可知以下結(jié)論：

（1）多倫站的年降水量序列可以分解為4個具有不同周期的波動分量和1個趨勢分量，反映了區(qū)域氣候系統(tǒng)變量變化所具有的復(fù)雜時域性。

（2）IMF1分量具有準(zhǔn)2～4年波動周期，其波動幅度在近60年來沒有明顯的趨勢變化。

（3）IMF2分量具有準(zhǔn)4～6年波動周期，其波動幅度在20世紀(jì)90年代之后較20世紀(jì)50—80年代為小。

（4）IMF3分量具有準(zhǔn)10年波動周期，其波動幅度在20世紀(jì)90年代之后較20世紀(jì)50—80年代為大。

（5）IMF4分量具有準(zhǔn)15～20年波動周期，其波動幅度在近60年來呈增加趨勢。

（6）Res分量顯示的是年降水量的整體變化趨勢，就整體而言，多倫站年降水量在近60年來呈衰減趨勢，降幅為15.63％。值得指出的是，多倫站年降水量序列的變化趨勢項(xiàng)可能屬于更長周期（更小頻率）波動的組成部分，而限于觀測時限，Res分量的波動周期和振幅目前還難以準(zhǔn)確測知，有待隨著觀測時限的延長逐步得以展現(xiàn)。

圖2　多倫站年降水量序列的IMF1分量

圖3　多倫站年降水量序列的IMF2分量

圖4　多倫站年降水量序列的IMF3分量

圖5　多倫站年降水量序列的IMF4分量

圖6　多倫站年降水量序列的Res分量

3　結(jié)語

筆者以多倫氣象站1956—2014年的年降水量序列為基礎(chǔ)，運(yùn)用EEMD分析了灤河上游流域年降水量的多時間尺度變化特性，揭示了其在不同時間尺度上所具有的不同周期和不同幅度的波動演化趨勢，年降水量值變化在整體上呈衰減趨勢。作為灤河的發(fā)源地和主要產(chǎn)流區(qū)之一，天然降水量衰減使得灤河上游流域產(chǎn)匯流系統(tǒng)的主要輸入項(xiàng)量值受到較大影響，加之人類活動對下墊面狀況的改變，使得河源區(qū)的地表與地下水資源量衰減更甚。在未來一段時期內(nèi)，隨著上游流域經(jīng)濟(jì)社會快速發(fā)展，工農(nóng)業(yè)生產(chǎn)和居民生活需水量將呈增加趨勢，流域水資源供需矛盾將逐漸突出。為此，需要切實(shí)加強(qiáng)節(jié)約用水工作、開展灤河上游冀蒙兩省省際水量分配、認(rèn)真貫徹落實(shí)最嚴(yán)格水資源管理制度，努力實(shí)現(xiàn)人水和諧的可持續(xù)發(fā)展。

［1］張少文，丁晶，廖杰，等.基于小波的黃河上游天然年徑流變化特性分析［J］.四川大學(xué)學(xué)報（工程科學(xué)版），2004，36 （3）：32-37.

［2］Tewfiki A H.On the optimal choice of a wavelet for signal representation［J］.IEEE Trans Information Theory，1992，38 （2）：747-765.

［3］Norden E H，Shen Z，Long S R，et a1.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sciences，1998，454：899-955.

［4］WU Z，HUANG N E.Ensemble empirical mode decomposition：a noise-assisted data analysis method［J］.Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sciences，1998，454：899-955.

［5］馮平，丁志宏，韓瑞光.基于EMD的洮河年徑流量變化多時間尺度分析［J］.干旱區(qū)資源與環(huán)境，2008，22 （12）：73-76.

TV125

A

1004-7328（2016）04-0037-03

10.3969/j.issn.1004-7328.2016.04.012

2016—03—18

張小鵬（1979—），男，助教，主要從事水利水電工程教學(xué)與研究工作。