環(huán)渤海地區(qū)GPS臺站時間序列分析

邱亞輝，薄志毅，郎 博，岳彩亞(. 北京工業(yè)職業(yè)技術(shù)學(xué)院，北京 0004； . 中國礦業(yè)大學(xué)(北京)，北京 00037)

環(huán)渤海區(qū)域主要包括燕山、中朝、魯東海和華北等4個二級板塊的部分區(qū)域，其西鄰鄂爾多斯塊體，南鄰寬闊穩(wěn)定的華南塊體，北鄰中蒙塊體，東鄰太平洋。由于該區(qū)域復(fù)雜多變的環(huán)境和氣候，以及不同地方的地質(zhì)、地貌差異較大，并且長時間受到周圍板塊的相互作用，因而使得該地區(qū)板塊較為活躍。而GPS基準(zhǔn)站坐標(biāo)時間序列為大地測量學(xué)及地球動力學(xué)研究提供了寶貴的基礎(chǔ)數(shù)據(jù)[1-2]。隨著GPS定位技術(shù)的不斷更新與發(fā)展，其定位精度目前可以達(dá)到10-8～10-9量級。

高精度GPS站點坐標(biāo)時間序列是對站點連續(xù)觀測數(shù)據(jù)經(jīng)過精密解算后得到的高精度位置信息序列，它既是對站點位置連續(xù)性的描述，也是對站點穩(wěn)定的反映。研究測站坐標(biāo)時間序列可以獲得臺站噪聲模型并將其剔除，有利于分析板塊水平方向運(yùn)動特性和板塊高程方向運(yùn)動規(guī)律，以及對地震的預(yù)測。前人已經(jīng)對該領(lǐng)域做了許多研究，Johnson等分析了隨機(jī)漫步噪聲給GPS站點速率估計可能帶來的影響[3]；符養(yǎng)等通過對GPS坐標(biāo)時間序列計算分析，獲取了全球高程季節(jié)變化和整體振蕩運(yùn)動的特性[4]；袁林果等對香港GPS基準(zhǔn)站時間序列特征進(jìn)行分析，證實了香港GPS連續(xù)運(yùn)行參考站網(wǎng)存在明顯的本地季節(jié)性周期變化[5]；張鵬等研究分析了中國地殼運(yùn)動觀測網(wǎng)絡(luò)基準(zhǔn)站的坐標(biāo)位置變化規(guī)律，發(fā)現(xiàn)位于中國境內(nèi)的基準(zhǔn)站的時間序列具有一定的周期性，在高程方向分量的周期性最為明顯，并且時間序列的時頻特性表現(xiàn)出明顯的區(qū)域性[6]。本文基于中國陸態(tài)網(wǎng)絡(luò)中環(huán)渤海區(qū)域連續(xù)觀測得到的臺站坐標(biāo)時間序列數(shù)據(jù)，分析了臺站主要受到的噪聲類型和噪聲頻率特性，研究了環(huán)渤海區(qū)域地殼運(yùn)動特性。

1 GPS坐標(biāo)時間序列噪聲譜類型分析原理

1.1 Lomb周期圖法

通常采用快速離散傅里葉變換、譜分析、Lomb周期圖法及最大似然估計對臺站坐標(biāo)序列進(jìn)行噪聲類型分析。其中快速傅里葉變換和譜分析要求觀測數(shù)據(jù)無間隔，在漫長的觀測時間段中，很難滿足這種要求，必須對缺失天數(shù)進(jìn)行插值，因而會帶來一些誤差。因此本文使用Lomb周期圖法[7-8]，它可以計算缺失天數(shù)的坐標(biāo)序列功率譜。

典型的Lomb周期圖法公式基于離散傅里葉變換

(1)

式中，fn=n/T為采樣頻率，T為采樣總長度；xi為GPS站點坐標(biāo)時間序列。典型Lomb周期圖法的使用和譜分析的方法類似，采樣必須是無間隔的。

簡化Lomb周期圖法公式為

(2)

式中，τ可通過式(3)求出。

(3)

式中，ω為采樣頻率；ti為所取觀測天數(shù)，時間單位為1 d。簡化Lomb周期圖法可用來計算缺失天數(shù)的GPS測站坐標(biāo)時間序列噪聲譜，除此法之外，近年來發(fā)展起來的小波分析也可以做到。正是由于這個原因，Lomb周期圖法和小波分析方法受到了該方向研究者的一致好評，避免了插值帶來的誤差影響。

1.2 GPS坐標(biāo)時間序列譜指數(shù)反演

地球物理現(xiàn)象功率譜函數(shù)通常表示為1/fα,而GPS坐標(biāo)時間序列功率譜可以表示為如下形式[8-9]

P(f)=P0f-α

(4)

式中，P0和α為待求量，P0為常量；α表示譜指數(shù)。自然界中現(xiàn)象多表現(xiàn)為-3<α<-1,此時為非平穩(wěn)過程，即分形布朗運(yùn)動，其中包括隨機(jī)漫步噪聲(RWN)。對于平穩(wěn)過程現(xiàn)象，當(dāng)α為-1時，稱為分形白噪聲，當(dāng)α為0，則為白噪聲(WN)；當(dāng)α為1時，則為閃爍噪聲(FN)；除此之外的冪指數(shù)噪聲系統(tǒng)稱為有色噪聲[9-11]。

在計算GPS站點坐標(biāo)時間序列譜時，通常的計算方法是將式(4)兩邊取對數(shù)，形式如下

lnP(f)=lnP0-αlnf

(5)

式中，P(f)可通過式(2)計算得到；P0和α可通過曲線擬合或最小二乘擬合得到；f為對應(yīng)時間序列頻率。

在GPS坐標(biāo)時間序列的功率譜分析中，還可以用到基于小波的多分辨率分析方法，這種方法不僅可以分析出序列的譜性質(zhì)，還可以分離出序列的頻率特性，可以一目了然地看出組成信號噪聲的頻率。除此之外，多分辨率分析方法可以提取坐標(biāo)時間序列高程方向的周期項，據(jù)以前學(xué)者研究，大部分臺站高程方向時間序列存在半周年項和周年項，少部分臺站存在兩周年項。

2 實例分析與討論

數(shù)據(jù)來源：由中國地震局與總參測繪局、中科院、原國家測繪地理信息局、中國氣象局、教育部合作共建的中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(以下簡稱陸態(tài)網(wǎng)絡(luò))已基本完成，其由260個連續(xù)GNSS基準(zhǔn)站、30個連續(xù)相對重力觀測站和100個絕對重力觀測站組成[12]。自2007年以來積累了將近10年的觀測數(shù)據(jù)，為監(jiān)測中國大陸地殼運(yùn)動、重力場形態(tài)及變化，以及現(xiàn)代大地測量基準(zhǔn)系統(tǒng)的建立和維持提供了充足的數(shù)據(jù)資源。取中國陸態(tài)網(wǎng)絡(luò)中位于環(huán)渤海區(qū)域自2010年至2014年的GPS臺站時間序列數(shù)據(jù)。

陸態(tài)網(wǎng)絡(luò)中的環(huán)渤海地區(qū)測站布設(shè)情況為：主要行政區(qū)域包括北京、天津、遼寧、河北和山東共三省兩市，共同構(gòu)成一個“C”形區(qū)域。在這個“C”形區(qū)域內(nèi)布設(shè)了密密麻麻的GPS測站點。陸態(tài)網(wǎng)絡(luò)的布設(shè)特點是緊緊圍繞著地震頻發(fā)地區(qū)和地殼相對不穩(wěn)定地區(qū)，這樣能夠精確地監(jiān)測中國大陸板塊整體運(yùn)動和區(qū)域地殼板塊運(yùn)動，為區(qū)域地塊的進(jìn)一步劃分提供了數(shù)據(jù)基礎(chǔ)，從而可更加精確地預(yù)測地震。

2.1 臺站時間序列分析

本文使用的是環(huán)渤海地區(qū)的GPS連續(xù)觀測站的觀測數(shù)據(jù)，使用GMT功能繪出臺站的N、E、U這3個方向的時間序列差值圖。限于篇幅限制，選擇了位于河北省的HETS臺站和遼寧省的LNYK臺站進(jìn)行分析，如圖1所示。

圖1 臺站HETS和臺站LNYK時間序列差值

圖1中N代表南北方向，北方向為正；E代表東西方向，東方向為正；U代表高程方向，向上的方向為正。橫軸是觀測的天數(shù)，以年積日標(biāo)記，并且年與年之間采用相疊加的方式，縱軸是時間段內(nèi)的單天解與第一天解的差值，圖1(a)為臺站HETS時間序列差值圖，(b)為臺站LNYK時間序列差值圖。從時序圖中可以看出臺站的垂直方向波動比水平方向波動要大，并且具有明顯的周期特性，類似于正弦函數(shù)，波峰多出現(xiàn)在每年的6、7月，波谷多出現(xiàn)在每年的12月與次年的1月，這說明環(huán)渤海區(qū)域的GPS臺站周年運(yùn)動較為規(guī)則，并且季節(jié)性是影響該區(qū)域臺站高程方向運(yùn)動的主要因素。據(jù)以前學(xué)者研究，臺站高程方向波動的周期性常表現(xiàn)出周年性和半周年性及季節(jié)性，有些臺站會表現(xiàn)出2年的周期性。這種波動特性是由地球周期性的膨脹和收縮，還是由周圍地塊間的相互作用，還是由該區(qū)域復(fù)雜環(huán)境和氣候造成，需要進(jìn)一步的研究。

臺站的水平方向同樣在做微弱的運(yùn)動，雖然海洋潮汐負(fù)荷與海洋非潮汐負(fù)荷、大氣潮汐負(fù)荷與大氣非潮汐負(fù)荷及固體潮負(fù)荷等對測站水平方向的運(yùn)動有影響，但影響較小，因此臺站坐標(biāo)的水平方向的變化仍主要以板塊運(yùn)動為主。從時序圖中可以看出臺站有向北和西的微弱運(yùn)動趨勢，運(yùn)動速度分別為1.10和3.15 mm/a，從而可得在2010—2014年，該地區(qū)板塊有往西北向的微弱運(yùn)動，但運(yùn)動量級較小，為3.34 mm/a，并且該區(qū)域內(nèi)所細(xì)分的次級塊體間相對運(yùn)動不是很活躍?？偠灾?，位于該區(qū)域的GPS臺站在水平方向比垂直方向穩(wěn)定。

2.2 臺站噪聲特性分析

對臺站HETS和臺站LNYK做lomb周期圖法計算并畫出頻譜和解算天數(shù)的對應(yīng)關(guān)系，如圖2所示。功率譜密度的計算結(jié)果一般會隨所使用的方法和計算中選定的參數(shù)不同而有差別，但總體上都可以反映臺站的頻譜關(guān)系。

圖3是使用Lomb周期圖法經(jīng)計算后繪出的頻譜圖。N、E是水平面在南北、東西兩個方向的分量，U是高程方向。從圖中可看出N、E兩個方向的坐標(biāo)序列表現(xiàn)出了微弱的周期特性，而U方向的坐標(biāo)序列則表現(xiàn)出了明顯的周期特性。正如上所述，臺站水平方向和高程方向的時間序列周期都存在一年和半年兩種形式，除此之外高程方向還呈現(xiàn)出了季節(jié)特性，這說明臺站垂直方向比水平方向容易受到季節(jié)性影響，這與張鵬等用小波分析法對IGS站坐標(biāo)序列頻譜分析所得出的結(jié)論一致[6]。

測站在水平方向的周期性弱于測站高程方向的周期性，是因為板塊的運(yùn)動基本可分解為兩個方向上的運(yùn)動：一是在水平方向的運(yùn)動，二是在垂直方向的運(yùn)動。板塊在水平方向的運(yùn)動主要表現(xiàn)為整體的移動，除此之外由于受到區(qū)域內(nèi)次級板塊的相互作用以及其他因素的影響還夾雜著微弱的周期振蕩運(yùn)動。板塊在垂直方向主要表現(xiàn)周期性的上升或者下沉，這種運(yùn)動與重力激發(fā)、熱效應(yīng)、水文動力學(xué)等有關(guān)，由于該區(qū)域測站位于渤海附近，容易受到海洋潮汐負(fù)荷的影響，因此，該區(qū)域臺站的周年和半周年運(yùn)動特性與海洋潮汐密切相關(guān)。

圖2 HETS頻譜分析圖

圖3 LNYK頻譜分析圖

從圖3中各測站振幅分析可得，臺站HETS水平方向分量在2010—2013年間振幅幾乎呈現(xiàn)正弦形式的遞增和遞減，并且周期約為1年，在2013年后的時間段內(nèi)臺站HETS在南北向上一直在做急劇的向北運(yùn)動，而東西向上在做向西的運(yùn)動，則反映出該臺站所在地塊可能受到了某些外力作用發(fā)生了明顯的運(yùn)動趨勢變化。對臺站LNYK在水平方向分析可得，2010—2011年時間段內(nèi)，該站所在的地塊在做不明顯的周期運(yùn)動，且整體表現(xiàn)出向西南方向的運(yùn)動趨勢，但在2011年后，該地區(qū)板塊呈現(xiàn)了明顯的東北向運(yùn)動。由此可得兩個測站雖都位于環(huán)渤海地區(qū)，但位于兩個不同的次級塊體上，充分反映了環(huán)渤海區(qū)域地質(zhì)、地貌復(fù)雜，斷裂帶較多，次級塊體的運(yùn)動差異較大的特性。據(jù)不完全統(tǒng)計，僅在北京、天津、遼寧、河北、山東等地出現(xiàn)了多達(dá)十幾處斷裂。

利用Lomb周期圖法計算出功率譜后，可通過式(5)對常數(shù)P0和噪聲類型參數(shù)α采取線性擬合，見表1。表1中列出了HETS和LNYK這2個測站擬合后的值。由表1可得出α的值在-1和1之間，則為平穩(wěn)過程，且α值越大，所代表的隨機(jī)過程越平滑。由上述可知環(huán)渤海地區(qū)的臺站主要受分形白噪聲、白噪聲、閃爍噪聲的影響，而隨機(jī)漫步噪聲(α=-2)的成分可能性不大，這與黃立人對中國地殼運(yùn)動觀測網(wǎng)絡(luò)研究得出的結(jié)論基本一致[13]。但是α值在0附近波動，因此可從三者的權(quán)重分析得出臺站受到的白噪聲影響較大，而分形白噪聲和閃爍噪聲所占的成分相對較小。因此可得，臺站主要受到白噪聲的影響。

3 結(jié) 論

文中分析了環(huán)渤海地區(qū)GPS臺站速度場及各臺站的坐標(biāo)時間序列信息，得出以下結(jié)論：

表1 坐標(biāo)時間序列擬合結(jié)果

(1) 該區(qū)域臺站高程方向呈現(xiàn)明顯的周期振蕩特性，周期約為1年，振幅最大可達(dá)2 cm，波峰和波谷分別出現(xiàn)在每年的6、7月和12月及下一年的1月；從U方向可以看出，在所觀測的時間段內(nèi)，該區(qū)域臺站高程方向的波動較為規(guī)則，即此處的板塊相對而言較為穩(wěn)定。

(2) 該處的地殼有向西北方向運(yùn)動的趨勢，但在2010—2014年間運(yùn)動量級較小，運(yùn)動值僅在3 mm左右。因此可判定這段時間內(nèi)該區(qū)域塊體較為穩(wěn)定。

(3) 使用Lomb周期圖法對環(huán)渤海地區(qū)的臺站坐標(biāo)時間序列進(jìn)行了分析，得出臺站水平方向和高程方向都表現(xiàn)出了一年、半年的周期性，但高程方向的周期性相對更加明顯，并且高程方向同時表現(xiàn)出了季節(jié)性變化，這是水平方向不具有的，進(jìn)一步證實了結(jié)論(1)。

(4) 通過對α和P0線性擬合，得出α主要在0附近波動，并對影響臺站的噪聲做了類型分析，結(jié)果表明臺站主要受到白噪聲、分形白噪聲、閃爍噪聲的影響，在這3種噪聲中，最主要的是白噪聲的影響。但此處所有的臺站幾乎不會受到隨機(jī)漫步噪聲的影響?？偟膩碚f，在求高精度測站坐標(biāo)和速度場信息時，為了消除季節(jié)性對測站造成的影響，需要選擇相同季節(jié)進(jìn)行GPS重復(fù)觀測和解算；考慮影響臺站高程方向的因素表現(xiàn)出周年性和半周年性，因此在對環(huán)渤海地區(qū)高精度地殼運(yùn)動監(jiān)測時，最佳方案是通過長期跟蹤觀測，從而可以消除周期性誤差。

參考文獻(xiàn)：

[1] 姜衛(wèi)平，夏傳義，李昭，等.環(huán)境負(fù)載對區(qū)域GPS基準(zhǔn)站時間序列的影響分析[J].測繪學(xué)報，2014,43(12)：1217-1223.

[2] 陳俊勇，張鵬，武軍酈，等.關(guān)于在中國構(gòu)建全球?qū)Ш叫l(wèi)星國家級連續(xù)運(yùn)行站系統(tǒng)的思考[J].測繪學(xué)報，2007,36(4)：366-369.

[3] 符養(yǎng).中國大陸現(xiàn)今地殼形變與GPS坐標(biāo)時間序列分析[D].上海：中國科學(xué)院上海天文臺，2002.

[4] JOHNSON H O，AGNEW D C.Monument Motion and Measurement of Crustal Velocities[J].Geophysical Research Letters，1995,22(21)：2905-2908.

[5] 袁林果，丁曉利，陳武，等.香港GPS基準(zhǔn)站坐標(biāo)序列特征分析[J].地球物理學(xué)報，2008,51(5)：1372-1384.

[6] 范朋飛.高精度GPS站點坐標(biāo)時間序列分析與應(yīng)用[D].西安：長安大學(xué)，2013：22-23.

[7] SCARGLE J D.Studies in Astronomical Time Series Analysis.Ⅱ-statistical Aspects of Spectral Analysis of Unevenly Spaced Data[J].The Astrophysical Journal,1982,263:835-853.

[8] 張鵬，蔣志浩，秘金鐘，等.我國GPS跟蹤站數(shù)據(jù)處理與時間序列特征分析[J].武漢大學(xué)學(xué)報，2007,32(3)：251-254.

[9] 朱文耀，符養(yǎng)，李彥.GPS高程導(dǎo)出的全球高程振蕩運(yùn)動及季節(jié)變化[J].中國科學(xué)，2003,33(5)：470-481.

[10] AGNEW D C.The Time-domain Behavior of Power-law Noises[J].Geophys Res Lett，1992,19:333-336.

[11] 田云鋒.GPS位置時間序列中的中長期誤差研究[D].武漢：中國地震局地質(zhì)研究所，2011.

[12] 趙斌，聶兆生，黃勇.大規(guī)模GPS揭示的華北地區(qū)現(xiàn)今垂直運(yùn)動[J].大地測量與地球動力學(xué)，2014，34(5)：35-39.

[13] 黃立人.GPS基準(zhǔn)站坐標(biāo)分量時間序列的噪聲特性分析[J].大地測量與地球動力學(xué)，2006,26(2)：31-38.