GPS坐標(biāo)時序共模誤差提取方法研究

嚴(yán)麗，羅正東，李萌，鄒小平

(1.東華理工大學(xué) 江西省數(shù)字國土重點實驗室,南昌 330013；2.東華理工大學(xué) 測繪工程學(xué)院,南昌 330013；3.江西省防震減災(zāi)與工程地質(zhì)災(zāi)害探測工程研究中心,南昌 330013))

0 引 言

GPS 坐標(biāo)時序中微小地殼形變信號，常被共模誤差(CME)掩蓋[1]，導(dǎo)致地殼構(gòu)造形變信息的信噪比(SNR)降低[2].CME 是區(qū)域GPS 測站坐標(biāo)時序中共同存在的一種或多種誤差集合[3].大區(qū)域內(nèi)衛(wèi)星軌道和天線相位中心(APC)誤差，小區(qū)域內(nèi)陸地水儲量等因素，是產(chǎn)生CME 的主要原因[4-5].區(qū)域內(nèi)，GPS坐標(biāo)時序間的相關(guān)性和距離呈負(fù)相關(guān)，測站間距離越大，相關(guān)性越小，當(dāng)測站間距離大于6 000 km 時，誤差間相關(guān)性趨于0.去除區(qū)域CME，對分析地表位移形變、構(gòu)造運動形變等具有重要的作用.因此，需采用有效方法提取CME，并予以去除.

CME 的提取方法，有堆棧濾波 (SF) 法、網(wǎng)絡(luò)反演濾波 (NIF) 法、主成分分析 (PCA) 法等.1997年，WDOWINSKI 等[6]在GPS 坐標(biāo)時序中發(fā)現(xiàn)CME，并提出通過SF 去除CME，減弱GPS 坐標(biāo)時序的空間相關(guān)性.同年，SEGALL 等[7]利用NIF 反演GPS 坐標(biāo)時序，可同時提取并去除CME.2006年，DONG 等[8]利用PCA 法提取CME，該方法適用于大空間跨度范圍CME 的提取.目前，仍有許多學(xué)者在研究GPS 坐標(biāo)時序及區(qū)域構(gòu)造運動時，進(jìn)行CME 分析[9-11].孫陽等[9]利用PCA 法提取2013—2018年環(huán)渤海區(qū)域27個GPS 測站的CME，研究表明，去除CME 后東(E)、北(N)、天頂(U)方向噪聲明顯減小，有效提高了GPS 坐標(biāo)時序的精度.常金龍等[10]利用PCA 和SF方法，有效提取華北地區(qū)CME，同時反映測站CME 空間分布特征.OZAWA 等[11]利用GPS 坐標(biāo)時序反演漫滑移時，利用PCA 和SF 方法提取CME，分析CME 對滑移反演的影響.

本文研究比較三種CME 處理方法：SF、NIF 和PCA 法.以2019—2021年日本房總半島區(qū)域GPS 坐標(biāo)時序為例：首先，進(jìn)行GPS 坐標(biāo)時序建模，提取殘差噪聲序列；其次，分別利用SF、NIF 和PCA 方法從殘差噪聲序列中提取CME；然后，比較不同GPS 站點空間分辨率條件下三種方法提取CME 的有效性和適用性；最后，研究CME 對日本房總半島2018年慢滑移的影響.

1 GPS 坐標(biāo)時序預(yù)處理

區(qū)域CME 隱含在各GPS 站點坐標(biāo)時序的噪聲殘差序列中，通過構(gòu)建GPS 坐標(biāo)時序模型，提取噪聲殘差時序x(t)[11-12]，方程為

式中：x(ti)為歷元i時的噪聲殘差時序；ti為時間；y(ti)為GPS 坐標(biāo)時序；a為常量項；vti為速度項；biC(ti-tj) 為歷元tj時由天線位置變動或同震等原因引起的階躍項；dln(1+(ti-te)/σ)為地震發(fā)生時刻te的震后弛豫項；τ為震后對數(shù)弛豫時間；Fi、Hi為年周期性運動系數(shù)；Mi、Ni為半年周期性運動系數(shù).

以日本房總半島為例，利用該區(qū)域62 個GPS 測站的數(shù)據(jù)，構(gòu)建上述坐標(biāo)時序模型，圖1顯示了站點的分布，并提取噪聲殘差時序.圖2展示了區(qū)域GPS站點坐標(biāo)時序預(yù)處理過程：首先，需進(jìn)行粗差剔除和階躍修復(fù)，剔除點如圖2(a)中標(biāo)注誤差棒所示，2019—2021年無明顯階躍；然后，進(jìn)行速度項(圖2(b))、周期項(圖2(c))建模，2019—2021年無強震弛豫項；最后，去除已建模的各項位移，獲取噪聲殘差時序，用于后續(xù)CME 的提取.

圖1 日本房總半島區(qū)域GPS 站空間分布

圖2 GPS 坐標(biāo)時序預(yù)處理

2 CME 的提取方法

2.1 SF 方法

假設(shè)CME 在區(qū)域空間上均勻分布，但隨著空間跨度范圍逐漸增大，CME 將越來越小直至消失.SF 法提取的CME 是測站殘差的加權(quán)平均值.因此，SF 法又稱加權(quán)平均法.區(qū)域GPS 網(wǎng)絡(luò)范圍內(nèi)存在n個測站，觀測歷元為m，可構(gòu)建一個m×n維矩陣X，X中元素x(ti,sj)表示第i天第j個測站的坐標(biāo)噪聲時序，CME 計算方程為[13]

式 中，σi,j為第i天第j個測站的中誤差.

2.2 NIF 方法

NIF 法是一種隨時間變化平滑斷層滑動時空分布的大地測量反演方法[7].NIF 不采用任何特定的函數(shù)形式，只要滑動在一定程度上暫時平滑，就可以恢復(fù)滑動的任意時間變化[14].通過NIF 不僅可以提取CME，還可以將微弱信號如慢滑移信號從CME 分離出來[15].利用NIF 反演模型，對GPS 測站坐標(biāo)時序u(x,t)建模為

式中：u(x,t)為GPS 坐標(biāo)時序；s(ξ,t)G(x,ξ)n(ξ)dA(ξ)為斷層面A(ξ) 上 ξ點處的累積滑移s(ξ,t)引發(fā)的地表位移；G(x,ξ)為彈性格林函數(shù)；L(x,t)為隨機基準(zhǔn)擺動；f(t)為CME；ε1(x,t)為隨機噪聲.

2.3 PCA 方法

PCA 法是一種通過線性變化，從多個變量中提取重要變量的多元統(tǒng)計分析方法.目前，利用PCA 法提取特征信息，已被廣泛應(yīng)用于地殼形變信號處理、遙感、地理信息系統(tǒng)和測繪等領(lǐng)域.利用PCA 法，可提取區(qū)域GPS 測站坐標(biāo)時序中的CME.設(shè)區(qū)域范圍內(nèi)存在n個測站，觀測歷元為m，構(gòu)成m×n維坐標(biāo)噪聲時序矩陣A，其中m＞n.矩陣A去趨勢項和均值化后的坐標(biāo)分量殘差時間序列矩陣為X(ti,xj)，其中i=1,2,···m；j=1,2,···n.將矩陣X(ti,xj)協(xié)方差矩陣定義為B，矩陣B中元素bi,j為[8,16]

將n×n維協(xié)方差矩陣B進(jìn)行正交分解

式中：VT為特征向量構(gòu)成的n×n維正交矩陣；Λ為k個非零對角元素構(gòu)成的特征值矩陣.矩陣B為滿秩矩陣(k=n)，則矩陣B的特征值為 λ1,λ2,λ3,···λn，特征向量為v1,v2,v3,···vn.X(ti,xj)可表示為

式中：j=1,2,···,n；ak為矩陣X(ti,xj)的第k個主成分，表達(dá)式為

式中：k=1,2,···,n；vk為第k個主成分ak的響應(yīng)特征矩陣；ak代表時間變化特征；vk代表空間響應(yīng)特征.

主成分分解的特征值越大，則主成分對原始序列的貢獻(xiàn)率越大.通常按貢獻(xiàn)率對特征值進(jìn)行降序排列，越靠前的主成分貢獻(xiàn)率越大，攜帶的信息越多，能夠反映區(qū)域的共同變化特征，反之越靠后的主成分只能反映自身變化特征.利用前p個主成分計算CME，表達(dá)式為

式中：ε(ti,xj)為CME；p為主成分個數(shù)；ak為第k個主成分；vk為ak對應(yīng)特征向量.

3 區(qū)域CME 分析

3.1 比較三種方法提取CME

研究利用SF、NIF 和PCA 方法提取日本房總半島區(qū)域2019—2021年GPS 坐標(biāo)時序中的CME，如圖3所示，三種方法提取的CME 均值近似為0，E、N 與U 方向的標(biāo)準(zhǔn)偏差均近似為1.4 mm、1.2 mm、3.5 mm，水平與高程方向CME 差異較大，水平方向振幅在5 mm 以內(nèi)，垂直方向振幅在15 mm 以內(nèi).本研究區(qū)域CME 提取的測試結(jié)果表明，SF、NIF 和PCA 方法提取CME 的性能基本相當(dāng).

圖3 三種方法提取CME 比較

3.2 GPS 站點空間分辨率對CME 的影響

GPS 站點在不同區(qū)域的空間分辨率差異較大，如我國西部大多數(shù)活動斷裂附近，GPS 站點的間距在80～150 km；日本陸地GPS 觀測站平均間距約為10 km；美國南加州采用重點地區(qū)密集觀測，GPS 站點間距可達(dá)2～5 km[17-18].謝樹明等[19]利用相關(guān)系數(shù)疊加濾波，計算較大空間(1 000 km 范圍內(nèi)9 個GPS站點與全球138 個站)的CME，結(jié)果表明不同空間尺度的GPS 站點分布得到的CME 存在差異.為了進(jìn)一步分析區(qū)域不同空間分辨率對提取CME 的影響，對研究區(qū)域進(jìn)行采樣，如圖1所示，GPS 站點間距約為10 km.采樣后的GPS 站點，如圖4所示，站間距分別約為20 km、40 km、80 km、120 km.在GPS 不同空間分辨率條件下，利用上述SF、NIF 和PCA 方法，提取日本房總半島區(qū)域2019—2021年GPS 坐標(biāo)時序中的CME，圖5為其均方根(RMS)值.圖5中，總體上隨著站點間距增加，CME 的RMS 值也隨之增加，表明隨著GPS 站點空間分辨率降低，所提取CME 的離散度將增大.

圖4 不同空間分辨率的GPS 站點

圖5 不同GPS 站點空間分辨率條件下提取CME 的RMS 值

3.3 CME 對慢滑移地表位移的影響

去除CME 一般是為了更好地分析GPS 坐標(biāo)時序中微弱的地殼形變信號[8].在GPS 坐標(biāo)時序中，微弱的慢滑移信號有時會隱藏在區(qū)CME 內(nèi)[15].NIF 法可從CME 中分離較小的慢滑移信號，在一定程度上提高慢滑移地表位移的準(zhǔn)確性[7,20].

以日本房總半島2018年慢滑移事件為例，研究CME 對慢滑移地表位移的影響.在CME 去除前后，GPS 站點J226 在E、N、U 方向的坐標(biāo)時序如圖6所示.圖6中，從2018～2019年，去除CME 后，坐標(biāo)時序的噪聲明顯減小，更加平滑.

圖6 去除CME 前后的慢滑移信號

房總半島2018年慢滑移事件的滑移中心區(qū)域位于東海岸附近(34.8°N～35.6°N，140.0°E～141.0°E)[11].統(tǒng)計研究區(qū)域內(nèi)所有站點去除基準(zhǔn)擺動和隨機噪聲后慢滑移的地表水平和高程位移，分別如圖7～8所示.其中，靠近滑移中心區(qū)域東海岸附近GPS 站點的水平和高程位移量明顯較大.圖7中，CME 去除前后，最大水平位移分別為3.5 cm 和3.8 cm，兩者相差0.3 cm；圖8中，CME 去除前后，最大高程位移為3.1 cm 和3.0 cm.高程方向，CME 的離散度較大，但其對慢滑移大小的影響并不十分明顯.未產(chǎn)生慢滑移的站點(緯度大于35.5°)，其位移主要為噪聲殘差，去除CME 后，噪聲明顯減??；而產(chǎn)生慢滑移的站點，其位移主要為慢滑移產(chǎn)生的地表運動，去除CME 后，水平位移略向東偏移.因此，去除CME，可減小因板塊運動、參考框架、大氣等產(chǎn)生的區(qū)域噪聲影響，提高慢滑移地表形變信號的SNR.

圖7 去除CME 前后的水平位移

圖8 去除CME 前后的高程位移

4 結(jié)束語

通過構(gòu)建GPS 坐標(biāo)時序模型，提取噪聲殘差時序，利用SF、NIF 和PCA 方法提取隱含噪聲殘差時序中的CME.比較三種方法提取CME 的性能，分析GPS 站點不同空間分辨率對CME 提取的影響，及共模誤差對慢滑移信號的影響：

1) 以日本房總半島區(qū)域2019—2021年GPS 坐標(biāo)時序為例，SF、NIF、PCA 三種方法提取CME 的性能基本相當(dāng).提取的CME 在不同方向上差異較大，E、N、U 方向的標(biāo)準(zhǔn)偏差分別為1.4 mm、1.2 mm、3.5 mm.水平方向振幅在5 mm 以內(nèi)，垂直方向振幅均在15 mm以內(nèi).

2) 區(qū)域GPS 站點的空間分辨率不同，提取的CME 略有差異.總體上，隨著站點間距的增加，CME的RMS 值也隨之增加，表明隨著GPS 站點空間分辨率的降低，所提取的CME 的離散度增大.

3) CME 會影響慢滑移信號的大小和方向，特別是對水平方向位移影響較為顯著.去除CME 可提高慢滑移地表形變信號的SNR.