張恒璟 程鵬飛

(1)遼寧工程技術(shù)大學測繪與地理科學學院，阜新 123000 2)中國測繪科學研究院，北京100830)

基于經(jīng)驗模式分解的CORS站高程時間序列分析*

張恒璟1，2)程鵬飛2)

(1)遼寧工程技術(shù)大學測繪與地理科學學院，阜新 123000 2)中國測繪科學研究院，北京100830)

提出基于經(jīng)驗模式分解技術(shù)的GPS高程時間序列時頻分析方法，并以兩個CORS基準站多年連續(xù)的高程時間序列數(shù)據(jù)為對象，獲取序列信號的各個本征模態(tài)函數(shù)分量，研究了信號周期運動的主要貢獻分量和序列趨勢項的合成與分離的方法。結(jié)果表明:在一定的序列區(qū)間長度上識別周期項與趨勢項才有意義;經(jīng)驗模式分解技術(shù)分離的序列趨勢項呈現(xiàn)大曲率的非線性形式，與傳統(tǒng)的周期函數(shù)高程序列擬合線性趨勢項存在較大的差異。

經(jīng)驗模式分解;本征模態(tài)函數(shù);趨勢項;GPS高程時間序列;方差貢獻率

1 引言

正確識別并分離序列趨勢項對研究GPS高程序列的頻譜特性和噪聲性質(zhì)具有十分重要的意義。傳統(tǒng)時頻信號分析前，須先識別并去除序列的趨勢項，例如基于傅里葉分析理論的功率譜估計和小波分解方法。經(jīng)典方法研究GPS高程時序信號的周期和噪聲特性時，事先給定一個具體的周期擬合模型［1，2］，并指定序列存在半年、年、兩年等的周期運動和長期的線性趨勢運動，通過函數(shù)擬合模型，求解各項的振幅相位等特性;接著在分析序列噪聲時，從原始序列中去除趨勢項和周期項后的殘差序列中，分析信號存在的各種噪聲特性，但實際上并不是每個GPS高程序列都具有事先設定的多種周期運動特性。因此有必要采用自適應的或者盡量不依賴事先給定特性的序列信號分解方法來識別序列的各種周期運動和趨勢項信息。本文引入處理非線性非平穩(wěn)信號的希爾伯特-黃變換方法［3］，采用經(jīng)驗模式分解(Empirical Mode Decomposition，EMD)和整體經(jīng)驗模式分解(Ensemble EMD，EEMD)，對我國CORS基準站高程序列數(shù)據(jù)進行分析，獲取了本征模態(tài)函數(shù)分量(Intrinsic Mode Function，IMF)，分析了各個IMF的方差貢獻率及對序列運動能量的貢獻大小，并研究了合并和提取序列高程運動趨勢項的方法。

2 經(jīng)驗模式分解技術(shù)

2.1 EMD分解原理

經(jīng)驗模式分解是基于非線性、非平穩(wěn)數(shù)據(jù)的自適應時頻信號分解方法，通過不斷地“篩選”，將原始時間序列信號分解為一系列頻率由高到低的本征模態(tài)函數(shù)分量和殘差項，分解產(chǎn)生的IMF函數(shù)需滿足兩個條件:一是待分解信號中極值點的個數(shù)與過零點的個數(shù)相等或至多相差1;二是在任意點上，由局部極大值定義的上包絡線與局部極小值定義的下包絡線的均值為零。分解的具體過程為［4，5］:

1)搜索原始信號所有的極大值點和極小值點，用三次樣條函數(shù)擬合所有極大值點得到上包絡線，同理得到下包絡線，取上下包絡的均值作為原始序列的平均包絡m1(t);從原始序列中扣除該平均包絡得到新的序列信號h1(t)=X(t)-m1(t)。如果h1(t)滿足IMF的條件，則作為第一個分解得到的本征模態(tài)函數(shù)分量，否則以h1(t)作為信號，繼續(xù)尋找其上下包絡線并求出平均包絡m2(t)，得到新信號h2(t)=h1(t)-m2(t)，迭代k次后直至hk(t)滿足IMF的條件，將其作為第一個IMF分量IMF1(t)。

2)計算殘余分量r1=X(t)-IMF1(t)，將r1作為信號，重復第一步的分解，得到第二個IMF分量IMF2(t)。

3)重復1)、2)步，直至殘余分量滿足分解中止的條件:即殘余分量為單調(diào)變化的殘差序列不能再分解出本征模態(tài)函數(shù)為止。原始信號經(jīng)過n次分解后，可描寫為:

這種分解基于數(shù)據(jù)的局部特征，與傅立葉頻譜分析分解成恒定振幅與頻率的一系列正弦函數(shù)和小波分解需預先給定基函數(shù)相比，EMD是一種自適應的分解方法。

2.2 整體經(jīng)驗模式分解

EMD分解的缺點之一是模式混疊現(xiàn)象:一個IMF分量包含著尺度互異的信號，或者尺度相似的信號存在于不同的IMF分量，這種現(xiàn)象主要是由于信號中斷引起的，不僅混淆了信號的時頻分布特性，更重要的使信號IMF分量的物理意義變得不清晰。為了減弱模式混疊現(xiàn)象的影響，文獻［6］在EMD分解基礎(chǔ)上提出了整體經(jīng)驗模式分解方法(EEMD)，一種噪聲輔助的數(shù)據(jù)分析方法(noise-assisted data analysis，NADA)，基本思路是在分解之前，將原始序列信號加入白噪聲，生成待分解信號，即

下標i表示第i次EMD分解。

白噪聲與原始信號的標準差之比例如可以設定為0.1，經(jīng)過EMD分解后生成一系列的IMF分量;由于每次隨機生成的白噪聲都是不同的，這樣可以通過重復上面的分解過程，得到信號的多次EMD分解的IMF分量，將對應分解尺度相同的IMF分量取平均，即得到了最終信號分解的IMF分量。文獻［6］建議整體平均的次數(shù)取100，而噪聲與信號的標準差之比可以取0.1、0.2、0.4等。本文最大分解層數(shù)n設定為:

其中N是序列長度，實驗采用的序列從1999—2010，共約4 000個數(shù)據(jù)點，因此得到的最大分解層數(shù)(IMF分量的個數(shù))不超過10。這樣就存在序列過度分解或分解不完全兩種情形，從實驗的幾個CORS站高程序列分解效果來看，剩余的殘差項均是單調(diào)的趨勢，符合EMD分解要求，不存在分解不完全的情形，而過度分解的情形需要對鄰近殘差項的IMF分量與殘差項合并，來判斷合并后的趨勢項是否仍然單調(diào)變化。

2.3 序列趨勢項的定義

傳統(tǒng)序列趨勢項是基于固定的數(shù)學模型定義，比如GPS時間序列常用的周期擬合模型為:

并由式(4)擬合出序列的常數(shù)項、線性速率項、周期項(年和半年周期)。式中，ti(i=1，…，N)是以年為單位的時間，以每日解構(gòu)成坐標時間序列，待求系數(shù)a和b為序列線性趨勢項，c、d和e、f分別為年周期和半年周期項的系數(shù)。

更復雜的情形，在上述模型中考慮各種原因引起的階躍式坐標突變(例如遠場大地震引起的同震位移、由于儀器或天線變更引起的位移、甚至由于某種不清楚原因引起的點位變化等)修正項，這是目前處理GPS時間序列常用的擬合模型。該模型扣除了線性速度項和常數(shù)項后的序列，是下一步時序信號頻譜和噪聲分析的基礎(chǔ)，比如基于傳統(tǒng)傅里葉理論的功率譜分析技術(shù)和目前常用的小波分解技術(shù)等。

文獻［6］指出，趨勢應該是序列信號的內(nèi)在特性之一，是序列的固有成分，內(nèi)在特性要求趨勢項的定義應該具有自適應性，應基于序列數(shù)據(jù)本身去分離，而不僅僅是用一個指定的函數(shù)模型去擬合和扣除。同時為區(qū)分趨勢項和長周期項，趨勢項的定義必須被限定在給定的序列尺度上，并且至多包含一個極值點?；谏厦娴目紤]，文獻［7］給出了趨勢項的定義，即:趨勢項是一個內(nèi)在的單調(diào)擬合函數(shù)，或在一個給定的序列數(shù)據(jù)區(qū)間上至多有一個極值點的函數(shù)，“給定的序列數(shù)據(jù)區(qū)間”可以是整個序列，也可以是序列數(shù)據(jù)的一部分。

3 實驗分析

3.1 EMD/EEMD分解IMF分量

對GPS高程序列數(shù)據(jù)分別采用EMD和EEMD兩種分解方法，獲取IMF分量。EEMD方法又設計了3種方案，整體平均的次數(shù)選用了文獻［6］的推薦值100，加入的白噪聲與原始序列信號的標準差之比分別取0.1、0.2、0.4，下面以EEMD01、EEMD02、EEMD04分別代表這3種EEMD分解方案。統(tǒng)計任一IMF分量的方差公式為:=E2(X)-(E (X))2，若X為有限長度的數(shù)據(jù)序列，則可統(tǒng)計出方差的估值。方差貢獻率是各個IMF分量的方差與分解得到的IMF分量方差之和的百分比。IMF分量方差貢獻率的大小反映了該頻率的分解信號在整個序列信號運動能量中的貢獻大小。

圖1(a)、(b)分別是BJFS(北京房山)和URUM (烏魯木齊)站采用 EMD/EEMD01/EEMD02/EEMD04 4種分解方案的IMF分量和殘差，限于篇幅未給出高頻的IMF1～5曲線。由公式(3)可知分解得到的IMF分量個數(shù)不超過10，從IMF1～10本征模態(tài)分量信號頻率依次降低，周期依次增大。

圖2是BJFS和URUM兩個基準站采用EMD/ EEMD01/EEMD02/EEMD04 4種方案統(tǒng)計得到的各個IMF分量中誤差和方差貢獻率。

BJFS站分解后的IMF1～3是高頻分量，IMF4和IMF5是近似月周期和雙月周期信號，但信號的振幅和方差較小，在整個序列信號中所占的比重較小。IMF6、7兩個信號分量的方差貢獻率之和達到50%以上，是序列的主要周期運動貢獻項，即BJFS站的周期運動以季節(jié)性變化和年周期變化為主。兩種分解方法得到的IMF6分量整體呈現(xiàn)季節(jié)性的周期變化，IMF7分量整體呈現(xiàn)明顯的年周期變化;但EMD分解與3種EEMD分解方案的方差貢獻率出現(xiàn)較大差異，EMD分解的IMF6、7方差貢獻率分別是25%與28%，EEMD分解分別是12%與40%左右，兩種技術(shù)得到的IMF6與IMF7方差貢獻率之和幾乎保持不變，其他本征模態(tài)分量的方差貢獻率基本相等。從圖1(a)明顯看出，EMD分解得到的IMF6季節(jié)性信號在橫軸2003—2006區(qū)間上出現(xiàn)了明顯的年周期信號，EMD分解得到的IMF7年周期信號在橫軸2003—2006區(qū)間上出現(xiàn)了明顯的兩年周期信號，說明EMD技術(shù)分解時部分IMF6信號被分解到IMF7分量上，這是由于模式混疊現(xiàn)象引起的，在圖1(a)中以橢圓部分表示。IMF8分量對應著明顯的兩周年運動，但其方差貢獻率和信號振幅相對較小。IMF9和IMF10分量振幅在5 mm之內(nèi)，是信號分解的長周期還是屬于過度分解的殘差信號，在趨勢項合并部分做進一步的分析。本文也設計了采用相同的EEMD方案對BJFS站多次分解的試驗，結(jié)果表明每次試驗之間低頻的本征模態(tài)分量IMF6、7的方差貢獻率出現(xiàn)了0.5～1.9%的微小變化，但是EMD的分解結(jié)果唯一，分析原因可能是EEMD分解時每次加入的白噪聲變化引起的，這也說明了EEMD技術(shù)仍存在著分解不唯一的缺點，需進一步的研究和完善EEMD分解算法。

圖1 IMF6～10分量與殘差項Fig.1 IMF6-10 components and residuals

圖2 IMF分量中誤差與方差貢獻率Fig.2 RMS and contribution rate of variance of each IMF

URUM站3種EEMD分解結(jié)果一致性很好，與EMD分解結(jié)果相比在中低頻IMF分量上有一定的差異。IMF1～3是分解得到的高頻信號，這3個分量方差貢獻率之和在20%左右，相對于BJFS站這3個高頻分量的方差貢獻率之和達到35%，GPS時間序列中高頻信號一般被認為是序列所含的噪聲，從這個角度上可認為URUM站高程序列信號中噪聲含量明顯低于BJFS站。IMF4分量周期大約是1個月，IMF5分量周期大約是兩個月。兩種分解技術(shù)得到的IMF6分量整體上呈現(xiàn)近似季節(jié)性和半年周期變化的信號，但在橫軸2002—2005區(qū)間上產(chǎn)生了明顯的年周期信號。兩種分解技術(shù)得到的IMF7分量均呈現(xiàn)明顯的年周期變化，但 EMD分解方法在IMF7序列的兩端產(chǎn)生了兩周年的信號。兩種分解技術(shù)得到的 IMF8整體上是近似2周年信號，但EMD分解在最右側(cè)產(chǎn)生了周期超過3年的分解信號。EMD分解在IMF6～8 3個分量均出現(xiàn)了模式混疊現(xiàn)象，雖然EEMD技術(shù)可以減弱模式混疊的影響，但從URUM站IMF6分量的分解結(jié)果可知，EEMD仍會出現(xiàn)頻譜混疊，需要對EEMD分解技術(shù)進一步的研究。從方差貢獻率指標定量分析:IMF6～8 3個分量EMD分解的方差貢獻率分別是32%、35%、4%，EEMD01和EEMD02方案分解結(jié)果幾乎完全相同方差貢獻率分別是21%、20%、21%，兩種分解技術(shù)得到的3個分量方差貢獻率之和在62～71%之間，說明URUM站的高程運動以季節(jié)性、年、兩年為主要周期。

3.2 趨勢項合成

圖3、圖4分別為BJFS和URUM站的趨勢項合成。

BJFS站IMF9、10分量的中誤差在0.8 mm左右，方差貢獻率為1.1～1.7%，這兩個本征模態(tài)分量在整個序列運動中所占的比重很小。從趨勢項合成圖3可知，最后的兩個本征模態(tài)分量IMF9、10一起與分解的殘差項合并后，曲線幾乎保持原有的單一變化趨勢，說明IMF9、10分量屬于過度分解，也有可能是以目前序列數(shù)據(jù)的尺度還無法識別的更長周期項，因此在研究的序列區(qū)間尺度上IMF9、10分量應歸于分解的殘差，一并作為序列趨勢項處理。

URUM站4種不同的分解方案對應的本征模態(tài)分量IMF10中誤差均在0.5mm左右，方差貢獻率在0.1～0.6%之間，在整個序列運動中所占比重幾乎可忽略不計。從趨勢項合成圖4可知，最后一個模態(tài)分量IMF10應歸于過度分解產(chǎn)生的分量，與殘差合并一起作為趨勢項處理。加入IMF9分量之后的趨勢項呈現(xiàn)明顯的周期振蕩，周期近似為序列長度的一半5.5年，加入IMF8合成后的趨勢項周期更小且振幅顯著增大，說明IMF8、9分量不能作為序列的趨勢項，而應當歸于序列運動較明顯的周期項分量，只是其振幅較小。

4 討論和總結(jié)

1)采用EMD/EEMD方法，獲取頻率從高到低的IMF分量，可明顯分解出GPS高程序列信號的趨勢項和周期項。從兩個基準站例子可知，分解后的低頻IMF分量為高程時間序列信號的周期運動給出了清晰的解釋，包含了明顯的1個月、2個月、季節(jié)性、年周期變化、兩周年周期變化與長周期變化因素，半周年項運動信號并不顯著，季節(jié)性和年周期變化是高程運動的主要貢獻項;EEMD分解相比EMD可以明顯減弱模式混疊對本征模態(tài)分量的影響，但并沒有完全消除模式混疊現(xiàn)象，同時受EEMD分解時加入白噪聲的影響，其分解結(jié)果也不唯一。EEMD比EMD方法有較大的優(yōu)點，但還需要進一步改進EEMD算法。

2)通過比較各個IMF分量的方差貢獻率及其中誤差大小，可以判斷對原始信號是否存在著過度分解的情形，進而合并過度分解的IMF分量與殘差項，獲得準確的序列趨勢項。合并后的趨勢項呈現(xiàn)大曲率的曲線變化，而非線性形式，傳統(tǒng)對序列趨勢項的線性擬合方法存在不足，序列趨勢項不僅包括了線性變化，也有暫時無法分離或識別的長周期信號，只有在一定的序列區(qū)間長度上探討并區(qū)分周期項和趨勢項才有意義。

圖3 BJFS站趨勢項合成Fig.3 Trend synthetic diagram of BJFS station

圖4 URUM站趨勢項合成Fig.4 Trend synthetic diagram of URUM station

致謝 感謝臺灣中央大學及Wu Zhaohua等人提供的EEMD分解算法，感謝美國麻省理工學院(MIT)和加州大學圣地亞哥分校 Scripps海洋研究所(SIO)及Robert B.King提供研究計算的GAMIT/ GLOBK軟件。

1 黃立人.GPS基準站坐標分量時間序列的噪聲特性分析［J］.大地測量與地球動力學，2006，(2):31-33.(Huang Liren.Niose properties in time series of coordinate component at GPS fiducial stations［J］.Journal of Geodesy and Geodynamics，2006，(2):31-33)

2 袁林果，等.香港GPS基準站坐標序列特征分析［J］.地球物理學報，2008，51(5):1 372-1 384.(Yuan Linguo，et al.Characteristics of daily position time series from the Hong Kong GPS fiducial network［J］.Chinese Journal of Geophysics，2008，51(5):1 372-1 384)

3 Huang Norden E，Shen Zheng and Long Steven R.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.The Roy-al Society，A(1998)454:903-995.

4 Huang Norden E and Shen Samuel S P.Hilbert-Huang transform and its applications［M］.Singapore:World Scientific，2005:2-14.

5 戴吾蛟，等.基于經(jīng)驗模式分解的濾波去噪法及其在GPS多路徑效應中的應用［J］.測繪學報，2006，35(4):321-327.(Dai Wujiao，et al.EMD filter method and its application in GPS multipath［J］.Acta Geodaetica at Cartographica Sinica，2006，35(4):321-327)

6 Wu Zhaohua and Huang Norden E.Ensemble empirical mode decomposition:a noise-assisted data analysis method［J］.Advances in Adaptive Data Analysis，2009，1(1):1-41.

7 Wu Zhaohua，et al.On the trend，detrending，and variability of nonlinear and nonstationary time series［J］.Proc.Natl.Acad.Sci.USA，2007，104(38):14 889-14 894.

8 徐世艷.經(jīng)驗模態(tài)分解的時頻分析方法及其應用［J］.吉林大學學報(信息科學版)，2009，27(5):487-492.(Xu Shiyan.Time-frequency analysis method and its application based on empirical mode decomposition［J］.Journal of Jilin University(Information Science Edition)，2009，27(5):487 -492)

9 蓋強，張海勇，徐曉剛.Hilbert-Huang變換的自適應頻率多分辨分析研究［J］.電子學報，2005，33(3):563-566.(Gai Qiang，Zhang Haiyong and Xu Xiaogang.Study of adaptive frequency multiresolution analysis of the Hilbert-Huang transform［J］.Acta Electronica Sinica，2005，33(3): 563-566)

10 徐佳，黃聲亨，麻鳳海.基于改進HHT理論的大型橋梁動態(tài)特性分析［J］.武漢大學學報(信息科學版)，2010，35(7):801-805.(Xu Jia，Huang Shengxiang and Ma Fenghai.The dynamic characteristics analysis for the large bridge based on the improved Hilbert-Huang transformation［J］.Geomatics and Information Science of Wuhan University，2010，35(7):801-805)

ANALYSIS ON TIME SERIES OF TWO CORS STATIONS’HEIGHT BASED ON EMD

Zhang Hengjing1，2)and Cheng Pengfei2)

(1)School of Geomatics of Liaoning Technical University，F(xiàn)uxin 123000 2)Chinese Academy of Surveyingamp;Mapping，Beijing100830)

The time-frequency analysis method base on EMD(empirical mode decomposition)for GPS height time series is put forward.Taking two national CORS stations’GPS height time series as an example，each intrinsic mode function is decomposed and the contribution rate of variance is also calculated.The problem of identifying father distinguishing the periodical signals and the trend is emphasizerd.The results indicate that it is only meaningful to distinguish the periodical and the trend signal on certain length of time series’section.The trend based on EMD or EEMD presents with large curvature while the traditional form of periodical fitting function only linear form.

empirical mode decomposition(EMD);intrinsic mode function;trend;GPS height time series;contribution rate of variance

1671-5942(2012)03-0129-06

2012-02-10

國家基礎(chǔ)測繪項目(B2551);中國測繪科學研究院基本科研業(yè)務費資助

張恒璟，男，1982年生，講師，博士生，主要從事地球動力學方面的研究.E-mail:zhhj1111@sohu.com

P207

A