徐創(chuàng)富,王恒,吉淵明,劉宇

(1.臺(tái)州市空間地理信息,浙江 臺(tái)州 318000;2.麗水市測(cè)繪中心,浙江 麗水 323000;3.浙江省測(cè)繪科學(xué)技術(shù)研究院, 浙江 杭州 310000;4.山東科技大學(xué)測(cè)繪科學(xué)與工程學(xué)院, 山東 青島 266000;5.中國(guó)測(cè)繪科學(xué)研究院, 北京 100830)

0 引 言

隨著連續(xù)運(yùn)行參考站(CORS)站點(diǎn)觀測(cè)精度的不斷提高,國(guó)內(nèi)外大量學(xué)者對(duì)CORS站時(shí)間序列的分析已經(jīng)進(jìn)行了深入的研究,其研究?jī)?nèi)容主要集中于對(duì)時(shí)間序列進(jìn)行擬合確定時(shí)間序列模型[1-5],通過(guò)濾波等方法研究其共模誤差分析其噪聲模型,通過(guò)功率譜、小波變換等方法研究其周期特征[6-8]．目前大多數(shù)研究都是對(duì)國(guó)際GNSS服務(wù)(IGS)站及一些小區(qū)域的CORS站進(jìn)行分析[9-11],對(duì)于浙江區(qū)域CORS站的分析還非常少．

浙江省衛(wèi)星導(dǎo)航定位基準(zhǔn)服務(wù)系統(tǒng)(ZJCORS)于2008年向全省用戶提供服務(wù),這些覆蓋全省的CORS站已持續(xù)運(yùn)行10年,一直在全天候連續(xù)監(jiān)測(cè)地面垂直形變．CORS站點(diǎn)能以毫米級(jí)精度全天候、連續(xù)監(jiān)測(cè)地表點(diǎn)位的三維位置變化．所以本文基于臺(tái)州及周邊地區(qū)15個(gè)CORS站2017-2019年的數(shù)據(jù),采用線性擬合、功率譜、小波變換等方法對(duì)時(shí)間序列進(jìn)行分析．

1 數(shù)據(jù)來(lái)源與預(yù)處理

1.1 數(shù)據(jù)來(lái)源

采用GAMIT/GLOBK軟件對(duì)數(shù)據(jù)進(jìn)行處理,首先使用GAMIT軟件進(jìn)行基線解算,解算過(guò)程中主要參數(shù)歷元間隔30 s,固體潮模型采用IERS2010,海潮模型采用FES2004,IGS站先驗(yàn)坐標(biāo)為SOPAC公布的ITRF2014下坐標(biāo)．然后利用GLOBK平差軟件進(jìn)行整網(wǎng)平差,得到ITRF2014框架下的站點(diǎn)坐標(biāo)變化時(shí)間序列及運(yùn)動(dòng)速度值,在進(jìn)行解算時(shí),選擇了11個(gè)IGS 站點(diǎn)參與計(jì)算,站點(diǎn)包括:AIRA、ARTU、BJFS、CHAN、DAEJ、LHAZ、PIMO、POL2、SHAO、URUM、YSSK,在計(jì)算過(guò)程中對(duì)于IGS站點(diǎn)E、N、U三個(gè)方向的約束量分別為5 mm、5 mm、10 mm,對(duì)于CORS站點(diǎn)三個(gè)方向的約束均為10 m,圖1為CORS站點(diǎn)分布圖．

圖1 CORS站分布

經(jīng)GAMIT/GLOBK軟件處理之后獲得15個(gè)CORS站2017年1月1日-2019年8月7日的時(shí)間序列,其大部分測(cè)站時(shí)間序列完整度在95%以上,僅有少部分站時(shí)間長(zhǎng)度在80%以下．

1.2 數(shù)據(jù)預(yù)處理

通過(guò)觀察CORS站點(diǎn)時(shí)間序列可以發(fā)現(xiàn)時(shí)間序列存在少量的間斷點(diǎn)以及突變點(diǎn),在對(duì)時(shí)間序列進(jìn)行分析時(shí)應(yīng)該首先將突變數(shù)據(jù)剔除,剔除標(biāo)準(zhǔn)是將大于三倍標(biāo)準(zhǔn)差的值刪掉,然后通過(guò)插值對(duì)數(shù)據(jù)進(jìn)行修復(fù),以獲得均勻采樣的時(shí)間序列,以此才能正確地提取時(shí)間序列中的周期項(xiàng)并對(duì)其分析．

常用的插值方法主要有線性插值和三次樣條插值,線性插值及三次樣條插值法適用于缺失值較少的時(shí)間序列[12],且缺失值越少插值的效果越好,目前有學(xué)者新提出正交多項(xiàng)式擬合[12]方法,對(duì)于數(shù)據(jù)連續(xù)缺失值較多的問(wèn)題提供了更好的結(jié)局方案,本文根據(jù)15個(gè)CORS站數(shù)據(jù)質(zhì)量情況,分別采用三次樣條插值以及正交多項(xiàng)式擬合的方法對(duì)數(shù)據(jù)進(jìn)行插值處理,

圖2為測(cè)站數(shù)據(jù)預(yù)處理前后的對(duì)比圖．圖2(a)為預(yù)處理前時(shí)間序列,圖2(b)為預(yù)處理后時(shí)間序列,因篇幅所限只展示DONT站U方向處理結(jié)果,可以發(fā)現(xiàn)時(shí)間序列中突變點(diǎn)有明顯減少,數(shù)據(jù)間斷也得到處理．

(a)DONT站U方向預(yù)處理前時(shí)間序列

(b)DONT站U方向預(yù)處理后時(shí)間序列圖2 預(yù)處理前后時(shí)間序列對(duì)比

1.3 線性項(xiàng)提取

使用譜分析進(jìn)行周期探測(cè)時(shí),要求時(shí)間序列的均值為零,即沒(méi)有趨勢(shì)項(xiàng),否則會(huì)影響功率譜分析的結(jié)果,根據(jù)式(1)對(duì)時(shí)間序列進(jìn)行趨勢(shì)項(xiàng)擬合．

y=a+bt+v，

(1)

式中：y為原始時(shí)間序列；a、b為趨勢(shì)項(xiàng)參數(shù)；b為線性速率；t為時(shí)間；v為時(shí)間序列殘差．圖3為CORS站時(shí)間序列插值與去線性項(xiàng)之后的結(jié)果圖,紅線為插值之后的時(shí)間序列,藍(lán)線為線性項(xiàng),綠線為殘差時(shí)間序列,由于篇幅所限,只展示DONT站的處理結(jié)果,由圖3可以看出臺(tái)州CORS站的線性運(yùn)動(dòng)主要存在于水平方向,垂直方向線性運(yùn)動(dòng)相較于水平方向較小．去線性項(xiàng)之后DONT站水平東方向時(shí)間序列殘差范圍在-10～10 mm,水平北方向時(shí)間序列殘差范圍在-5～5 mm,垂直方向時(shí)間序列殘差范圍在-20～20 mm．

(a)DONT站E方向時(shí)間序列

(c)DONT站U方向時(shí)間序列圖3 去趨勢(shì)項(xiàng)時(shí)間序列圖

2 三維速度場(chǎng)

利用15座CORS站線性項(xiàng),以時(shí)序長(zhǎng)度為權(quán),變差函數(shù)選擇高斯模型,按加權(quán)克里金進(jìn)行格網(wǎng)化,再對(duì)格網(wǎng)采用滑動(dòng)平均濾波,生成的大地高年變率格網(wǎng)模型．如圖4所示，為了更好地顯示水平速度場(chǎng)特征，繪圖過(guò)程中，移去了臺(tái)州市整體水平速度，由圖可以看出臺(tái)州地區(qū)地面沉降呈明顯的空間分布不均勻特點(diǎn),從總體上來(lái)看,臺(tái)州市西北部天臺(tái)臨海地區(qū)沉降明顯,臺(tái)州東南部地區(qū)隆升明顯．

圖4 大地高年變率格網(wǎng)模型

表1為臺(tái)州地區(qū)附近部分CORS站運(yùn)動(dòng)速度統(tǒng)計(jì)表,根據(jù)統(tǒng)計(jì)2017-2019年三年時(shí)間中臺(tái)州及其附近地區(qū)各基準(zhǔn)站平均運(yùn)動(dòng)速度值為北方向-12.9 mm/a,東方向31 mm/a,總體運(yùn)動(dòng)趨勢(shì)為東南方向,這與前人所進(jìn)行的研究也是相一致的．CORS站高程方向呈不同的運(yùn)動(dòng)趨勢(shì),LINH站沉降最為明顯其沉降速率超過(guò)1 mm/a,臺(tái)州站抬升最為明顯,隆升速率超過(guò)3 mm/a．

表1 臺(tái)州CORS站運(yùn)動(dòng)速度統(tǒng)計(jì)

3 時(shí)序周期特征分析

3.1 功率譜分析

功率譜是功率譜密度函數(shù)的簡(jiǎn)稱,它定義為單位頻帶內(nèi)的信號(hào)功率,功率譜表示了信號(hào)功率隨著頻率的變化關(guān)系[13-14]．本文使用周期圖法對(duì)時(shí)間序列進(jìn)行功率譜分析,對(duì)于隨機(jī)信號(hào)s(ω)其相應(yīng)的功率譜表達(dá)式如式(2)所示:

(2)

由公式(2)可知,功率譜分析只是在頻域?qū)π盘?hào)進(jìn)行分析,并不包含任何時(shí)頻的信息,因此通過(guò)功率譜分析僅能獲得CORS站時(shí)間序列的周期信息,而不能得到其各個(gè)周期具體的持續(xù)時(shí)間以及振幅等信息．圖5為DNOT站E、N、U三個(gè)方向的功率譜密度圖,橫坐標(biāo)為周期,縱坐標(biāo)為功率譜能量密度．

根據(jù)計(jì)算,臺(tái)州地區(qū)CORS站E、N、U三個(gè)方向均具有一定的周期性,但是U方向上的功率譜能量明顯高于E、N方向,約為水平方向的十倍,這說(shuō)明高程方向上的周期性變化明顯大于水平方向．根據(jù)圖中兩個(gè)CORS站的功率譜密度可以看出U方向主要的年周期大概在320天左右,半年周期在200天左右,季節(jié)周期在90～120天．

圖5 功率譜密度

3.2 小波變換

對(duì)CORS站進(jìn)行功率譜分析只能得到其周期信號(hào)的頻率范圍,如果想要得到其精確的周期、振幅及持續(xù)時(shí)間,我們應(yīng)該使用小波分析對(duì)基準(zhǔn)站的數(shù)據(jù)進(jìn)行處理,小波分析克服了功率譜分析的缺點(diǎn),具有多分辨率分析的特點(diǎn),在時(shí)域和頻域都有表征信號(hào)局部信息的能力,其可以根據(jù)信號(hào)的頻率改變時(shí)間分辨率,從而探測(cè)正常信號(hào)中的瞬態(tài),并展示其頻率成分．小波變換的含義是把某一被稱為基本小波的函數(shù)作位移b后,再在不同尺度a下,與待分析信號(hào)f(x)作內(nèi)積,其原理如式(3)所示:

f(x)∈L2(R).

(3)

上述積分變換為f(x)以ψ(x)為基的積分連續(xù)小波變換(CWT),a為尺度因子,表示與頻率相關(guān)的伸縮,b為時(shí)間平移因子．

在實(shí)際應(yīng)用中,連續(xù)小波變換需要離散化,相應(yīng)的離散化小波變換系數(shù)為.

(4)

式中,ψj,k為離散小波函數(shù)．

根據(jù)之前功率譜密度得到的結(jié)果,本文將主要的周期信號(hào)分為年周期、半年周期以及季節(jié)周期,對(duì)于其余的短周期以及不規(guī)則周期信號(hào)不進(jìn)行處理,圖6為通過(guò)小波變換所獲得的DONT站三個(gè)方向各周期振幅及時(shí)間信息結(jié)果圖．

通過(guò)小波分析結(jié)果圖6可以看出,臺(tái)州地區(qū)及其附近CORS站U方向年周期振幅在5～10 mm,波峰出現(xiàn)在每年8月份左右,波谷出現(xiàn)在每年3月份左右,半年周期振幅在0～5 mm,波峰出現(xiàn)在每年6月份和12月份,波谷出現(xiàn)在每年3月份和9月份,季節(jié)周期振幅在5～10 mm,波峰出現(xiàn)在每年的3、6、9、12月份,波谷出現(xiàn)在每年的2、5、8、11月份,CORS站水平方向的年周期振幅在0～2 mm,半年周期以及季節(jié)周期振幅均在0～1 mm,由此可見(jiàn),CORS站水平方向上的周期性明顯弱于垂直方向,垂直方向相較于水平方向受到地球物理因素的影響更加明顯．

根據(jù)研究表明,引起CORS站周期性變化的原因非常復(fù)雜,通常來(lái)說(shuō)主要分為以下三類,第一類是由日月引力攝動(dòng)引起的,包括極移、固體潮、海潮以及大氣負(fù)荷等,第二類是由流體引起的,包括大氣壓的季節(jié)性變化、非潮汐負(fù)荷的波動(dòng)．第三類是由各種誤差引起的,包括衛(wèi)星軌道模型、電離層模型,相位中心的變化等[15-16]．

臺(tái)州地處浙江東南部沿海地區(qū),屬于亞熱帶季風(fēng)氣候,四季分明,夏季氣溫炎熱,降雨量豐富,冬季氣溫低,降雨量少,因此,臺(tái)州CORS站受到的非構(gòu)造負(fù)荷的影響十分明顯,CORS站周期性變化更加突出．

圖6 小波分析結(jié)果圖

4 結(jié) 論

本文利用線性擬合法、功率譜分析和小波變換技術(shù),對(duì)臺(tái)州及其附近地區(qū)15個(gè)CORS站2017年1月-2019年8月的連續(xù)觀測(cè)數(shù)據(jù)進(jìn)行時(shí)間序列分析,結(jié)果表明,該地區(qū)CORS站的水平運(yùn)動(dòng)主要以東南方向的線性運(yùn)動(dòng)為主,平均運(yùn)動(dòng)速度值為東方向31 mm/a、南方向12.9 mm/a,總體運(yùn)動(dòng)速度為33.58 mm/a．在高程方向,臺(tái)州地區(qū)沉降具有明顯的區(qū)域性特征,具體表現(xiàn)為臺(tái)州西北部CORS站沉降明顯,最大沉降速率超過(guò)1 mm/a,臺(tái)州東南部CORS站抬升明顯,最大隆升速率超過(guò)3 mm/a,位于臺(tái)州市區(qū)內(nèi)．