李經(jīng)緯，占 偉，梁洪寶，王 勇，郭南男

(1.中國(guó)地震局第一監(jiān)測(cè)中心，天津 300180；2.天津城建大學(xué) 地質(zhì)與測(cè)繪學(xué)院，天津 300384)

0 引言

近年來(lái)，我國(guó)許多區(qū)域(如川滇地區(qū))GPS垂向存在顯著的周期性運(yùn)動(dòng)(其中年周期運(yùn)動(dòng)最為突出)(張飛鵬等，2002；劉任莉等，2013；梁洪寶等，2015；王巖等，2018)，盛傳貞等(2014)認(rèn)為滇西地區(qū)GPS坐標(biāo)變化時(shí)間序列的垂向分量中，普遍包含有明顯的年周期非構(gòu)造形變波動(dòng)，高值可達(dá)12 mm；Zhan等(2017)基于2010—2015年云南27個(gè)GPS連續(xù)站垂向時(shí)間序列得出該時(shí)段云南地區(qū)GPS連續(xù)站垂向年周期振幅平均值為9.7 mm，認(rèn)為云南地區(qū)季節(jié)性的降雨與GPS垂向年周期運(yùn)動(dòng)相關(guān)性較好。GPS垂向周期性運(yùn)動(dòng)研究有重要的研究意義和實(shí)際應(yīng)用價(jià)值，目前我國(guó)GPS流動(dòng)觀測(cè)的觀測(cè)數(shù)據(jù)量太少、無(wú)法消除周期性運(yùn)動(dòng)的影響，導(dǎo)致得到的垂向線性運(yùn)動(dòng)速率結(jié)果精度不高，GPS垂向速度場(chǎng)應(yīng)用較少(方穎等，2014；朱爽等，2017；李長(zhǎng)軍等，2018)，因此精確獲取垂向周期性運(yùn)動(dòng)信息能夠提高GPS流動(dòng)觀測(cè)的垂向線性運(yùn)動(dòng)速率精度，極大地促進(jìn)GPS垂向速度場(chǎng)的應(yīng)用。

一些學(xué)者使用了不同的方法獲取GPS流動(dòng)站的垂向周期性運(yùn)動(dòng)參數(shù)(Liangetal，2013；占偉等，2016)，但是對(duì)這些方法的效果評(píng)估、分析較少，尤其是中國(guó)大陸幅員廣闊，不同區(qū)域影響垂向周期性運(yùn)動(dòng)的主要因素不一，不同方法的適用性也不同，因此對(duì)GPS垂向周期性運(yùn)動(dòng)參數(shù)獲取方法進(jìn)行全面的評(píng)估十分必要。本文以云南為研究區(qū)域，選用4種垂向年周期運(yùn)動(dòng)參數(shù)獲取方法，以研究區(qū)域內(nèi)2010—2015年陸態(tài)網(wǎng)絡(luò)26個(gè)GPS連續(xù)站觀測(cè)數(shù)據(jù)得到的垂向年周期運(yùn)動(dòng)結(jié)果，評(píng)估這4種方法的效果，綜合評(píng)估各方法的優(yōu)缺點(diǎn)，對(duì)不同的區(qū)域如何獲取垂向年周期運(yùn)動(dòng)參數(shù)給出建議。

1 方法介紹

常用的垂向年周期運(yùn)動(dòng)獲取方法分為2類，一類是根據(jù)一些物理模型(如負(fù)荷模型或者GRACE)直接求解某一空間位置的垂向年周期運(yùn)動(dòng)參數(shù)；另一類是對(duì)區(qū)域內(nèi)GPS連續(xù)站實(shí)測(cè)的垂向年周期運(yùn)動(dòng)參數(shù)進(jìn)行空間內(nèi)插。本文選用4種GPS垂向年周期運(yùn)動(dòng)獲取方法，分別為：

(1)方法1：負(fù)荷模型方法。地球表面各種負(fù)荷變化是GPS垂向周期性運(yùn)動(dòng)的主要影響因素，而在云南地區(qū)，陸地水的質(zhì)量變化是該區(qū)域GPS垂向周期性運(yùn)動(dòng)的主要影響因素(姜衛(wèi)平等，2013；Jiangetal，2017)，因此筆者從EOST Loading Service(1)http：//loading.-ustrasbg.fr/ditrf_form.php.下載了2010—2015年全球陸地水文負(fù)荷模型(GLDAS)解算的垂向位移格網(wǎng)數(shù)據(jù)(CM)，數(shù)據(jù)空間分辨率為0.5°×0.5°，時(shí)間分辨率為3 h?；谠颇系貐^(qū)各格網(wǎng)點(diǎn)垂向位移時(shí)間序列，采用線性+年周期+半年周期方法擬合得到各格網(wǎng)點(diǎn)年周期運(yùn)動(dòng)正弦波振幅A和余弦波振幅B，為了與GPS連續(xù)站結(jié)果比對(duì)，采用多核函數(shù)法空間插值獲取云南地區(qū)26個(gè)GPS連續(xù)站所在空間位置的年周期運(yùn)動(dòng)正弦波振幅A1和余弦波振幅B1。

(2)方法2：GRACE模型方法。利用GRACE數(shù)據(jù)給出的地球重力場(chǎng)球諧系數(shù)(Bruinsmaetal，2010)和負(fù)荷勒夫數(shù)(Farrell，1972)，可以根據(jù)彈性負(fù)荷形變理論獲取負(fù)荷引起的形變信息(Wahretal，1998；Kusche，2005；Fuetal，2012)。筆者按照Hao等(2016)的方法處理了CSR提供的GRACE-RL05數(shù)據(jù)：使用GCM05S模型去掉平均重力場(chǎng)；用SLR觀測(cè)的C20項(xiàng)系數(shù)項(xiàng)和一階項(xiàng)系數(shù)替換模型中相應(yīng)值；采用高斯濾波方法去除高階模型噪聲，最終得到了26個(gè)GPS連續(xù)站所在空間位置的垂向位移結(jié)果(采樣率為每月)，采用線性+年周期+半年周期方法擬合得到各GPS連續(xù)站的年周期運(yùn)動(dòng)正弦波振幅A2和余弦波振幅B2。

(3)方法3：空間加權(quán)法?；赯han等(2017)給出的2010—2015年云南地區(qū)26個(gè)GPS連續(xù)站垂向位移時(shí)間序列(進(jìn)行了預(yù)處理，包括去除了儀器更換等原因產(chǎn)生的突跳)，采用線性+年周期+半年周期方法擬合得到各連續(xù)站年周期運(yùn)動(dòng)正弦波振幅A0和余弦波振幅B0。對(duì)某一個(gè)GPS連續(xù)站(記為插值點(diǎn))，根據(jù)其余25個(gè)連續(xù)站的年周期運(yùn)動(dòng)正弦波振幅A0和余弦波振幅B0，采用空間加權(quán)法插值得到插值點(diǎn)的年周期運(yùn)動(dòng)正弦波振幅A3和余弦波振幅B3：

(1)

式中：dj為插值點(diǎn)與第j個(gè)GPS連續(xù)站的水平距離；A0j和B0j為第j個(gè)GPS連續(xù)站的年周期運(yùn)動(dòng)正弦波和余弦波振幅。

(4)方法4：多核函數(shù)法。該方法計(jì)算過(guò)程與方法3類似，對(duì)于每一個(gè)GPS連續(xù)站，都是其它25個(gè)GPS連續(xù)站垂向年周期運(yùn)動(dòng)參數(shù)(A0和B0)空間插值得到插值點(diǎn)的年周期運(yùn)動(dòng)參數(shù)(A4和B4)，唯一的區(qū)別就是空間插值采用多核函數(shù)法。多核函數(shù)法由美國(guó)Hardy(1971)提出，其空間插值計(jì)算過(guò)程如下：

(2)

式中：f(λ，φ)為多核函數(shù)；sj(λ，φ，λj，φj)為多核函數(shù)；dj1.1為球面上2點(diǎn)間的大地線長(zhǎng)度(單位為km)；cj為待定系數(shù)(楊國(guó)華等，2012)。將所有可用的數(shù)據(jù)點(diǎn)作為核函數(shù)點(diǎn)，可以獲得數(shù)學(xué)解析式，計(jì)算任一空間位置的結(jié)果。

2 效果評(píng)估

本文依據(jù)4種方法與GPS連續(xù)觀測(cè)得到的垂向年周期運(yùn)動(dòng)振幅差異和一致性評(píng)估各方法的效果。

2.1 振幅差異

圖1給出了GPS連續(xù)觀測(cè)得到的垂向年周期運(yùn)動(dòng)與4種方法結(jié)果之間的差異，從中可以看出：①方法1和方法2與GPS連續(xù)觀測(cè)的垂向年周期運(yùn)動(dòng)振幅差異較大，且方向上大體一致，說(shuō)明這2種方法得到的垂向年周期運(yùn)動(dòng)與GPS連續(xù)觀測(cè)結(jié)果存在系統(tǒng)的時(shí)間差；②方法3和方法4與GPS連續(xù)觀測(cè)的垂向年周期運(yùn)動(dòng)振幅差異比前2種方法偏小，且差異的方向較為雜亂、沒(méi)有表現(xiàn)出較為明顯的規(guī)律，說(shuō)明后2種方法更為接近GPS連續(xù)站得到的垂向年周期運(yùn)動(dòng)。表1給出了4種方法與GPS實(shí)測(cè)結(jié)果的振幅差異統(tǒng)計(jì)結(jié)果，4種方法與GPS連續(xù)站實(shí)測(cè)的振幅差異平均值分別為4.97，3.31，1.38，1.43 mm，而26個(gè)GPS連續(xù)站實(shí)測(cè)的垂向年周期運(yùn)動(dòng)振幅平均值為9.87 mm，可見(jiàn)方法1和方法2與GPS實(shí)測(cè)結(jié)果的差異較為顯著(分別達(dá)到實(shí)測(cè)振幅的50.35%和34.95%)。

圖1 4種方法與GPS連續(xù)觀測(cè)的年周期運(yùn)動(dòng)差異Fig.1 Difference of crustal annual motions between continuous GPS and four methods

表1 4種方法與GPS連續(xù)觀測(cè)的年周期運(yùn)動(dòng)振幅差異Tab.1 Difference of crustal annual motion amplitude between continuous GPS and four methods (單位：mm)

2.2 一致性評(píng)估

本文用WRMS減速比(Fu，F(xiàn)reymueller，2012)來(lái)評(píng)估各方法得到的垂向年周期參數(shù)與GPS連續(xù)站時(shí)間序列的一致性，定義為：

(3)

式中：RatioW表示W(wǎng)RMS減速比；WRMSG是2010—2015年26個(gè)GPS連續(xù)站垂向時(shí)間序列去除線性后WRMS值；WRMSm是GPS垂向時(shí)間序列去除線性后再去除年周期運(yùn)動(dòng)后的殘差WRMS值；WRMSf是GPS垂向時(shí)間序列去除線性和年周期運(yùn)動(dòng)后的殘差WRMS值。

RatioW能夠反映GPS連續(xù)站的垂向時(shí)間序列和各方法求得的年周期運(yùn)動(dòng)相位和振幅的一致性。圖2給出了4種方法的RatioW結(jié)果，從中可以看出后2種方法的絕大多數(shù)測(cè)站RatioW大于0.9，遠(yuǎn)遠(yuǎn)優(yōu)于前2種方法，但方法3和4中YNZD站RatioW相對(duì)較小(圖2c，d)。表2給出了4種方法的RatioW統(tǒng)計(jì)結(jié)果，其平均值分別為0.70，0.84，0.95，0.94，表明方法3和4效果大體相當(dāng)，優(yōu)于方法1和2，相對(duì)而言，方法2的效果優(yōu)于方法1。

表2 4種方法的RatioW統(tǒng)計(jì)結(jié)果Tab.2 Static results of WRMS reduction ratio of the four methods

圖2 4種方法的RatioW
Fig.2 WRMS reduction ratio of four methods

3 討論

從4種方法的整體評(píng)估結(jié)果可以看出，方法1和方法2得到的垂向年周期運(yùn)動(dòng)無(wú)論是振幅大小還是運(yùn)動(dòng)方向，均與GPS連續(xù)觀測(cè)的結(jié)果存在較大的差異，而方法3和方法4利用GPS連續(xù)觀測(cè)進(jìn)行空間內(nèi)插，得到的垂向年周期運(yùn)動(dòng)與GPS連續(xù)觀測(cè)結(jié)果一致性較好。

由于全球負(fù)荷模型是基于全球觀測(cè)數(shù)據(jù)建立的模型，模型建立時(shí)使用的中國(guó)大陸范圍內(nèi)的實(shí)測(cè)數(shù)據(jù)較少，因此在中國(guó)大陸尤其是云南這樣的局部區(qū)域與實(shí)測(cè)數(shù)據(jù)一致性較差，使得全球負(fù)荷模型得到的垂向年周期運(yùn)動(dòng)與GPS實(shí)測(cè)結(jié)果差異較大。

方法2的GRACE模型得到的垂向周期性信號(hào)反映的是區(qū)域空間尺度達(dá)300～400 km的整體形變，GPS垂向周期性信號(hào)不僅反映了區(qū)域整體形變，還可能疊加一些局部效應(yīng)(姜衛(wèi)平等，2013；Yanetal，2016)，因此GRACE得到垂向年周期運(yùn)動(dòng)振幅與GPS實(shí)測(cè)結(jié)果存在一定差異。

需要說(shuō)明的是：不同的全球負(fù)荷模型(方法)在中國(guó)大陸的符合程度存在差異，同時(shí)GRACE數(shù)據(jù)處理方法(方法2)也存在一定的差異，這2種方法在中國(guó)大陸的應(yīng)用優(yōu)化值得更為深入的研究。

云南地區(qū)GPS連續(xù)站較多，且區(qū)域垂向周期性運(yùn)動(dòng)空間一致性較好，各點(diǎn)的垂向年周期運(yùn)動(dòng)無(wú)論是振幅大小還是初相均較為一致(Zhanetal，2017)，因此2種空間插值的方法效果要優(yōu)于前2種方法，而且方法3和方法4評(píng)估結(jié)果大體相當(dāng)，說(shuō)明只要區(qū)域內(nèi)垂向周期性運(yùn)動(dòng)空間一致性較好，即使選用不同的空間內(nèi)插方法，對(duì)結(jié)果影響不大。值得注意的是，在26個(gè)GPS連續(xù)站中，YNZD站空間位置最為突出(最北端)，周邊的GPS連續(xù)站與之距離較遠(yuǎn)，因此在方法3和方法4中，該站利用其余25個(gè)測(cè)站空間內(nèi)插得到的垂向年周期運(yùn)動(dòng)參數(shù)與實(shí)測(cè)結(jié)果差異最大，這說(shuō)明2種空間插值方法要求GPS連續(xù)站的分布盡量均勻。

上述4種方法各有優(yōu)劣，雖然在本文的實(shí)驗(yàn)中，方法1和方法2的效果不如方法3和4，但前2種方法可以直接利用已有的全球負(fù)荷模型和GRACE給出的任意一點(diǎn)垂向位求得任意一點(diǎn)的垂向年周期運(yùn)動(dòng)。方法3和方法4雖然效果較好，但是對(duì)區(qū)域內(nèi)GPS連續(xù)站的數(shù)量和分布有一定的要求，在實(shí)際應(yīng)用中可能受到限制。事實(shí)上，目前中國(guó)大陸仍有一些區(qū)域GPS連續(xù)站較少，以陸態(tài)網(wǎng)絡(luò)為例，目前中國(guó)大陸范圍內(nèi)只有260個(gè)GPS連續(xù)站，平均站間距為130 km，僅少數(shù)區(qū)域如首都圈、南北地震帶GPS連續(xù)站密度相對(duì)較高(云南地區(qū)GPS連續(xù)站平均站間距79 km)，在西藏、東北、華南等區(qū)域GPS連續(xù)站較少，因此方法1和2在這些GPS連續(xù)站較少的區(qū)域仍然有一定的實(shí)際使用價(jià)值。

4 結(jié)論

本文選用4種獲取GPS垂向年周期運(yùn)動(dòng)參數(shù)的方法，以GPS連續(xù)觀測(cè)得到的垂向年周期運(yùn)動(dòng)結(jié)果評(píng)估了這4種方法的效果，結(jié)果表明：由于云南地區(qū)GPS連續(xù)站較多且垂向年周期運(yùn)動(dòng)空間一致性較好，利用GPS連續(xù)觀測(cè)的垂向年周期運(yùn)動(dòng)參數(shù)進(jìn)行空間插值(方法3和方法4)都能夠取得較好的效果，GRACE模型方法(方法2)效果稍差，負(fù)荷模型方法(方法1)效果最差，因此對(duì)于獲取區(qū)域內(nèi)(如GPS流動(dòng)站)的垂向年周期運(yùn)動(dòng)參數(shù)，在GPS連續(xù)觀測(cè)站較多的區(qū)域，建議采用空間加權(quán)法(方法3)或者多核函數(shù)法(方法4)；在GPS連續(xù)站較少的區(qū)域，建議采用GRACE模型方法。

由于本文的主要目的是驗(yàn)證垂向年周期運(yùn)動(dòng)信息的構(gòu)建模型有效性，因此只考慮了季節(jié)性運(yùn)動(dòng)中最主要的年周期運(yùn)動(dòng)，但本方法也適用于其它周期性運(yùn)動(dòng)(如半年周期運(yùn)動(dòng))。

感謝NASA提供GRACE數(shù)據(jù)、EOST提供全球陸地水文負(fù)荷模型垂向位移數(shù)據(jù)。