張雙成 李 軍 安寧康 馮智杰 呂佳明 王 杰 葉志磊

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院,西安市雁塔路126號, 710054

2 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室,西安市雁塔路中段1號,710054

高精度GNSS坐標(biāo)序列中蘊(yùn)含的構(gòu)造運(yùn)動和非構(gòu)造運(yùn)動信息已廣泛應(yīng)用于地殼運(yùn)動研究、地下水儲量反演和位移探測等領(lǐng)域[1-2],但GNSS坐標(biāo)序列中存在著一種與時空特性相關(guān)的共模誤差(CME),會使測站位置和速度估計產(chǎn)生偏差,導(dǎo)致解譯出錯誤的地球物理信息。

目前,CME提取方法可分為區(qū)域?yàn)V波法、參考框架法和統(tǒng)計信號分解法3類,其中常用的統(tǒng)計信號分解法對CME提取效果最佳。區(qū)域?yàn)V波法前提需要假設(shè)CME空間分布均勻,但此假設(shè)與實(shí)際情況相悖。Wang等[3]基于PCA對CME提取進(jìn)行研究,結(jié)果表明,PCA相比于區(qū)域?yàn)V波法有更好的效果。然而,PCA只能對符合高斯分布的信號進(jìn)行提取,無法對具有非高斯特性的有色噪聲進(jìn)行識別。

ICA是一種顧及高階統(tǒng)計信號的盲源分離方法。劉斌等[4]將ICA與PCA應(yīng)用于垂向GNSS坐標(biāo)序列的時空濾波,結(jié)果表明,二者提取的共模分量有著不同的時空分布特征;李斐等[5]基于ICA對GNSS坐標(biāo)序列進(jìn)行矩陣特征分解,重構(gòu)共模分量,有效提取了坐標(biāo)序列中的CME,但I(xiàn)CA分量順序及其空間特性表現(xiàn)出隨機(jī)性,并且ICA只能對一種高斯信號進(jìn)行識別,這對共模分量確定和CME空間特征分析帶來嚴(yán)重困擾。

vbICA 是一種使用變分貝葉斯推理執(zhí)行ICA的信號分解方法。目前,該方法在大地測量時間序列中的瞬態(tài)信號和季節(jié)性信號提取方面得到初步應(yīng)用;Gao等[6]也初次使用該方法對云南、四川GNSS測站坐標(biāo)序列進(jìn)行有效濾波。針對PCA與ICA在CME提取方面的不足,本文采用vbICA對實(shí)驗(yàn)區(qū)20個GNSS測站坐標(biāo)序列進(jìn)行CME提取,并與PCA、ICA結(jié)果對比,分析vbICA方法的有效性和先進(jìn)性,最后基于vbICA濾波估計速度場變化。

1 原理與方法

1.1 vbICA原理

針對GNSS坐標(biāo)序列中的有色噪聲,建立帶噪聲的標(biāo)準(zhǔn)ICA混疊模型[7]:

x=As+e

(1)

式中,x為觀測坐標(biāo)序列矩陣,s=[s1,s2,…,sL]為信源向量,A為混合矩陣,e為高斯噪聲向量。

(2)

(3)

式中,m為獨(dú)立分量個數(shù),N為觀測歷元總數(shù),λk介于0～1區(qū)間。PCA的分量貢獻(xiàn)率由特征值計算。

對所有獨(dú)立分量貢獻(xiàn)率按降序排列,選取貢獻(xiàn)率最大的前p個獨(dú)立分量表示最顯著信號,最終CME可表示為:

(4)

1.2 性能評估指標(biāo)

1.2.1 距離相關(guān)系數(shù)

由于環(huán)境負(fù)載會造成基準(zhǔn)站周期性位移,致使 GNSS坐標(biāo)序列存在非線性信號。進(jìn)行站間相關(guān)性分析的常用方法為Pearson相關(guān)系數(shù)法,其在分析殘留有非線性信號的坐標(biāo)殘差序列時會導(dǎo)致相關(guān)性誤差。因此,為準(zhǔn)確反映扣除CME前后站間相關(guān)性變化,本文將采用能夠衡量非線性變量之間相關(guān)程度的距離相關(guān)系數(shù)法[8]。

1.2.2 均方根變化

(5)

2 數(shù)據(jù)來源與預(yù)處理

在PBO(plate boundary observatory)網(wǎng)絡(luò)中選取均勻分布于俄勒岡州、內(nèi)華達(dá)州、加利福利亞州內(nèi)20個連續(xù)運(yùn)行GNSS觀測站,對其2009-01～2019-01期間的坐標(biāo)序列進(jìn)行處理。原始數(shù)據(jù)來源于Nevada Geodetic Laboratory,該數(shù)據(jù)采用GipsyX軟件處理,同時加入極潮、海潮及固體潮等模型改正,并對電離層、對流層延遲、天線相位中心偏差進(jìn)行校正,最終得到IGS2014框架下GNSS測站的三維坐標(biāo)序列。經(jīng)統(tǒng)計,20個測站數(shù)據(jù)完整率均在98%以上。

2.1 原始數(shù)據(jù)預(yù)處理

由于受儀器故障和外界因素干擾,致使GNSS坐標(biāo)序列不可避免地存在隨機(jī)的數(shù)據(jù)缺失和粗差。首先,根據(jù)E、N、U方向上坐標(biāo)序列的中誤差設(shè)置閾值分別為20 mm、20 mm、30 mm進(jìn)行初次粗差剔除;然后,利用3倍四分位數(shù)法進(jìn)行二次粗差剔除;最終得到3個方向坐標(biāo)序列最大缺失率分別為1.48%、1.51%、1.51%,且均發(fā)生在P059測站。為獲得等時間間隔的坐標(biāo)序列,采用三次樣條插值對空缺位置數(shù)據(jù)進(jìn)行恢復(fù)。此外,由于強(qiáng)震同震位移或設(shè)備更換等原因,可能會導(dǎo)致原始坐標(biāo)序列中存在階躍現(xiàn)象。根據(jù)Nevada Geodetic Laboratory提供的站點(diǎn)log文件,P386測站點(diǎn)在2016-10-13由于更換天線導(dǎo)致階躍,其中N方向階躍明顯,約為10 mm。具體修正流程為:1)根據(jù)站點(diǎn)log文件確定站點(diǎn)階躍時刻;2)利用Hector軟件估計階躍大小;3)將階躍從原始坐標(biāo)序列中扣除。對于其他存在階躍的坐標(biāo)序列,采用同樣流程處理。

2.2 坐標(biāo)殘差序列獲取

經(jīng)過預(yù)處理后,原始坐標(biāo)序列將不含階躍項(xiàng)和粗差。此時,任意t時刻,單站單方向的GNSS坐標(biāo)序列y(t)可簡化表示為:

y(t)=a0+a1(t-t0)+a2(t)+a3(t)+ε(t)

(6)

式中,t0為起始?xì)v元,a0為初始位置,a1為線性趨勢,a2為周期運(yùn)動,a3為CME,ε為有色噪聲。根據(jù)式(6),對原始坐標(biāo)序列采用加權(quán)最小二乘擬合并扣除a1、a2,得到坐標(biāo)殘差序列。

3 結(jié)果與分析

3.1 CME提取分析

Dong等[9]基于PCA/KLE濾波方法對CME進(jìn)行研究,提出利用共模模式進(jìn)行CME判斷,即被作為共模分量的主分量(principal component,PC)具有最高的分量占比,可認(rèn)為是共模分量,貢獻(xiàn)率最大的分量綜合了原始坐標(biāo)序列最多的信息。圖1為3個方向所有分量的貢獻(xiàn)率統(tǒng)計結(jié)果。由圖可知,vbICA、PCA和ICA分解后,3個方向第一主分量PC1貢獻(xiàn)率均表現(xiàn)為最大,其中,E方向PC1貢獻(xiàn)率分別為94.06%、68.92%、73.52%;N方向PC1貢獻(xiàn)率為90.73%、61.21%、72.09%;U方向PC1貢獻(xiàn)率為88.62%、56.18%、67.47%,分量PC2～PC20貢獻(xiàn)率依次減小,且均小于10%。因此,本文以PC1作為共模分量重構(gòu)CME,結(jié)果如圖2所示。

圖1 E、N、U方向主分量貢獻(xiàn)率Fig.1 Contribution rate of PCS in E,Nand U directions

圖2 vbICA、PCA、ICA重構(gòu)的CME Fig.2 CME reconstructed by vbICA, PCA and ICA

為對比vbICA與PCA和ICA提取的CEM差異,分別計算3個方向上CME之間的Pearson相似度。統(tǒng)計結(jié)果如表1所示,可以看出,PCA與ICA的相似度較低,在E、U方向上相似度均低于0.6,根據(jù)ICA的高階統(tǒng)計特性,說明PCA所得共模分量未考慮CME的高階信息;vbICA重構(gòu)的CME與PCA和ICA的結(jié)果具有較強(qiáng)的一致性,其最大相似度為0.80,最低為0.62,說明相同歷元時刻的CME相近,且在整體表現(xiàn)為周期上的一致性。然而,vbICA與PCA在E方向的相似度較高,由于PCA存在潛在的“聚類”問題可能導(dǎo)致E方向PC1包含不同的物理模式,進(jìn)而導(dǎo)致CME存在虛假信息,對于該現(xiàn)象后文將根據(jù)相關(guān)評價指標(biāo)進(jìn)一步分析。整體而言,以上結(jié)果初步證實(shí)了vbICA在提取CME方面的有效性。

表1 CME 相似度

3.2 時空特征分析

對3個方向20個主分量的空間響應(yīng)(spatial response, SR)進(jìn)行標(biāo)準(zhǔn)化。圖3分別為vbICA、PCA和ICA三種方法在E方向上的共模分量PC1及其空間響應(yīng)SR1,箭頭向上表示SR為正,向下為表示SR為負(fù)。由圖可知,3種方法獲取的共模分量PC1具有共同的周期特征,且所有站點(diǎn)的SR1表現(xiàn)一致,均為正響應(yīng)或負(fù)響應(yīng);箭頭大小表示SR大小,不同站點(diǎn)的SR1大小不同,表明所有測站的CME在空間分布上的不均勻性。

圖3 E方向PC1及其SR1Fig.3 PC1 and SR1 in the direction E

對比圖3發(fā)現(xiàn),vbICA和PCA在E方向上空間響應(yīng)為正,ICA為負(fù),且vbICA和PCA第一分量PC1方向相反。該現(xiàn)象的產(chǎn)生與算法自身特性密切相關(guān):對于vbICA來說,信源模型選擇是利用vbICA進(jìn)行GNSS坐標(biāo)序列分解的核心,而不同信源模型初始化先驗(yàn)信息具有唯一性,最終表現(xiàn)為當(dāng)選定信源模型后,針對同一觀測數(shù)據(jù)無論vbICA分解多少次,主分量及其空間響應(yīng)符號均不會發(fā)生變化;而在ICA分解中,由于權(quán)重矩陣初始化采用的是隨機(jī)矩陣,故而導(dǎo)致循環(huán)收斂后的混合矩陣具有隨機(jī)性,最終重構(gòu)的PC分量順序也是隨機(jī)的;PCA基于正交分解方差最大原則實(shí)現(xiàn)降維,在進(jìn)行坐標(biāo)序列特征分解時,并沒有引入隨機(jī)變量,其結(jié)果與vbICA表現(xiàn)一致。結(jié)合圖2,雖然不同方法所得空間響應(yīng)SR符號和對應(yīng)PC符號存在差異,但最終并不會對CME重構(gòu)造成影響。

3.3 濾波效果分析

為分析vbICA與PCA、ICA的CME提取效果差異,本文從原始坐標(biāo)序列CME濾波前后均方根變化和測站間距離相關(guān)系數(shù)變化角度進(jìn)行分析。

根據(jù)式(5)計算E、N、U方向CME濾波后各測站坐標(biāo)殘差序列RMS減少百分比變化,由于篇幅限制,僅對任意10個測站結(jié)果進(jìn)行統(tǒng)計,如表2(單位%)所示。由表可知,3種方法均可有效降低各測站E、N、U方向坐標(biāo)殘差序列的RMS,說明CME濾波后測站坐標(biāo)序列離散程度降低、振幅減小。通過濾波前后3個方向RMS平均減少百分比可以發(fā)現(xiàn),所有測站E、N、U方向上的坐標(biāo)序列經(jīng)vbICA濾波后RMS平均減小百分比分別為36.57%、31.63%、10.97%;經(jīng)PCA濾波后RMS平均減小百分比分別為28.57%、22.77%、10.25%;經(jīng)ICA濾波后RMS平均減小百分比分別為14.90%、23.34%、10.44%,顯然vbICA結(jié)果明顯好于PCA和ICA。針對上述vbICA與PCA相似度高的問題,由RMS結(jié)果可知,PCA提取的CME存在虛假信號。

表2 濾波后各站RMS減少百分比

進(jìn)一步分析發(fā)現(xiàn),P730經(jīng)ICA濾波后E方向RMS變化為-0.16%,相比濾波前RMS增大,結(jié)合圖3(c)發(fā)現(xiàn),P730測站空間特性明顯;同理,vbICA和ICA對P166測站RMS改善也并不明顯,說明測站的局部效應(yīng)對于CME時空濾波具有不同程度的影響。垂直方向RMS變化明顯低于水平方向,表明同一測站CME對該測站垂直和水平方向影響并不一致,且水平方向影響較大,主要原因?yàn)?垂直方向除CME以外,還受到非潮汐海洋負(fù)荷、環(huán)境負(fù)載以及基巖熱膨脹效應(yīng)等非線性變化因素的影響。

根據(jù)各測站距離的遠(yuǎn)近,以較遠(yuǎn)的PABH測站為基準(zhǔn),計算并統(tǒng)計其余19個測站相對于基準(zhǔn)站在E、N、U方向上濾波前后的R變化百分比。整體而言,vbICA濾波后,E、N、U方向上R平均減小率分別為60.53%、56.84%、25.80%;ICA平均減小率為11.58%、33.54%、12.11%;PCA平均減小率為25.94%、22.79%、11.57%。由此可知,扣除vbICA提取的CME后,測站間相關(guān)性減弱程度明顯高于ICA和PCA。此外,R在水平方向上的變化趨勢較為明顯,而在垂直方向上變化差異較小,一方面說明CME的影響因素在水平方向和垂直方向上存在差異,另一方面說明時間跨度長的GNSS坐標(biāo)序列在水平方向受板塊運(yùn)動影響較大。

3.4 CME對速度場影響

GNSS坐標(biāo)序列CME濾波不僅可以獲得更加精細(xì)的形變信息,還可以優(yōu)化GNSS速度場,對于全球板塊運(yùn)動、mm級動態(tài)參考框架建立與維持等地學(xué)研究具有重要意義[1]。

根據(jù)黃立人等[10]相關(guān)學(xué)者研究結(jié)果,假設(shè)水平方向和垂向分別以FN+WN和 FN+PL為最優(yōu)噪聲模型,采用極大似然估計法分別對測站3個方向原始坐標(biāo)序列進(jìn)行趨勢項(xiàng)估計。表3統(tǒng)計任意10個測站經(jīng)vbICA濾波前后速度及其不確定度變化;圖4(a)、圖4(b)為濾波前后水平和垂向GNSS速度場及其不確定度變化,箭頭長短表示速度的相對大小,圓圈大小表示不確定度的相對大小。對比分析可知,水平方向上,該區(qū)域整體向西南方向運(yùn)動,且有渦旋趨勢;在垂向上整體表現(xiàn)為地表下沉。剔除CME后,20個測站3個方向速度場估計的標(biāo)準(zhǔn)差平均降低52.71%、49.88%、18.04%。其中,水平速度場精度改善明顯,不考慮個別特殊站,E方向測站速度平均變化不高于0.04 mm/a,N方向不高于0.01 mm/a。此外,濾波后P276測站的速度變化達(dá)到4.47 mm/a,P730變化達(dá)到了1.57 mm/a,且由圖4(b)可見,P305、P365測站運(yùn)動速率為正,即垂直向上運(yùn)動,說明垂向變化不僅與區(qū)域整體運(yùn)動有關(guān),局部形變因素也不容忽視。CME濾波可凸顯局部效應(yīng)較強(qiáng)的部分測站所隱藏的信號,提高GNSS速度場估計的可靠性,有助于地球物理現(xiàn)象的清晰解釋。

表3 CME濾波前后測站速度及其不確定度變化

圖4 基于vbICA的GNSS速度場Fig.4 GNSS velocity field based on vbICA

4 結(jié) 語

本文針對PCA、ICA在提取CME方面的不足,研究vbICA對 CME的提取。通過對濾波前后坐標(biāo)序列RMS值、距離相關(guān)性以及相似度等指標(biāo)進(jìn)行濾波效果分析,結(jié)果表明,vbICA方法濾波效果明顯優(yōu)于PCA和ICA;濾波后GNSS坐標(biāo)序列的精度明顯提高,且vbICA比ICA表現(xiàn)出更強(qiáng)的魯棒性。

致謝：感謝Nevada Geodetic Laboratory提供實(shí)驗(yàn)數(shù)據(jù)。