基于臺州及周邊地區(qū)CORS站網(wǎng)地表形變時空特征分析

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

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

0 引 言

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

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

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

1.1 數(shù)據(jù)來源

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

圖1 CORS站分布

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

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

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

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

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

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

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

1.3 線性項提取

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

y=a+bt+v，

(1)

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

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

(c)DONT站U方向時間序列圖3 去趨勢項時間序列圖

2 三維速度場

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

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

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

表1 臺州CORS站運動速度統(tǒng)計

3 時序周期特征分析

3.1 功率譜分析

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

(2)

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

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

圖5 功率譜密度

3.2 小波變換

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

f(x)∈L2(R).

(3)

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

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

(4)

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

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

通過小波分析結(jié)果圖6可以看出,臺州地區(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,由此可見,CORS站水平方向上的周期性明顯弱于垂直方向,垂直方向相較于水平方向受到地球物理因素的影響更加明顯．

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

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

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

4 結(jié) 論

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