張恒璟 程鵬飛
(1)遼寧工程技術(shù)大學(xué)測繪與地理科學(xué)學(xué)院,阜新 123000 2)中國測繪科學(xué)研究院,北京100830)
基于經(jīng)驗(yàn)?zāi)J椒纸獾腃ORS站高程時(shí)間序列分析*
張恒璟1,2)程鵬飛2)
(1)遼寧工程技術(shù)大學(xué)測繪與地理科學(xué)學(xué)院,阜新 123000 2)中國測繪科學(xué)研究院,北京100830)
提出基于經(jīng)驗(yàn)?zāi)J椒纸饧夹g(shù)的GPS高程時(shí)間序列時(shí)頻分析方法,并以兩個(gè)CORS基準(zhǔn)站多年連續(xù)的高程時(shí)間序列數(shù)據(jù)為對象,獲取序列信號的各個(gè)本征模態(tài)函數(shù)分量,研究了信號周期運(yùn)動的主要貢獻(xiàn)分量和序列趨勢項(xiàng)的合成與分離的方法。結(jié)果表明:在一定的序列區(qū)間長度上識別周期項(xiàng)與趨勢項(xiàng)才有意義;經(jīng)驗(yàn)?zāi)J椒纸饧夹g(shù)分離的序列趨勢項(xiàng)呈現(xiàn)大曲率的非線性形式,與傳統(tǒng)的周期函數(shù)高程序列擬合線性趨勢項(xiàng)存在較大的差異。
經(jīng)驗(yàn)?zāi)J椒纸?本征模態(tài)函數(shù);趨勢項(xiàng);GPS高程時(shí)間序列;方差貢獻(xiàn)率
正確識別并分離序列趨勢項(xiàng)對研究GPS高程序列的頻譜特性和噪聲性質(zhì)具有十分重要的意義。傳統(tǒng)時(shí)頻信號分析前,須先識別并去除序列的趨勢項(xiàng),例如基于傅里葉分析理論的功率譜估計(jì)和小波分解方法。經(jīng)典方法研究GPS高程時(shí)序信號的周期和噪聲特性時(shí),事先給定一個(gè)具體的周期擬合模型[1,2],并指定序列存在半年、年、兩年等的周期運(yùn)動和長期的線性趨勢運(yùn)動,通過函數(shù)擬合模型,求解各項(xiàng)的振幅相位等特性;接著在分析序列噪聲時(shí),從原始序列中去除趨勢項(xiàng)和周期項(xiàng)后的殘差序列中,分析信號存在的各種噪聲特性,但實(shí)際上并不是每個(gè)GPS高程序列都具有事先設(shè)定的多種周期運(yùn)動特性。因此有必要采用自適應(yīng)的或者盡量不依賴事先給定特性的序列信號分解方法來識別序列的各種周期運(yùn)動和趨勢項(xiàng)信息。本文引入處理非線性非平穩(wěn)信號的希爾伯特-黃變換方法[3],采用經(jīng)驗(yàn)?zāi)J椒纸?Empirical Mode Decomposition,EMD)和整體經(jīng)驗(yàn)?zāi)J椒纸?Ensemble EMD,EEMD),對我國CORS基準(zhǔn)站高程序列數(shù)據(jù)進(jìn)行分析,獲取了本征模態(tài)函數(shù)分量(Intrinsic Mode Function,IMF),分析了各個(gè)IMF的方差貢獻(xiàn)率及對序列運(yùn)動能量的貢獻(xiàn)大小,并研究了合并和提取序列高程運(yùn)動趨勢項(xiàng)的方法。
經(jīng)驗(yàn)?zāi)J椒纸馐腔诜蔷€性、非平穩(wěn)數(shù)據(jù)的自適應(yīng)時(shí)頻信號分解方法,通過不斷地“篩選”,將原始時(shí)間序列信號分解為一系列頻率由高到低的本征模態(tài)函數(shù)分量和殘差項(xiàng),分解產(chǎn)生的IMF函數(shù)需滿足兩個(gè)條件:一是待分解信號中極值點(diǎn)的個(gè)數(shù)與過零點(diǎn)的個(gè)數(shù)相等或至多相差1;二是在任意點(diǎn)上,由局部極大值定義的上包絡(luò)線與局部極小值定義的下包絡(luò)線的均值為零。分解的具體過程為[4,5]:
1)搜索原始信號所有的極大值點(diǎn)和極小值點(diǎn),用三次樣條函數(shù)擬合所有極大值點(diǎn)得到上包絡(luò)線,同理得到下包絡(luò)線,取上下包絡(luò)的均值作為原始序列的平均包絡(luò)m1(t);從原始序列中扣除該平均包絡(luò)得到新的序列信號h1(t)=X(t)-m1(t)。如果h1(t)滿足IMF的條件,則作為第一個(gè)分解得到的本征模態(tài)函數(shù)分量,否則以h1(t)作為信號,繼續(xù)尋找其上下包絡(luò)線并求出平均包絡(luò)m2(t),得到新信號h2(t)=h1(t)-m2(t),迭代k次后直至hk(t)滿足IMF的條件,將其作為第一個(gè)IMF分量IMF1(t)。
2)計(jì)算殘余分量r1=X(t)-IMF1(t),將r1作為信號,重復(fù)第一步的分解,得到第二個(gè)IMF分量IMF2(t)。
3)重復(fù)1)、2)步,直至殘余分量滿足分解中止的條件:即殘余分量為單調(diào)變化的殘差序列不能再分解出本征模態(tài)函數(shù)為止。原始信號經(jīng)過n次分解后,可描寫為:
這種分解基于數(shù)據(jù)的局部特征,與傅立葉頻譜分析分解成恒定振幅與頻率的一系列正弦函數(shù)和小波分解需預(yù)先給定基函數(shù)相比,EMD是一種自適應(yīng)的分解方法。
EMD分解的缺點(diǎn)之一是模式混疊現(xiàn)象:一個(gè)IMF分量包含著尺度互異的信號,或者尺度相似的信號存在于不同的IMF分量,這種現(xiàn)象主要是由于信號中斷引起的,不僅混淆了信號的時(shí)頻分布特性,更重要的使信號IMF分量的物理意義變得不清晰。為了減弱模式混疊現(xiàn)象的影響,文獻(xiàn)[6]在EMD分解基礎(chǔ)上提出了整體經(jīng)驗(yàn)?zāi)J椒纸夥椒?EEMD),一種噪聲輔助的數(shù)據(jù)分析方法(noise-assisted data analysis,NADA),基本思路是在分解之前,將原始序列信號加入白噪聲,生成待分解信號,即
下標(biāo)i表示第i次EMD分解。
白噪聲與原始信號的標(biāo)準(zhǔn)差之比例如可以設(shè)定為0.1,經(jīng)過EMD分解后生成一系列的IMF分量;由于每次隨機(jī)生成的白噪聲都是不同的,這樣可以通過重復(fù)上面的分解過程,得到信號的多次EMD分解的IMF分量,將對應(yīng)分解尺度相同的IMF分量取平均,即得到了最終信號分解的IMF分量。文獻(xiàn)[6]建議整體平均的次數(shù)取100,而噪聲與信號的標(biāo)準(zhǔn)差之比可以取0.1、0.2、0.4等。本文最大分解層數(shù)n設(shè)定為:
其中N是序列長度,實(shí)驗(yàn)采用的序列從1999—2010,共約4 000個(gè)數(shù)據(jù)點(diǎn),因此得到的最大分解層數(shù)(IMF分量的個(gè)數(shù))不超過10。這樣就存在序列過度分解或分解不完全兩種情形,從實(shí)驗(yàn)的幾個(gè)CORS站高程序列分解效果來看,剩余的殘差項(xiàng)均是單調(diào)的趨勢,符合EMD分解要求,不存在分解不完全的情形,而過度分解的情形需要對鄰近殘差項(xiàng)的IMF分量與殘差項(xiàng)合并,來判斷合并后的趨勢項(xiàng)是否仍然單調(diào)變化。
傳統(tǒng)序列趨勢項(xiàng)是基于固定的數(shù)學(xué)模型定義,比如GPS時(shí)間序列常用的周期擬合模型為:
并由式(4)擬合出序列的常數(shù)項(xiàng)、線性速率項(xiàng)、周期項(xiàng)(年和半年周期)。式中,ti(i=1,…,N)是以年為單位的時(shí)間,以每日解構(gòu)成坐標(biāo)時(shí)間序列,待求系數(shù)a和b為序列線性趨勢項(xiàng),c、d和e、f分別為年周期和半年周期項(xiàng)的系數(shù)。
更復(fù)雜的情形,在上述模型中考慮各種原因引起的階躍式坐標(biāo)突變(例如遠(yuǎn)場大地震引起的同震位移、由于儀器或天線變更引起的位移、甚至由于某種不清楚原因引起的點(diǎn)位變化等)修正項(xiàng),這是目前處理GPS時(shí)間序列常用的擬合模型。該模型扣除了線性速度項(xiàng)和常數(shù)項(xiàng)后的序列,是下一步時(shí)序信號頻譜和噪聲分析的基礎(chǔ),比如基于傳統(tǒng)傅里葉理論的功率譜分析技術(shù)和目前常用的小波分解技術(shù)等。
文獻(xiàn)[6]指出,趨勢應(yīng)該是序列信號的內(nèi)在特性之一,是序列的固有成分,內(nèi)在特性要求趨勢項(xiàng)的定義應(yīng)該具有自適應(yīng)性,應(yīng)基于序列數(shù)據(jù)本身去分離,而不僅僅是用一個(gè)指定的函數(shù)模型去擬合和扣除。同時(shí)為區(qū)分趨勢項(xiàng)和長周期項(xiàng),趨勢項(xiàng)的定義必須被限定在給定的序列尺度上,并且至多包含一個(gè)極值點(diǎn)?;谏厦娴目紤],文獻(xiàn)[7]給出了趨勢項(xiàng)的定義,即:趨勢項(xiàng)是一個(gè)內(nèi)在的單調(diào)擬合函數(shù),或在一個(gè)給定的序列數(shù)據(jù)區(qū)間上至多有一個(gè)極值點(diǎn)的函數(shù),“給定的序列數(shù)據(jù)區(qū)間”可以是整個(gè)序列,也可以是序列數(shù)據(jù)的一部分。
對GPS高程序列數(shù)據(jù)分別采用EMD和EEMD兩種分解方法,獲取IMF分量。EEMD方法又設(shè)計(jì)了3種方案,整體平均的次數(shù)選用了文獻(xiàn)[6]的推薦值100,加入的白噪聲與原始序列信號的標(biāo)準(zhǔn)差之比分別取0.1、0.2、0.4,下面以EEMD01、EEMD02、EEMD04分別代表這3種EEMD分解方案。統(tǒng)計(jì)任一IMF分量的方差公式為:=E2(X)-(E (X))2,若X為有限長度的數(shù)據(jù)序列,則可統(tǒng)計(jì)出方差的估值。方差貢獻(xiàn)率是各個(gè)IMF分量的方差與分解得到的IMF分量方差之和的百分比。IMF分量方差貢獻(xiàn)率的大小反映了該頻率的分解信號在整個(gè)序列信號運(yùn)動能量中的貢獻(xiàn)大小。
圖1(a)、(b)分別是BJFS(北京房山)和URUM (烏魯木齊)站采用 EMD/EEMD01/EEMD02/EEMD04 4種分解方案的IMF分量和殘差,限于篇幅未給出高頻的IMF1~5曲線。由公式(3)可知分解得到的IMF分量個(gè)數(shù)不超過10,從IMF1~10本征模態(tài)分量信號頻率依次降低,周期依次增大。
圖2是BJFS和URUM兩個(gè)基準(zhǔn)站采用EMD/ EEMD01/EEMD02/EEMD04 4種方案統(tǒng)計(jì)得到的各個(gè)IMF分量中誤差和方差貢獻(xiàn)率。
BJFS站分解后的IMF1~3是高頻分量,IMF4和IMF5是近似月周期和雙月周期信號,但信號的振幅和方差較小,在整個(gè)序列信號中所占的比重較小。IMF6、7兩個(gè)信號分量的方差貢獻(xiàn)率之和達(dá)到50%以上,是序列的主要周期運(yùn)動貢獻(xiàn)項(xiàng),即BJFS站的周期運(yùn)動以季節(jié)性變化和年周期變化為主。兩種分解方法得到的IMF6分量整體呈現(xiàn)季節(jié)性的周期變化,IMF7分量整體呈現(xiàn)明顯的年周期變化;但EMD分解與3種EEMD分解方案的方差貢獻(xiàn)率出現(xiàn)較大差異,EMD分解的IMF6、7方差貢獻(xiàn)率分別是25%與28%,EEMD分解分別是12%與40%左右,兩種技術(shù)得到的IMF6與IMF7方差貢獻(xiàn)率之和幾乎保持不變,其他本征模態(tài)分量的方差貢獻(xiàn)率基本相等。從圖1(a)明顯看出,EMD分解得到的IMF6季節(jié)性信號在橫軸2003—2006區(qū)間上出現(xiàn)了明顯的年周期信號,EMD分解得到的IMF7年周期信號在橫軸2003—2006區(qū)間上出現(xiàn)了明顯的兩年周期信號,說明EMD技術(shù)分解時(shí)部分IMF6信號被分解到IMF7分量上,這是由于模式混疊現(xiàn)象引起的,在圖1(a)中以橢圓部分表示。IMF8分量對應(yīng)著明顯的兩周年運(yùn)動,但其方差貢獻(xiàn)率和信號振幅相對較小。IMF9和IMF10分量振幅在5 mm之內(nèi),是信號分解的長周期還是屬于過度分解的殘差信號,在趨勢項(xiàng)合并部分做進(jìn)一步的分析。本文也設(shè)計(jì)了采用相同的EEMD方案對BJFS站多次分解的試驗(yàn),結(jié)果表明每次試驗(yàn)之間低頻的本征模態(tài)分量IMF6、7的方差貢獻(xiàn)率出現(xiàn)了0.5~1.9%的微小變化,但是EMD的分解結(jié)果唯一,分析原因可能是EEMD分解時(shí)每次加入的白噪聲變化引起的,這也說明了EEMD技術(shù)仍存在著分解不唯一的缺點(diǎn),需進(jìn)一步的研究和完善EEMD分解算法。
圖1 IMF6~10分量與殘差項(xiàng)Fig.1 IMF6-10 components and residuals
圖2 IMF分量中誤差與方差貢獻(xiàn)率Fig.2 RMS and contribution rate of variance of each IMF
URUM站3種EEMD分解結(jié)果一致性很好,與EMD分解結(jié)果相比在中低頻IMF分量上有一定的差異。IMF1~3是分解得到的高頻信號,這3個(gè)分量方差貢獻(xiàn)率之和在20%左右,相對于BJFS站這3個(gè)高頻分量的方差貢獻(xiàn)率之和達(dá)到35%,GPS時(shí)間序列中高頻信號一般被認(rèn)為是序列所含的噪聲,從這個(gè)角度上可認(rèn)為URUM站高程序列信號中噪聲含量明顯低于BJFS站。IMF4分量周期大約是1個(gè)月,IMF5分量周期大約是兩個(gè)月。兩種分解技術(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個(gè)分量均出現(xiàn)了模式混疊現(xiàn)象,雖然EEMD技術(shù)可以減弱模式混疊的影響,但從URUM站IMF6分量的分解結(jié)果可知,EEMD仍會出現(xiàn)頻譜混疊,需要對EEMD分解技術(shù)進(jìn)一步的研究。從方差貢獻(xiàn)率指標(biāo)定量分析:IMF6~8 3個(gè)分量EMD分解的方差貢獻(xiàn)率分別是32%、35%、4%,EEMD01和EEMD02方案分解結(jié)果幾乎完全相同方差貢獻(xiàn)率分別是21%、20%、21%,兩種分解技術(shù)得到的3個(gè)分量方差貢獻(xiàn)率之和在62~71%之間,說明URUM站的高程運(yùn)動以季節(jié)性、年、兩年為主要周期。
圖3、圖4分別為BJFS和URUM站的趨勢項(xiàng)合成。
BJFS站IMF9、10分量的中誤差在0.8 mm左右,方差貢獻(xiàn)率為1.1~1.7%,這兩個(gè)本征模態(tài)分量在整個(gè)序列運(yùn)動中所占的比重很小。從趨勢項(xiàng)合成圖3可知,最后的兩個(gè)本征模態(tài)分量IMF9、10一起與分解的殘差項(xiàng)合并后,曲線幾乎保持原有的單一變化趨勢,說明IMF9、10分量屬于過度分解,也有可能是以目前序列數(shù)據(jù)的尺度還無法識別的更長周期項(xiàng),因此在研究的序列區(qū)間尺度上IMF9、10分量應(yīng)歸于分解的殘差,一并作為序列趨勢項(xiàng)處理。
URUM站4種不同的分解方案對應(yīng)的本征模態(tài)分量IMF10中誤差均在0.5mm左右,方差貢獻(xiàn)率在0.1~0.6%之間,在整個(gè)序列運(yùn)動中所占比重幾乎可忽略不計(jì)。從趨勢項(xiàng)合成圖4可知,最后一個(gè)模態(tài)分量IMF10應(yīng)歸于過度分解產(chǎn)生的分量,與殘差合并一起作為趨勢項(xiàng)處理。加入IMF9分量之后的趨勢項(xiàng)呈現(xiàn)明顯的周期振蕩,周期近似為序列長度的一半5.5年,加入IMF8合成后的趨勢項(xiàng)周期更小且振幅顯著增大,說明IMF8、9分量不能作為序列的趨勢項(xiàng),而應(yīng)當(dāng)歸于序列運(yùn)動較明顯的周期項(xiàng)分量,只是其振幅較小。
1)采用EMD/EEMD方法,獲取頻率從高到低的IMF分量,可明顯分解出GPS高程序列信號的趨勢項(xiàng)和周期項(xiàng)。從兩個(gè)基準(zhǔn)站例子可知,分解后的低頻IMF分量為高程時(shí)間序列信號的周期運(yùn)動給出了清晰的解釋,包含了明顯的1個(gè)月、2個(gè)月、季節(jié)性、年周期變化、兩周年周期變化與長周期變化因素,半周年項(xiàng)運(yùn)動信號并不顯著,季節(jié)性和年周期變化是高程運(yùn)動的主要貢獻(xiàn)項(xiàng);EEMD分解相比EMD可以明顯減弱模式混疊對本征模態(tài)分量的影響,但并沒有完全消除模式混疊現(xiàn)象,同時(shí)受EEMD分解時(shí)加入白噪聲的影響,其分解結(jié)果也不唯一。EEMD比EMD方法有較大的優(yōu)點(diǎn),但還需要進(jìn)一步改進(jìn)EEMD算法。
2)通過比較各個(gè)IMF分量的方差貢獻(xiàn)率及其中誤差大小,可以判斷對原始信號是否存在著過度分解的情形,進(jìn)而合并過度分解的IMF分量與殘差項(xiàng),獲得準(zhǔn)確的序列趨勢項(xiàng)。合并后的趨勢項(xiàng)呈現(xiàn)大曲率的曲線變化,而非線性形式,傳統(tǒng)對序列趨勢項(xiàng)的線性擬合方法存在不足,序列趨勢項(xiàng)不僅包括了線性變化,也有暫時(shí)無法分離或識別的長周期信號,只有在一定的序列區(qū)間長度上探討并區(qū)分周期項(xiàng)和趨勢項(xiàng)才有意義。
圖3 BJFS站趨勢項(xiàng)合成Fig.3 Trend synthetic diagram of BJFS station
圖4 URUM站趨勢項(xiàng)合成Fig.4 Trend synthetic diagram of URUM station
致謝 感謝臺灣中央大學(xué)及Wu Zhaohua等人提供的EEMD分解算法,感謝美國麻省理工學(xué)院(MIT)和加州大學(xué)圣地亞哥分校 Scripps海洋研究所(SIO)及Robert B.King提供研究計(jì)算的GAMIT/ GLOBK軟件。
1 黃立人.GPS基準(zhǔn)站坐標(biāo)分量時(shí)間序列的噪聲特性分析[J].大地測量與地球動力學(xué),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基準(zhǔn)站坐標(biāo)序列特征分析[J].地球物理學(xué)報(bào),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)驗(yàn)?zāi)J椒纸獾臑V波去噪法及其在GPS多路徑效應(yīng)中的應(yīng)用[J].測繪學(xué)報(bào),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)驗(yàn)?zāi)B(tài)分解的時(shí)頻分析方法及其應(yīng)用[J].吉林大學(xué)學(xué)報(bào)(信息科學(xué)版),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 蓋強(qiáng),張海勇,徐曉剛.Hilbert-Huang變換的自適應(yīng)頻率多分辨分析研究[J].電子學(xué)報(bào),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 徐佳,黃聲亨,麻鳳海.基于改進(jìn)HHT理論的大型橋梁動態(tài)特性分析[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),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ǔ)測繪項(xiàng)目(B2551);中國測繪科學(xué)研究院基本科研業(yè)務(wù)費(fèi)資助
張恒璟,男,1982年生,講師,博士生,主要從事地球動力學(xué)方面的研究.E-mail:zhhj1111@sohu.com
P207
A