香港CORS站坐標時間序列分析研究

陳晨,魏冠軍,高志鈺,寇瑞雄

(1.蘭州交通大學 測繪與地理信息學院,甘肅 蘭州 730070;2.地理國情監(jiān)測技術應用國家地方聯(lián)合工程研究中心，甘肅 蘭州 730070；3.甘肅省地理國情監(jiān)測工程實驗室,甘肅 蘭州 730070)

0 引 言

全球定位系統(tǒng)(GPS)基準站坐標時間序列既可以反映出測站的線性變化,又可以反映出測站存在的非線性變化[1]．建立合適的區(qū)域連續(xù)運行參考站(CORS)噪聲模型對研究CORS站的穩(wěn)定性及區(qū)域地表形變規(guī)律具有重要的意義．Zhang等[2]首次將噪聲分析引入到GPS坐標時間序列分析中來,研究結果表明坐標時間序列中具有“白噪聲(WN)+閃爍噪聲(FN)”的特征．Mao等[3]采用譜分析和極大似然估計方法分析了23個全球IGS站的噪聲特性,得出最佳噪聲模型為WN加FN．Langbein[4]分析了加利福尼亞南部及內華達南部地區(qū)236個連續(xù)GPS運行站的噪聲模型,認為最佳噪聲模型以FN或者隨機游走噪聲RW為主．蔣志浩等[5]對我國國家CORS網(wǎng)進行了數(shù)據(jù)處理,認為國家CORS網(wǎng)基準站坐標時間序列噪聲表現(xiàn)為WN、FN及RW．李昭等[6]對中國區(qū)域11個IGS基準站坐標時間序列進行地表質量負荷改正后,得到各坐標分量上的最優(yōu)噪聲模型主要以WN加FN為主,同時也存在冪律噪聲(PL)、帶通噪聲(BP)及一階高斯馬爾科夫噪聲(FOGM)．李斐等[7]對南極半島地區(qū)8個IGS基準站坐標時間序列進行空間濾波后認為各站點坐標時間序列中不僅存在WN,還存在有色噪聲(CN),在有些站點的E方向可能還存在WN．袁林果等[8]對香港12個GPS連續(xù)運行參考站坐標時間序列分析后,認為最優(yōu)噪聲模型為可變自噪聲(VW)加FN．

由此可見,不同地區(qū)噪聲特性存在著差異,針對此項問題,本文以香港衛(wèi)星參考站網(wǎng)(SatRef)為研究對象,來分析測站坐標時間序列的周期特性與噪聲特性．在袁林果研究的基礎上,選擇更多的噪聲模型來分析確定最優(yōu)噪聲模型．考慮到香港地區(qū)區(qū)域較小(1106.34平方千米),最長基線HKNP-HKWS僅為49.89 km,參考站之間具有很強的空間相關性．因此,本文決定選取石碑山(HKOH)、昂坪(HKNP)、錦田(HKKT)、藍地(HKLT)和南丫島(HKLM)這五個數(shù)據(jù)質量較好的站點來進行分析．利用功率譜分析的方法來獲取各個參考站點的周期項．對單天解坐標分量時間序列建立擬合參數(shù)模型,采用6種噪聲類型,將其組合成14種噪聲模型,運用極大似然法估計參考站坐標時間序列在各個噪聲模型下的噪聲量級、周期性振幅、測站速度及不確定度等信息．選用Langbein提出的保守估計準則判斷不同噪聲模型的優(yōu)劣,確定測站各分量上的最優(yōu)噪聲模型．比較分析在僅考慮WN和在顧及CN的影響下,參考站的線性速度以及周期振幅和相位等參數(shù)估計值的差異性,建立合適的區(qū)域CORS網(wǎng)速度場模型．

1 區(qū)域CORS網(wǎng)數(shù)據(jù)預處理

香港衛(wèi)星參考站網(wǎng)(SatRef)截止到2018年10月1日共有18個參考站,其中T430站和大埔滘(HKTP)站的觀測時間較短,所以本次數(shù)據(jù)處理不采用這兩個點．選取數(shù)據(jù)長度為2014-05-01到2017-08-31共40個月的參考站網(wǎng)觀測數(shù)據(jù),同時引入ITRF2014參考框架下BIFS、CUSV、LHAZ、PIMO、SHAO、TWTF、TCMS和PBRI共8個IGS站的同步觀測數(shù)據(jù),起到一個約束的作用．基線解算時海潮改正模型采用FES2004模型,對流層模型采用Saastamoinen改正模型,地潮改正模型采用IERS03．檢驗基線解算的成果,其最短、最長基線重復率分別為0.4 mm和1.9 mm,最短基線的相對重復率為2.24×10-8量級,最長基線相對重復率為5.87×10-10量級．基線解得標準均方根誤差(NRMS)值大多處于0.15～0.18之間,最大不超過0.2,表明基線解算精度很高,符合網(wǎng)平差要求．GAMIT基線解算完成后,利用GLOBK軟件進行平差,得到各參考站N、E、U三個方向上的平均NRMS值分別為1.9 mm、1.8 mm和1.3 mm,均優(yōu)于3 mm,達到了較高的解算精度．

經(jīng)過GAMIT/GLOBK軟件處理后得到各站的原始坐標時間序列,由于篇幅限制,此處只列出 HKOH站三個坐標分量上的原始坐標時間序列圖,并對其進行了簡單的線性擬合,如圖1～3所示．

圖1 HKOH站N方向原始坐標時間序列圖

圖2 HKOH站E方向原始坐標時間序列圖

圖3 HKOH站U方向原始坐標時間序列圖

從圖1～3可以得出以下結論:

1)HKOH站在N、E兩個水平方向主要表現(xiàn)為線性變化,且有整體向東南方向運動的趨勢,N方向的運行速度約為-12 mm/a,E方向上的運行速度約為32 mm/a．

2)U方向上主要呈現(xiàn)為周期性變化,周期性振幅具有夏季大,冬季小的季節(jié)性特征,振幅約為10 mm左右．

3)該站在U方向上要比在N、E方向上的離散程度大,這是因為高程方向上的精度要低于水平方向精度,即在水平方向上該點較穩(wěn)定,在高程方向上波動較大．

2 參考站坐標時間序列功率譜分析

功率譜分析法是分析時間序列的常用方法之一．功率譜密度表示了各個頻率上的能量強度.觀察離散數(shù)據(jù)系列的周期性,如果某個信號具有周期特性,則其周期運動所對應的功率在全部功率中占有比重較大,在功率譜圖上表現(xiàn)為峰值[9-10]．

地球物理現(xiàn)象的功率譜可以用下式表示:

P(f)=P0f-α,

(1)

其對數(shù)形式為

InP(f)=InP0-αInf,

(2)

式中:P(f)為功率譜密度;f為頻率;P0為系數(shù);α為譜指數(shù),P0和α都為待定參數(shù)．

許多研究表明GPS坐標時間序列中存在著周期信號,表現(xiàn)為年周期與半年周期特性[11-13]．本文選擇了SatRef中HKOH、HKKT、HKNP、HKLM和HKLT五個具有代表性的站點來探測香港地區(qū)的衛(wèi)星參考站的坐標時間序列是否也具有周期特性．主要方法是利用快速傅里葉變換(FFT)來獲取各個站點的功率譜密度P(f),通過分析它們的功率譜圖,能夠直觀地揭示各參考站點N、E、U三個坐標分量上的周期特性．需要注意的是,在進行FFT的時候,原始坐標時間序列需要滿足均勻采樣和零均值的特性．在不滿足上述條件的情況下,也可以利用周期圖法進行譜分析,弊端是計算時間較長[14]．

去除掉原始坐標時間序列線性趨勢,經(jīng)過FFT得到各參考站的功率譜圖,如圖4～6所示．

圖4 各基準站N方向功率譜圖

圖5 各基準站E方向功率譜圖

圖6 各基準站U方向功率譜圖

圖4～6中橫坐標表示年周期數(shù),縱坐標表示單位周期內的能量．從圖中可以直觀地看出各個參考站點在N、E、U三個坐標分量方向上均表現(xiàn)出一定的周期特性,且以年周期項最為明顯．各個分量上還存在較為明顯的半年周期項．HKNP和HKLT兩個站的U方向上還存在1.5年周期項．經(jīng)過統(tǒng)計發(fā)現(xiàn),所有參考站點在U分量方向上的功率譜能量高于N、E分量方向2～5倍,這說明高程方向上的周期性變化要比水平方向變化更明顯．經(jīng)過統(tǒng)計發(fā)現(xiàn),除HKOH站以外,其余各站N分量方向上年周期所對應的功率譜能量普遍大于E方向,說明在N方向的年周期變化要比E方向明顯．

3 參考站坐標時間序列噪聲模型建立與分析

3.1 極大似然估計法

GPS坐標時間序列噪聲分析可以用極大似然估計(MLE)的方法來完成．MLE可以同時估計噪聲類型、周期性振幅、測站速率及不確定度.

MLE可以估計參考站坐標時間序列殘差中包含的WN與相關噪聲的振幅．主要原理是對單天解坐標分量時間序列建立如下的參數(shù)模型[6,15]:

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

esin(4πti)+fcos(4πti)+

(3)

NIn(2π)],

(4)

式中:N為時間序列長度;C為協(xié)方差陣;det表示矩陣C的行列式;lik為似然值．

3.2 噪聲類型

大量的研究成果表明,GPS坐標時間序列中同時存在WN和CN．本文主要分析WN與CN的組合噪聲模型．在實驗中選取的6種噪聲分別為:WN、VW、FN、RW、PL和一階高斯馬爾科夫噪聲(GM)．將6種噪聲分成兩組組合成14種噪聲模型,第一組噪聲模型分別為:WN、WN+FN、WN+RW、WN+FN+RW、WN+GM、WN+RW+GM和WN+PL．第二組噪聲模型分別為:VW、VW+FN、VW+RW、VW+FN+RW、VW+GM、VW+RW+GM和VW+PL．

3.3 最優(yōu)噪聲模型評價準則

根據(jù)極大似然估計法的原理可知,不同的噪聲模型將得到不同的MLE值,一般認為,極大似然值越大,結果越可靠[3-4]．但是噪聲模型中包含的未知參數(shù)越多,其MLE值也越大．因此,為了確保結果的可靠性,不能只以MLE值的大小來判定最優(yōu)噪聲模型．對于不同的噪聲模型,蒙特·卡洛實驗表明:95%的顯著水平下,當兩種噪聲模型的MLE之差大于3.0時,兩種模型具有可區(qū)分性[15]．

3.4 香港衛(wèi)星參考站坐標時間序列噪聲分析

3.4.1 WN模型的選取

針對SatRef中選取的五個參考站點,運用極大似然估計法得到其在14種噪聲模型下的MLE值[16]．對五個站點N、E、U三個分量上在WN模型和VW模型下的平均MLE值進行比較,結果如表1所示．

表1 WN模型與VW模型MLE值對比

從表1中可以看出,簡單WN模型的MLE值大于VW模型的MLE值,根據(jù)最優(yōu)噪聲模型評價準則可知,就WN模型而言,WN模型比VW模型更適合于SatRef．此外,利用MLE計算得到的五個參考站分量VW+FN模型的MLE值明顯大于其他組合模型,這與袁林果等對SatRef的噪聲模型研究結果一致[8]．但此后的研究結果認為VW模型僅能反映測站分量的質量好壞,并不能作為參考站最優(yōu)噪聲模型[6]．根據(jù)上述兩點原因,本文決定采用WN與其組合模型來分析SatRef的最優(yōu)噪聲模型．

3.4.2 最優(yōu)噪聲組合模型的確定

將WN+FN、WN+RW、WN+FN+RW、WN+GM、WN+RW+GM和WN+PL六種噪聲模型的MLE值與WN的MLE值作差,得到MLE值的差值圖如圖7所示.

從圖7中可以看出:6種噪聲組合模型與WN模型MLE值的差值都大于0,表明SatRef的坐標時間序列中不僅含有WN,還存在CN．WN+FN、WN+GM、WN+RW+GM和WN+PL四種噪聲模型與WN模型的差值大于其他兩種噪聲模型,說明這4種噪聲模型較其它模型更優(yōu)越,且這四種噪聲模型的MLE值相差不大．

通過3.3節(jié)的介紹可知,不能只以MLE值的大小來判定參考站點最優(yōu)噪聲模型,為了確保結果的可靠性,本文選用Langbein提出的保守估計準則判斷不同噪聲模型的優(yōu)劣．首先,選取WN+FN和WN+RW兩種組合模型MLE值較大者作為零假設,然后將WN+FN+RW、WN+GM和WN+PL噪聲模型的MLE值與零假設作比較,若MLE差值大于2.6,則拒絕零假設,認為該模型更優(yōu),否則接受零假設,認為所選的復雜模型無效．若WN+FN+RW、WN+GM和WN+PL均優(yōu)于零假設,則選擇MLE較大的作為最優(yōu)模型,根據(jù)此方法來尋找各個參考站點N、E、U分量上的最優(yōu)噪聲模型．本文將接受WN+RW+GM組合模型的閾值設為5.2,接受其余噪聲組合模型的閾值設為2.6．表2示出參考站各分量的最優(yōu)噪聲模型．

圖7 五個參考站N/E/U方向MLE值差值

參考站點N分量E分量U分量 HKOHWN+FNWN+FNWN+PL HKKTWN+FNWN+FNWN+GM HKNPWN+GMWN+GMWN+PL HKLMWN+PLWN+FNWN+FN HKLTWN+GMWN+FNWN+FN

從表2中可以看出SatRef坐標時間序列的噪聲特性呈現(xiàn)多樣化的特征,15個坐標分量上WN+FN組合模型所占比例最大,為53.3%．WN+GM和WN+PL組合模型所占比例分別為26.7%和20%．由此可見,盡管香港地區(qū)區(qū)域較小,參考站點之間相距較近,但不同的站點的坐標時間序列仍表現(xiàn)出不同的噪聲特性．目前對產(chǎn)生這種差異性的原因尚不明確,推測可能是與框架的定義以及地球物理效應的空間相關性有關,有的測站可能存在較大的系統(tǒng)誤差[5-6]．

4 顧及有色噪聲對參數(shù)估計的影響

利用CATS軟件進行極大似然估計時,會得出公式(3)中的各個參數(shù)．為顧及CN對參數(shù)估計的影響,在ITRF2014框架下,本文對SatRef中選取的五個參考站點采用WN+FN噪聲模型估計其線性速率、年周期振幅和相位三個參數(shù),并與WN模型進行比較[7,17-18]．

表3所示為參數(shù)估計結果,可以得出以下兩個結論:

1)在WN+FN和WN兩種不同的噪聲模型下所估計的參數(shù)具有明顯的差異．三個參數(shù)中,速率最大差異在HKLM站的N方向,為2.64 mm/a;年周期振幅最大差異在HKNP站的N方向,為0.54 mm;相位最大差異在HKLT站的U方向,為2.18 rad．

2)比較在WN+FN和WN兩種噪聲模型下所得到的線性速率、年周期振幅和相位三個參數(shù)估計的最大不確定度,前者約為后者的10倍、5倍和8倍．由此說明,僅考慮WN所獲取的參考站線性速率精度并不能反映實際的SatRef速度場的精度,同時也在很大程度上低估了年周期振幅和相位的參數(shù)估計不確定度．

結合結論1)、2)可知,在進行參考站坐標時間序列分析時,CN對參數(shù)估計的影響是不容忽視的．

5 顧及CN的水平速度場分析

由第4章可以得出在WN+FN噪聲模型下HKOH、HKKT、HKNP、HKLM和HKLT五個站點的線性速率估值及其精度,其中N方向的平均速率約為-11.46 mm/a,E方向的平均速率約為31.45 mm/a,U方向的平均速率為-0.18 mm/a,水平速度場圖如圖8所示．由文獻[19]和[20]可知位于華南塊體的SHAO、WUHNS、XIAM的平均速度分別為為35.50 mm/a、35.82 mm/a和34.54 mm/a[19-20],結合SatRef的水平速度為33.47±1.5 mm/a可知,華南塊體整體穩(wěn)定,運動速度與趨勢基本一致．

圖8 顧及CN的SatRef站速度場

6 結束語

運用GAMIT/GLOBK數(shù)據(jù)處理軟件對SatRef中的16個站點進行基線解算與網(wǎng)平差,得到了高精度的參考站原始坐標時間序列．通過對選取的HKOH、HKKT、HKNP、HKLM和HKLT五個站點進行周期特性與噪聲特性分析后可得到以下結論:

1)各個站點的三個坐標分量上都存在年周期與半年周期的特性,在某些站點的U方向還存在1.5年周期特性,且高程方向上的周期性變化比水平方向更明顯．

2)采用極大似然估計法來分析噪聲特性后,發(fā)現(xiàn)SatRef站點的坐標時間序列噪聲特性呈現(xiàn)多樣化的特征,其中WN+FN組合噪聲模型所占比例最大．

3)在進行測站速度,年周期振幅和相位等參數(shù)估計時,不能僅僅考慮WN,還必須考慮CN的影響．

4)根據(jù)最優(yōu)噪聲模型WN+FN估計得出SatRef站速度場在水平方向上有整體向東南方向運動的趨勢,與周邊IGS站點以及華南塊體的運動方向基本一致．