vbICA方法用于GNSS坐標序列共模誤差提取研究

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

1 長安大學地質工程與測繪學院,西安市雁塔路126號, 710054

2 地理信息工程國家重點實驗室,西安市雁塔路中段1號,710054

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

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

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

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

1 原理與方法

1.1 vbICA原理

針對GNSS坐標序列中的有色噪聲,建立帶噪聲的標準ICA混疊模型[7]:

x=As+e

(1)

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

(2)

(3)

式中,m為獨立分量個數,N為觀測歷元總數,λk介于0～1區間。PCA的分量貢獻率由特征值計算。

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

(4)

1.2 性能評估指標

1.2.1 距離相關系數

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

1.2.2 均方根變化

(5)

2 數據來源與預處理

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

2.1 原始數據預處理

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

2.2 坐標殘差序列獲取

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

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

(6)

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

3 結果與分析

3.1 CME提取分析

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

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

圖2 vbICA、PCA、ICA重構的CME Fig.2 CME reconstructed by vbICA, PCA and ICA

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

表1 CME 相似度

3.2 時空特征分析

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

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

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

3.3 濾波效果分析

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

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

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

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

根據各測站距離的遠近,以較遠的PABH測站為基準,計算并統計其余19個測站相對于基準站在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后,測站間相關性減弱程度明顯高于ICA和PCA。此外,R在水平方向上的變化趨勢較為明顯,而在垂直方向上變化差異較小,一方面說明CME的影響因素在水平方向和垂直方向上存在差異,另一方面說明時間跨度長的GNSS坐標序列在水平方向受板塊運動影響較大。

3.4 CME對速度場影響

GNSS坐標序列CME濾波不僅可以獲得更加精細的形變信息,還可以優化GNSS速度場,對于全球板塊運動、mm級動態參考框架建立與維持等地學研究具有重要意義[1]。

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

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

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

4 結 語

本文針對PCA、ICA在提取CME方面的不足,研究vbICA對 CME的提取。通過對濾波前后坐標序列RMS值、距離相關性以及相似度等指標進行濾波效果分析,結果表明,vbICA方法濾波效果明顯優于PCA和ICA;濾波后GNSS坐標序列的精度明顯提高,且vbICA比ICA表現出更強的魯棒性。

致謝：感謝Nevada Geodetic Laboratory提供實驗數據。