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

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

0 引 言

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

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

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

香港衛(wèi)星參考站網(wǎng)(SatRef)截止到2018年10月1日共有18個參考站,其中T430站和大埔滘(HKTP)站的觀測時間較短,所以本次數(shù)據(jù)處理不采用這兩個點(diǎn)．選取數(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．檢驗基線解算的成果,其最短、最長基線重復(fù)率分別為0.4 mm和1.9 mm,最短基線的相對重復(fù)率為2.24×10-8量級,最長基線相對重復(fù)率為5.87×10-10量級．基線解得標(biāo)準(zhǔn)均方根誤差(NRMS)值大多處于0.15～0.18之間,最大不超過0.2,表明基線解算精度很高,符合網(wǎng)平差要求．GAMIT基線解算完成后,利用GLOBK軟件進(jìn)行平差,得到各參考站N、E、U三個方向上的平均NRMS值分別為1.9 mm、1.8 mm和1.3 mm,均優(yōu)于3 mm,達(dá)到了較高的解算精度．

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

圖1 HKOH站N方向原始坐標(biāo)時間序列圖

圖2 HKOH站E方向原始坐標(biāo)時間序列圖

圖3 HKOH站U方向原始坐標(biāo)時間序列圖

從圖1～3可以得出以下結(jié)論:

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

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

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

2 參考站坐標(biāo)時間序列功率譜分析

功率譜分析法是分析時間序列的常用方法之一．功率譜密度表示了各個頻率上的能量強(qiáng)度.觀察離散數(shù)據(jù)系列的周期性,如果某個信號具有周期特性,則其周期運(yùn)動所對應(yīng)的功率在全部功率中占有比重較大,在功率譜圖上表現(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坐標(biāo)時間序列中存在著周期信號,表現(xiàn)為年周期與半年周期特性[11-13]．本文選擇了SatRef中HKOH、HKKT、HKNP、HKLM和HKLT五個具有代表性的站點(diǎn)來探測香港地區(qū)的衛(wèi)星參考站的坐標(biāo)時間序列是否也具有周期特性．主要方法是利用快速傅里葉變換(FFT)來獲取各個站點(diǎn)的功率譜密度P(f),通過分析它們的功率譜圖,能夠直觀地揭示各參考站點(diǎn)N、E、U三個坐標(biāo)分量上的周期特性．需要注意的是,在進(jìn)行FFT的時候,原始坐標(biāo)時間序列需要滿足均勻采樣和零均值的特性．在不滿足上述條件的情況下,也可以利用周期圖法進(jìn)行譜分析,弊端是計算時間較長[14]．

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

圖4 各基準(zhǔn)站N方向功率譜圖

圖5 各基準(zhǔn)站E方向功率譜圖

圖6 各基準(zhǔn)站U方向功率譜圖

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

3 參考站坐標(biāo)時間序列噪聲模型建立與分析

3.1 極大似然估計法

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

MLE可以估計參考站坐標(biāo)時間序列殘差中包含的WN與相關(guān)噪聲的振幅．主要原理是對單天解坐標(biāo)分量時間序列建立如下的參數(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坐標(biāo)時間序列中同時存在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)噪聲模型評價準(zhǔn)則

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

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

3.4.1 WN模型的選取

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

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

從表1中可以看出,簡單WN模型的MLE值大于VW模型的MLE值,根據(jù)最優(yōu)噪聲模型評價準(zhǔn)則可知,就WN模型而言,WN模型比VW模型更適合于SatRef．此外,利用MLE計算得到的五個參考站分量VW+FN模型的MLE值明顯大于其他組合模型,這與袁林果等對SatRef的噪聲模型研究結(jié)果一致[8]．但此后的研究結(jié)果認(rèn)為VW模型僅能反映測站分量的質(zhì)量好壞,并不能作為參考站最優(yōu)噪聲模型[6]．根據(jù)上述兩點(diǎn)原因,本文決定采用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的坐標(biāo)時間序列中不僅含有WN,還存在CN．WN+FN、WN+GM、WN+RW+GM和WN+PL四種噪聲模型與WN模型的差值大于其他兩種噪聲模型,說明這4種噪聲模型較其它模型更優(yōu)越,且這四種噪聲模型的MLE值相差不大．

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

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

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

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

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

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

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ù)估計不確定度．

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

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

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

圖8 顧及CN的SatRef站速度場

6 結(jié)束語

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

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

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

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

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