基于EMD方法的地心運(yùn)動(dòng)時(shí)間序列分析

喬靈娜，趙春梅，馬天明,

(1. 山東科技大學(xué)測(cè)繪科學(xué)與工程學(xué)院，山東 青島 266590； 2. 中國(guó)測(cè)繪科學(xué)研究院，北京 100830； 3. 遼寧工程技術(shù)大學(xué)，遼寧 阜新 123000)

地球質(zhì)心(center of mass，CM)是整個(gè)地球的質(zhì)量中心，地表物質(zhì)變遷引起CM相對(duì)于地球參考框架原點(diǎn)的位移稱作地心運(yùn)動(dòng)[1]。CM是地球參考系的原點(diǎn)，建立地心運(yùn)動(dòng)的觀測(cè)模型用于參考框架原點(diǎn)的改正，將會(huì)進(jìn)一步提高空間大地測(cè)量精度。同時(shí)地心運(yùn)動(dòng)蘊(yùn)含了豐富的地球物理信息，在研究時(shí)變重力場(chǎng)等領(lǐng)域中也發(fā)揮著重要作用[2]。故測(cè)定分析地心運(yùn)動(dòng)成為一項(xiàng)必要的工作。

許多學(xué)者利用GPS和SLR等空間大地測(cè)量技術(shù)在測(cè)定地心運(yùn)動(dòng)方面作了大量相關(guān)研究。目前常用的地心運(yùn)動(dòng)反演方法有網(wǎng)平移法[3-4]、動(dòng)力法[5-8]和一階形變法[9-10]等。研究均發(fā)現(xiàn)了地心運(yùn)動(dòng)存在周年和近半年的周期分量。但不同方法求得的振幅存在一定差異，對(duì)于其他周期項(xiàng)的估計(jì)也不完全相同，信號(hào)混淆是造成這種差異的原因之一。由于相對(duì)較低頻的周期信息會(huì)湮沒在噪聲中，影響對(duì)周期的分析及確定[6]，因此也有學(xué)者在分析地心運(yùn)動(dòng)周期性變化之前，利用低通濾波器或小波變換剔除高頻信息的影響，再利用“干凈”的信號(hào)分析地心運(yùn)動(dòng)規(guī)律[5-6]。但這一類方法不具有較強(qiáng)的自適應(yīng)性，會(huì)對(duì)地心運(yùn)動(dòng)的周期分析造成一定影響。

針對(duì)信號(hào)混淆問題，如果可以采用一種自適應(yīng)的盲源信號(hào)分解方法，剔除地心運(yùn)動(dòng)時(shí)間序列中的高頻信息，分析重構(gòu)后的干凈時(shí)序，則可能會(huì)有更好的效果。經(jīng)驗(yàn)?zāi)B(tài)分解方法(empirical mode decomposition，EMD)以其較強(qiáng)的自適應(yīng)性，目前已被廣泛應(yīng)用于時(shí)間序列的信號(hào)分離等領(lǐng)域，并取得了良好的效果。如文獻(xiàn)[11]運(yùn)用EMD方法對(duì)微震信號(hào)降噪，效果明顯；文獻(xiàn)[12]利用EMD方法從GPS時(shí)間序列中分離出降水因素導(dǎo)致的準(zhǔn)兩年周期地表垂直形變；文獻(xiàn)[13]證明EMD方法能有效分離出GPS信號(hào)中混雜的噪聲信號(hào)，從而削弱噪聲對(duì)結(jié)果的影響；文獻(xiàn)[14]運(yùn)用EMD和WD聯(lián)合算法，分析了GPS水汽時(shí)間序列，并探測(cè)到各個(gè)GPS站均存在周年、半周年的周期震蕩。

地心運(yùn)動(dòng)時(shí)序與GPS時(shí)間序列等類似，均混有高頻信息的干擾，因此也可以利用EMD方法對(duì)地心運(yùn)動(dòng)時(shí)序進(jìn)行預(yù)處理，對(duì)相對(duì)“干凈”的時(shí)序進(jìn)行分析。本文在現(xiàn)有研究的基礎(chǔ)上，將EMD方法應(yīng)用于地心運(yùn)動(dòng)時(shí)間序列分析領(lǐng)域；采用網(wǎng)平移法對(duì)IGS提供的GNSS周解進(jìn)行解算得到2007—2017年間地心運(yùn)動(dòng)時(shí)間序列，對(duì)其進(jìn)行分解重構(gòu)，剔除高頻項(xiàng)，并驗(yàn)證該方法的有效性；最后，利用重構(gòu)后的時(shí)間序列對(duì)地心運(yùn)動(dòng)的周期和振幅作進(jìn)一步的分析探討。

1 地心運(yùn)動(dòng)的確定

1.1 網(wǎng)平移法

利用IGS提供的2007—2017年GNSS周解數(shù)據(jù)，采用網(wǎng)平移法計(jì)算其與ITRF2014框架原點(diǎn)的平移量，得到地心運(yùn)動(dòng)時(shí)間序列。本文所使用的是Helmert七參數(shù)轉(zhuǎn)換法，計(jì)算公式如下[7]

(1)

式中，(x，y，z)為測(cè)站在ITRF2014框架下的坐標(biāo)；(X，Y，Z)為測(cè)站坐標(biāo)的GNSS周解；εx、εy、εz為旋轉(zhuǎn)參數(shù)；s為尺度參數(shù)；Tx、Ty、Tz為平移參數(shù)，即地心位移量[15]。為提高計(jì)算精度，本文盡可能均勻地選取核心測(cè)站數(shù)據(jù)進(jìn)行地心運(yùn)動(dòng)反演，并對(duì)粗差進(jìn)行修復(fù)，對(duì)于時(shí)間序列中的缺失部分，采用三次樣條插值方法補(bǔ)全。

1.2 地心運(yùn)動(dòng)解算結(jié)果精度分析

本文求解的地心運(yùn)動(dòng)時(shí)間序列如圖1所示。對(duì)所得地心運(yùn)動(dòng)結(jié)果的精度評(píng)價(jià)主要包括兩方面：一方面Tx、Ty和Tz的標(biāo)準(zhǔn)偏差可作為其內(nèi)符合精度[16]；另一方面美國(guó)得克薩斯大學(xué)空間研究中心(CSR)提供的利用SLR數(shù)據(jù)(跟蹤lageos1/2兩顆衛(wèi)星)進(jìn)行動(dòng)力法解算的地心運(yùn)動(dòng)序列是國(guó)際公認(rèn)的最佳結(jié)果，可以參考此時(shí)間序列評(píng)定計(jì)算結(jié)果的外符合精度。表1給出了地心運(yùn)動(dòng)時(shí)序主要信息，以及內(nèi)、外符合精度統(tǒng)計(jì)。

mm

從表1中可以看出，地心運(yùn)動(dòng)的量級(jí)在3個(gè)方向均為毫米級(jí)，平均值分別為0.27、-2.75和-2.39 mm，內(nèi)符合精度在2.5 mm之內(nèi)，外符合精度在5 mm之內(nèi)，精度較好，并且Tx、Ty方向的內(nèi)、外符合精度均優(yōu)于Tz方向。本文所計(jì)算的地心運(yùn)動(dòng)平均值和內(nèi)外符合精度與CGS(centro de geodasia spaziale)采用幾何法和動(dòng)力法的計(jì)算結(jié)果[17]及文獻(xiàn)[6—7]給出的統(tǒng)計(jì)結(jié)果相比較具有很好的一致性。

2 EMD方法重構(gòu)地心運(yùn)動(dòng)時(shí)序

高頻信息會(huì)影響地心運(yùn)動(dòng)周期的分析及確定，在分析地心運(yùn)動(dòng)之前，應(yīng)抑制噪聲(即高頻部分)的影響[6]。EMD方法是一種處理非線性非平穩(wěn)信號(hào)的時(shí)頻分析方法[18]，它基于信號(hào)本身自適應(yīng)地從高頻到低頻逐次分解，不需要任何的先驗(yàn)值，即可獲得一組固有模態(tài)函數(shù)(intrinsic mode function，IMF)分量，這些模態(tài)函數(shù)能夠很好地反映信號(hào)在任何時(shí)間局部的頻率特征。EMD的基本思路是把時(shí)間序列x(t)根據(jù)原序列自身特征分解成若干個(gè)模態(tài)分量IMF和一個(gè)趨勢(shì)項(xiàng)，表達(dá)式為

(2)

式中，Tk(t)為趨勢(shì)項(xiàng)；IMFk為模態(tài)分量。

為了得到相對(duì)“干凈”的時(shí)間序列，本文利用EMD方法按如下步驟對(duì)時(shí)序進(jìn)行分解和重構(gòu)[11，13]：

(1) 將原始時(shí)序分解為k個(gè)IMF分量，并計(jì)算各IMF分量與原始序列的相關(guān)系數(shù)。

(2) 利用相關(guān)系數(shù)局部最小值原則確定高頻項(xiàng)分界值，即一組相關(guān)系數(shù)中第一個(gè)取得局部最小值對(duì)應(yīng)的IMF，記為IMFn。

(3) IMFn及IMFn之前的分量為高頻項(xiàng)，將剩下的低頻項(xiàng)(即周期項(xiàng)和趨勢(shì)項(xiàng)之和)通過重構(gòu)方法獲得一組“干凈”的時(shí)間序列。

這種序列重構(gòu)方法無需設(shè)定經(jīng)驗(yàn)值判定高頻項(xiàng)的分界值，避免了人為因素造成的誤差，具有自適應(yīng)性。各個(gè)方向?qū)?yīng)各IMF分量與原始序列的相關(guān)系數(shù)如圖2所示，從圖中可以看出Tx、Ty和Tz方向的分界層均集中在第2個(gè)IMF分量。Tx方向上與原序列一致性較好的是第4個(gè)IMF分量，相關(guān)系數(shù)接近0.9，Ty和Tz方向上與原序列一致性較好的是第4、7個(gè)IMF分量，相關(guān)系數(shù)在0.6以上。同時(shí)，每一個(gè)方向上的IMF分量的個(gè)數(shù)不同，Ty和Tz方向的IMF分量個(gè)數(shù)為7，而Tx方向有8個(gè)IMF分量，這說明EMD方法是根據(jù)原始序列的自身特征進(jìn)行分解。

根據(jù)分界層及分界層之前的IMF為高頻項(xiàng)的原則，對(duì)IMF3及之后的分量進(jìn)行重構(gòu)。Tx、Ty和Tz3個(gè)方向重構(gòu)序列分別如圖3(a)—(c)所示，其中圖Ⅰ為重構(gòu)前的時(shí)間序列和重構(gòu)的低頻項(xiàng)，圖Ⅱ?yàn)榉蛛x出的高頻項(xiàng)。

由圖3可知，原始序列與低頻項(xiàng)具有良好的一致性，其包含的周期信息更為明顯，這會(huì)使得周期振幅的分析結(jié)果更加可靠。去除高頻項(xiàng)后Tx、Ty和Tz方向的振幅范圍分別為：-4～4 mm、-7～1 mm、-8～2 mm。分離出的高頻項(xiàng)(即IMF1和IMF2分量之和)均在0附近隨機(jī)波動(dòng)，周期性特征并不明顯。

為了驗(yàn)證EMD方法抑制噪聲的效果，可以通過計(jì)算重構(gòu)時(shí)間序列占原序列的能量百分比[13]，以及兩者之間的相關(guān)系數(shù)這兩個(gè)參數(shù)作為評(píng)價(jià)利用EMD方法重構(gòu)地心運(yùn)動(dòng)時(shí)間序列效果的標(biāo)準(zhǔn)。表2給出了3個(gè)方向EMD重構(gòu)時(shí)間序列對(duì)應(yīng)的相關(guān)系數(shù)和剩余能量百分比統(tǒng)計(jì)結(jié)果。

表2 EMD重構(gòu)時(shí)間序列的評(píng)價(jià)參數(shù)統(tǒng)計(jì) (%)

由表2可以看出，重構(gòu)時(shí)序與原序列的相關(guān)系數(shù)在95%左右，說明兩者的符合度較好，保留了序列中較多的有效信息。重構(gòu)后時(shí)間序列的剩余能量百分比均在90%以上，說明被剔除掉的基本上是能量較低的高頻信息。綜上所述，在利用基于EMD方法重構(gòu)得到的時(shí)間序列進(jìn)行地心運(yùn)動(dòng)結(jié)果分析時(shí)，可以在一定程度上減小高頻信息的混入對(duì)結(jié)果的影響。

3 地心運(yùn)動(dòng)時(shí)間序列分析

對(duì)重構(gòu)的2007—2017年地心運(yùn)動(dòng)時(shí)序進(jìn)行快速傅里葉變換，計(jì)算功率譜密度，由此識(shí)別周期，并計(jì)算每個(gè)周期的貢獻(xiàn)率大小。圖4給出了Tx、Ty和Tz3個(gè)方向的功率譜密度圖(采樣頻率為7 d)。表3給出了使用EMD方法重構(gòu)地心運(yùn)動(dòng)變化時(shí)間序列前后各個(gè)周期的貢獻(xiàn)率統(tǒng)計(jì)結(jié)果。

從圖4中可以看到，3個(gè)方向上在0.019 Hz上(對(duì)應(yīng)一年的周期)均有較高的能量出現(xiàn)，說明3個(gè)方向存在明顯的周年變化；Ty和Tz方向在低頻部分也有較大的能量，這說明Ty和Tz還存在較明顯的長(zhǎng)期變化項(xiàng)，并且Tz方向更為明顯。從表3可以看出，在利用EMD方法重構(gòu)時(shí)間序列后，各主要周期的貢獻(xiàn)率均有所提高，使得周期性特征的識(shí)別更加明顯準(zhǔn)確。經(jīng)計(jì)算，Tx、Ty和Tz方向分別提高了12.3%、16.7%、6.3%。另外，Tz方向上的長(zhǎng)周期變化趨勢(shì)較Tx、Ty方向更為明顯，其趨勢(shì)項(xiàng)的貢獻(xiàn)率在各周期貢獻(xiàn)率中最高。

表3 利用EMD重構(gòu)時(shí)間序列前后的周期貢獻(xiàn)率統(tǒng)計(jì)

為了得到各個(gè)周期對(duì)應(yīng)振幅和相位，結(jié)合最小二乘方法對(duì)重構(gòu)后的時(shí)間序列進(jìn)行擬合[5]

(3)

式中，x(t)為t時(shí)刻的地心運(yùn)動(dòng)；a為常數(shù)項(xiàng)；b為趨勢(shì)項(xiàng)；t0為初始時(shí)間；A為振幅；φ為相位；m為周期個(gè)數(shù)；T為周期；A和φ為調(diào)和系數(shù)。地心運(yùn)動(dòng)周期性變化振幅和相位的擬合結(jié)果見表4。

表4 地心運(yùn)動(dòng)周期性變化的振幅和相位

現(xiàn)有研究表明，地心運(yùn)動(dòng)存在明顯的周年項(xiàng)和相對(duì)較弱的半年項(xiàng)[9]，這種變化主要與大氣、海洋等質(zhì)量再分布有關(guān)[1]。從表4中可以看出，周年項(xiàng)與半周年項(xiàng)對(duì)應(yīng)的振幅均為毫米級(jí)，并且周年項(xiàng)的振幅大于半年振幅，這與已有研究成果基本一致[2，3,6,19]。Tx、Ty和Tz方向的周年項(xiàng)振幅均為各周期對(duì)應(yīng)振幅的最大值，分別為2.32、1.89和2.07 mm。但Tz方向中周期貢獻(xiàn)率最大的是4004 d，接近于本文的研究時(shí)間跨度，這說明了Tz方向的趨勢(shì)項(xiàng)較明顯，但Tz方向的周年變化貢獻(xiàn)率與長(zhǎng)期趨勢(shì)的貢獻(xiàn)率僅相差1.5%左右，因此，可以認(rèn)為在Tz方向上也存在比較明顯的周年變化。另外，從表3中可以看出，3個(gè)方向的半年項(xiàng)均較弱，其中Tx和Ty方向上還存在4～6個(gè)月的震蕩。研究表明，半年項(xiàng)周期具有時(shí)變性[6]，這一特性的存在可能導(dǎo)致半年項(xiàng)的能量分布較為離散，從而導(dǎo)致這種周期震蕩。

表5列出了本文及相關(guān)學(xué)者求得的地心運(yùn)動(dòng)長(zhǎng)期變化。相比于周期性運(yùn)動(dòng)，長(zhǎng)期運(yùn)動(dòng)速度較小，Tx、Ty和Tz方向分別為-0.02、0.13和-0.27 mm/a。本文結(jié)果與現(xiàn)有研究基本一致，3個(gè)方向上的趨勢(shì)性變化基本都在1 mm/a之內(nèi)。另外，趨勢(shì)項(xiàng)的絕對(duì)值在Tz方向上最大。

表5 地心運(yùn)動(dòng)長(zhǎng)期變化 (mm/a)

表6為學(xué)者根據(jù)不同研究方法得出的振幅和相位與本文所得結(jié)果的比較。需要說明的一點(diǎn)是，由于Tx和Ty方向上的半年項(xiàng)都具有時(shí)變性，因此在統(tǒng)計(jì)這兩個(gè)方向上的半年項(xiàng)時(shí)，本文根據(jù)表4選取4～8個(gè)月長(zhǎng)度的周期中振幅最大的一項(xiàng)，作為半年項(xiàng)的統(tǒng)計(jì)結(jié)果，即Tx方向?yàn)?21 d，Ty方向?yàn)?82 d。

表6 不同研究得出地心運(yùn)動(dòng)周期性變化的振幅和相位

除了周年項(xiàng)和半周年項(xiàng)的變化，以及長(zhǎng)期的運(yùn)動(dòng)趨勢(shì)外，本文還探測(cè)到了其他的周期信息。Tx和Ty方向有800 d的周期，但Tz方向不存在；Tx和Tz方向存在1334 d的周期性變化，而Ty方向則不明顯。這類周期變化被稱為年際變化[6]，其可能與地表水、大氣、海平面和地幔對(duì)流等的不規(guī)則變化有關(guān)。

4 結(jié) 語

本文針對(duì)IGS提供的GNSS周解數(shù)據(jù)，采用Helmert七參數(shù)法計(jì)算了2007—2017年的地心運(yùn)動(dòng)時(shí)間序列，其內(nèi)外符合精度均優(yōu)于5 mm，與已有研究具有良好的一致性。采用EMD方法對(duì)時(shí)間序列進(jìn)行分解重構(gòu)，發(fā)現(xiàn)重構(gòu)時(shí)序的高頻信息得到有效抑制，周期貢獻(xiàn)率有所提高。利用該序列求得的地心運(yùn)動(dòng)在Tx、Ty和Tz3個(gè)方向的周年振幅分別為2.32、1.89和2.07 mm；半年項(xiàng)較小，且在Tx和Ty方向上具有時(shí)變性；此外探測(cè)發(fā)現(xiàn)還存在一些其他年際變化及小于1 mm/a的長(zhǎng)期運(yùn)動(dòng)趨勢(shì)。