利用GPS和水文負載模型研究云南地區(qū)垂向季節(jié)性波動變化和構造變形

胡順強, 王坦, 管雅慧, 楊振宇*

1 首都師范大學資源環(huán)境與旅游學院， 北京 100048 2 中國地震臺網(wǎng)中心， 北京 100045

0 引言

GPS技術在地殼運動監(jiān)測中已經(jīng)得到了成功的應用(Wang et al., 2001;Gan et al., 2007; Zheng et al., 2017; Wang and Shen, 2020).可靠的GPS垂向速度場主要是通過GPS連續(xù)站觀測數(shù)據(jù)解算得到，由于目前GPS連續(xù)站較少，且解算得到的垂向位移誤差大致是水平方向位移誤差的2～3倍以上(Liang et al., 2013)，因此，GPS垂向位移在地殼形變中應用較少.GPS垂向位移中除了真實的地殼構造變形之外，還受以周年和半周年的季節(jié)性波動變化的影響(Dong et al., 1998; Nikolaidis, 2002; 田云鋒和沈正康, 2009).大氣、水文及非潮汐海洋等地表質量負載是引起GPS垂向位移中季節(jié)性變化的主要原因(Van Dam et al., 1994, 1997; Zerbini et al., 2004; 王敏等, 2005; Wu et al., 2019).在一些陸地水變化較大的區(qū)域，水文負載引起的季節(jié)性變化能達到厘米級(王林松等, 2014).因此，如何去除GPS垂向位移中因大氣、水文及非潮汐海洋負載引起的垂向季節(jié)性變化，最大限度地獲得真實的地殼構造運動引起的變形量，成為地殼動力學研究的熱點(Pan et al., 2018;Zhan et al., 2020; Li et al., 2020).

分析大氣、水文和非潮汐海洋等地表質量負載引起的垂向季節(jié)性變化方法主要有地表負載模型和GRACE模型等.對大氣和非潮汐海洋負載的研究，非潮汐海洋負載對沿海地區(qū)的GPS垂向位移改正效果更好(Van Dam et al., 2012).Petrov 和Boy(2004)，Tregoning和Van Dam(2005)分別研究了大氣負載引起的垂向季節(jié)性變化可達20 mm和18 mm.對水文負載的研究，由地表負載模型得到的水文負載形變對GPS垂向季節(jié)性變化的影響為毫米到厘米不等，且改正GPS垂向位移后的均方根值都會減少(Van Dam et al., 2001; Dill and Dobslaw, 2013; Xiang et al., 2018).多位學者使用GRACE模型數(shù)據(jù)所得的水文負載形變和GPS垂向季節(jié)性變化進行對比研究，大部分結果表明兩者在總體上具有較好的相關性和一致性(Davis et al., 2004;Nahmani et al., 2012;Fu and Freymueller, 2012;Pan et al., 2016;丁一航等, 2018; Li et al., 2019a).Dong等(2002)的研究結果表明大氣、水文和非潮汐海洋負載引起的GPS連續(xù)站最大周年振幅分別可達4 mm、7～8 mm和2～3 mm，可以解釋約為40%的GPS垂向季節(jié)性變化.綜合上述研究可知，不同地區(qū)不同地表質量負載(大氣、水文、非潮汐海洋)對GPS垂向季節(jié)性變化影響均不一樣.

云南地區(qū)(21°N—29°N，97°E—106°E)位于青藏高原東南側，是多塊體組成的重要地質研究區(qū)，內(nèi)部具有瀾滄江、小江和紅河等多條大型且復雜的斷裂帶，是中國大陸現(xiàn)今構造活動最為劇烈的區(qū)域之一(Gao et al., 2017).部分學者分別使用時間跨度在2010—2012年(盛傳貞等(2014)) , 2010—2015年(Hao等(2016), Zhan等(2017))的GRACE模型和GPS數(shù)據(jù)研究云南地區(qū)的GPS垂向位移與水文負載形變情況，認為水文負載是引起云南地區(qū)GPS垂向運動季節(jié)性變化的主要因素之一.然而，GRACE只能有效分辨大約400 km范圍內(nèi)水文負載變化，對于GPS連續(xù)站局部小尺度范圍的水文負載影響不能有效辨別，因此，本文使用時間跨度更長(2011—2020年)且空間分辨率為0.5°×0.5°的格網(wǎng)LSDM模型和GPS數(shù)據(jù)來探討云南地區(qū)垂向運動的季節(jié)性變化和現(xiàn)今構造變形.前人分析水文負載對云南地區(qū)GPS垂向運動的季節(jié)性變化影響時，在季節(jié)性信號提取方面關注比較少.因此，本文使用奇異譜分析(Singular Spectrum Analysis，SSA)方法提取GPS垂向位移和LSDM形變的季節(jié)性信號.由于GPS垂向位移和LSDM形變中含有其他非周年信號且相位差可能隨時間變化，傳統(tǒng)的皮爾遜相關系數(shù)只能簡單的從單一時間尺度上衡量兩者的相關性，而忽略了兩者在多時間尺度上的相關性.小波分析可以分析不同時間序列在不同時段和頻率尺度上的相關性，能夠揭示兩者在時頻空間中的相位關系(Xu, 2016).因此，本文使用連續(xù)小波變換、交叉小波變換(Cross Wavelet Transform, XWT)和小波相干性在時頻空間下多時間尺度研究GPS垂向位移與LSDM形變的周期特性.

1 數(shù)據(jù)與方法

1.1 GPS數(shù)據(jù)

本文選取的云南地區(qū)27個GPS連續(xù)站觀測數(shù)軟件解算得到，首先利用GAMIT進行單日解解算，估計站位置、衛(wèi)星軌道參數(shù)及方差-協(xié)方差矩陣等；然后使用GLOBK進行網(wǎng)平差，從而獲得連續(xù)站的單日解坐標時間序列，詳細的解算策略參考文獻(Zhao et al., 2015).在對GPS垂向位移進行分析之前需要進行粗差剔除、階躍改正、缺失數(shù)據(jù)插值等預處理.預處理主要方法如下：

1)采用四分位距法(Inter Quartile Range, IQR)對GPS垂向位移進行粗差探測和剔除，IQR判別準則原理如下：

IQR=Q2-Q1.

(1)

異常值探測區(qū)間為

[Q1-1.5·IQR,Q2+1.5·IQR],

(2)

式(2)中，Q1和Q2分別表示為最靠近1/4和3/4處的下分位值和上分位值.

圖1 云南地區(qū)27個GPS連續(xù)站分布和數(shù)據(jù)缺失圖Fig.1 Location and Data availability of 27 GPS stations in Yunnan region

2)使用最小二乘線性擬合方法對GPS垂向位移的階躍項進行改正.

3) 使用Schneider (2001)提出的正則期望最大化(Regularized EM, RegEM)方法對云南地區(qū)27個連續(xù)站的GPS垂向位移中缺失的數(shù)據(jù)進行插值，該方法顧及了27個GPS連續(xù)站之間的相關性和物理背景，不需要先驗信息和依賴數(shù)據(jù)模型，只根據(jù)數(shù)據(jù)自身特性進行插值，在GPS坐標時間序列插值中已經(jīng)得到了廣泛應用(明鋒, 2019; Li et al., 2019b).其主要原理如下:

由p個連續(xù)站和m個觀測天數(shù)的GPS垂向位移組成m×p維矩陣X，對于矩陣X中的任一個連續(xù)站的GPS垂向位移，非缺失與缺失的GPS垂向位移可使用線性回歸模型來描述：

xm=um+(xa-ua)B+e,

(3)

式(3)中，矩陣B∈Rpa×pm為回歸系數(shù)，殘差e的均值為0，協(xié)方差陣C∈Rpa×pm為未知的隨機變量，非缺失與缺失的GPS垂向位移組成的向量分別為xa和xm，均值分別為ua和um.給定的均值u和協(xié)方差陣通過條件最大似然估計法計算X中每一行包含數(shù)據(jù)缺失的GPS垂向位移的回歸系數(shù)B和殘差協(xié)方差陣C，然后再使用式(4)對缺失的GPS垂向位移進行插值.

(4)

圖2為KMIN、XIAG、YNMH和YNTC連續(xù)站經(jīng)過預處理后的GPS垂向位移，從圖2中可看出使用RegEM方法的插值效果與整體運動趨勢一致.

圖2 預處理后的KMIN (a)、XIAG (b)、YNMH (c)和YNTC (d)連續(xù)站的GPS垂向位移Fig.2 The GPS vertical displacement of KMIN (a)，XIAG (b)，YNMH (c) and YNTC (d) stations after pretreatment

1.2 地表負載模型數(shù)據(jù)

地表負載模型使用德國波斯坦地學中心提供的由于大氣負載、水文負載和非潮汐海洋負載形變的規(guī)則全球格網(wǎng)數(shù)據(jù)(Dill and Dobslaw, 2013)，空間分辨率為0.5°×0.5°.其中，水文負載形變由時間分辨率24 h的LSDM模型數(shù)據(jù)(Dill, 2008)計算得到，大氣和非潮汐海洋負載形變分別由時間分辨率3 h的ECMWF模型和MPIOM模型數(shù)據(jù)計算得到(Marsland et al., 2003).各個GPS連續(xù)站相對應的大氣負載，非潮汐海洋負載和水文負載形變通過雙線性法對模型數(shù)據(jù)進行插值得到.為了和GPS垂向位移的日分辨率統(tǒng)一，分別將大氣和非潮汐海洋負載形變的3 h時間分辨率統(tǒng)一為每日時間分辨率.圖3為經(jīng)過處理后KMIN連續(xù)站對應的水文負載、大氣負載、非潮汐海洋負載形變.

圖3 KMIN連續(xù)站對應的水文(a)、大氣(b) 、非潮汐海洋負載(c)形變Fig.3 hydrological (a), Atmospheric (b), non-tidal ocean (c) load deformation of KMIN station

1.3 奇異譜分析

提取GPS垂向位移中季節(jié)性信號的方法主要有小波分解法(Bogusz and Klos, 2016)、卡爾曼濾波法(Davis et al., 2012)、基于最小二乘擬合的諧波模型法(Blewitt and Lavallée, 2002)、SSA方法(Chen et al., 2013)等，但這些方法各具優(yōu)缺點，如小波分解法中需要根據(jù)經(jīng)驗提前定義小波基函數(shù)和分解層數(shù)，不同的小波基函數(shù)和分解層數(shù)對季節(jié)性信號提取影響較大；卡爾曼濾波法的隨機過程需要提前估計，程序運行耗時較長；基于最小二乘擬合的諧波模型法最為常用，然而，該方法提取出的周年、半周年項振幅和相位均為常數(shù),不符合GPS垂向位移中實際的季節(jié)性變化.SSA方法是一種針對時間序列中信號分組與重構的方法，它的優(yōu)勢在于不需要任何外在的假設條件和先驗信息，不僅能夠提取時間序列數(shù)據(jù)中的非線性趨勢信號，且不受正弦波假定的約束，能夠穩(wěn)定可靠的識別與強化周期信號等，因此，特別適合分析和提取有周期振蕩的時間序列數(shù)據(jù).

將時間跨度為2011年1月—2020年9月云南地區(qū)27個連續(xù)站的GPS垂向位移作為原始時間序列{y1,y2,y3,…,yn}，設定窗口長度M(M

(5)

1)計算滯后矩陣X的自協(xié)方差矩陣C及特征值λ1≥λ2≥…λM≥0和特征向量Ekj，如下式(6):

(6)

2)計算滯后矩陣X在特征向量Ekj上的投影，即為時間主成分，如式(7)所示:

(7)

3)時間序列的重建(Reconstruction Component, RC).由第k個時間主成分和特征向量按照式(8)進行重建(Vautard et al., 1992).

(8)

(9)

1.4 小波分析

連續(xù)小波變換將小波基函數(shù)作為帶通濾波器運用于分析具有時間序列的數(shù)據(jù)(Grinsted et al.,2004)，可以很好地揭示單一時間序列數(shù)據(jù)的多時間-尺度變換特征的振蕩周期.對于某個連續(xù)站的GPS垂向位移xn(n=1,2,…,N)，將其連續(xù)小波變換定義為

(10)

WXY=WXWY*,

(11)

式(11)中，*表示復共軛，小波功率受邊緣效應影響較大的區(qū)域為影響堆(Cone of Influence, COI).為了準確描述xn,yn的相位關系，選擇COI區(qū)域之外一系列相位角的圓域平均值來衡量：

(12)

交叉小波的相似性定義為

ρi=cos(am),

(13)

式(13)中,ρ取值范圍為-1到1，0代表不相關.

交叉小波變換能夠很好地揭示GPS垂向位移和LSDM形變共同的高能量區(qū)及時頻域中的相位關系.小波相干性則能夠很好地彌補交叉小波變換在低能量區(qū)的不足，用來度量時頻域中兩者的局部相關密切程度，即對應交叉小波功率譜中低能量值區(qū).連續(xù)站的GPS垂向位移和LSDM形變xn,yn的小波相干譜定義為

(14)

式(14)中，S是平滑窗口，s是小波尺度，Rn(s)是局部相干系數(shù).

2 結果

2.1 GPS垂向位移與水文負載形變比較

為了更好的對比分析GPS垂向位移與LSDM形變的季節(jié)性變化關系，從GPS垂向位移中去除由1.2節(jié)中介紹的方法得到的大氣和非潮汐海洋負載形變,然后使用最小二乘法(Least Squares Fitting, LSF)對27個連續(xù)站的GPS垂向位移進行去線性趨勢處理.27個連續(xù)站的GPS垂向位移和LSDM形變均出現(xiàn)季節(jié)性變化且整體運動趨勢較為一致，大部分連續(xù)站的GPS垂向位移變化范圍在-30 mm至30 mm，相對應的LSDM形變位移值范圍在-10 mm至10 mm，說明水文負載形變并不能完全解釋GPS垂向位移的季節(jié)性變化，水文負載形變是引起云南地區(qū)GPS垂向季節(jié)性變化的主要影響因素之一(Tesmer et al., 2011).圖4為KMIN、XIAG、YNMH和YNTC連續(xù)站的GPS垂向位移和LSDN形變的比較.圖4中YNMH、YNTC和XIAG連續(xù)站的相位和振幅具有較好的一致性，YNMH、YNTC和XIAG連續(xù)站位于云南地區(qū)陸地水變化較大的區(qū)域，說明在陸地水變化較大的區(qū)域，水文負載形變能夠很好地解釋GPS垂向位移的季節(jié)性變化.KMIN連續(xù)站在2011—2014年的相位吻合度較差，說明該時段水文負載形變不能很好地解釋GPS垂向位移的季節(jié)性變化，可能受站點局部環(huán)境、系統(tǒng)誤差及噪聲影響較大.如果GPS垂向位移和LSDM形變的周期性具有一定的相關性，則有兩種情況：1)GPS垂向位移和LSDM形變的周期項為同一周期，他們都是由同一個物理因素造成的，只是因他們的分辨率或者獲取的方法不同，從而顯現(xiàn)出兩個周期；2)如果GPS垂向位移和LSDM形變周期項的物理因素之間具有相關性.那么，可以使用皮爾遜相關系數(shù)(如式(15))來判斷GPS垂向位移與LSDM形變是否為同一周期項或者由同一物理因素產(chǎn)生的周期振動.27個連續(xù)站的GPS垂向位移與LSDM形變的相關系數(shù)計算結果如表1所示，兩者的相關性平均為0.59，相關性最小的為YNGM連續(xù)站，值為0.15；相關性最大的為YNLJ連續(xù)站，值為0.76，說明大部分連續(xù)站的GPS垂向位移與LSDM形變具有較強的相關性.為了進一步說明GPS垂向位移與LSDM形變季節(jié)性變化的一致性，本文使用從GPS垂向位移中扣除LSDM形變后的RMS減少量(如式(16))來定量評估LSDM形變能否有效改正GPS垂向位移的非構造形變，若RMS減少量值為正值，說明LSDM形變能夠有效的改正GPS垂向位移中的非構造形變，計算的RMS減少量結果如表1所示.從表1可知，除了YNGM連續(xù)站的RMS減少量值為負值之外，其他連續(xù)站的值都為正值，說明LSDM形變都能夠有效去除這些連續(xù)站中的非構造形變；所有連續(xù)站的RMS減少量平均值為17.13%，說明GPS垂向位移中所包含的非構造形變，平均約17.13%來源于LSDM形變.

圖4 KMIN (a)、XIAG (b)、YNMH (c)和YNTC (d)連續(xù)站的GPS垂向位移與相對應的LSDM形變Fig.4 GPS vertical displacement and corresponding LSDM deformation for KMIN (a)，XIAG (b)，YNMH (c) and YNTC (d)

表1 GPS與LSDM相關系數(shù)和RMS減少量Table 1 Correlation coefficient and RMS reduction between GPS and LSDM

(15)

(16)

2.2 GPS共模誤差與水文負載形變比較

Wdowinsk 等(1997)研究表明，區(qū)域GPS網(wǎng)中不同連續(xù)站的GPS垂向位移中會存在時空相關的非構造形變噪聲，稱為共模誤差，這類非構造形變噪聲來源不明確，有可能來源于觀測環(huán)境的好壞，地表質量負載和GPS數(shù)據(jù)解算的系統(tǒng)誤差等.通常情況下在GPS垂向位移中去除速度項和周期項后的殘差時間序列中進行共模誤差提取，但是周期項信號中可能包含部分共模誤差，因此，對去除速度項之后的GPS垂向位移殘差時間序列進行共模誤差提取.對云南地區(qū)的27個GPS連續(xù)站使用主成分分析計算出該地區(qū)的共模誤差.圖5為KMIN、XIAG、YNMH和YNTC連續(xù)站的GPS共模誤差和LSDM形變的對比，從圖5中可知，GPS共模誤差和LSDM形變的整體運動趨勢更為一致，使用式(15)定量計算兩者的相關性，結果如表2所示，所有連續(xù)站的GPS共模誤差與LSDM形變的平均相關系數(shù)為0.75.若從GPS共模誤差中去除LSDM形變后，所有連續(xù)站的RMS減少量平均為35.72%，所以在云南地區(qū)的GPS共模誤差與LSDM形變具有很好的相關性和一致性.因此，水文負載是該地區(qū)GPS共模誤差的主要成份，這一結果與盛傳貞等(2014)的結果一致.

圖5 KMIN (a)、XIAG (b)、YNMH (c)和YNTC (d)連續(xù)站的GPS共模誤差與LSDM形變對比Fig.5 Comparison of GPS common mode error and LSDM deformation for KMIN (a)，XIAG (b)，YNMH (c) and YNTC (d)

表2 GPS共模誤差與LSDM形變的相關系數(shù)和RMS減少量Table 2 Correlation coefficient and RMS reduction between GPS common mode error and LSDM deformation

2.3 小波譜分析

為了進一步分析GPS垂向位移與LSDM形變的關系，分別使用連續(xù)小波變換、交叉小波變換和小波相干性對兩者進行周期分析,選取Morlet小波作為母函數(shù).由于連續(xù)站數(shù)目較多，本文以YNTC連續(xù)站作為實例進行詳細介紹.圖7為YNTC連續(xù)站的GPS垂向位移與LSDM形變的連續(xù)小波、交叉小波和小波相干性的功率譜.圖6a和6b的連續(xù)小波功率譜中黃色與藍色分別表示能量密度的峰值和谷值，從藍色到黃色，表示小波能量譜依次遞增.黑色粗實線封閉區(qū)域表示通過了95%置信水平的顯著性檢驗，黑色細實線下方為COI區(qū)域.從圖6a中YNTC連續(xù)站的連續(xù)小波功率譜可知，GPS垂向位移在全時域范圍內(nèi)存在256～512天(0.7～1.4a)的主振蕩周期且通過了顯著性檢驗，在2013—2016和2018—2020時間范圍內(nèi)存在128～256天(0.4～0.7a)的主振蕩周期且通過了顯著性檢驗，圖6a中小于128天的GPS垂向位移中高頻部分沒有顯著的振蕩周期.從圖6b中LSDM形變的連續(xù)小波功率譜可知，LSDM形變在全時域范圍內(nèi)存在128～256天(0.4～0.7a)，256～512天(0.7～1.4a)的主振蕩周期且通過了顯著性檢驗.由圖6a和6b的結果表明，YNTC連續(xù)站的GPS垂向位移與LSDM形變在全時間域內(nèi)均具有顯著的近年周期(256～512天)變化.與此同時，使用連續(xù)小波變換對其他連續(xù)站進行周期分析，其他連續(xù)站都呈現(xiàn)出相似的結果.

連續(xù)小波變換只是反映了GPS垂向位移與LSDM形變自身的時間-尺度變換特征，為了進一步分析GPS垂向位移與LSDM形變的相互周期特性，使用交叉小波變換和小波相干性來反映GPS垂向位移與LSDM形變在高能量區(qū)和低能量區(qū)及其相位關系.如圖6c所示，黑色箭頭反映了參與交叉小波變換分析的GPS垂向位移與LSDM形變的相位關系: 若規(guī)定箭頭方“→”表示二者相位相同，則“←”是相位相反;若箭頭“↑”表示GPS提前LSDM 1/4 周期，則“↓”為LSDM提前GPS 1/4 周期.

從圖6c中交叉小波功率譜可知，YNTC連續(xù)站的GPS垂向位移與LSDM形變在全時域內(nèi)存在128～256天(0.4～0.7a)，256～512天(0.7～1.4a)的共振周期，且通過了顯著性檢驗.從圖6d中小波相干功率譜低能量區(qū)，GPS垂向位移與LSDM形變存在部分時域內(nèi)小于128天的共振周期，說明LSDM形變在高頻部分不是造成GPS垂向位移主要的驅動力，可能是由GPS中系統(tǒng)誤差引起的(Ray et al., 2008).其他連續(xù)站的GPS垂向位移與LSDM形變在全時域內(nèi)均存在256～512天(0.7～1.4a)的共振周期，在部分時域內(nèi)存在128～256天(0.4～0.7a)的共振周期.因此，本文主要關注GPS垂向位移和LSDM形變共同存在的256～512天的近年周期變化，從圖6c可知，在256～512天的共振周期內(nèi)的相位比較穩(wěn)定，所以本文使用年周期的平均相位來表示GPS垂向位移與LSDM形變在256～512天的平均相位關系.表3為所有連續(xù)站在年周期變化的GPS垂向位移與LSDM形變的平均相位關系.除此之外，交叉小波的平均相位相似性更能反映GPS垂向位移與LSDM形變在不同時域上的相關性，所以計算GPS垂向位移與LSDM形變在年周期變化的平均相位相似性，結果如圖7所示：YNWS、KMIN、YNMZ、YNHZ、YNML、YNTH、YNDC和YNJP連續(xù)站年周期變化的平均相位相似性大小低于0.9.這些平均相位相似性低于0.9的連續(xù)站位于云南地區(qū)東部小江斷裂帶附近，年降雨量較少.云南地區(qū)降雨量在時間上具有顯著的年周期性特征，季節(jié)性降雨是導致該地區(qū)地殼垂向周期性變化的主要因素.云南地區(qū)的降雨量主要受孟加拉灣和南海的暖濕氣流影響，形成了從南到北逐漸減少，東部小于西部的特點(趙寧坤等，2009; Zhan et al., 2017).因此，在降雨量較少，陸地水變化較小的云南東部地區(qū)的GPS連續(xù)站，水文負載引起的形變不能全部解釋GPS年周期變化，其他因素(如其他地球物理因素、LSDM模型誤差等)和水文負載形變共同作用引起了這幾個連續(xù)站的年周期項變化，其他大部分連續(xù)站的年周期平均相位相似性大小在0.9～1之間，說明大部分連續(xù)站的GPS垂向位移與LSDM形變的年周期變化是物理相關的，水文負載形變是引起GPS年周期變化的主要原因.

圖7 27個連續(xù)站的GPS垂向位移與LSDM的交叉小波平均相似性Fig.7 The mean semblance of cross wavelet between 27 GPS vertical displacement and corresponding LSDM deformation

表3 27個連續(xù)站的GPS垂向位移與LSDM形變的平均相位差Table 3 The average phase angle between 27 GPS vertical displacement and corresponding LSDM deformation

為了驗證交叉小波變換方法的可靠性，圖8展示了云南地區(qū)27個連續(xù)站通過XWT與LSF兩種方法計算GPS垂向位移與LSDM形變年周期變化的相位差值結果：差異大部分在±6°以內(nèi)；大部分連續(xù)站相吻合，這表明交叉小波變換可以在時頻空間下有效檢驗GPS垂向位移與LSDM形變的相位關系.

圖8 XWT與LSF兩種方法獲取27個GPS連續(xù)站與LSDM的年周期相位差對比Fig.8 The comparison of annual period phase difference between 27 GPS vertical displacement and corresponding LSDM deformation for XWT and LSF method

2.4 季節(jié)性信號提取

通過小波分析結果可知，不同連續(xù)站的GPS垂向位移與LSDM形變在不同時域內(nèi)存在不同周期的季節(jié)性變化，因此，有必要使用一種高精度的季節(jié)性信號提取方法.在使用SSA方法對GPS垂向位移與LSDM形變進行季節(jié)性信號提取時，窗口大小和截取的主分量個數(shù)對于能否精確提取季節(jié)性信號顯得尤為重要.多位學者使用SSA方法提取季節(jié)性信號時，通過多次實驗得到窗口大小取兩年最為合適(Chen et al., 2013; Li et al., 2017; Xiang et al., 2018),所以本文設定的窗口大小M=730，截取主分量個數(shù)則依據(jù)ω-correlation方法和統(tǒng)計重建后RC成分的方差貢獻率大小來進行綜合判斷.

本文同樣以YNTC連續(xù)站作為實例進行分析，使用ω-correlation方法分別對YNTC連續(xù)站的GPS垂向位移與LSDM形變的前20階RC成分進行分析，從圖9可知，GPS垂向位移從第5個RC成分之后，LSDM形變從第10個RC成分之后， RC成分之間的周期性分離較差，說明它們中噪聲占主要部分.圖9a中，GPS垂向位移中RC1-RC2、RC3-RC4的ω-correlation系數(shù)分別為0.98，0.99；圖9b中，LSDM形變中RC1-RC2、RC3-RC4、RC6-RC7和RC8-RC9的ω-correlation系數(shù)分別為0.99、0.99、0.83、0.85,ω-correlation系數(shù)都大于0.6，可以看作是同一周期成分進行合并.同時分別統(tǒng)計前20階RC成分的方差貢獻率，從圖10a可知，RC1和RC2的方差貢獻率最大,分別為45%和32.4%；RC3和RC4分別為5.7%和4.8%，其他RC成分的貢獻率較低,前4階的方差貢獻率總和為88%.從圖10b可知，LSDM形變中RC1和RC2的方差貢獻率最大，分別為51.4%和38%；RC3、RC4、RC5的貢獻率分別為3.3%，3%，2.4%,其他RC成分方差貢獻率較低，前9階的方差貢獻率總和為99%.綜合方差貢獻率和ω-correlation周期性分析，YNTC連續(xù)站的GPS垂向位移與LSDM形變的RC主分量截斷數(shù)分別為4和9，所以利用SSA方法提取YNTC連續(xù)站的GPS季節(jié)性信號為RC1-RC4成分之和，LSDM季節(jié)性信號為RC1-RC9成分之和.為了對比SSA方法提取GPS垂向位移與LSDM形變的季節(jié)性信號效果，使用LSF方法分別對YNTC連續(xù)站的GPS垂向位移與LSDM形變進行季節(jié)性信號提取，提取結果如圖11所示， SSA與LSF方法提取的季節(jié)性信號與原始數(shù)據(jù)整體運動趨勢一致；相比于LSF方法，SSA方法提取的GPS季節(jié)性信號能夠更加反映原始數(shù)據(jù)的季節(jié)性變化趨勢.同理分別使用SSA和LSF方法對其他26個GPS連續(xù)站的GPS垂向位移與LSDM形變進行季節(jié)性信號提取.為了定量評估SSA和LSF方法提取季節(jié)性信號結果，分別計算SSA和LSF方法提取的季節(jié)性信號和GPS垂向位移和LSDM形變的相關系數(shù)和RMS減少量.計算結果如表4所示，SSA和LSF方法提取的季節(jié)性信號與GPS垂向位移的相關系數(shù)平均為0.73和0.64；RMS減少量平均為32.49%和24.05%；SSA和LSF方法提取的季節(jié)性信號與LSDM形變的相關系數(shù)平均為0.98和0.95，RMS減少量平均為82.77%和69.47%.通過對比可知，SSA方法提取的季節(jié)性信號要整體優(yōu)于LSF方法.

圖9 YNTC連續(xù)站的前20階GPS(a)與LSDM(b)時間序列的w-correlation分析結果Fig.9 w-correlation analysis results of the first 20 order between GPS vertical displacement (a) and corresponding LSDM deformation (b) for YNTC station

圖10 YNTC連續(xù)站的GPS(a)與LSDM(b)時間序列的前20階的RC方差貢獻率統(tǒng)計Fig.10 The first 20 order RC variance contribution rate statistics between GPS vertical displacement (a) and corresponding LSDM deformation (b) for YNTC station

圖11 SSA與LSF方法提取YNTC連續(xù)站的GPS(a)與LSDM(b)的季節(jié)性信號Fig.11 Seasonal signals of GPS vertical displacement (a) and corresponding LSDM deformation (b) extracted by SSA and LSF for YNTC station

表4 SSA與LSF方法提取的GPS與LSDM季節(jié)性信號的相關系數(shù)和RMS減少量Table 4 Correlation coefficient and RMS reduction between GPS and LSDM seasonal signals extracted by SSA and LSF methods

2.5 云南地區(qū)垂向地殼形變

通過表1可知，除YNGM連續(xù)站的RMS減少量為負值之外，其他連續(xù)站的RMS減少量均為正值.因此，為了獲得高精度的云南地區(qū)27個GPS垂向速度場，除了YNGM連續(xù)站之外，其他GPS連續(xù)站都使用水文負載模型數(shù)據(jù)進行改正處理，然后對所有連續(xù)站再使用主成分分析方法(Dong et al., 2006)進行共模誤差剔除，國內(nèi)外學者普遍認為白噪聲和閃爍噪聲組合的噪聲模型能夠代表大部分GPS垂向位移的噪聲特性(Williams et al., 2004;Didova et al., 2016;He et al., 2019)，因此，在經(jīng)過水文負載模型和共模誤差預處理后，最后使用Hector軟件(Bos et al., 2013)通過白噪聲和閃爍噪聲組合模型估計所有連續(xù)站的速率及不確定度.圖12為最后得到的2011—2020年云南地區(qū)的垂向速度場，從圖12中可知，云南地區(qū)以紅河斷裂帶為界劃分滇西南塊體，以紅河斷裂帶、小江斷裂帶和麗江-小金河斷裂組成川滇塊體南部.滇西南塊體中除YNMJ連續(xù)站以0.89 mm·a-1速率抬升之外，其他連續(xù)站以0.01～1.9 mm·a-1的速率沉降；而川滇塊體南部除YNYS和YNYM連續(xù)站速率分別以0.1 mm·a-1和0.05 mm·a-1沉降之外，其他連續(xù)站以0.13～2 mm·a-1速率抬升.

圖12 云南地區(qū)2011—2020年的GPS垂向速度場(紅色和藍色箭頭分別代表抬升和沉降速率)Fig.12 GPS vertical velocity field in Yunnan region during 2011—2020(Red and blue respectively represent the up and down vertical rates)

川滇塊體南部的整體抬升與小江斷裂的左旋剪切運動密切相關，由空間大地測量的研究結果表明小江斷裂帶的運動速率為7～13 mm·a-1(Shen et al., 2005; Zheng et al.,2017; 李長軍等, 2019)，導致川滇塊體整體南向運動，并受阻于南部的紅河斷裂帶，使得川滇塊體南部發(fā)生抬升.而滇西南塊體沿著高黎貢右旋走滑斷裂和南汀河左旋斷裂與勐興左旋斷裂大規(guī)模旋轉運動(Tong et al., 2021)，在滇西南塊體的中東部形成拉張而沉降.

3 討論

3.1 云南地區(qū)垂向速度場差異分析

Hao等(2016)使用GRACE模型數(shù)據(jù)去除GPS垂向位移中非構造形變之后得到的云南地區(qū)2010—2015年GPS垂向速度場整體以0.1～3.4 mm·a-1的速率抬升為主(圖13a).但是，本文得到的2011—2020年GPS垂向速度場結果與之(Hao 等(2016))相差較大，卻與Zhan 等(2017)同樣使用GRACE模型數(shù)據(jù)改正非構造形變之后得到的2010—2015年速度場(圖13b)在運動趨勢上大體一致，都是以紅河斷裂帶為邊界，滇西南塊體以沉降為主，川滇塊體南部以抬升為主，然而在具體數(shù)值上仍存在差異.本文得到的2011—2020年GPS垂向速度場與前人(Hao 等(2016)和Zhan等(2017))的結果有差異，原因可能如下：1)使用的GPS觀測數(shù)據(jù)時間跨度不同，本文數(shù)據(jù)源是2011—2020年GPS連續(xù)站觀測數(shù)據(jù)，Zhan 等(2017)和Hao 等(2016)所使用的數(shù)據(jù)源是2010—2015年GPS連續(xù)站觀測數(shù)據(jù)，不同時間跨度內(nèi)的構造形變信息存在不一致性，時間跨度越長對獲取可靠準確的速度場更有益；2)使用的水文負載模型數(shù)據(jù)不同.本文所使用的是LSDM模型數(shù)據(jù)，Zhan 等(2017)和Hao 等(2016)所使用的是GRACE模型數(shù)據(jù)，地球物理數(shù)據(jù)源不同會導致垂直位移數(shù)值不相同，改正效果也會存在差異；3)本文對GPS垂向位移進行共模誤差去除和噪聲模型估計，雖然Hao等(2016)也做了共模誤差和噪聲模型估計處理，但是噪聲模型和空間濾波方法不一樣；而Zhan等(2017)并未做共模誤差和噪聲模型估計，因此導致部分連續(xù)站速度場數(shù)值存在差異.

圖13 云南地區(qū)2010—2015年的GPS垂向速度場(a) 根據(jù)Hao等, 2016中速度場數(shù)據(jù)繪制; (b) Zhan等,2017文獻中速度場圖.Fig.13 GPS vertical velocity field in Yunnan region during 2010—2015(a) The velocity field is plotted according to the data from (Hao et al., 2016); (b) The velocity field diagram from (Zhan et al., 2017).

3.2 異常的YNGM連續(xù)站分析

從表1中可知，YNGM連續(xù)站的GPS垂向位移與LSDM形變的相關性最低，且RMS減少量為負值，說明水文負載形變不是引起YNGM連續(xù)站GPS垂向位移中季節(jié)性變化的原因.本文分別從YNGM連續(xù)站的觀測環(huán)境和建設質量情況以及周圍氣象站的降雨量和溫度等多方面進行綜合分析，判斷造成YNGM連續(xù)站異常的原因.從陸態(tài)網(wǎng)絡收集的資料可知建設的YNGM連續(xù)站為基巖墩標，且連續(xù)站附近沒有大型湖泊、河流.通過對收集到2011年1月—2020年4月時間段YNGM連續(xù)站觀測的MP1和MP2多路徑效應值及數(shù)據(jù)觀測有效率來判斷周圍環(huán)境的觀測質量，如圖14所示，MP1和MP2值大部分在0.5 m以下，且數(shù)據(jù)有效率大部分在90%以上，由此可知YNGM連續(xù)站觀測環(huán)境質量不是造成異常的原因.為了進一步分析YNGM連續(xù)站是否受降雨量與溫度的影響，本文收集到了YNGM連續(xù)站距離最近約為29 km氣象站在2011年1月—2020年6月間觀測的溫度與降雨量數(shù)據(jù)，收集的月降雨量和溫度如圖15所示，降雨量與溫度變化均存在著季節(jié)性變化，月平均溫度，月平均最高溫度和月平均最低溫度范圍分別在11.6°～25.7°、20.1°～33.2°、5.2°～21.6°，因此可以排除基巖的熱脹冷縮，冰川覆蓋、冰雪消融等情況造成YNGM異常的原因.YNGM連續(xù)站附近氣象站觀測的月累積降雨量比較大，由于LSDM模型數(shù)據(jù)得到的水文負載形變不包含地下水變化引起的形變.綜合因素考慮，地下水變化、LSDM模型和GPS解算的系統(tǒng)誤差可能是造成YNGM連續(xù)站中GPS垂向位移與LSDM形變相關性和一致性較差的主要原因.

圖14 YNGM連續(xù)站的數(shù)據(jù)觀測有效率和多路徑效應誤差，數(shù)據(jù)來源：中國地震局GNSS數(shù)據(jù)產(chǎn)品服務平臺(ftp:∥ftp.cgps.ac.cn)Fig.14 Data quality inspection results of YNGM station from 2011 to 2020, Data source: GNSS data product service platform of China Earthquake Administration (ftp:∥ftp.cgps.ac.cn)

圖15 YNGM連續(xù)站在2011年1月—2020年6月的月降水量和最高最低平均溫度，數(shù)據(jù)來源:中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn)Fig.15 Monthly precipitation and max, mean, min temperature changes at YNGM station from January 2011—June 2020. Data source: China Meteorological Data Network (http:∥data.cma.cn)

4 結論

本文使用時間跨度在2011—2020年的27個連續(xù)站的GPS垂向位移和LSDM形變來研究云南地區(qū)的垂向運動的季節(jié)性變化，對比分析了GPS垂向位移及GPS共模誤差與LSDM形變的定量關系和變化特征，并使用了SSA方法提取GPS垂向位移和LSDM形變的季節(jié)性信號和使用連續(xù)小波變換、交叉小波變換和小波相干性分析GPS垂向位移和LSDM形變的周期關系，最后得到了云南地區(qū)2011—2020年的GPS垂向速度場，本文得出的主要結論如下：

(1)GPS垂向位移和LSDM形變的整體運動趨勢一致，但是，LSDM形變的位移值范圍要整體小于GPS的，說明LSDM形變能解釋部分GPS垂向運動季節(jié)性變化.GPS共模誤差與LSDM形變的平均相關系數(shù)為0.75，若從GPS共模誤差中去除LSDM形變后， RMS減少量平均為35.72%，說明GPS共模誤差物理成因中大約有35.72%來源于LSDM形變.

(2)通過小波分析結果可知，大部分GPS連續(xù)站的周年平均相位相似性大小為0.9～1，說明LSDM形變與GPS垂向位移在年周期項變化是物理相關的，進一步說明水文負載形變是GPS年周期變化的主要驅動力，少部分GPS連續(xù)站是由其他因素和水文負載形變共同作用造成了GPS年周期項變化.相比于LSF方法提取的GPS垂向位移和LSDM季節(jié)性信號，SSA方法提取的季節(jié)性信號更加符合GPS垂向位移與LSDM實際的季節(jié)性運動.

(3)2011—2020年云南地區(qū)垂向速度場顯示滇西南塊體主要以0.01～1.9 mm·a-1的速率沉降，滇中塊體主要以0.13～2 mm·a-1速率抬升.

致謝“中國大陸構造環(huán)境監(jiān)測網(wǎng)絡”提供的GNSS數(shù)據(jù)產(chǎn)品，德國波斯坦地學中心提供的地表質量負載數(shù)據(jù)，在此一并表示感謝！