席換，曲國慶，張建霞，王暉

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

坐標時間序列的噪聲信息會影響坐標的解算精度，造成測站非線性運動模型一定的偏差，在研究連續(xù)運行參考站(continuously operating reference station,CORS)單日解坐標時間序列的噪聲時，有學者認為坐標時間序列中僅存在白噪聲。隨著對噪聲進一步地研究發(fā)現(xiàn)，GNSS坐標時間序列不僅存在白噪聲還存在有色噪聲，目前最優(yōu)隨機模型描述為白噪聲+閃爍噪聲已達成共識[1-8]。唐江森等[9]研究了山東部分CORS站兩年坐標時間序列的噪聲特征，確定其最佳噪聲模型為白噪聲+閃爍噪聲。由于研究數(shù)據(jù)時間跨度短，并不能完全反映站點的噪聲特征。通常認為時間序列跨度大于2.5 a，能夠準確估計時間序列的線性項、周期項及其精度[10]。

因此，本文在已有研究的基礎上，以山東CORS站2015年1月1日至2018年12月31日連續(xù)4年的坐標時間序列為研究對象，分析CORS站坐標時間序列的噪聲特征，以確定山東CORS各站點的最佳噪聲模型。

1 噪聲特征分析

1.1 噪聲序列的獲取

GPS基準站的擬合模型[11]可用式(1)表示：

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

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

(1)

式中：y(ti)為基準站各分量坐標時間序列；a為觀測序列的初始位置；b為觀測序列線性速度；ti為時間；c、d和e、f分別為年、半年周期項系數(shù)；gi為tgj時刻的階躍；H為一階梯函數(shù)；vi為噪聲序列。若顧及更復雜的情況，如某時刻測站運動速率的改變或某一事件發(fā)生后測站運動速率呈指數(shù)衰減的情況等，可采用文獻[12]提供的坐標分量每日解觀測序列參數(shù)模型。通過對坐標時間序列進行趨勢項擬合及去除、小波分析方法探測和提取周期項、修正含有階躍的站點等步驟得到噪聲序列，作為噪聲特征分析的數(shù)據(jù)基礎。

1.2 粗差的剔除與插補

坐標時間序列受到外界觀測條件、傳輸信號干擾或電氣元器件的突然變化的影響，會存在粗差。粗差的存在會影響噪聲特征的分析，應將其剔除，以免造成結果的偏差。采用拉依達準則(3σ)對坐標時間序列的粗差進行探測與剔除，該方法的標準為

(2)

連續(xù)均勻采樣的坐標時間序列是數(shù)據(jù)分析的基礎，而實際的觀測序列并不滿足此條件；加之剔除粗差的歷元，時間序列總會缺失某些天的數(shù)據(jù)，因此需利用插值方法將缺失的數(shù)據(jù)補全。通過對多種插值方法進行對比分析，最終選用最符合數(shù)據(jù)大致趨勢的三次多項式插值法進行插補。

線性趨勢項擬合實際上是對坐標時間序列進行零均值化處理，去除趨勢項的坐標時間序列稱為殘差序列。此序列在零附近連續(xù)均勻波動，是用于探測與提取周期項的數(shù)據(jù)序列。

1.3 周期項的探測與提取

小波分析方法是一種可以同時在時間域和頻率域分析信號的方法，實現(xiàn)坐標時間序列周期項的探測與提取。其基本原理如下：

(3)

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

AK=H*AK+1+G*DK+1，K=N,…,2,1，

(4)

式(4)為小波重構公式。式中H*、G*為H、G的共軛轉置。

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

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

P(f)∝fk，

(5)

稱k為譜指數(shù)。將式兩邊取對數(shù)為

lnP(f)∝klnf，

(6)

k在雙對數(shù)坐標系中為擬合直線的斜率。不同的譜指數(shù)對應不同的噪聲類型，k=0為白噪聲(wn)，k=-1為閃爍噪聲(fn)，k=-2為隨機漫步噪聲(rwn)。

最大似然估計是確定不同噪聲模型下噪聲分量的方法，它是按最大似然準則使得在此噪聲模型下噪聲序列與其協(xié)方差的概率密度最大，其模型[13]為

式中：X為時間序列列向量；σw,σrw,σf分別為白噪聲、隨機漫步噪聲和閃爍噪聲；N為序列長度；QXX為X的協(xié)方差。

2 山東CORS站噪聲特征分析

2.1 山東CORS站數(shù)據(jù)分析

以山東27個CORS站(圖1)三個方向的坐標時間序列為研究對象，在進行了數(shù)據(jù)預處理后，如剔除粗差與插值、階躍項修正和去除趨勢項得到殘差序列。在殘差序列的基礎上進行周期項的探測與提取，得到噪聲時間序列，用于確定噪聲類型與估計各噪聲分量值，分析山東CORS站的噪聲特征。

圖1 山東CORS站分布

以QUFU站為例，為顯示坐標軸刻度，將N、E和U三個方向分別減去參考值(3 962 977 m、10 587 692 m、59 m)，圖2、圖3分別給出了原始坐標時間序列及趨勢項和殘差序列。由圖2、圖3可知，QUFU站坐標時間序列水平方向存在著明顯的線性趨勢。垂直方向不僅存在著線性變化，還存在著周期性變化。殘差序列在零附近均勻的波動，說明已具有零均值的特性，可進行坐標時間序列的周期分析。

圖2 原始坐標序列及趨勢項

對QUFU站垂直方向進行小波分析，可得到坐標時間序列中存在的周期項(年、半年、季節(jié)周期項)，如圖4所示。

由此可得到QUFU站垂直方向的噪聲序列，如圖5所示。按照上述方法可得到山東CORS站的全部噪聲序列，為噪聲特征分析提供數(shù)據(jù)基礎。

2.2 山東CORS站噪聲特征分析

根據(jù)式(6)譜指數(shù)的定義，計算求得基準站各坐標分量的譜指數(shù)，見表1。

表1 基準站各方向譜指數(shù)

由表1可知，基準站水平方向的譜指數(shù)均在-1～0之間，大部分垂直分量譜指數(shù)也在-1～0之間，ZAZH站垂直方向譜指數(shù)為-1.011，說明基準站各坐標分量噪聲均不具有純白噪聲的特性，且同一測站的不同坐標分量的噪聲模型有所不同。因此需根據(jù)譜指數(shù)逐個確定基準站的噪聲模型。為確定基準站坐標分量的最佳噪聲模型，根據(jù)最大似然估計原理，不同噪聲模型的極大似然估計值越大，結果越可靠。本文以ZAZH站為例，計算了五種噪聲模型：白噪聲(wn)、白噪聲+閃爍噪聲(wn+fn)、白噪聲+隨機漫步噪聲(wn+rwn)、閃爍噪聲+隨機漫步噪聲(fn+rwn)和白噪聲+閃爍噪聲+隨機漫步噪聲(wn+fn+rwn)的最大似然估計值，結果如圖6所示。

由圖6可知，不同噪聲模型對應的最大似然估計值不同，白噪聲模型的最大似然估計值最小，說明坐標分量中不僅僅含有白噪聲，還含有有色噪聲。水平方向上，白噪聲+閃爍噪聲和白噪聲+閃爍噪聲+隨機漫步噪聲的最大似然估計值最大且相等，說明在沒有隨機漫步噪聲存在的情況下，白噪聲+閃爍噪聲與白噪聲+閃爍噪聲+隨機漫步噪聲所估計出來的噪聲分量是相同的，因此，ZAZH站水平方向的最佳噪聲模型為白噪聲+閃爍噪聲或白噪聲+閃爍噪聲+隨機漫步噪聲；垂直方向上，最大似然估計值為-4 202.15，所對應的最佳噪聲模型為白噪聲+閃爍噪聲+隨機漫步噪聲。由此可以看出同一基準站的不同方向可以有不同的最佳噪聲模型。假設基準站的噪聲模型為白噪聲+閃爍噪聲+隨機漫步噪聲，對山東地區(qū)27個基準站坐標分量在此噪聲模型下進行的噪聲分量估計統(tǒng)計見表2。

圖6 ZAZH站不同噪聲模型的最大似然估計值

表2 山東CORS站在wn+fn+rwn下的噪聲分量值

由表2可以看出，觀測時間序列坐標分量具有不同噪聲特征。在N方向上，有25個基準站僅含有白噪聲和閃爍噪聲，2個站含有白噪聲、閃爍噪聲和隨機漫步噪聲；E方向上，有20個站點含有白噪聲和閃爍噪聲，7個站點含有白噪聲、閃爍噪聲和隨機漫步噪聲；U方向上有16個站點含有白噪聲和閃爍噪聲，11個站點含有白噪聲、閃爍噪聲和隨機漫步噪聲。因此山東CORS站的最佳噪聲模型為：N方向上采用wn+fn模型，E和U方向上采用wn+fn+rwn模型。同時，U方向的噪聲分量估計值要明顯高于其它兩個方向，這與高程方向精度低于水平方向精度的結論是一致的。

3 結論

通過對坐標時間序列預處理得到連續(xù)均勻、不含階躍項的殘差序列，再利用小波分析方法提取其周期項，得到不含明顯周期項的噪聲時間序列，用于分析山東CORS的噪聲特征,最終得到以下結論：

1)計算出的譜指數(shù)顯示坐標時間序列中不僅存在白噪聲，還存在有色噪聲。在計算基準站的運動速度時，應顧及有色噪聲對其的影響。

2) 山東CORS各坐標分量具有不同的噪聲特征，可以有不同的噪聲模型。N方向的最佳噪聲模型為wn+fn，E、U方向的最佳噪聲模型為wn+fn+rwn。