陶國(guó)強(qiáng)

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

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

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

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

(1)

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

y=Ax+e

(2)

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

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

(3)

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

2 基于SSA的GNSS坐標(biāo)時(shí)間序列粗差探測(cè)與噪聲估計(jì)

2.1 SSA基本原理

對(duì)于時(shí)間序列{y1,y2,…,yN}，首先選擇一個(gè)合適的窗口L，將其進(jìn)行滯后排列：

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

(4)

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

C=VΛVT

(5)

式中，Λ∈RL×L為對(duì)角陣，其對(duì)角元λk(1≤k≤L)為C的奇異值，V=(v1,v2,…,vL)為正交陣，vk為第k個(gè)特征向量。可通過(guò)V和Y確定主成分矩陣：

A=VY

(6)

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

(7)

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

(8)

則原序列yi可表示成：

(9)

2.2 粗差探測(cè)

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

原坐標(biāo)時(shí)間序列的粗差大都反映在ε中，將其中的元素進(jìn)行升序排列?？紤]到GNSS坐標(biāo)時(shí)間序列的精度隨歷元而變化，采用開(kāi)窗IQR準(zhǔn)則進(jìn)行粗差探測(cè)，其判定準(zhǔn)則為：

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

(10)

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

(11)

式中，d=c/3。對(duì)殘差序列ε進(jìn)行遍歷，以當(dāng)前值為中心，構(gòu)建窗口序列，當(dāng)滿足準(zhǔn)則(11)時(shí)，即認(rèn)為當(dāng)前歷元觀測(cè)值為粗差。

2.3 噪聲分量估計(jì)

(12)

式中，y″表示生成的新序列。考慮到式(12)嚴(yán)格滿足Gauss-Markov模型，因此可采用LS_VCE方法估計(jì)噪聲分量：

(13)

(14)

3 算例分析

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

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

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

(15)

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

圖1 模擬的坐標(biāo)時(shí)間序列Fig.1 The simulated coordinate time series

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

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

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

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

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

圖6 500次模擬實(shí)驗(yàn)2種算法估計(jì)的噪聲分量估值絕對(duì)誤差Fig.6 The absolute error of noise components estimationsderived by two algorithms for each of 500 simulations

3.2 實(shí)測(cè)數(shù)據(jù)

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

表2 站點(diǎn)信息

圖7 LUZH站坐標(biāo)時(shí)間序列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種算法估計(jì)的信號(hào)Fig.9 Signals estimated by two algorithmsat LUZH station

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

表3 LUZH站噪聲分量估值

圖10 2種算法探測(cè)剔除粗差后的坐標(biāo)時(shí)間序列Fig.10 The coordinate time series after removinggross error by two algorithms

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

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

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

4 結(jié) 語(yǔ)

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