中國區域內IGS站時間序列的非線性變化分析

戴海亮,孫付平,朱新慧

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

0 引 言

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

本文以ITRF網站上提供的 GPS 臺站時間序列為研究對象,在對時間序列數據預處理的基礎上,對中國區域內8個IGS站時間序列的非線性變化規律進行頻譜分析,提取了臺站時間序列可能存在的周期并根據三角函數法進行了建模研究.

1 數據獲取

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

2 測站時間序列分析

2.1 預處理

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

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

2.2 線性擬合去趨勢項

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

yt=a+bt+xt.

(1)

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

2.3 傅里葉變換提取周期

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

(2)

其反變換為

(3)

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

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

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

表1 中國區域內GPS站主要周期項及振幅統計表

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

2.4 三角函數擬合

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

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

(3)

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

(4)

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

3 結束語

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

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

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

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