基于EEMD的鹽城大豐區(qū)降水量時間序列多尺度分析

陳 超，李鵬飛，陳尚云，陳 凱，陳則煌

(1.江蘇省鹽城市大豐區(qū)氣象局，江蘇鹽城 224100；2.中國氣象局氣溶膠與云降水重點開放實驗室，江蘇南京 210044；3.南京信息工程大學(xué)大氣物理學(xué)院，江蘇南京 210044)

?

基于EEMD的鹽城大豐區(qū)降水量時間序列多尺度分析

陳 超1，李鵬飛2，3，陳尚云1，陳 凱1，陳則煌2，3

(1.江蘇省鹽城市大豐區(qū)氣象局，江蘇鹽城 224100；2.中國氣象局氣溶膠與云降水重點開放實驗室，江蘇南京 210044；3.南京信息工程大學(xué)大氣物理學(xué)院，江蘇南京 210044)

摘要[目的] 研究鹽城大豐地區(qū)年降水量時間序列規(guī)律。 [方法] 基于集合經(jīng)驗?zāi)B(tài)分解(EEMD)的數(shù)據(jù)處理方法，通過Hilbert變換對降水量數(shù)據(jù)進行多尺度特征研究，并對研究結(jié)果進行顯著性和誤差檢驗。 [結(jié)果] 近55 a鹽城大豐地區(qū)年降水量序列可以分解成4個IMF分量和1個趨勢分量；主要存在0～0.05和0.45～0.50 Hz 2個頻率區(qū)間，其間主周期為2.12 a年際周期和27.50 a年代際周期，二者均通過95%顯著性檢驗。EEMD可以有效地重構(gòu)原始降水量數(shù)據(jù)變化曲線，與原始數(shù)據(jù)對比誤差百分比始終控制在±0.6%范圍內(nèi)，由白噪聲帶來的干擾誤差幾乎可以完全忽略。 [結(jié)論] 該研究可為合理開發(fā)與充分利用水文資源、防汛抗旱等提供科學(xué)依據(jù)。

關(guān)鍵詞年降水量；時間序列；EEMD；Hilbert變換；顯著性檢驗；鹽城大豐地區(qū)

近年來，全球溫度升高和氣候變化對水文氣象的影響受到各國越來越多的關(guān)注[1]，其中關(guān)于降水量多尺度變化的周期性規(guī)律研究已成為當(dāng)今學(xué)者們研究的熱點之一[2]。隨著統(tǒng)計學(xué)與計算數(shù)學(xué)等領(lǐng)域的飛速發(fā)展，基于多種研究方法和輔助計算工具對水文歷史累積數(shù)據(jù)深入地再分析有了很大的推進作用。吳亞楠等[3]利用滑動平均、小波分析、M-K突變檢驗等方法對第二松花江流域65 a間的降水量序列變化特征進行了系統(tǒng)的分析。李光強等[4]利用小波變換方法分析了近41 a來河北地區(qū)降水序列的周期變化，得出河北地區(qū)中東部、西部以及北部不同地區(qū)的周期性和差異性。以上主要采用小波分析法等[3-5]對降水量的周期性進行分析。然而小波分析依然沒脫離海森堡測不準原理的束縛，某種尺度下不能在時間和頻率上同時具有較高的精度，以及小波是非適應(yīng)性的導(dǎo)致周期計算不夠準確。HHT(Hilbert Huang Transform)技術(shù)，作為一種基于Hilbert變換的新的時域分析方法，由于具有較好的自適應(yīng)性、可以對非線性非平穩(wěn)的信號進行分析等優(yōu)點，被許多學(xué)者所推崇，已成為取代小波分析方法研究水文氣象的周期及變化的主流方法。然而，隨著對HHT-Huang的核心經(jīng)驗?zāi)Ｊ椒纸馑惴?Empirical Mode Decomposition,EMD)的深入研究，發(fā)現(xiàn)其存在一定的模態(tài)混疊問題，導(dǎo)致進行的Hlibert變換所得的瞬時頻率與Hlibert譜的準確度降低[6]。近年，Huang等[7]與法國的Flaridrin研究小組嘗試利用白噪聲來對信號中斷所缺失的尺度進行彌補，創(chuàng)造性地提出了基于集合經(jīng)驗?zāi)B(tài)分解(EEMD)改進HHT的全新方法，這很大程度上克服了模態(tài)混疊的問題，保留了原序列的信息。因此，EEMD方法勢必將在水文氣象領(lǐng)域中得到廣泛的運用，為合理開發(fā)與充分利用水文資源、防汛抗旱等提供科學(xué)依據(jù)。筆者基于改進型經(jīng)驗?zāi)＝?整體經(jīng)驗?zāi)B(tài)分解法(EEMD)，并以1960～2014年鹽城大豐地區(qū)年降水量時間序列為例進行應(yīng)用，探討近55 a鹽城大豐地區(qū)年降水量時間序列的變化規(guī)律與周期特征，并對所得結(jié)果進行置信檢驗和誤差分析。

1資料與方法

1.1資料選取選取鹽城大豐地區(qū)1960～2014年月降水量時間序列數(shù)據(jù)作為基礎(chǔ)研究資料，數(shù)據(jù)來源于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)。

1.2分析方法經(jīng)驗?zāi)B(tài)分解EMD是 HHT 變換的主要步驟，作為HHT變換進行分析的前提條件。先通過自適應(yīng)EMD求取IMF，然后利用Hilbert變換來求取各IMF的瞬時頻率以及Hilbert譜。然而由于現(xiàn)實的信號存在諸如脈沖、噪聲以及間斷信號等其他無關(guān)因素的干擾，影響在求取包絡(luò)過程中對于極值點選取的準確度，一旦所求取的包絡(luò)不是真實信號的極值點時，就會導(dǎo)致篩選得到的IMF不準確，因而產(chǎn)生模態(tài)混疊現(xiàn)象，使IMF失去真實的物理意義[8]。

鑒于受原始信號頻率因素影響和EMD自身算法的缺陷等原因產(chǎn)生的模態(tài)混疊問題，Huang等[7]曾提出使用中斷檢測法來處理該問題，但由于此方法需要人為檢驗而存在效率不高等缺點未被廣泛使用。譚善文等[9]提出利用多分辨率思想來實現(xiàn)EMD的分解，但這種方法是以EMD的自適應(yīng)為代價的。最后由Wu等[10]提出使用添加白噪聲的EEMD算法來進行IMF分量篩選，該算法有效地解決了EMD方法中存在的模態(tài)混疊等問題，簡單有效地改善了HHT變換的準確度問題。

EEMD實質(zhì)上是一種引用白噪聲對信號進行輔助處理的方法，即噪聲輔助信號處理。與降噪技術(shù)相反，通過將白噪聲加入原始信號，然后通過平滑技術(shù)去除脈沖干擾。對于時變信號x(t)，可以將x(t)表述為[7]：

X(t)=x(t)+ω(t)

(1)

式中，X(t)是加入白噪聲后的信號；ω(t)是符合正態(tài)分布的小幅度高斯白噪聲，服從下式：

(2)

式中,N為噪聲總數(shù)；ε為高斯白噪聲的幅值；εn為根據(jù)IMF重構(gòu)的信號與原始信號之間的誤差值。基于式(1)，利用EMD對X(t)進行分解，得到各階IMF分量，則X(t)可以表示為：

(3)

向X(t)加入符合(2)式的隨機白噪聲并重復(fù)上述過程，可得：

(4)

由于高斯白噪聲具有零均值特性，當(dāng)進行多次平均后，噪聲將相互抵消，并不會影響到原信號的變化情況，則原始信號所對應(yīng)的IMF分量imfn(t)為：

(5)

因此，x(t)可以分解為：

(6)

式中，imfn(t)為各個IMF分量；rm(t)為趨勢項，表征整個事件序列信號的變化趨勢。對各個IMF分量分別進行Hilbert變換，可得到其頻率、周期以及Hilbert譜H(w，t)。將H(w，t)對時間求積分，即可得Hilbert邊際譜h(w)[8]：

(7)

式中，w為瞬時頻率；t為時間；T為所研究序列的時間長度。

2結(jié)果與分析

圖1　1960～2014年鹽城大豐地區(qū)降水量時間序列變化圖Fig.1　Changes of precipitation time series in Dafeng District of Yancheng City from 1960 to 2014

2.1鹽城大豐地區(qū)降水變化特征由圖1可見，鹽城大豐地區(qū)最大年降水量為1 718.60mm，發(fā)生在1991年；最小年降水量為494.00mm，發(fā)生在1978年；55a年平均降水量為1 046.15mm。1960～2014年鹽城大豐地區(qū)相鄰年降水量變化波動較大，尤其是在1972、1987、1990年，相鄰的兩年之間降水量落差情況比較明顯。從總體趨勢來看，55 a間鹽城大豐地區(qū)年降水量在波動中增加。對該地區(qū)55 a年降水量監(jiān)測數(shù)據(jù)運用最小二乘法計算年降水量趨勢定量變化一次線性方程為：y= 0.099 8x+847.93，其中氣候傾向率為0.998 mm/10 a，表示鹽城大豐地區(qū)從20世紀60年代開始，年降水量總體上呈遞增趨勢，增幅為1 mm/10 a。

2.2EEMD在鹽城大豐地區(qū)降水時間序列值中的應(yīng)用使用EEMD方法對鹽城大豐地區(qū)降水量進行IMF分解時，在55 a年降水量值的序列中加入的高斯白噪聲序列的參數(shù)如下：標準差為0.01，100組。同時，通過平滑技術(shù)將導(dǎo)致模態(tài)混疊的奇異點進行濾除，而后將得到55 a鹽城大豐地區(qū)降水量時間序列分解后的真實固有模態(tài)函數(shù)IMFs。這里的IMFs可以分為4個固有模態(tài)函數(shù)和1個趨勢項，說明55 a鹽城大豐地區(qū)年降水量序列包含多個時間尺度變化特征(圖2)。根據(jù)圖2對分解得到的IMF各分量的變化情況進行分析，其中IMF1～IMF4的頻率值逐漸減小，平均波動幅值逐漸變小，波動周期依次變長。另外，趨勢分量r代表了整個時間序列上鹽城大豐地區(qū)年降水量的變化情況，即減少-增大-減小的勺形變化趨勢。而圖1展示的是年降水量最小二乘法處理后的結(jié)果，代表的是以整個時間序列作為一個整體的降水趨勢。圖2是對圖1更為細致的趨勢分解，可以理解為將整個時間序列進行多次劃分后展示的降水量趨勢。因此，二者并不沖突。

圖2　基于EEMD對1960～2014年鹽城大豐地區(qū)年降水量變化IMFs分解Fig.2　IMFs decomposition of precipitation changes in Dafeng District of Yancheng City based on EEMD from 1960 to 2014

根據(jù)IMF平均周期計算公式[6]，利用EEMD法分解后求得IMFs中心頻率與平均周期。由表1可知，IMF1～IMF4中心頻率分別是0.470、0.160、0.055和0.036 Hz，表明鹽城大豐地區(qū)年降水量的變化分別存在2.12、6.10 a的年際變化和18.30、27.50 a的年代際變化。

表1　各階IMF分量統(tǒng)計特征值

根據(jù)年降水量數(shù)據(jù)進行HHT-Huang變換后得到的55 a鹽城大豐地區(qū)年降水量的Hilbert譜圖(圖3)。由圖3可以看出，色標中紅色標記值最大，圖中紅色標記線主要存在的頻率范圍為0～0.05和0.45～0.50 Hz。雖然在頻率0.05～0.45 Hz也存在標記值較大的標記線，但其范圍分散，時間上跨度較大，因此不具有統(tǒng)計顯著性。由在HHT-Huang變換的基礎(chǔ)上頻率對時間積分后得到的結(jié)果(圖4)可見，信號幅值在整個0～0.50 Hz內(nèi)有分布，但幅值較大的主要分布在0～0.05和0.45～0.50 Hz，也表明在這2個區(qū)間中信號的能量最大，振蕩最為劇烈，是整個時間序列中頻率顯著性最強的2個范圍。圖4進一步驗證了圖3中的分析結(jié)果，相應(yīng)地，55 a鹽城大豐地區(qū)年降水量的變化周期主要存在于2個區(qū)間，即20.00 a以上的代際變化和2.00～2.20 a的年際變化。與55 a鹽城大豐地區(qū)年降水量變化趨勢分解圖中IMF分量的波動情況(圖2)對比，邊際譜和IMF分量變化均能給出序列的周期變化，但由于邊際譜更能反映各個頻率在整個周期上的能量幅值強度，進而更容易辨別整個頻率序列上主變化頻率所在的范圍。因此，綜合以上分析，55 a鹽城大豐地區(qū)年降水量的變化周期主要為2.12 a的年際周期和27.50 a的年代際周期。

圖3　1960～2014年鹽城大豐地區(qū)年降水量時間序列的Hilbert譜圖Fig.3　Hilbert spectrogram of precipitation time series in Dafeng District of Yancheng City from 1960 to 2014

圖4　1960～2014年鹽城大豐地區(qū)年降水量時間序列的邊際譜圖Fig.4　Marginal spectrogram of precipitation time series in Dafeng District of Yancheng City from 1960 to 2014

2.3EEMD降水時間序列值研究結(jié)果誤差分析EEMD是在EMD的基礎(chǔ)上發(fā)展起來的，它繼承了EMD的自適性，但為避免尺度混合問題，使最終分解的IMFs保持物理上的唯一性，由于從開始之初引入了白噪聲，那么對于分解到的IMF判斷顯得十分重要，即判斷信號分解后IMF分量是否真實有效。將基于EEMD變換后得到的IMF分量按照時間序列進行疊加，將疊加后的重構(gòu)信號進行繪圖得到鹽城大豐地區(qū)55 a年降水量的邊際譜圖。從圖5可看出，重構(gòu)信號與原始信號幾乎重疊。為進一步分析重構(gòu)信號與原始信號的失真度，將重構(gòu)信號與原始信號進行數(shù)學(xué)計算并引入誤差分析后進行繪圖，得到基于IMF重構(gòu)波形與原始波形的對比誤差變化(圖6)，在整個時間序列上，對比誤差百分比始終控制在±0.6%范圍以內(nèi)，處于10-3的量級范圍之內(nèi)，對原始信號的研究幾乎沒有影響。由此可見，基于EEMD的HHT法可以有效地重構(gòu)原始降水量時間序列值數(shù)據(jù)變化曲線。

圖5　基于IMF重構(gòu)波形與原始波形的對比Fig.5　Comparison between reconstruction waveform and original waveform based on IMF

2.4EEMD降水時間序列值研究結(jié)果顯著性檢驗分析從圖7可看出，E點處于95%與99%置信區(qū)間，B點較C、D點更接近95%的置信區(qū)間，可信水平約為95%。因此，綜上B點與E點的可信度最高，即IMF1分量和IMF4分量具有顯

圖6　基于IMF重構(gòu)波形與原始波形的對比誤差變化Fig.6　Changes of comparative error between reconstruction waveform and original waveform based on IMF

著性意義，這與根據(jù)IMF分量中心頻率得到的結(jié)果一致。因此 IMF1 與 IMF4分量顯著性較強，則其對應(yīng)的鹽城市大豐地區(qū)年降水量時間序列值的變化周期為主周期，具有統(tǒng)計意義。

注：A點為原始數(shù)據(jù)；B點為IMF1分量；C點為IMF2分量；D點為IMF3分量；E點為IMF4分量。Note：A was original data；B was MF1 component；C was IMF2 component；D was IMF3 component; E was IMF4 component. 圖7　各階IMF趨勢的顯著性檢驗Fig.7　Significance test of IMF trend in each order

3結(jié)論

該研究基于EEMD分解的HHT變換分析法對1960～2014年鹽城市大豐地區(qū)年降水量時間序列值變化進行了系統(tǒng)的研究，得出以下結(jié)論：

(1)1960～2014年鹽城大豐地區(qū)相鄰年降水量變化波動較大，從總體趨勢來看，55 a間鹽城大豐地區(qū)年降水量時間序列值在波動中增加，其中氣候傾向率為0.998 mm/10 a。

(2)通過EEMD可以將55 a鹽城大豐地區(qū)年降水量時間序列數(shù)據(jù)分解成4個IMF分量和1個趨勢分量。結(jié)合HHT變換邊際譜，55 a鹽城大豐地區(qū)年降水量時間序列值的變化主要有2個周期區(qū)間，即2.00～2.20 a的年際變化和20.00 a以上的代際變化，其間主周期主要為2.12 a的年際周期和27.50 a的年代際周期，并通過95%顯著性檢驗。

(3)EEMD可以有效地重構(gòu)鹽城大豐地區(qū)原始降水量時間序列數(shù)據(jù)的變化曲線，與原始數(shù)據(jù)的對比誤差百分比始終控制在±0.6%范圍以內(nèi)，由白噪聲帶來的干擾誤差幾乎可以完全忽略。

參考文獻

[1] CAI W J，SANTOSO A，WANG G J，et al.Increased frequency of extreme Indian Ocean Dipole events due to greenhouse warming[J].Nature，2014，510(7504)：254-258.

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

[3] 吳亞楠，李鴻雁，王玉新.第二松花江流域近75年降水量變化特征分析[J].安徽農(nóng)業(yè)科學(xué)，2012，40(4)：2257-2259.

[4] 李春強，杜毅光，李保國.1965-2005年河北省降水量變化的小波分析[J].地理科學(xué)進展，2010，29(11)：1340-1344.

[5] 辛渝，陳洪武，張廣興，等.新疆年降水量的時空變化特征[J].高原氣象，2008，27(5)：993-1003.

[6] 時世晨.EEMD時頻分析方法研究和仿真系統(tǒng)設(shè)計[D].上海：華東師范大學(xué)，2011：36-45.

[7] HUANG N E，SHEN Z，LONG S R.A new view of nonlinear water waves：The Hilbert spectrum[J].Annual review of fluid mechanics，2003，31(1)：417-457.

[8] 陳則煌，張云峰，謝菲，等.EEMD在雷暴日趨勢特征分析中的應(yīng)用[J].熱帶地理，2015(4)：601-606.

[9] 譚善文，秦樹人，湯寶平.Hilbert-Huang變換的濾波特性及其應(yīng)用[J].重慶大學(xué)學(xué)報(自然科學(xué)版)，2004，27(2)：9-12.

[10] WU Z，HUANG N E.Ensemble empirical mode decomposition：A noise-assisted data analysis method[J].Advances in adaptive data analysis，2011，1(1)：1-41.

Multiscale Analysis of the Precipitation Time Series in Dafeng District of Yancheng City Based on EEMD

CHEN Chao1, LI Peng-fei2,3, CHEN Shang-yun1et a

(1. Dafeng District Meteorological Bureau, Yancheng, Jiangsu 224100; 2. Key Laboratory for Aerosol-Cloud-Precipitation of China Meteorological Administration, Nanjing, Jiangsu 210044; 3. College of Atmospheric Physics, Nanjing University of Information Science & Technology, Nanjing, Jiangsu 210044)

Abstract[Objective] To research the precipitation time series in Dafeng District of Yancheng City. [Method] Based on the data treatment method of ensemble empirical mode decomposition (EEMD), multi-scale characteristics research of precipitation data were researched by Hilbert transform. Significance and error tests on research results were carried out. [Result] The sequence of annual precipitation of Dafeng District in recent 55 years could be broken down into four IMF (Intrinsic Mode Function) components and a trend component. There were mainly two frequency regions, which were 0 - 0.05 Hz and 0.45- 0.50 Hz. The annual cycle was 2.12 years and the decadal cycle was 27.50 years during the main period. Both of them passed the 95% significance test. EEMD precipitation could effectively reconstruct the curve of original data. Compared with the original data, error percentage was always controlled within ± 0.6%. Then the interference error by the white noise could be almost completely ignored. [Conclusion] This research provides scientific basis for the rational development and fully utilization of hydrologic resources, flood control and drought relief.

Key wordsAnnual precipitation; Time series; EEMD; Hilbert transform; Significance test; Dafeng District in Yancheng City

基金項目國家自然科學(xué)基金項目(41175003)；江蘇高校優(yōu)勢學(xué)科建設(shè)工程資助項目。

作者簡介陳超(1987- )，男，江蘇鹽城人，工程師，在讀碩士，從事雷電過電壓防護技術(shù)與氣象數(shù)據(jù)統(tǒng)計研究。

收稿日期2016-03-08

中圖分類號S 161.6；P 448

文獻標識碼A

文章編號0517-6611(2016)10-209-03