于亞杰，唐江森

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

1 引 言

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

2 小波譜分析

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

(1)

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

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

(2)

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

3 小波分析

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

小波分解表達(dá)式為：

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

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

(3)

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

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ù)預(yù)處理

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

圖1 所選站點(diǎn)分布圖

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

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

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

圖2 保留階躍項(xiàng)的坐標(biāo)時(shí)間序列及其擬合直線

圖3 去除階躍項(xiàng)的坐標(biāo)時(shí)間序列及其擬合直線

采用最小二乘線性擬合對SDCORS基準(zhǔn)站坐標(biāo)時(shí)間序列的線性趨勢項(xiàng)進(jìn)行擬合，以SDCORS中GAOM站U方向坐標(biāo)時(shí)間序列為例，其坐標(biāo)時(shí)間序列及其擬合直線如圖4所示。

圖4 GAOM站U方向原始坐標(biāo)時(shí)間序列及其擬合直線

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

5 SDCORS部分站點(diǎn)頻譜分析

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

圖5GAOM站三個(gè)坐標(biāo)分量的小波譜圖及各周期信號(hào)圖

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

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

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

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

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

6 結(jié) 語

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