基于奇異譜分析的GNSS坐標時間序列粗差探測與噪聲估計

陶國強

1 東華理工大學(xué)長江學(xué)院，江西省撫州市學(xué)府路56號，344000

1 傳統(tǒng)諧波模型

GNSS基準站連續(xù)坐標時間序列中蘊含著豐富的信號與噪聲。其中信號包含線性信號與非線性信號，線性信號代表基準站由構(gòu)造應(yīng)力導(dǎo)致的構(gòu)造運動，一般以速度場形式表示；非線性信號主要表現(xiàn)為由環(huán)境負載導(dǎo)致的測站季節(jié)性位移。噪聲則包含白噪聲與有色噪聲。一般以諧波模型來描述GNSS坐標時間序列：

dkcos(2kπti)+e(ti)

(1)

式中，y(ti)為ti歷元時刻的基準站坐標，a、b為測站初始位置和速度，ck、dk(k=1,2)為年和半年周期振幅，e(ti)為觀測誤差。若顧及所有觀測歷元ti(i=1,2,…,N)，則式(1)可寫為：

y=Ax+e

(2)

式中，y、e∈RN×1為觀測值及其誤差向量，A∈RN×6為系數(shù)矩陣，x∈R6×1為未知參數(shù)。各變量的構(gòu)成為：

則參數(shù)的最小二乘估值為：

(3)

式中，Σy為觀測值的方差-協(xié)方差陣。由式(3)可知，要求得最優(yōu)估值，需已知觀測值的協(xié)方差陣。目前公認的最優(yōu)噪聲模型為白噪聲+閃爍噪聲的組合，因此需要估計各個噪聲的噪聲分量，可采用方差分量估計方法確定。

2 基于SSA的GNSS坐標時間序列粗差探測與噪聲估計

2.1 SSA基本原理

對于時間序列{y1,y2,…,yN}，首先選擇一個合適的窗口L，將其進行滯后排列：

模板質(zhì)量控制之一就是要求底面模板的曲線精度符合內(nèi)拱底緣坐標要求，而且具有足夠的強度和剛度，保證在安裝鋼筋及澆筑混凝土過程中不會變形。為此，首先在模架系統(tǒng)上方鋪設(shè)10 cm×10 cm方木，在其上方安裝底模。底模安裝完成后，由測量員對底模高程、曲度進行校核定位，將調(diào)整后的拱肋底模進行加固，防止變形。

(4)

式中，K=N-L+1，Y∈RK×L為滯后矩陣。其方差-協(xié)方差陣為C=YTY，對其進行奇異值分解：

C=VΛVT

(5)

式中，Λ∈RL×L為對角陣，其對角元λk(1≤k≤L)為C的奇異值，V=(v1,v2,…,vL)為正交陣，vk為第k個特征向量?？赏ㄟ^V和Y確定主成分矩陣：

A=VY

(6)

A的第k個行向量被稱為第k個主成分，其第i個元素可由式(7)計算：

(7)

式中，vk,j為vk的第j個元素。由第k個主成分可重構(gòu)第k個分量：

(8)

則原序列yi可表示成：

(9)

2.2 粗差探測

已有大量研究表明，GNSS坐標時間序列中的粗差主要反映在其噪聲(殘差)中，因此可直接對噪聲序列進行粗差探測。為避免殘留信號對粗差探測的干擾，需將信號與噪聲進行有效分離。經(jīng)SSA處理后，GNSS坐標時間序列已被分解為若干能量不等的分量，其中能量較高(軌跡矩陣對應(yīng)分量的特征值較大)的表示低頻信號，能量較低(軌跡矩陣對應(yīng)分量的特征值較小)的表示高頻噪聲。本文采用特征值比例法進行信噪分離，取總能量占比90%以上的特征值對應(yīng)的重構(gòu)分量作為信號s，剩下的分量作為噪聲ε。

原坐標時間序列的粗差大都反映在ε中，將其中的元素進行升序排列。考慮到GNSS坐標時間序列的精度隨歷元而變化，采用開窗IQR準則進行粗差探測，其判定準則為：

|εi-med(εi-L/2,εi+L/2)|>c·IQR(εi-L/2,εi+L/2)

(10)

式中，εi為ε的第i個分量，med(·)為求中位數(shù)算子，IQR(·)為求四分位距算子，L為窗口長度(一般取182)，c為比例因子(一般取3)。因基于IQR準則的粗差探測有較強的邊緣效應(yīng)，式(10)僅適用于序列中部的粗差探測，對于序列兩端的粗差探測效果較差。本文對式(10)進行改進，對于邊緣數(shù)據(jù)適當降低c值，將邊緣數(shù)據(jù)非邊緣化，則有：

(11)

式中，d=c/3。對殘差序列ε進行遍歷，以當前值為中心，構(gòu)建窗口序列，當滿足準則(11)時，即認為當前歷元觀測值為粗差。

2.3 噪聲分量估計

(12)

式中，y″表示生成的新序列?？紤]到式(12)嚴格滿足Gauss-Markov模型，因此可采用LS_VCE方法估計噪聲分量：

(13)

(14)

3 算例分析

3.1 模擬數(shù)據(jù)

為了驗證本文算法的正確性，采用模擬的坐標時間序列數(shù)據(jù)進行驗證。模擬策略同文獻[9]，即認為GNSS坐標時間序列由3個部分組成：

y(ti)=s1(ti)+s2(ti)+ε(ti)

(15)

表1 模擬時間序列參數(shù)

圖1 模擬的坐標時間序列Fig.1 The simulated coordinate time series

采用SSA對模擬的坐標時間序列進行分析，窗口長度取365。根據(jù)SSA原理構(gòu)建軌跡矩陣，并求其協(xié)方差陣的特征值，結(jié)果見圖2。由圖可見，前4個特征值占總特征值的91%，因此可選擇前4個特征值對應(yīng)的RC作為信號，剩余RC作為噪聲。圖3和圖4分別為SSA和諧波模型分離出的信號與噪聲，顯然SSA提取出的信號更能反映出時間序列的真實非線性變化，其殘差也相對較為平穩(wěn)；而諧波模型僅能粗略地描述時間序列，其殘差中仍然含有部分信號。此外，原始序列中的粗差大都傳遞至殘差中。分別采用傳統(tǒng)算法和本文算法探測殘差中的粗差，結(jié)果見圖5?？梢钥闯觯瑐鹘y(tǒng)算法僅能探測到83%的粗差，而本文算法可探測到94%的粗差。

圖2 軌跡矩陣協(xié)方差陣的特征值及其占總特征值之比Fig.2 The eigenvalues of the covariance matrix of the trajectory matrix and its ratio to the total eigenvalues

圖3 諧波模型和SSA提取的信號Fig.3 Signals extracted by harmonic model and SSA

圖4 諧波模型和SSA提取信號后剩余的殘差(噪聲)Fig.4 Residual (noise) after filtering out signals using harmonic model and SSA

圖5 2種算法探測出的粗差與模擬粗差的對比Fig.5 Comparison of gross error detected by two algorithms and simulated gross error

圖6 500次模擬實驗2種算法估計的噪聲分量估值絕對誤差Fig.6 The absolute error of noise components estimationsderived by two algorithms for each of 500 simulations

3.2 實測數(shù)據(jù)

為進一步驗證本文算法的有效性，選取陸態(tài)網(wǎng)絡(luò)提供的10個基準站觀測數(shù)據(jù)進行分析。所有數(shù)據(jù)可通過中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(http://www.cgps.ac.cn/)獲取，表2為各站點的相關(guān)信息。以LUZH站為例，圖7為其1999～2019年的連續(xù)觀測坐標序列，其中含有大量的粗差，必須剔除。分別采用SSA和諧波模型對去趨勢后的坐標序列進行信噪分離，圖8為3個方向軌跡矩陣協(xié)方差陣的特征值及其占總特征值之比。由圖可知，N方向前3個特征值之和占總特征值之比達90.88%，E方向前5個特征值之和占總特征值之比達90.05%，U方向前10個特征值之和占總特征值之比達90.03%。由此可確定，信號分量和噪聲分量的分界層分別為3(N)、5(E)和10(U)。2種算法估計的信號如圖9所示。可以看出，傳統(tǒng)諧波模型僅能粗略地描述GNSS坐標時間序列，特別是對于水平方向的坐標序列，其擬合度較差；而SSA不受正弦波的假定約束，它從時間序列的動力重構(gòu)出發(fā)，能自適應(yīng)地提取出時間序列的信號，因此其能較好地描述GNSS坐標時間序列的非線性變化。

表2 站點信息

圖7 LUZH站坐標時間序列Fig.7 The coordinate time series at LUZH station

圖8 軌跡矩陣協(xié)方差陣的特征值及其占總特征值之比Fig.8 The eigenvalues of the covariance matrixof the trajectory matrix and its ratio to the total eigenvalues

圖9 LUZH站2種算法估計的信號Fig.9 Signals estimated by two algorithmsat LUZH station

將信號從原坐標序列中去除，分別采用本文算法和傳統(tǒng)算法對殘差序列進行粗差探測，圖10為2種算法探測粗差剔除后的坐標序列。由圖可見，采用本文算法處理后的坐標時間序列中幾乎沒有粗差，序列整體較為干凈；而采用傳統(tǒng)算法處理后的坐標時間序列仍存在一些小粗差(如1999～2003年)，這主要是因為諧波模型對該時段的擬合效果較差。分別采用2種算法估計坐標時間序列中的噪聲分量，結(jié)果如表3所示。從結(jié)果來看，傳統(tǒng)算法估計的噪聲分量略大于本文算法。

表3 LUZH站噪聲分量估值

圖10 2種算法探測剔除粗差后的坐標時間序列Fig.10 The coordinate time series after removinggross error by two algorithms

圖11 2種算法探測到的粗差比例Fig.11 The proportions of gross error detected by two algorithms

圖12 2種算法估計的白噪聲Fig.12 White noise estimated by two algorithms

圖13 2種算法估計的閃爍噪聲Fig.13 Flicker noise estimated by two algorithms

4 結(jié) 語

由于多種因素的影響，GNSS坐標時間序列受到各類噪聲和粗差的污染。為了能夠更加精確地估計地殼形變信號，需將其中混有的粗差進行有效地識別、剔除并且對噪聲進行準確的建模。以往進行粗差探測和噪聲估計的算法大都是基于諧波模型進行構(gòu)建的，然而諧波模型難以精確地描述GNSS坐標時間序列的非線性變化，進一步影響算法的性能。本文提出一種基于SSA的GNSS坐標時間序列粗差探測和噪聲估計算法，并采用模擬實驗對算法進行驗證。實驗結(jié)果表明，相比于傳統(tǒng)算法，本文算法的粗差探測成功率更高，且估計的噪聲分量與真值更為接近。陸態(tài)網(wǎng)絡(luò)實測數(shù)據(jù)分析結(jié)果表明，高程方向時間序列噪聲強度要高于水平方向，且閃爍噪聲的振幅要高于白噪聲。此外，由于諧波模型不正確的模型假設(shè)導(dǎo)致所獲得的殘差序列中仍含有部分信號，進一步導(dǎo)致其估計的噪聲振幅略高于本文算法。