戴海亮,孫付平,朱新慧
(信息工程大學(xué) 地理空間信息學(xué)院,河南 鄭州 450001)
隨著 GPS 觀測數(shù)據(jù)的不斷豐富,分析測站時(shí)間序列中的非線性變化可以獲得各種地球物理現(xiàn)象和季節(jié)性變化規(guī)律,如果對(duì)這部分殘差規(guī)律進(jìn)行建模擬合從而消除其影響,那么可以進(jìn)一步提高GPS臺(tái)站坐標(biāo)的精度[1-2].時(shí)間序列數(shù)據(jù)不僅可以反映出臺(tái)站的趨勢變化,而且也可以反映出臺(tái)站存在著的非線性變化.研究表明,幾乎所有IGS站(尤其是高程方向)都呈現(xiàn)顯著的非線性運(yùn)動(dòng)趨勢,且振幅可達(dá)1~2 cm[3-6].國內(nèi)外學(xué)者針對(duì)測站時(shí)間序列中的季節(jié)性變化規(guī)律進(jìn)行了廣泛的研究,取得了豐富的成果.文獻(xiàn)[7]對(duì)加利福尼亞南部和內(nèi)華達(dá)州南部GPS位移測量中的時(shí)間序列研究發(fā)現(xiàn),時(shí)間序列中的季節(jié)性信號(hào)通常表現(xiàn)為周年或半周年形式,并能夠以常振幅、常相位的諧波模型來描述.文獻(xiàn)[8]和文獻(xiàn)[9]對(duì)季節(jié)性信號(hào)建模時(shí),證實(shí)用周年加半周年的諧波模型并不能很好地描述GPS測站時(shí)間序列中的季節(jié)性變化,還要考慮一些其他長周期和短周期季節(jié)性信號(hào)的影響.文獻(xiàn)[10]通過不同的測站時(shí)間序列數(shù)據(jù)實(shí)驗(yàn)比較發(fā)現(xiàn),三角函數(shù)法能更加直觀地體現(xiàn)出非線性變化的規(guī)律, 對(duì)于進(jìn)一步分析非線性變化規(guī)律的性質(zhì)與機(jī)制具有重要的參考價(jià)值.通過對(duì)ITRF2008測站時(shí)間序列的分析發(fā)現(xiàn)時(shí)間序列中包含了大量的周期性規(guī)律,所以,僅以周年加半周年或者其他幾個(gè)主要周期項(xiàng)對(duì)時(shí)間序列中的周期性規(guī)律進(jìn)行描述顯然是不夠的.
本文以ITRF網(wǎng)站上提供的 GPS 臺(tái)站時(shí)間序列為研究對(duì)象,在對(duì)時(shí)間序列數(shù)據(jù)預(yù)處理的基礎(chǔ)上,對(duì)中國區(qū)域內(nèi)8個(gè)IGS站時(shí)間序列的非線性變化規(guī)律進(jìn)行頻譜分析,提取了臺(tái)站時(shí)間序列可能存在的周期并根據(jù)三角函數(shù)法進(jìn)行了建模研究.
實(shí)驗(yàn)采用的數(shù)據(jù)來源于國際地球參考框架ITRF2008的GPS坐標(biāo)殘差數(shù)據(jù),可在官網(wǎng)上(http://itrf.ensg.ign.fr/ITRF_solutions/2008/)下載.絕大部分GPS時(shí)間序列包含了1997年至2009年約12 a的殘差數(shù)據(jù),采樣間隔為7 d,而且殘差序列已經(jīng)剔除了固體潮、極潮等部分環(huán)境負(fù)荷影響.在中國區(qū)域內(nèi),由于南北跨緯度廣,而且氣候復(fù)雜多樣,測站所處的環(huán)境差異較大.本文根據(jù)測站分布狀況和穩(wěn)定程度,選取了國內(nèi)不同區(qū)域的8個(gè)GPS站,通過分析不同地區(qū)的臺(tái)站坐標(biāo)殘差序列的規(guī)律,在頻譜分析的基礎(chǔ)上提取時(shí)間序列的周期項(xiàng)和振幅,然后對(duì)非線性變化的規(guī)律進(jìn)行分析.鑒于篇幅,下面以BJFS站為例,其測站東向、垂向和北向原始?xì)埐顣r(shí)間序列如圖1所示.
多路徑效應(yīng)、軌道異常以及測站相關(guān)的誤差等因素都會(huì)導(dǎo)致GPS坐標(biāo)時(shí)間序列存在粗差[11].通常采用兩種方法來剔除粗差:一種是基于最小二乘(LS)的標(biāo)準(zhǔn)化殘差,采用假設(shè)檢驗(yàn)進(jìn)行粗差探測;另一種是根據(jù)“3σ”準(zhǔn)則對(duì)粗差進(jìn)行識(shí)別并剔除[12].本文采用后者進(jìn)行粗差的探測和剔除,然后利用鄰近點(diǎn)取平均值的方法來擬合產(chǎn)生粗差的歷元.
在一般情況下,ITRF網(wǎng)站提供的GPS時(shí)間序列正常的采樣間隔是7 d,但是由于一些客觀因素的影響導(dǎo)致時(shí)間序列不能均勻采樣,有的甚至存在間隔數(shù)星期甚至更長的情況,使得坐標(biāo)時(shí)間序列存在間斷點(diǎn).對(duì)于測站時(shí)間序列存在間斷點(diǎn)的情況,文獻(xiàn)[1]和文獻(xiàn)[13]在進(jìn)行頻譜分析時(shí)都作了插值處理.同樣,本文考慮到處理效率,以及頻譜分析時(shí)要求數(shù)據(jù)均勻采樣且具備零均值特性,所以在實(shí)驗(yàn)時(shí)采用了三次埃爾米特多項(xiàng)式的方法對(duì)測站時(shí)間序列中的間斷點(diǎn)進(jìn)行了插值處理.這種插值方法中插值函數(shù)及其一階導(dǎo)數(shù)都是連續(xù)的,所以插值結(jié)果比較光滑,而且速度也比三次樣條插值以及分段線性插值要快.
對(duì)測站時(shí)間序列進(jìn)行頻譜分析時(shí)要求原始信號(hào)為平穩(wěn)信號(hào),即要求原始數(shù)據(jù)為零均值.因此,需要對(duì)時(shí)間序列數(shù)據(jù)進(jìn)行零均值化處理,扣除其趨勢項(xiàng),得到變化平穩(wěn)的時(shí)間序列,擬合模型為
yt=a+bt+xt.
(1)
式中:yt為原始時(shí)間序列;xt為去除常數(shù)項(xiàng)和線性趨勢項(xiàng)之后的時(shí)間序列.線性擬合求得常數(shù)項(xiàng)和一次項(xiàng)系數(shù)后,在原始時(shí)間序列中減去擬合值,剩下的殘差就是可以用于功率譜分析的新時(shí)間序列[14].從數(shù)據(jù)中刪除趨勢可以將分析集中在數(shù)據(jù)趨勢本身的波動(dòng)上.
傅里葉變換是一種研究信號(hào)的頻譜分析方法,把時(shí)域和頻域聯(lián)系在了一起.它可以將一個(gè)時(shí)域信號(hào)分解為多個(gè)頻域信號(hào),而眾多的頻域信號(hào)又可以準(zhǔn)確無誤地重構(gòu)原來的時(shí)域信號(hào).這種變換是可逆的且能量保持不變[15].1965 年,美國工程師 Cooley 和 Tukey提出了快速傅里葉變換(FFT)的概念.快速傅里葉變換可以有效地從給定的信號(hào)中找出頻率與能量的關(guān)系,而這種關(guān)系是通過轉(zhuǎn)換信號(hào)時(shí)間域的計(jì)算關(guān)系得到的,其定義如下式:
(2)
其反變換為
(3)
式中:f(t)為時(shí)間域的信號(hào);F(jω)為快速傅里葉變換的結(jié)果,稱為函數(shù)f(t)的頻譜密度函數(shù).
如圖2所示, BJFS站垂直方向坐標(biāo)時(shí)間序列的頻譜圖中,振幅最高點(diǎn)的頻率是3.255e-8 Hz,振幅是4.415 mm,可以得出所對(duì)應(yīng)的周期為355.6天,近似為一周年.另外,在周年項(xiàng)附近還存在著振幅為2.0 mm,周期為0.5年以及振幅為1.7 mm,周期為10年的周期項(xiàng).
本文通過對(duì)選定的中國區(qū)域內(nèi)8個(gè)GPS站的殘差時(shí)間序列進(jìn)行分析發(fā)現(xiàn),測站時(shí)間序列中很少存在嚴(yán)格的周年項(xiàng)、半周年項(xiàng),主要是與之相近的類周年項(xiàng)、類半年項(xiàng).所以本文在提取周期項(xiàng)時(shí),將時(shí)間介于350~380 天之間的周期項(xiàng)都視為周年項(xiàng),將時(shí)間介于170~200天之間的周期項(xiàng)都視為半周年項(xiàng),以此類推,提取時(shí)間序列中可能存在的周期.鑒于篇幅,將振幅影響較大的幾個(gè)主要周期項(xiàng)及其振幅統(tǒng)計(jì)如表1所示.
表1 中國區(qū)域內(nèi)GPS站主要周期項(xiàng)及振幅統(tǒng)計(jì)表
從表1可以看出,中國區(qū)域內(nèi)大部分GPS測站東向、垂向和北向時(shí)間序列的周期規(guī)律中,周年項(xiàng)較為普遍,表現(xiàn)出來的振幅影響也最大,因此是其主要周期項(xiàng);其次是半年周期項(xiàng),然后是2~3 a周期項(xiàng),其他周期項(xiàng)規(guī)律如10 a、13 a等則不具有普遍性;個(gè)別測站如XIAN站的3個(gè)方向上還表現(xiàn)出12 a以上的長周期項(xiàng)為主要周期規(guī)律.同時(shí),還發(fā)現(xiàn)每個(gè)測站各個(gè)分量上表現(xiàn)出來的周期規(guī)律不盡相同,而且不同測站的周期規(guī)律也存在差異,究其原因,與它們所處的地理位置和環(huán)境氣候等因素有著不可分割的關(guān)系.
根據(jù)諧波分析的基本思想,包含非線性變化的測站坐標(biāo)時(shí)間序列可以看成是由一系列正余弦(sinkx、coskx,k=0,1,2,…)函數(shù)構(gòu)成的.若將去趨勢項(xiàng)后的時(shí)間序列X(t)看作是由不同頻率fi的余弦波疊加而成,則公式可以表示為
②扎賚特旗水文地質(zhì)條件相對(duì)較好,含水層顆粒粗大,厚度較大,水量相對(duì)較豐富,水質(zhì)較好,具備實(shí)施節(jié)水改造、推廣高效節(jié)水灌溉的技術(shù)條件。
(3)
式中:k、Ai、fi(i=0,1,2,…,k)都為常數(shù);k為時(shí)間序列中周期項(xiàng)的個(gè)數(shù);Ai、fi、φi分別為第i個(gè)周期波的振幅、頻率與初始相位;φi為在(-π,π)內(nèi)均勻分布的獨(dú)立隨機(jī)變量,ε(t)為三角函數(shù)擬合后的殘差,對(duì)式(3)系數(shù)化可得:
(4)
此外,在進(jìn)行數(shù)據(jù)預(yù)處理時(shí),剔除粗差采用的是“3σ”準(zhǔn)則,當(dāng)粗差量級(jí)較大時(shí),殘差的加權(quán)平方和也往往偏大,導(dǎo)致“漏判”量級(jí)較小的粗差.實(shí)驗(yàn)證明使用該方法時(shí),粗差剔除率僅為0.1%左右,這可能是影響周期項(xiàng)提取和擬合精度的因素之一.而且實(shí)驗(yàn)在數(shù)據(jù)處理時(shí)默認(rèn)的采樣間隔是7 d,實(shí)際上在對(duì)采樣時(shí)間進(jìn)行差分時(shí),還存在一些6 d、8 d的采樣間隔,并不是完全意義上的7 d,有的還存在間隔數(shù)星期甚至更長的情況,使得坐標(biāo)時(shí)間序列存在間斷點(diǎn).如在對(duì)BJFS站的采樣時(shí)間進(jìn)行差分時(shí),采樣間隔大部分是7 d以及在7 d上下波動(dòng),而不是嚴(yán)格意義上的7 d,并且在210周左右時(shí)出現(xiàn)了30多天的間斷,如圖4所示.有的測站殘差時(shí)間序列間斷時(shí)間更長,因此考慮到處理效率,需要對(duì)原始?xì)埐顣r(shí)間序列數(shù)據(jù)進(jìn)行預(yù)處理.另外,在對(duì)間斷點(diǎn)進(jìn)行插值時(shí),選擇的插值方法等都會(huì)對(duì)周期項(xiàng)的提取以及擬合精度產(chǎn)生一定的影響.
本文通過對(duì)國際地球參考框架ITRF2008在中國區(qū)域內(nèi)的8個(gè)IGS站的時(shí)間序列進(jìn)行分析,提取了時(shí)間序列中可能存在的周期項(xiàng),并根據(jù)三角函數(shù)法對(duì)周期項(xiàng)進(jìn)行建模分析,得出了如下結(jié)論:
1)周年項(xiàng)具有普遍性,半年項(xiàng)也具有一定的普遍性,影響比周年項(xiàng)稍小,其他周期規(guī)律如3 a、10 a、13 a等周期項(xiàng)則不具有普遍性,這與測站所處的地理位置和環(huán)境氣候等因素有很大的關(guān)系.
2)在測站東向、北向和垂向3個(gè)方向的周期規(guī)律中,東向和北向周期規(guī)律振幅較小,一般都在cm級(jí)以內(nèi),而垂向上周期規(guī)律最為明顯,振幅也最大,主要原因是測站垂直方向受到多種地球物理因素(如未模型化的固體潮、海洋潮汐以及大氣潮引起的荷載等周期性變化的影響)疊加的結(jié)果.
3)每個(gè)測站各個(gè)分量的周期規(guī)律不盡相同,而且各個(gè)測站的周期項(xiàng)也存在差異,因此,難以建立一個(gè)統(tǒng)一的周期模型來對(duì)所有測站的時(shí)間序列進(jìn)行擬合.而且使用三角函數(shù)只對(duì)時(shí)間序列中的季節(jié)性因素進(jìn)行了擬合,還有一些非季節(jié)性因素(如地震、火山爆發(fā)等較大的自然現(xiàn)象)也不容忽視,因此研究如何建立一個(gè)最佳的擬合模型對(duì)于建立和維持高精度地球參考框架是個(gè)關(guān)鍵問題.