徐良葉，邵德盛，吳學(xué)群，牛 甜，韓 逍

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

0 引言

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

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

1 PCA 和KLE 方法

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

如果把所有站點(diǎn)的時(shí)間序列排列起來(lái)，1個(gè)站點(diǎn)時(shí)間序列為一列，可表為m×n的矩陣（通常m＞n）：

它是一個(gè)n×n的矩陣，這時(shí)X的方差矩陣可表為：

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

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

式中的ak（t）表示第k個(gè)主成分，vk（x）為對(duì)應(yīng)主成分的響應(yīng)特征矩陣，ak（t）代表時(shí)間特征，vk（x）代表空間響應(yīng)。將B 矩陣的特征值λ 按從大到小的排列，特征值越大代表對(duì)殘余時(shí)間序列方差的貢獻(xiàn)越大[10]?？梢杂锰卣髦道鄯e貢獻(xiàn)率mk來(lái)代表主成分對(duì)矩陣的貢獻(xiàn)率，mk可由下式表示：

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

這個(gè)過(guò)程就稱為Karhunen-Loeve 展開（KLE）[11]。

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

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

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

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

2.1 獲取GNSS 時(shí)間序列殘差

GNSS 臺(tái)站坐標(biāo)每日觀測(cè)時(shí)間序列擬合函數(shù)為[12]：

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

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

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

2.2 異常站點(diǎn)剔除與共模誤差提取

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

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

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

圖3 PCA 解算的前3 個(gè)主成分對(duì)應(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 之間，說(shuō)明云南GNSS 坐標(biāo)時(shí)間序列中存在比較明顯的共模誤差。

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

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

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

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

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

2.3 異常站點(diǎn)共模誤差提取及形變分析

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

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

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

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

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

圖9 YNXP 站去噪后的坐標(biāo)時(shí)間序列Fig.9 Coordinate time series of YNXP station after noise remova

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

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

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

圖11 YNXP 站去噪后的坐標(biāo)時(shí)間序列Fig.11 Coordinate time series of YNXP station afternoise removal

3 結(jié)論

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

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