CORS基準站坐標時間序列周期性分析

于亞杰，唐江森

(河北省第二測繪院，河北 石家莊 050031)

1 引 言

山東省衛(wèi)星定位連續(xù)運行綜合應用服務系統(tǒng)(Shandong Continuous Operational Reference System，簡稱SDCORS)是“數(shù)字山東”的基礎設施之一，也是山東省一項重要的空間數(shù)據(jù)基礎建設工程。SDCORS已經(jīng)連續(xù)運行多年，積累的大量坐標時間序列可以反映出測站的線性與非線性的運動特征，為山東地區(qū)地殼形變監(jiān)測及一些地球物理現(xiàn)象的解釋提供了寶貴的基礎數(shù)據(jù)。國內(nèi)外大量學者已對CORS系統(tǒng)坐標時間序列進行了深入研究，其方法主要集中在通過公式擬合來確定坐標時間序列模型[1，2]、通過濾波等方法研究其共模誤差[3～7]、研究坐標時間序列噪聲模型[8～12]等方面，對坐標時間序列周期信號特征和產(chǎn)生原因的分析仍不充分。本文擬對SDCORS部分測站坐標時間序列進行頻譜分析，通過研究坐標時間序列的周期性變化特征進一步探索周期信號的成因。

2 小波譜分析

由于基準站受到多種地球物理現(xiàn)象的共同影響，在其坐標時間序列中會表現(xiàn)出一定的周期性。小波譜可以探測信號的所有頻率成分，反映這些頻率成分在整個觀測區(qū)間內(nèi)的幅值變化情況[13]。信號x(t)關于小波函數(shù)ψ(t)的小波變換為：

(1)

信號x(t)的小波變換幅度的平方即為小波譜：

PWx(s，u)=|Wx(s，u)|2

(2)

小波譜度量了信號x(t)在小波時-頻窗中的能量，反映了其局部時頻能量密度大小[14]。

3 小波分析

小波分解可以將基準站坐標時間序列按頻帶分解，將利用小波譜分析得到的各周期信號分離在一個或者幾個頻帶之內(nèi)，再利用小波重構(gòu)將各周期項提取出來。小波分析相當于一個時頻信號濾波器，它能夠?qū)⑿盘柗纸庠谧涌臻g中，并在子空間里分析信號[15～17]。

小波分解表達式為：

AK+1=HAK，K=1，2，…，N

DK+1=GDK，K=1，2，…，N

(3)

小波重構(gòu)表達式為：

AK+1=H*AK+1+G*DK+1，K=N，…，1

(4)

式(4)和式(5)中，H和G分別是小波低通濾波器和高通濾波器；H*和G*分別是H和G的對偶算子；AK，DK分別表示第K層的低頻部分和高頻部分。

4 數(shù)據(jù)預處理

為了研究SDCORS基準站坐標時間序列的周期性變化規(guī)律，選取SDCORS中GAOM(高密)、LAIY(萊陽)、LAIZ(萊州)、LOKO(龍口)、RUSH(乳山)、WEND(文登)、ZHUC(諸城)，DEZH(德州)、DOYI(東營)、HEZE(菏澤)、JYRS(濟陽)、LICH(聊城)、QUFU(曲阜)、SDTA(泰安)、TEZH(滕州)、YAXI(陽信)、ZHCU(周村)，采用2016年1月1日～2017年12月31日單日解坐標時間序列進行分析。所選測站基本覆蓋整個山東區(qū)域，且分布較為均勻，能夠反映SDCORS的整體情況，其站點分布如圖1所示。

圖1 所選站點分布圖

GPS基準站單日解坐標時間序列中存在部分數(shù)據(jù)缺失，必須進行插值以獲得完整的時間序列，才能提取出正確的周期項。根據(jù)SDCORS坐標時間序列的數(shù)據(jù)缺失情況，本次實驗采用三次多項式插值對缺失數(shù)據(jù)進行插補。

由于如更換天線或地震等原因，坐標時間序列可能會從某個時間開始產(chǎn)生階躍式點位突變，去除這種階躍式點突變的方法是從突變時刻開始，對突變時間以后的坐標時間序列加上突變項改正[2]。由于GAOM站E方向原始坐標時間序列中存在明顯的點位突變，因此階躍項的去除以GAOM站E方向坐標時間序列為例。圖1和圖2分別是SDCORS中GAOM站E方向保留階躍項和加階躍項改正之后的坐標時間序列及其擬合直線。

圖2和圖3中，折線表示保留階躍項的GAOM站E方向原始坐標時間序列和加階躍改正后的坐標時間序列，直線表示其擬合直線。從圖1可以看出GAOM站E方向坐標時間序列在第434天發(fā)生了明顯的坐標突變，坐標下降了約 7 mm，圖3是加階躍改正后的坐標時間序列，經(jīng)過改正后的坐標時間序列呈現(xiàn)出較好的一致性。

圖2 保留階躍項的坐標時間序列及其擬合直線

圖3 去除階躍項的坐標時間序列及其擬合直線

采用最小二乘線性擬合對SDCORS基準站坐標時間序列的線性趨勢項進行擬合，以SDCORS中GAOM站U方向坐標時間序列為例，其坐標時間序列及其擬合直線如圖4所示。

圖4 GAOM站U方向原始坐標時間序列及其擬合直線

圖4中折線代表GAOM站U方向插值后坐標時間序列，直線是其擬合直線，即線性趨勢項。從圖4中可以看出坐標時間序列存在明顯的周期性，周期約為 1年。為了確定坐標時間序列中所包含的周期項，并在時域研究周期項的振幅等特征，需要對插值和去除線性趨勢項后的坐標時間序列進行小波譜分析和小波分析。

5 SDCORS部分站點頻譜分析

為了探測SDCORS站點坐標時間序列中的周期項對站點坐標時間序列進行小波譜分析，在時間域研究小波譜探測得到的周期項，對時間序列進行小波分析，將各周期信號轉(zhuǎn)化到時域。為了提高小波譜圖的分辨率，得到時間序列中的所有頻率成分，采用0天～600天數(shù)據(jù)識別長周期項，0天～150天用于識別半年周期以下的短周期項。以SDCORS中GAOM站為例，其N、E、U三個坐標分量的小波譜圖和小波分析如圖5所示：

圖5GAOM站三個坐標分量的小波譜圖及各周期信號圖

圖5中(a)(b)(d)(e)(g)(h)分別是N、E、U三個坐標分量的小波譜圖，(e)(g)(h)是三個坐標分量提高分辨率后的小波譜圖，(i)(c)(h)分別是利用小波分析提取的三個坐標分量的各個周期信號圖。從圖5中可以明顯看出GAOM站在N、E、U三個方向上都存在明顯的年周期、半年周期、季節(jié)周期和月周期信號。各個周期信號的強度和出現(xiàn)的時間不盡相同，其中年周期信號和半年周期信號出現(xiàn)在整個時間段內(nèi)，季節(jié)周期在2016年的11月～12月以及2017年的1月～5月較弱，N方向和U方向月周期在每年的6月下旬～9月底較強，E方向和U方向第一年的1月下旬～4月上旬也出現(xiàn)月周期信號。

對上述17個基準站的小波譜分析結(jié)果表明，其呈現(xiàn)出一定的規(guī)律，總體來說垂直方向上的各周期信號的強度要比水平方向強，在同一周期下的振幅大于水平方向，這和GPS在垂直方向上的測量精度比水平方向上精度低的結(jié)論相符。表明垂直方向受到多種復雜地球物理現(xiàn)象的影響比水平方向嚴重，這也是GPS在垂直方向上的應用受到限制的原因之一。

E方向上的年周期比N方向和U方向上的年周期要長，具體原因尚不明確，可能是多種復雜現(xiàn)象疊加之后在東西方向形成了較長的周期。推測海潮的影響比較大，山東半島東部沿海，海洋潮汐自東向西拍打海岸，使海底負荷發(fā)生變化，海底負荷產(chǎn)生的變化自東向西從沿海傳遞到內(nèi)陸測站，形成了E方向較長的年周期。

在N方向和U方向上，季節(jié)周期在2016年的11月～12月以及2017年的1月～5月較弱，而E方向相對較強。這段時間山東地區(qū)溫度較低，上方的電離層活動緩慢，電子含量較低，對流層濕度相對穩(wěn)定并且較小，GPS信號在經(jīng)過大氣層時產(chǎn)生的大氣延遲在這段時間內(nèi)比較穩(wěn)定，不會出現(xiàn)明顯變化，此外這段時間內(nèi)山東地區(qū)的降水量偏小，農(nóng)作物灌溉用水量較少，地下水位在這段時間呈現(xiàn)均勻的緩慢下降的趨勢，對GPS定位的影響緩慢，因此在坐標時間序列中沒有出現(xiàn)顯著的周期性變化。而在E方向上，海洋半年潮汐造成海底負荷的變化，海底負荷的變化自東向西傳遞到內(nèi)陸，在E方向上產(chǎn)生季節(jié)性周期信號?？偟膩碚f，造成GPS時間序列呈現(xiàn)出季節(jié)性周期的因素主要是海潮、大氣延遲以及地下水的變化，其中海潮對E方向的影響較為嚴重。

N方向和U方向的月周期在每年的6月下旬～9月較明顯，這段時間山東地區(qū)氣溫偏高，降水豐富，大氣濕度變化較快，電離層電子濃度增大，山東地區(qū)地下水位維持在較高的水平，這些因素相互作用，引起了GPS坐標時間序列在這段時間呈現(xiàn)月周期信號。而E方向的月周期則在2016年的6月下旬～9月底，2017年的2月～5月較強，這段時間比N方向和U方向提前100天左右，此外2016年的1月下旬～4月上旬E方向上的月周期比其他兩個方向明顯，具體成因需要進一步研究。

6 結(jié) 語

小波譜能夠探測到CORS基準站時間序列信號中所隱含的周期項。而小波分析能將小波譜探測到的周期項通過小波分解和小波重構(gòu)從原始坐標時間序列中分離并提取出來，進而得到了SDCORS基準站在N、E、U三方向的各周期信號。通過對各周期信號的分析，得到海潮負荷、大氣延遲和地下水位變化等原因是這些周期項成因的結(jié)論。