胡順強(qiáng), 王坦, 管雅慧, 楊振宇*

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

0 引言

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

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

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

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

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

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

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

IQR=Q2-Q1.

(1)

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

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

(2)

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

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

2)使用最小二乘線性擬合方法對(duì)GPS垂向位移的階躍項(xiàng)進(jìn)行改正.

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

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

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

(3)

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

(4)

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

圖2 預(yù)處理后的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 地表負(fù)載模型數(shù)據(jù)

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

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

1.3 奇異譜分析

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

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

(5)

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

(6)

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

(7)

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

(8)

(9)

1.4 小波分析

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

(10)

WXY=WXWY*,

(11)

式(11)中，*表示復(fù)共軛，小波功率受邊緣效應(yīng)影響較大的區(qū)域?yàn)橛绊懚?Cone of Influence, COI).為了準(zhǔn)確描述xn,yn的相位關(guān)系，選擇COI區(qū)域之外一系列相位角的圓域平均值來衡量：

(12)

交叉小波的相似性定義為

ρi=cos(am),

(13)

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

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

(14)

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

2 結(jié)果

2.1 GPS垂向位移與水文負(fù)載形變比較

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

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

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

(15)

(16)

2.2 GPS共模誤差與水文負(fù)載形變比較

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

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

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

2.3 小波譜分析

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

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

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

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

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

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

圖8 XWT與LSF兩種方法獲取27個(gè)GPS連續(xù)站與LSDM的年周期相位差對(duì)比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é)性信號(hào)提取

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

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

圖9 YNTC連續(xù)站的前20階GPS(a)與LSDM(b)時(shí)間序列的w-correlation分析結(jié)果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)時(shí)間序列的前20階的RC方差貢獻(xiàn)率統(tǒng)計(jì)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é)性信號(hào)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é)性信號(hào)的相關(guān)系數(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減少量為負(fù)值之外，其他連續(xù)站的RMS減少量均為正值.因此，為了獲得高精度的云南地區(qū)27個(gè)GPS垂向速度場(chǎng)，除了YNGM連續(xù)站之外，其他GPS連續(xù)站都使用水文負(fù)載模型數(shù)據(jù)進(jìn)行改正處理，然后對(duì)所有連續(xù)站再使用主成分分析方法(Dong et al., 2006)進(jìn)行共模誤差剔除，國(guó)內(nèi)外學(xué)者普遍認(rèn)為白噪聲和閃爍噪聲組合的噪聲模型能夠代表大部分GPS垂向位移的噪聲特性(Williams et al., 2004;Didova et al., 2016;He et al., 2019)，因此，在經(jīng)過水文負(fù)載模型和共模誤差預(yù)處理后，最后使用Hector軟件(Bos et al., 2013)通過白噪聲和閃爍噪聲組合模型估計(jì)所有連續(xù)站的速率及不確定度.圖12為最后得到的2011—2020年云南地區(qū)的垂向速度場(chǎng)，從圖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垂向速度場(chǎng)(紅色和藍(lán)色箭頭分別代表抬升和沉降速率)Fig.12 GPS vertical velocity field in Yunnan region during 2011—2020(Red and blue respectively represent the up and down vertical rates)

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

3 討論

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

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

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

圖14 YNGM連續(xù)站的數(shù)據(jù)觀測(cè)有效率和多路徑效應(yīng)誤差，數(shù)據(jù)來源：中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)(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ù)來源:中國(guó)氣象數(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 結(jié)論

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

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

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

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

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