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

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

0 引 言

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

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

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

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

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

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

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

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

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

2 CME 的提取方法

2.1 SF 方法

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

式 中，σi,j為第i天第j個(gè)測(cè)站的中誤差.

2.2 NIF 方法

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

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

2.3 PCA 方法

PCA 法是一種通過(guò)線性變化，從多個(gè)變量中提取重要變量的多元統(tǒng)計(jì)分析方法.目前，利用PCA 法提取特征信息，已被廣泛應(yīng)用于地殼形變信號(hào)處理、遙感、地理信息系統(tǒng)和測(cè)繪等領(lǐng)域.利用PCA 法，可提取區(qū)域GPS 測(cè)站坐標(biāo)時(shí)序中的CME.設(shè)區(qū)域范圍內(nèi)存在n個(gè)測(cè)站，觀測(cè)歷元為m，構(gòu)成m×n維坐標(biāo)噪聲時(shí)序矩陣A，其中m＞n.矩陣A去趨勢(shì)項(xiàng)和均值化后的坐標(biāo)分量殘差時(shí)間序列矩陣為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è)非零對(duì)角元素構(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個(gè)主成分，表達(dá)式為

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

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

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

3 區(qū)域CME 分析

3.1 比較三種方法提取CME

研究利用SF、NIF 和PCA 方法提取日本房總半島區(qū)域2019—2021年GPS 坐標(biāo)時(shí)序中的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 提取的測(cè)試結(jié)果表明，SF、NIF 和PCA 方法提取CME 的性能基本相當(dāng).

圖3 三種方法提取CME 比較

3.2 GPS 站點(diǎn)空間分辨率對(duì)CME 的影響

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

圖4 不同空間分辨率的GPS 站點(diǎn)

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

3.3 CME 對(duì)慢滑移地表位移的影響

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

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

圖6 去除CME 前后的慢滑移信號(hào)

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

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

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

4 結(jié)束語(yǔ)

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

1) 以日本房總半島區(qū)域2019—2021年GPS 坐標(biāo)時(shí)序?yàn)槔?，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 站點(diǎn)的空間分辨率不同，提取的CME 略有差異.總體上，隨著站點(diǎn)間距的增加，CME的RMS 值也隨之增加，表明隨著GPS 站點(diǎn)空間分辨率的降低，所提取的CME 的離散度增大.

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