基于EEMD和Husrt指數的GNSS基準站的垂向速率估計

彭　葳，戴吾蛟

(中南大學 地球科學與信息物理學院,湖南 長沙 410083)

?

基于EEMD和Husrt指數的GNSS基準站的垂向速率估計

彭葳，戴吾蛟

(中南大學 地球科學與信息物理學院,湖南 長沙 410083)

摘要：連續的全球衛星導航系統(GNSS)基準站的坐標時間序列中包含了復雜的噪聲信號、非構造形變及其它因素的影響，尤其在垂直方向，對GNSS基準站在國際地球參考框架(ITRF)下的運動速率估計產生了較大的干擾。為進一步提高速率精度，文中采用整體模態分解(EEMD)方法對GNSS基準站的垂向觀測時間序列進行分解，并根據各種信號的Hurst值進行分類及重構為噪聲信號、季節性信號和長期趨勢信號，采用最小二乘方法擬合長期趨勢信號得到垂向速率。通過對中國大陸構造環境監測網絡(CMONOC)的GNSS臺站從2001—2013年近13 a的垂向坐標時間序列的實例分析，采用基于EEMD和Husrt指數的最小二乘法能夠準確地估計GNSS基準站的垂向速率。

關鍵詞：整體經驗模態分解；Hurst指數；負荷改正；速率估計；最小二乘法

GNSS臺站的坐標和運動速率的準確估計可以用于參考框架的維持及地殼板塊運動的研究等方面。在國際參考框架(ITRF)下，GNSS時間序列在東方向(E)和北方向(N)的年變化量為厘米級，而噪聲和季節性變化為毫米級，故在N、E方向的速率估計較為準確；但垂向的形變較小(毫米級)，而噪聲和季節性變化影響卻達到厘米級。Nikolaidis根據GNSS時間序列中的噪聲、周期變化、趨勢變化的特性，建立了GNSS坐標時間序列的幾何模型，并采用加權最小二乘法求解，可以估計出站點運動速率[1]。但是GNSS觀測序列中的噪聲和周期變化會影響速率的估計，為了提高估計幾何模型中速率的估計精度，需要準確地分析觀測時間序列的噪聲和周期形變。GNSS坐標序列中的噪聲的估計方法有很多種，包括功率譜、極大似然估計、最小二乘方差分量估計等方法[2-3]。而GNSS坐標時間序列中的季節性變化主要由大氣負荷、土壤水負荷、非潮汐海洋負荷、積雪負荷等因素的影響[4]，許多研究者將季節性變化歸結于GNSS臺站之間的共模誤差，并采用主成分分析等方法將共模誤差削弱[5]。在削弱共模誤差之后，顧及有色噪聲擬合GNSS臺站的幾何模型，可以估計出GNSS站的速率[6]。但這些方法估計GNSS站的垂向速率仍需要顧及較多的多余參數，影響速率的估計。

GNSS時間序列中的隱含信號包含著特有的物理特征，采用時間序列分析方法可以分析GNSS時間序列中的隱含信息[7]。本文根據GNSS時間序列中包含的隱含信號(白噪聲、有色噪聲、周期形變和長期趨勢形變等)存在著不同的時頻特征，采用整體經驗模態分解(EEMD)根據信號的時頻特征進行分解[8]；同時，GNSS時間序列中包含的隱含信號又存在著不同的分形特征，采用Hurst指數可以對EEMD方法分解的子信號進行歸類。本文從單站的單維GNSS時間序列分析的角度出發，采用EEMD方法分解單維的GNSS垂向的觀測時間序列，得到一組固有模態分量(IMF)，然后通過估計EEMD分解的IMF分量的Hurst指數，可以較為準確地提取噪聲序列[9]。而參考去噪后的負荷改正模型的Hurst指數，能夠準確提取GNSS的坐標時間序列中的非構造形變。而準確地分析噪聲和非構造形變能夠減少GNSS幾何模型中的參數，從而提高模型的精度，對于形變較小的GNSS垂向時間序列的速率項的估計有著重要的作用。

1GNSS垂向時序速率估計模型

根據GNSS垂向時間序列的特征，本文建立的GNSS垂向時間序列的速率估計流程如圖1所示。

圖1　GNSS垂向時間序列的速率估計流程

1.1經驗模態分解

經驗模態分解(EMD)是一種廣泛用于從含有噪聲的非線性、非平穩的數據中提取信號的自適應的時頻分析方法[10]。假設一組時間序列x(t)，EMD方法它是根據信號的時頻特征將時間序列x(t)分解成若干個子序列和殘差序列，而這些不同頻率的子序列也被稱為固有模態函數(IMF)。

EMD的一般公式為

(1)

其中:di(t)為第i個IMF分量，rn(t)為殘差序列。

EMD分解的主要步驟如下：

1)尋找信號中的所有極大值和極小值，擬合成上包絡線emax(t)和下包絡線emin(t)。

3)計算時間序列與上下包絡線均值的差值d(t)=x(t)-m(t)。

4)重復1)～3)，迭代殘余值m(t)。

EMD可以從信號中分解頻率從高到低的IMF分量，高頻的IMF分量一般認為是噪聲信號，中低頻的IMF分量需要物理因素來進行分析解釋，但是EMD方法存在著模態混疊的問題，這會使得IMF分量的物理意義不明確。為了解決EMD方法的不足，Wu等提出了整體經驗模態分解(EEMD)方法[8]。

首先，在時間序列x(t)加入不同的白噪聲序列。

(2)

其中，wi(t)為第i次加入的白噪聲，且每次隨機加入的白噪聲的方差相同。

然后，將得到的xi(t)進行EMD分解，得到N組IMF分量，取N組中對應尺度的IMF分量的均值，得到EEMD分解的最終IMF分量。EEMD分解的信號可以表示為

(3)

1.2噪聲和非構造形變的Hurst指數分析

Hurst指數會發生在許多應用數學的領域，包括分形和混沌理論、長記憶過程和譜分析。隨機過程也可以產生自相似分形，其中一個重要的例子是布朗運動(Brown motion)。分數布朗運動(fractional Brown motion,fBm) 是非平穩的隨機過程，而分數高斯噪聲(fractional Gaussian noise,fGn) 是一個零均值平穩過程，可看作為分數布朗運動的增量。

GNSS時間序列的分數高斯噪聲序列的頻率和功率譜密度的關系可以描述為[11]

(4)

其中，α為譜指數。

消除趨勢波動分析(DFA)方法是一種分形標度指數計算方法，通常用于檢測非平穩時間序列的自相似性，可用于時間序列的物理特征的分析。DFA方法計算的Hurst指數(H)可以描述噪聲序列的分形特征，也作為噪聲分類及提取的依據。Hurst指數與噪聲的關系可以描述為

(5)

在上式的定義下，H=0.5為白噪聲，H=1為閃爍噪聲，H=1.5為隨機游走信號[9]。白噪聲和閃爍噪聲與環境負荷效應的影響相關性不強，且白噪聲加閃爍噪聲模型(WN+FN)符合大多數GNSS基準站數據中的噪聲特征[12]，故將H≤1的IMF分量作為噪聲信號提取[13]；此公式定義下的隨機游走信號可能與負荷效應有關，需要進一步分析[14]。

在中國大陸地區，環境負荷(大氣負荷、土壤水負荷、積雪負荷等)對地殼的影響在幾年至十幾年的時間范圍內是較為平穩的，負荷影響的地殼季節性變化也具有特有的物理特征(頻率和功率譜密度之間的關系)。通過Hurst指數進行信號特征分析，重新組合EEMD分解的IMF分量，可以快速提取GNSS時間序列中的非構造形變。

1.3最小二乘速率估計

GNSS單測站的任一分量的幾何模型可以表示為[1]

其中：a為初始偏差；b為速率；c(t)為年周期和半年周期項；d(t)為其他項，包括偏移項、震后速率變化等等；v為殘差項，且在殘差項中還包含有白噪聲和有色噪聲[15]。

在準確估計出噪聲序列和非構造形變序列之后，GNSS單測站的某一分量的坐標時間序列y′(t)可以表示為

其中：a′為初始偏差；b′為速率；v′為殘差，且v′中不包含噪聲信號。采用最小二乘方法對上式求解，求出速率項系數b′。

2非構造形變分析

本文選取了中國大陸構造環境監測網絡的17個站的垂向時間序列，從2001年至2013年的垂向時間序列作為實例分析。GNSS垂向時間序列的變化都表現為含有白噪聲、有色噪聲、周期波動和長期趨勢運動[16]。采用EEMD方法對GNSS垂向時間序列進行分解，加入的白噪聲與信號的標準差的比例為0.2，平均次數為200次。圖2為作為示例的EEMD方法分解的BJFS、JIXN、KMIN 3個站的垂向時間序列的IMF和殘差項，其余14個站的EEMD分解的情況與之相似。

圖2　EEMD方法分解的GNSS站(BJFS、JIXN、KMIN)的垂向時間序列的IMF分量及殘差項(Resi)

計算EEMD分解的IMF分量的Hurst值，將H≤1的IMF分量重構為噪聲信號，然后將剩余IMF重新編號(1,2,3,…)，如圖3所示。

圖3　去噪后的剩余IMF分量的Hurst值

從圖3中可以發現各個站的剩余的IMF分量中，第2個IMF分量和第3個IMF分量會產生明顯的分離，而這個分離的界限的Hurst值約為1.7，故可以選取H≤1.7作為區分非構造形變和趨勢形變的界限值。

為了驗證信號區分的準確性，本文引入環境負荷的改正模型，用于對比分析EEMD方法和Hurst指數提取的GNSS垂向時間序列中的非構造形變信號。QOCA軟件是一款GPS數據的后處理軟件，它能夠較為準確地計算各種負荷效應的改正值。大氣負荷、土壤水負荷、積雪負荷的數據從NCEP(National Center for Environmental Prediction) 獲得。

1) 大氣地表壓力數據的采樣時間間隔為6 h，空間的分辨率為2.5°×2.5°；

2) 土壤水數據的采樣間隔為1 d，空間的分辨率為1.875°×1.875°；

3) 積雪數據的采樣間隔為1 d，空間的分辨率為1.875°×1.875°。

通過QOCA求出各種負荷每天的均值，得到時間分辨率為1 d的負荷改正值的時間序列。同樣采用EEMD方法和Hurst指數濾除信號中的噪聲，計算17個GNSS監測站的剩余分量之和的Hurst值，其均值約為1.4，與H=1.7界限存在著差異的可能原因如下：

1) EEMD將非構造形變分解成多個IMF分量；

2) DFA方法計算Hurst指數存在計算誤差(δH≤0.1)[17]。

通過對比圖4中環境負荷改正值總和與使用本文方法提取的非構造形變信號，可以看出通過EEMD和Hurst指數的方法提取的非構造形變信號與環境負荷的改正值信號在波形上高度相似，說明該方法能有效地提取非構造形變。

圖4　環境負荷改正值總和(紅色線)與EEMD+Hurst方法提取的季節性信號(黑色線)的對比圖

3GNSS時間序列分析及速率估計

3.1模擬實驗

根據GNSS基準站的垂向時間序列的信號特征[16]，本文模擬了白噪聲、閃爍噪聲、季節性形變和趨勢信號。其中白噪聲由MATLAB產生，閃爍噪聲由文獻[13]中的遞歸方法產生，季節性信號是由QOCA軟件(http://gipsy.jpl.nasa.gov/qoca/)解算的長春站(CHUN)的去噪后的負荷改正值，趨勢信號由直線模擬。模擬的兩個信號中，一個加入的是較小的趨勢信號，另一個加入較大的趨勢信號，模擬的時間長為4 747 d，模擬的兩個不同速率的信號如圖5所示。

EEMD方法對兩個模擬信號進行分解，加入的白噪聲與信號的標準差的比例為0.2，平均次數為100次。模擬信號1被分解成11個IMF分量和1個殘余項，模擬信號2被分解成9個IMF分量和1個殘余項。計算得到的IMF分量的Hurst值，提取H≤1的噪聲信號，因為是模擬數據是模擬中國地區的GNSS監測站,故直接提取H≤1.7的IMF分量為非構造形變信號；剩余的IMF分量重構為趨勢信號，并采用最小二乘法進行擬合。將本文方法計算結果與CATS軟件計算結果進行對比，CATS (The Create and Analyze Time Series) 是一款采用最小二乘擬合時間序列的多參數模型，同時顧及了殘差時間序列的噪聲序列[18]，圖6給出了對比的結果。

圖5　模擬的兩個不同速率的模擬信號

圖6　基于EEMD和Hurst指數的最小二乘擬合的直線(紅線)、CATS擬合的直線(藍線)與模擬直線(黑線)的對比

圖6顯示的是基于EEMD和Hurst指數的最小二乘擬合的直線(紅線)、CATS擬合的直線(藍線)與模擬直線(黑線)的對比，可以發現紅線和黑線的斜率(速率)更為接近，它們的值如表1所示。

表1　基于EEMD和Hurst指數最小二乘擬合方法與CATS得到的垂向速率及擬合精度的對比

通過模擬實驗可以說明EEMD+Hurst方法估計得到的速率更趨近于真值，且該方法的擬合精度要高于CATS的精度。

3.2GNSS垂向時間序列的速率估計

在準確地估計了GNSS垂向時間序列中的噪聲信號和非構造形變信號后，采用最小二乘估計速率。圖7為作為示例的BJFS、JIXN、KMIN原始坐標序列、長期趨勢信號和擬合直線的對比。

圖7　原始時間序列(灰線)、長期趨勢信號(紅線)和擬合直線(藍線)的對比

從圖7中可以看出，提取的長期趨勢信號能較好地反映原始時間序列的趨勢變化的特征，且擬合的直線與長期趨勢信號的相似度也很高，說明使用基于EEMD和Hurst指數的垂向速率估計是較為準確的。

為了驗證EEMD和Hurst指數的垂向速率估計方法得到的GNSS時間序列的垂向速率的有效性和準確性，將本文方法計算結果與CATS軟件計算結果進行對比，對比結果如表2所示。

表2　基于EEMD和Hurst指數最小二乘擬合方法與CATS得到的垂向速率及擬合精度的對比

從表2中可以看出，采用基于EEMD和Hurst指數的最小二乘法估計的垂向速率和CATS軟件估計的垂向速率很接近。但是，當站臺的運動速率變化較小時，CATS軟件的擬合精度相對垂向速率較大，影響了垂向速率分析的準確性；而本文算法由于準確提取了趨勢信號，擬合精度較高，能較為準確地估計趨勢變化垂向速率。

4結論

通過基于EEMD和最小二乘方法模型可以準確地分析GNSS垂向時間序列及估計速率項。通過對中國大陸構造環境監測網絡(CMONOC)的17個站的垂向時間序列的分析及模擬試驗，可以得到以下兩個結論：①Hurst指數及環境負荷可以有效地對EEMD分解的IMF分量中的噪聲和非構造形變信號進行識別和提取；②基于EEMD+Hurst的最小二乘法得到的垂向能夠較好地反映地殼趨勢運動變化，得到的垂向速率的擬合精度較高，這對于有些垂向形變較小的GNSS參考站的速率估計有非常大的意義。

參考文獻：

[1]NIKOLAIDIS R.Observation of geodetic and seismic deformation with the Global Positioning System[D].San Diego：University of California,2002.

[2]ZHANG J,BOCK Y,JOHNSON H,et al.Southern California Permanent GPS Geodetic Array：Error analysis of daily position estimates and site velocities[J].Journal of Geophysical Research：Solid Earth (1978～2012),1997,102(B8)：18035-18055.

[3]TEUNISSEN P J G,AMIRI-SIMKOOEI A R.Least-squares variance component estimation[J].Journal of Geodesy,2008,82(2)：65-82.

[4]DONG D,FANG P,BOCK Y,et al.Anatomy of apparent seasonal variations from GPS‐derived site position time series[J].Journal of Geophysical Research：Solid Earth (1978-2012),2002,107(B4)：ETG 9-1-ETG 9-16.

[5]袁林果,丁曉利,陳武,等.香港 GPS 基準站坐標序列特征分析[J].地球物理學報,2008,51(5)：1372-1384.

[6]蔣志浩,張鵬,秘金鐘,等.顧及有色噪聲影響的CGCS2000下我國CORS站速度估計[J].測繪學報,2010,39(4)：355.

[7]章浙濤,朱建軍,盧駿,等.小波變換在時間序列特征提取中的應用[J].測繪工程,2014,23(6)：21-26.

[8]WU Z,HUANG N E.Ensemble empirical mode decomposition：a noise-assisted data analysis method[J].Advances in adaptive data analysis,2009,1(1)：1-41.

[9]MONTILLET J P,TREGONING P,MCCLUSKY S,et al.Extracting white noise statistics in GPS coordinate time series[J].Geoscience and Remote Sensing Letters,IEEE,2013,10(3)：563-567.

[10] HUANG N E,SHEN Z,LONG S R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society of London.Series A：Mathematical,Physical and Engineering Sciences,1998,454(1971)：903-995.

[11] AGNEW D C.The time-domain behavior of power‐law noises[J].Geophysical research letters,1992,19(4)：333-336.

[12] 黃立人.GPS 基準站坐標分量時間序列的噪聲特性分析[J].大地測量與地球動力學,2006,26(2)：31-33.

[13] RODRIGUEZ E,ECHEVERRIA J C,ALVAREZ-RAMIREZ J.1/fαfractal noise generation from Grünwald-Letnikov formula[J].Chaos,Solitons & Fractals,2009,39(2)：882-888.

[14] 田云鋒,沈正康.GPS 坐標時間序列中非構造噪聲的剔除方法研究進展[J].地震學報,2009,31(1)：68-81.

[15] WILLIAMS S D P.The effect of coloured noise on the uncertainties of rates estimated from geodetic time series[J].Journal of Geodesy,2003,76(9-10)：483-494.

[16] 榮敏,孫付平,賈小林,等.利用GPS基準站數據淺析我國地殼垂直運動[J].測繪工程,2009,18(4)：7-9.

[17] GRECH D,MAZUR Z.On the scaling ranges of detrended fluctuation analysis for long-term memory correlated short series of data[J].Physica A：Statistical Mechanics and its Applications,2013,392(10)：2384-2397

[18] WILLIAMS S D P.CATS：GPS coordinate time series analysis software[J].GPS solutions,2008,12(2)：147-153.

[責任編輯：劉文霞]

Vertical velocity estimation of GNSS reference stationbased on EEMD method and Husrt exponent

PENG Wei,DAI Wujiao

(School of Geosciences and Info-Physics,Central South University,Changsha 410083,China)

Abstract：The daily coordinate time series of Global satellite navigation system (GNSS) reference stations contain the complex effects of noise signals,non-tectonic deformation and other factors,particularly in the vertical direction,of which the factors have a greater impact on the velocity estimation of GNSS reference station on the International Terrestrial Reference Frame (ITRF).In order to improve the accuracy of the velocity estimation,this paper uses the Ensemble Empirical Mode Decomposition (EEMD) method to decompose the vertical time series of GNSS reference station,and the noise signal,seasonal signal and long-term trend signal can be reconstructed respectively based on the Hurst value of signals.Then the vertical velocity can be estimated by using the least square method to fit the long-term trend signals.By analyzing the vertical time series of GNSS stations of Crustal Movement Observation Network of China (CMONOC) from January 2001 to December 2013,the EEMD method and Husrt exponent based on the least square method can accurately estimate the vertical velocity of GNSS reference station.

Key words：EEMD；Hurst exponent；load correction；velocity estimation；least square method

中圖分類號：P228

文獻標識碼：A

文章編號：1006-7949(2016)04-0060-06

作者簡介：彭葳(1989-)，男，碩士研究生.

基金項目：國家973計劃資助項目(2013CB733303)

收稿日期：2015-02-04；修回日期：2015-03-08