席換，曲國慶

(山東理工大學(xué)建筑工程學(xué)院，山東淄博 255049)

0 引言

GPS坐標(biāo)時間序列不僅反映坐標(biāo)分量的整體趨勢，還反映出基準(zhǔn)站的非線性變化.相關(guān)研究表明，由于全球衛(wèi)星導(dǎo)航系統(tǒng)（GNSS）基準(zhǔn)站坐標(biāo)時間序列受到地球物理因素及系統(tǒng)誤差的影響，不僅存在線性變化，還存在周期變化，尤其在高程方向表現(xiàn)得最為明顯[1-6].分析坐標(biāo)時間序列的非線性變化特征，有利于研究基準(zhǔn)站的穩(wěn)定性和基于非線性運(yùn)動模型獲得基準(zhǔn)站的瞬時位置.坐標(biāo)時間序列的噪聲信息會影響坐標(biāo)解算精度，造成測站非線性運(yùn)動模型有一定的偏差.GNSS坐標(biāo)時間序列不僅存在白噪聲（WN）還存在有色噪聲（CN），目前最優(yōu)隨機(jī)模型廣泛地被描述為WN+閃爍噪聲（FN）[7-14].就山東部分連續(xù)運(yùn)行參考站（CORS）而言，唐江森等[15]研究了兩年坐標(biāo)時間序列的噪聲特征，表明其最佳噪聲模型也為WN+FN.通常認(rèn)為時間序列跨度大于2.5年，才能準(zhǔn)確估計時間序列的線性項、周期項及其精度[16].

本文在已有研究的基礎(chǔ)上，以山東CORS 2015-01-01—2018-12-31連續(xù)四年的坐標(biāo)時間序列為研究對象，深入分析CORS站坐標(biāo)時間序列的周期與噪聲特征，以確定山東CORS站點(diǎn)的最佳噪聲模型.

1 數(shù)據(jù)處理策略

選取了較為均勻分布在山東地區(qū)的CORS站，如圖1所示.

圖1 站點(diǎn)分布

受到外界條件、多路徑效應(yīng)或傳輸信號干擾等因素的影響，觀測時間序列存在粗差，不僅影響數(shù)據(jù)質(zhì)量，而且影響站點(diǎn)速度的計算及噪聲特性的分析，有必要將其剔除.采用拉依達(dá)準(zhǔn)則法（3σ準(zhǔn)則）探測并剔除觀測序列的粗差.外界觀測條件、站點(diǎn)受損或網(wǎng)絡(luò)故障等都會影響時間序列的連續(xù)性和完整性，加之剔除粗差的歷元，坐標(biāo)時間序列會出現(xiàn)在某天或某幾天的觀測數(shù)據(jù)缺失的現(xiàn)象，經(jīng)過對多種插值方法反復(fù)實(shí)驗(yàn)對比，采用三次多項式插值補(bǔ)全數(shù)據(jù).

坐標(biāo)時間序列在較長時間內(nèi)受某些因素的影響，形成一種整體變動的趨勢，該趨勢使得坐標(biāo)時間序列為非平穩(wěn)序列.對坐標(biāo)時間序列進(jìn)行零均值化處理，實(shí)際上就是平穩(wěn)化的過程，其擬合模型如下：

式中:yt為坐標(biāo)時間序列；10?8分別為坐標(biāo)時間序列擬合的常數(shù)值和一次項系數(shù)值；10?9為去除趨勢項的殘差序列.

2 CORS站坐標(biāo)時間序列分析

2.1 周期特征分析

小波分析方法是一種同時在時、頻兩域內(nèi)分析時間序列的方法，實(shí)現(xiàn)坐標(biāo)時間序列周期項的探測與提取.其基本原理如下：

式(2)為小波分解公式.式中：AK，DK分別表示分解后第K層的低頻和高頻部分；H、G分別表示小波低通和高通濾波器；N表示信號長度.

式(3)為小波重構(gòu)公式.式中：H*、G*為H、G的共軛轉(zhuǎn)置.

通過小波分析提取出QUFU站的年、半年和季節(jié)周期項，如圖2所示.

由圖2可知，QUFU站N、E、U三個坐標(biāo)方向均表現(xiàn)出一定的周期特征，具有明顯的年周期項和半年周期項變化，整體周期性變化較為穩(wěn)定.對27個CORS站的周期特征進(jìn)行分析可得：1）U方向的半年周期變化特征最為明顯，其次是年周期特征，最后是3～4個月的季節(jié)周期特征；2）N、E方向具有明顯的年周期或半年周期，不同站點(diǎn)，差異較大；3）坐標(biāo)時間序列的周期信號并不是嚴(yán)格的年周期、半年周期等周期項，而是接近于年周期、半年周期等周期項.不同測站所表現(xiàn)出來的周期特征存在著差異，同一測站的不同方向坐標(biāo)分量的周期規(guī)律也有所不同，這是因?yàn)闇y站所處的地理位置不同，環(huán)境因素成為影響它們周期規(guī)律不同的主要因素.

圖2 QUFU站周期項

2.2 噪聲特征分析

2.2.1 噪聲序列的獲取

GPS基準(zhǔn)站的擬合模型[17]可用下式表示

式中：y(ti)為CORS站原始坐標(biāo)時間序列；a為觀測序列的初始位置；b為線性速度；ti為時間；c、d、e、f分別為年、半年周期項系數(shù)；gi為tgj時刻的階躍；H為階梯函數(shù)；vi為噪聲序列.若顧及更復(fù)雜的情況，如某時刻測站運(yùn)動速率的改變或某一事件發(fā)生后測站運(yùn)動速率呈指數(shù)衰減的情況等，可采用文獻(xiàn)[8]提供的坐標(biāo)分量每日解觀測序列參數(shù)模型.在數(shù)據(jù)處理策略中可實(shí)現(xiàn)坐標(biāo)時間序列線性趨勢項、階躍的去除，再利用小波分析提取坐標(biāo)時間序列的周期項，由此獲得坐標(biāo)分量的噪聲序列.

2.2.2 譜指數(shù)與最大似然估計

噪聲的功率譜密度P(f)與噪聲頻率f之間存在著某種冪次關(guān)系，即

式中，α為譜指數(shù).式(5)兩邊取對數(shù)為

明顯地，α在雙對數(shù)直角坐標(biāo)系中為擬合直線的斜率.不同的譜指數(shù)對應(yīng)不同的噪聲類型，α=0對應(yīng)WN，α=?1對應(yīng)FN，α=?2對應(yīng)RWN.

最大似然估計（MLE）是確定不同噪聲模型下噪聲分量的方法，使得在相應(yīng)噪聲模型下噪聲序列與其協(xié)方差的概率密度最大[18]，即

式中:X為時間序列列向量；σw,σrw,σf分別為WN、RWN和FN；N為序列長度；QXX為X的協(xié)方差.

根據(jù)式(6)譜指數(shù)的定義，計算求得CORS站各坐標(biāo)分量的譜指數(shù)，如表1所示.

由表1可知，CORS站水平方向的譜指數(shù)均在?1～0之間，垂直方向除ZAZH站外，譜指數(shù)均在?1～0之間，說明CORS站各坐標(biāo)分量噪聲均不具有純WN的特性，且同一測站不同方向的噪聲模型有所不同.為確定CORS站坐標(biāo)分量的最佳噪聲模型，根據(jù)最大似然估計準(zhǔn)則，以ZAZH站為例，計算了五種噪聲模型：WN、WN+FN、WN+RWN、FN+RWN和WN+FN+RWN的MLE值，結(jié)果如表2所示.

由表2可知，不同噪聲模型對應(yīng)的MLE值不同，WN模型MLE值最小，說明坐標(biāo)分量中不僅含有WN，還含有CN.水平方向上，WN+FN模型和WN+FN+RWN模型的MLE值最大且相等，說明在沒有RWN模型存在的情況下，WN+FN模型與WN+FN+RWN模型所估計出來的噪聲分量是相同的.因此，ZAZH站水平方向的最佳噪聲模型為WN+FN或WN+FN+RWN；垂直方向上，MLE值為?4 202.15所對應(yīng)的最佳噪聲模型為WN+FN+RWN.由此可以看出同一CORS站的不同方向可以有不同的最佳噪聲模型.假設(shè)CORS站的噪聲模型為WN+FN+RWN，在此模型下計算了山東地區(qū)27個CORS站的噪聲分量，如表3所示.

表2 ZAZH噪聲模型M LE 值

表3 WN+FN+RWN噪聲模型下的噪聲分量估計值

由表3可知，觀測時間序列坐標(biāo)分量具有不同的噪聲特征.N方向上，有25個站可用模型WN+FN描述，2個站可用模型WN+FN+RWN描述；E方向上，20個站可用模型WN+FN描述，7個站可用模型WN+FN+RWN描述；U方向上，16個站可用模型WN+FN描述，11個站可用模型WN+FN+RWN描述.因此山東CORS站的最佳噪聲模型為：N方向上采用WN+FN模型，E和U方向上采用WN+FN+RWN模型.同時，U方向的噪聲分量估計值要明顯高于其它兩個方向，這與高程方向精度低于水平方向精度的結(jié)論是一致的.

3 結(jié)束語

本文以山東27個CORS站的坐標(biāo)時間序列為研究對象，通過小波分析方法提取周期項進(jìn)行分析，得出了坐標(biāo)時間序列中并不嚴(yán)格地存在年周期項、半年周期項等周期信號，而是與之接近的周年、半周年項等周期信號；另外，在觀測序列的坐標(biāo)分量中，U方向的周期信號較水平方向明顯.

CORS站的譜指數(shù)表明坐標(biāo)時間序列中不僅存在WN，還存在CN.根據(jù)MLE方法，估計了CORS站各坐標(biāo)分量在WN+FN+RWN模型下的噪聲分量，得出了CORS站各坐標(biāo)分量可以有不同的噪聲模型，N方向的最佳噪聲模型為WN+FN，E、U方向的最佳噪聲模型為WN+FN+RWN的結(jié)論.