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

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

?

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

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

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

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

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

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

1資料與方法

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

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

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

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

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

式中，w為瞬時(shí)頻率；t為時(shí)間；T為所研究序列的時(shí)間長(zhǎng)度。

2結(jié)果與分析

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

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

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

圖2　基于EEMD對(duì)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平均周期計(jì)算公式[6]，利用EEMD法分解后求得IMFs中心頻率與平均周期。由表1可知，IMF1～I(xiàn)MF4中心頻率分別是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ì)特征值

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

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

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

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

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

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

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

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

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

3結(jié)論

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

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

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

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

參考文獻(xiàn)

[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é)報(bào)(自然科學(xué)版)，2011，47(4)：627-635.

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

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

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

[6] 時(shí)世晨.EEMD時(shí)頻分析方法研究和仿真系統(tǒng)設(shè)計(jì)[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] 陳則煌，張?jiān)品?，謝菲，等.EEMD在雷暴日趨勢(shì)特征分析中的應(yīng)用[J].熱帶地理，2015(4)：601-606.

[9] 譚善文，秦樹(shù)人，湯寶平.Hilbert-Huang變換的濾波特性及其應(yīng)用[J].重慶大學(xué)學(xué)報(bào)(自然科學(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

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

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

收稿日期2016-03-08

中圖分類(lèi)號(hào)S 161.6；P 448

文獻(xiàn)標(biāo)識(shí)碼A

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