中國區(qū)域內(nèi)IGS站時間序列的非線性變化分析

戴海亮,孫付平,朱新慧

(信息工程大學(xué) 地理空間信息學(xué)院,河南 鄭州 450001)

0 引 言

隨著 GPS 觀測數(shù)據(jù)的不斷豐富,分析測站時間序列中的非線性變化可以獲得各種地球物理現(xiàn)象和季節(jié)性變化規(guī)律,如果對這部分殘差規(guī)律進行建模擬合從而消除其影響,那么可以進一步提高GPS臺站坐標(biāo)的精度[1-2].時間序列數(shù)據(jù)不僅可以反映出臺站的趨勢變化,而且也可以反映出臺站存在著的非線性變化.研究表明,幾乎所有IGS站(尤其是高程方向)都呈現(xiàn)顯著的非線性運動趨勢,且振幅可達1～2 cm[3-6].國內(nèi)外學(xué)者針對測站時間序列中的季節(jié)性變化規(guī)律進行了廣泛的研究,取得了豐富的成果.文獻[7]對加利福尼亞南部和內(nèi)華達州南部GPS位移測量中的時間序列研究發(fā)現(xiàn),時間序列中的季節(jié)性信號通常表現(xiàn)為周年或半周年形式,并能夠以常振幅、常相位的諧波模型來描述.文獻[8]和文獻[9]對季節(jié)性信號建模時,證實用周年加半周年的諧波模型并不能很好地描述GPS測站時間序列中的季節(jié)性變化,還要考慮一些其他長周期和短周期季節(jié)性信號的影響.文獻[10]通過不同的測站時間序列數(shù)據(jù)實驗比較發(fā)現(xiàn),三角函數(shù)法能更加直觀地體現(xiàn)出非線性變化的規(guī)律, 對于進一步分析非線性變化規(guī)律的性質(zhì)與機制具有重要的參考價值.通過對ITRF2008測站時間序列的分析發(fā)現(xiàn)時間序列中包含了大量的周期性規(guī)律,所以,僅以周年加半周年或者其他幾個主要周期項對時間序列中的周期性規(guī)律進行描述顯然是不夠的.

本文以ITRF網(wǎng)站上提供的 GPS 臺站時間序列為研究對象,在對時間序列數(shù)據(jù)預(yù)處理的基礎(chǔ)上,對中國區(qū)域內(nèi)8個IGS站時間序列的非線性變化規(guī)律進行頻譜分析,提取了臺站時間序列可能存在的周期并根據(jù)三角函數(shù)法進行了建模研究.

1 數(shù)據(jù)獲取

實驗采用的數(shù)據(jù)來源于國際地球參考框架ITRF2008的GPS坐標(biāo)殘差數(shù)據(jù),可在官網(wǎng)上(http://itrf.ensg.ign.fr/ITRF_solutions/2008/)下載.絕大部分GPS時間序列包含了1997年至2009年約12 a的殘差數(shù)據(jù),采樣間隔為7 d,而且殘差序列已經(jīng)剔除了固體潮、極潮等部分環(huán)境負荷影響.在中國區(qū)域內(nèi),由于南北跨緯度廣,而且氣候復(fù)雜多樣,測站所處的環(huán)境差異較大.本文根據(jù)測站分布狀況和穩(wěn)定程度,選取了國內(nèi)不同區(qū)域的8個GPS站,通過分析不同地區(qū)的臺站坐標(biāo)殘差序列的規(guī)律,在頻譜分析的基礎(chǔ)上提取時間序列的周期項和振幅,然后對非線性變化的規(guī)律進行分析.鑒于篇幅,下面以BJFS站為例,其測站東向、垂向和北向原始殘差時間序列如圖1所示.

2 測站時間序列分析

2.1 預(yù)處理

多路徑效應(yīng)、軌道異常以及測站相關(guān)的誤差等因素都會導(dǎo)致GPS坐標(biāo)時間序列存在粗差[11].通常采用兩種方法來剔除粗差:一種是基于最小二乘(LS)的標(biāo)準(zhǔn)化殘差,采用假設(shè)檢驗進行粗差探測;另一種是根據(jù)“3σ”準(zhǔn)則對粗差進行識別并剔除[12].本文采用后者進行粗差的探測和剔除,然后利用鄰近點取平均值的方法來擬合產(chǎn)生粗差的歷元.

在一般情況下,ITRF網(wǎng)站提供的GPS時間序列正常的采樣間隔是7 d,但是由于一些客觀因素的影響導(dǎo)致時間序列不能均勻采樣,有的甚至存在間隔數(shù)星期甚至更長的情況,使得坐標(biāo)時間序列存在間斷點.對于測站時間序列存在間斷點的情況,文獻[1]和文獻[13]在進行頻譜分析時都作了插值處理.同樣,本文考慮到處理效率,以及頻譜分析時要求數(shù)據(jù)均勻采樣且具備零均值特性,所以在實驗時采用了三次埃爾米特多項式的方法對測站時間序列中的間斷點進行了插值處理.這種插值方法中插值函數(shù)及其一階導(dǎo)數(shù)都是連續(xù)的,所以插值結(jié)果比較光滑,而且速度也比三次樣條插值以及分段線性插值要快.

2.2 線性擬合去趨勢項

對測站時間序列進行頻譜分析時要求原始信號為平穩(wěn)信號,即要求原始數(shù)據(jù)為零均值.因此,需要對時間序列數(shù)據(jù)進行零均值化處理,扣除其趨勢項,得到變化平穩(wěn)的時間序列,擬合模型為

yt=a+bt+xt.

(1)

式中:yt為原始時間序列;xt為去除常數(shù)項和線性趨勢項之后的時間序列.線性擬合求得常數(shù)項和一次項系數(shù)后,在原始時間序列中減去擬合值,剩下的殘差就是可以用于功率譜分析的新時間序列[14].從數(shù)據(jù)中刪除趨勢可以將分析集中在數(shù)據(jù)趨勢本身的波動上.

2.3 傅里葉變換提取周期

傅里葉變換是一種研究信號的頻譜分析方法,把時域和頻域聯(lián)系在了一起.它可以將一個時域信號分解為多個頻域信號,而眾多的頻域信號又可以準(zhǔn)確無誤地重構(gòu)原來的時域信號.這種變換是可逆的且能量保持不變[15].1965 年,美國工程師 Cooley 和 Tukey提出了快速傅里葉變換(FFT)的概念.快速傅里葉變換可以有效地從給定的信號中找出頻率與能量的關(guān)系,而這種關(guān)系是通過轉(zhuǎn)換信號時間域的計算關(guān)系得到的,其定義如下式:

(2)

其反變換為

(3)

式中:f(t)為時間域的信號;F(jω)為快速傅里葉變換的結(jié)果,稱為函數(shù)f(t)的頻譜密度函數(shù).

如圖2所示, BJFS站垂直方向坐標(biāo)時間序列的頻譜圖中,振幅最高點的頻率是3.255e-8 Hz,振幅是4.415 mm,可以得出所對應(yīng)的周期為355.6天,近似為一周年.另外,在周年項附近還存在著振幅為2.0 mm,周期為0.5年以及振幅為1.7 mm,周期為10年的周期項.

本文通過對選定的中國區(qū)域內(nèi)8個GPS站的殘差時間序列進行分析發(fā)現(xiàn),測站時間序列中很少存在嚴(yán)格的周年項、半周年項,主要是與之相近的類周年項、類半年項.所以本文在提取周期項時,將時間介于350～380 天之間的周期項都視為周年項,將時間介于170～200天之間的周期項都視為半周年項,以此類推,提取時間序列中可能存在的周期.鑒于篇幅,將振幅影響較大的幾個主要周期項及其振幅統(tǒng)計如表1所示.

表1 中國區(qū)域內(nèi)GPS站主要周期項及振幅統(tǒng)計表

從表1可以看出,中國區(qū)域內(nèi)大部分GPS測站東向、垂向和北向時間序列的周期規(guī)律中,周年項較為普遍,表現(xiàn)出來的振幅影響也最大,因此是其主要周期項;其次是半年周期項,然后是2～3 a周期項,其他周期項規(guī)律如10 a、13 a等則不具有普遍性;個別測站如XIAN站的3個方向上還表現(xiàn)出12 a以上的長周期項為主要周期規(guī)律.同時,還發(fā)現(xiàn)每個測站各個分量上表現(xiàn)出來的周期規(guī)律不盡相同,而且不同測站的周期規(guī)律也存在差異,究其原因,與它們所處的地理位置和環(huán)境氣候等因素有著不可分割的關(guān)系.

2.4 三角函數(shù)擬合

根據(jù)諧波分析的基本思想,包含非線性變化的測站坐標(biāo)時間序列可以看成是由一系列正余弦(sinkx、coskx,k=0,1,2,…)函數(shù)構(gòu)成的.若將去趨勢項后的時間序列X(t)看作是由不同頻率fi的余弦波疊加而成,則公式可以表示為

②扎賚特旗水文地質(zhì)條件相對較好，含水層顆粒粗大，厚度較大，水量相對較豐富，水質(zhì)較好，具備實施節(jié)水改造、推廣高效節(jié)水灌溉的技術(shù)條件。

(3)

式中:k、Ai、fi(i=0,1,2,…，k)都為常數(shù);k為時間序列中周期項的個數(shù);Ai、fi、φi分別為第i個周期波的振幅、頻率與初始相位;φi為在(-π,π)內(nèi)均勻分布的獨立隨機變量,ε(t)為三角函數(shù)擬合后的殘差,對式(3)系數(shù)化可得:

(4)

此外,在進行數(shù)據(jù)預(yù)處理時,剔除粗差采用的是“3σ”準(zhǔn)則,當(dāng)粗差量級較大時,殘差的加權(quán)平方和也往往偏大,導(dǎo)致“漏判”量級較小的粗差.實驗證明使用該方法時,粗差剔除率僅為0.1%左右,這可能是影響周期項提取和擬合精度的因素之一.而且實驗在數(shù)據(jù)處理時默認的采樣間隔是7 d,實際上在對采樣時間進行差分時,還存在一些6 d、8 d的采樣間隔,并不是完全意義上的7 d,有的還存在間隔數(shù)星期甚至更長的情況,使得坐標(biāo)時間序列存在間斷點.如在對BJFS站的采樣時間進行差分時,采樣間隔大部分是7 d以及在7 d上下波動,而不是嚴(yán)格意義上的7 d,并且在210周左右時出現(xiàn)了30多天的間斷,如圖4所示.有的測站殘差時間序列間斷時間更長,因此考慮到處理效率,需要對原始殘差時間序列數(shù)據(jù)進行預(yù)處理.另外,在對間斷點進行插值時,選擇的插值方法等都會對周期項的提取以及擬合精度產(chǎn)生一定的影響.

3 結(jié)束語

本文通過對國際地球參考框架ITRF2008在中國區(qū)域內(nèi)的8個IGS站的時間序列進行分析,提取了時間序列中可能存在的周期項,并根據(jù)三角函數(shù)法對周期項進行建模分析,得出了如下結(jié)論:

1)周年項具有普遍性,半年項也具有一定的普遍性,影響比周年項稍小,其他周期規(guī)律如3 a、10 a、13 a等周期項則不具有普遍性,這與測站所處的地理位置和環(huán)境氣候等因素有很大的關(guān)系.

2)在測站東向、北向和垂向3個方向的周期規(guī)律中,東向和北向周期規(guī)律振幅較小,一般都在cm級以內(nèi),而垂向上周期規(guī)律最為明顯,振幅也最大,主要原因是測站垂直方向受到多種地球物理因素(如未模型化的固體潮、海洋潮汐以及大氣潮引起的荷載等周期性變化的影響)疊加的結(jié)果.

3)每個測站各個分量的周期規(guī)律不盡相同,而且各個測站的周期項也存在差異,因此,難以建立一個統(tǒng)一的周期模型來對所有測站的時間序列進行擬合.而且使用三角函數(shù)只對時間序列中的季節(jié)性因素進行了擬合,還有一些非季節(jié)性因素(如地震、火山爆發(fā)等較大的自然現(xiàn)象)也不容忽視,因此研究如何建立一個最佳的擬合模型對于建立和維持高精度地球參考框架是個關(guān)鍵問題.