徐良葉，邵德盛，吳學群，牛 甜，韓 逍

（1. 云南省地震局，云南 昆明 650200；2. 昆明理工大學 國土資源工程學院，云南 昆明 650093）

0 引言

連續(xù)GNSS 時間坐標序列中存在一種空間相關(guān)的誤差稱之為共模誤差（commno Mode Error，CME），共模誤差被普遍認為是GNSS 坐標時間序列的主要來誤差源之一，一些弱構(gòu)造或瞬態(tài)構(gòu)造信號會被共模誤差所掩蓋，它會影響區(qū)域網(wǎng)中站點形變特征獲取的準確性[1-3]。關(guān)于共模誤差的起源，還有待進一步研究，已有研究表明它可能是由以下幾個方面引起的[4-5]：軟件、算法及數(shù)據(jù)處理策略等不完善而引起的系統(tǒng)誤差；GNSS 衛(wèi)星軌道誤差，海潮改正殘差，地球旋轉(zhuǎn)參數(shù)誤差等模型誤差；環(huán)境負載，小規(guī)模（區(qū)域性） 的地殼變形影響。

坐標時間序列中的共模誤差可以通過PCA/KLE 相結(jié)合的濾波方法有效地提取出來[6-9]，為了更準確的提取共模誤差，需要剔除本地效應(yīng)強的站點，而被剔除的站點時間序列里仍然包含有用的構(gòu)造運動信息，相關(guān)方面的研究較少，所以本文基于云南2015—2019 年GNSS 坐標時間序列，首先根據(jù)2 倍標準差的方法剔除本地效應(yīng)強的站點，利用PCA/KLE 法來提取云南GNSS 時間序列的共模誤差；接著對本地效應(yīng)強的站點的時間序列進行去噪處理，結(jié)合地震信息，分析震前站點的形變信息。

1 PCA 和KLE 方法

主成份分析（Principal Component Analysis）PCA 也叫經(jīng)驗正交函數(shù)分析（Empirical Ortho—gonal function） 即EOF 或特征向量分析。它把區(qū)域網(wǎng)站點時間序列分解成時間域的主分量和空間域的特征向量，目的在于通過利用降維的思想把多指標降為少數(shù)幾個綜合指標[6]。

如果把所有站點的時間序列排列起來，1個站點時間序列為一列，可表為m×n的矩陣（通常m＞n）：

它是一個n×n的矩陣，這時X的方差矩陣可表為：

矩陣VT是由特征向量構(gòu)成的n×n的正交矩陣，式中Λ 矩陣是由k個特征值構(gòu)成的對角矩陣，這時V為X的正交基底，X可以由此正交基底表示為：

式中的展開系數(shù)ak（ti）的序列即為k階主分量，它可以由下式求出：

式中的ak（t）表示第k個主成分，vk（x）為對應(yīng)主成分的響應(yīng)特征矩陣，ak（t）代表時間特征，vk（x）代表空間響應(yīng)。將B 矩陣的特征值λ 按從大到小的排列，特征值越大代表對殘余時間序列方差的貢獻越大[10]。可以用特征值累積貢獻率mk來代表主成分對矩陣的貢獻率，mk可由下式表示：

相關(guān)矩陣C可由協(xié)方差矩陣B標準得到，即dij=bij/（σiσj），其中σ 的定義為：

這個過程就稱為Karhunen-Loeve 展開（KLE）[11]。

PCA 和KLE 的唯一區(qū)別就是用來計算正交向量的矩陣不同。PCA 算法中使用的是協(xié)方差矩陣，KLE 算法中使用的是標準化后的相關(guān)系數(shù)矩陣。由于PCA 方法能得到精度較高的主成分值和空間響應(yīng)，但如果本地噪聲較大的話，則無法準確地分辨出CME，而KLE 方法正好彌補了這個缺點，它可以抑制本地噪聲的影響，并提取出CME，它的缺點是得不到準確的空間相應(yīng)[6，8]。

2 數(shù)據(jù)處理

本文數(shù)據(jù)來自于中國地震局第一監(jiān)測中心解算并免費提供的GNSS 數(shù)據(jù)產(chǎn)品（http://www.eqdsc.com/)。選取云南及周邊2015—2019 年數(shù)據(jù)比較完整的33 個GNSS 基準站的觀測資料做分析（數(shù)據(jù)缺失率小于10%）。站點分布如圖1所示。

圖1 云南及周邊站點分布圖Fig.1 Distribution of stations in Yunnan and surrounding areas

2.1 獲取GNSS 時間序列殘差

GNSS 臺站坐標每日觀測時間序列擬合函數(shù)為[12]：

式中ti表示歷元（單位：年），a為初始位置，b為速度，c、d為半周年系數(shù)，e、f為年數(shù)為Tj時刻的階躍式偏移量，vi為殘差序列。

使用QOCA 軟件去除原始坐標時間序列中的階躍、速度、周年、半周年項，得到坐標殘差時間序列vi，限于篇幅本文以彌勒站（YNML）為例，如圖2 所示。

圖2 YNML 站原始坐標時間序列（左圖） 及殘差坐標時間序列（右圖）Fig.2 Original coordinate time series of YNML Station（left） and residual coordinate time series（right）

2.2 異常站點剔除與共模誤差提取

如果沒有剔除本地效應(yīng)較強的站點，會導致PCA 結(jié)果有一定程度的偏差，文獻[6]研究表明，本地效應(yīng)較強的站點具有顯著的空間響應(yīng)，且分布不均勻，此時PCA結(jié)果的PC1（第一主成份） 對應(yīng)的空間響應(yīng)與KLE 結(jié)果與大多數(shù)站相比有明顯的差異。所以本文根據(jù)PCA結(jié)果與KLE 結(jié)果的PC1 對應(yīng)的空間響應(yīng)差值的2 倍標準差來剔除異常站點，如下式所示：

根據(jù)2.1 節(jié)獲取的殘差時間序列作為PCA/KLE 的輸入數(shù)據(jù)，由（13）式來剔除異常站點，再重復進行PCA/KLE 計算來提取PC1，直到?jīng)]有異常站點。本文循環(huán)了4 次，共剔除了9 個站點，剩余24 個站點，將剩余站點通過PCA方法提取共模誤差。

PCA提取前3 個主成分在N、E、U 方向的貢獻率分別為67.25%、4.07%、3.81%；43.04%、8.55%、5.90%；41.69%、6.74%、5.59%，在N、E、U 方向的累計貢獻率分別為75.13%、57.49%、53.95%，可見主要的方差集中在PC1上。由PCA 得到的前3 個主成分空間響應(yīng)如圖3 所示，PC1 的空間響應(yīng)分布均勻且都為正，在N、E、U 方向的均值為0.8867、0.8798、0.8622，PC2、PC3 的空間響應(yīng)分布不均勻，有正有負，由于共模誤差表示空間相關(guān)的誤差，PC1 的空間響應(yīng)一致性較好，所以本文選取PC1 作為共模誤差。

圖3 PCA 解算的前3 個主成分對應(yīng)的空間響應(yīng)Fig.3 Spatial responses of the first three principal components solved by PCA

PCA 提取的共模誤差如圖4 所示（以YNML 站為例），E、N 方向的共模誤差大多數(shù)在-8～8mm 之間，U 方向的共模誤差大多數(shù)在-14～14mm 之間，說明云南GNSS 坐標時間序列中存在比較明顯的共模誤差。

圖4 PCA 提取的共模誤差Fig.4 Common mode error extracted by PCA

對比PCA 去除共模誤差前與去除共模誤差后的坐標時間序列（圖5） 可以看出，去除共模誤差后，N、E、U 方向的坐標時間序列波動幅度減小了，并趨于平穩(wěn)。

圖5 YNML 站PCA 濾波前和濾波后的坐標時間序列Fig.5 Coordinate time series of YNML station before and after PCA filtering

如圖6（藍色小三角形表示濾波前的RMS，紅色小三角形表示濾波后的RMS） 所示，PCA濾波前的坐標時間序列在在N、E、U 方向的均方差（RMS） 為2.27mm、2.25mm、6.14mm，濾波后為1.07mm、1.39mm、3.98mm，說明該方法很好地剔除了共模誤差，提高了坐標精度。

圖6 PCA 濾波前和濾波后的坐標時間序列均方差Fig.6 Mean square deviation of coordinate time series before and after PCA filtering

2.3 異常站點共模誤差提取及形變分析

使用KLE 濾波能夠有效地減小本地噪聲的影響，并提取出CME，所以本文基于KLE 方法提取本地效應(yīng)較強的站點（2.2 節(jié)中被剔除的異常站點） 坐標時間序列中的CME，如圖7 所示，以新平站（YNXP）為例，圖8 為去除共模誤差后的坐標時間序列，N、E、U 方向的坐標時間序列比之前穩(wěn)定。KLE 濾波前YNXP 站坐標時間序列在N、E、U 方向的均方差（RMS）為3.22mm、4.12mm、6.31mm，濾波后為2.10mm、2.76mm、4.73mm，說明剔除CME 后站點的坐標精度和可靠性得到提高。

圖7 YNXP 站原始坐標時間序列（左圖） 與殘差坐標時間序列（右圖）Fig.7 Original coordinate time series（left） and residual coordinate time series of YNXP Station（right）

研究表明[13]，去除共模誤差后的坐標時間序列依舊含有噪聲，包括白噪聲、閃爍噪聲、隨機漫步噪聲，其中的兩種或三種混合噪聲等。而自適應(yīng)卡爾曼濾波（Kalman） 的前提要求是數(shù)據(jù)中的噪聲為高斯分布的白噪聲[14]，小波變換被譽為“數(shù)學顯微鏡”，它能從時頻域的局部信號中有效地提取信息，具有較好的去噪特點[15]。

所以本文利用文獻[16]中的小波變換和自適應(yīng)Kalman 濾波算法，將2.3 節(jié)中去除共模誤差后的坐標時間序列進行小波變換處理，再將處理后的結(jié)果作為自適應(yīng)Kalman 濾波處理的輸入數(shù)據(jù)，其目的是更科學、有效地降低坐標時間序列中的噪聲，最終處理后的結(jié)果（圖9）與沒有進行去噪處理的坐標時間序列（圖8 中右圖）相比，可以看出，去噪后坐標時間序列曲線更光滑，更能直觀地凸顯地殼的形變信息。

圖8 YNXP 站KLE 濾波前坐標時間序列（左圖） 與濾波后坐標時間序列（右圖）Fig.8 Coordinate time series before YNXP Station KLE filtering（left） and after filtering coordinate time series（right）

圖9 YNXP 站去噪后的坐標時間序列Fig.9 Coordinate time series of YNXP station after noise remova

為了驗證該方法的合理性，本文對濾除的噪聲（即去噪之前的坐標時間序列與去噪后的時間序列的差值） 進行統(tǒng)計分析，如圖10 所示，濾除的噪聲符合正態(tài)分布，說明本文去除噪聲方法的可靠性。

圖10 去除的噪聲統(tǒng)計直方圖Fig.10 Histogram of removed noise statistics

本文結(jié)合2018 年8 月13 日通海5.0 級地震，分析地震前YNXP 站點的異常形變信息，如圖11 所示，圖中的紅線表示坐標時間序列的均值加減2 倍標準差所對應(yīng)的限差值，明顯可以看出從2018 年5 月份開始，YNXP 站點在N、E 方向的運動趨勢嚴重偏離了均值線，并超過了2 倍的標準差，N 方向的位移在2018 年5月30 日達到最大值-4.82mm，E 方向的位移在2018 年6 月17 日達到最大值-5.54mm；同樣從圖7 中也可以看出YNXP 站點在N、E 方向原始坐標時間序列和殘差坐標時間序列均出現(xiàn)了凹槽，地震后站點恢復到原有的運動趨勢。相關(guān)工作人員在2018 年5 月對站點進行異常核實，排除了數(shù)據(jù)質(zhì)量問題和周邊環(huán)境影響，認為YNXP 站點的變化為真實的構(gòu)造運動引起。說明此方法能凸顯地殼形變和構(gòu)造運動信息，為提取GNSS 坐標時間序列中的信息提供了參考。

圖11 YNXP 站去噪后的坐標時間序列Fig.11 Coordinate time series of YNXP station afternoise removal

3 結(jié)論

本文采用PCA 和KLE 方法提取云南GNSS坐標時間序列的共模誤差，然后對坐標時間序列進行去噪處理，結(jié)合2018 年通海5.0 級地震，分析地震前站點異常形變信息。

結(jié)果表明，本文提出的定量去除本地效應(yīng)較強的站點的方法能更準確地從GNSS 時間序列中去除共模誤差，有效地提高了坐標精度；對減小本地效應(yīng)誤差后的時間序列進行去噪處理后，坐標時間序列曲線更加光滑，通海5.0級地震前YNXP 站點在N、E 方向的運動趨勢嚴重偏離了均值線，加速向西南方向運動，并超過了2 倍的標準差，說明此方法能凸顯地殼形變和構(gòu)造運動信息，為提取GNSS 坐標時間序列中的信息提供了參考和新的思路。