环渤海地区GPS台站时间序列分析

邱亚辉，薄志毅，郎 博，岳彩亚(. 北京工业职业技术学院，北京 0004； . 中国矿业大学(北京)，北京 00037)

环渤海区域主要包括燕山、中朝、鲁东海和华北等4个二级板块的部分区域，其西邻鄂尔多斯块体，南邻宽阔稳定的华南块体，北邻中蒙块体，东邻太平洋。由于该区域复杂多变的环境和气候，以及不同地方的地质、地貌差异较大，并且长时间受到周围板块的相互作用，因而使得该地区板块较为活跃。而GPS基准站坐标时间序列为大地测量学及地球动力学研究提供了宝贵的基础数据[1-2]。随着GPS定位技术的不断更新与发展，其定位精度目前可以达到10-8～10-9量级。

高精度GPS站点坐标时间序列是对站点连续观测数据经过精密解算后得到的高精度位置信息序列，它既是对站点位置连续性的描述，也是对站点稳定的反映。研究测站坐标时间序列可以获得台站噪声模型并将其剔除，有利于分析板块水平方向运动特性和板块高程方向运动规律，以及对地震的预测。前人已经对该领域做了许多研究，Johnson等分析了随机漫步噪声给GPS站点速率估计可能带来的影响[3]；符养等通过对GPS坐标时间序列计算分析，获取了全球高程季节变化和整体振荡运动的特性[4]；袁林果等对香港GPS基准站时间序列特征进行分析，证实了香港GPS连续运行参考站网存在明显的本地季节性周期变化[5]；张鹏等研究分析了中国地壳运动观测网络基准站的坐标位置变化规律，发现位于中国境内的基准站的时间序列具有一定的周期性，在高程方向分量的周期性最为明显，并且时间序列的时频特性表现出明显的区域性[6]。本文基于中国陆态网络中环渤海区域连续观测得到的台站坐标时间序列数据，分析了台站主要受到的噪声类型和噪声频率特性，研究了环渤海区域地壳运动特性。

1 GPS坐标时间序列噪声谱类型分析原理

1.1 Lomb周期图法

通常采用快速离散傅里叶变换、谱分析、Lomb周期图法及最大似然估计对台站坐标序列进行噪声类型分析。其中快速傅里叶变换和谱分析要求观测数据无间隔，在漫长的观测时间段中，很难满足这种要求，必须对缺失天数进行插值，因而会带来一些误差。因此本文使用Lomb周期图法[7-8]，它可以计算缺失天数的坐标序列功率谱。

典型的Lomb周期图法公式基于离散傅里叶变换

(1)

式中，fn=n/T为采样频率，T为采样总长度；xi为GPS站点坐标时间序列。典型Lomb周期图法的使用和谱分析的方法类似，采样必须是无间隔的。

简化Lomb周期图法公式为

(2)

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

(3)

式中，ω为采样频率；ti为所取观测天数，时间单位为1 d。简化Lomb周期图法可用来计算缺失天数的GPS测站坐标时间序列噪声谱，除此法之外，近年来发展起来的小波分析也可以做到。正是由于这个原因，Lomb周期图法和小波分析方法受到了该方向研究者的一致好评，避免了插值带来的误差影响。

1.2 GPS坐标时间序列谱指数反演

地球物理现象功率谱函数通常表示为1/fα,而GPS坐标时间序列功率谱可以表示为如下形式[8-9]

P(f)=P0f-α

(4)

式中，P0和α为待求量，P0为常量；α表示谱指数。自然界中现象多表现为-3<α<-1,此时为非平稳过程，即分形布朗运动，其中包括随机漫步噪声(RWN)。对于平稳过程现象，当α为-1时，称为分形白噪声，当α为0，则为白噪声(WN)；当α为1时，则为闪烁噪声(FN)；除此之外的幂指数噪声系统称为有色噪声[9-11]。

在计算GPS站点坐标时间序列谱时，通常的计算方法是将式(4)两边取对数，形式如下

lnP(f)=lnP0-αlnf

(5)

式中，P(f)可通过式(2)计算得到；P0和α可通过曲线拟合或最小二乘拟合得到；f为对应时间序列频率。

在GPS坐标时间序列的功率谱分析中，还可以用到基于小波的多分辨率分析方法，这种方法不仅可以分析出序列的谱性质，还可以分离出序列的频率特性，可以一目了然地看出组成信号噪声的频率。除此之外，多分辨率分析方法可以提取坐标时间序列高程方向的周期项，据以前学者研究，大部分台站高程方向时间序列存在半周年项和周年项，少部分台站存在两周年项。

2 实例分析与讨论

数据来源：由中国地震局与总参测绘局、中科院、原国家测绘地理信息局、中国气象局、教育部合作共建的中国大陆构造环境监测网络(以下简称陆态网络)已基本完成，其由260个连续GNSS基准站、30个连续相对重力观测站和100个绝对重力观测站组成[12]。自2007年以来积累了将近10年的观测数据，为监测中国大陆地壳运动、重力场形态及变化，以及现代大地测量基准系统的建立和维持提供了充足的数据资源。取中国陆态网络中位于环渤海区域自2010年至2014年的GPS台站时间序列数据。

陆态网络中的环渤海地区测站布设情况为：主要行政区域包括北京、天津、辽宁、河北和山东共三省两市，共同构成一个“C”形区域。在这个“C”形区域内布设了密密麻麻的GPS测站点。陆态网络的布设特点是紧紧围绕着地震频发地区和地壳相对不稳定地区，这样能够精确地监测中国大陆板块整体运动和区域地壳板块运动，为区域地块的进一步划分提供了数据基础，从而可更加精确地预测地震。

2.1 台站时间序列分析

本文使用的是环渤海地区的GPS连续观测站的观测数据，使用GMT功能绘出台站的N、E、U这3个方向的时间序列差值图。限于篇幅限制，选择了位于河北省的HETS台站和辽宁省的LNYK台站进行分析，如图1所示。

图1 台站HETS和台站LNYK时间序列差值

图1中N代表南北方向，北方向为正；E代表东西方向，东方向为正；U代表高程方向，向上的方向为正。横轴是观测的天数，以年积日标记，并且年与年之间采用相叠加的方式，纵轴是时间段内的单天解与第一天解的差值，图1(a)为台站HETS时间序列差值图，(b)为台站LNYK时间序列差值图。从时序图中可以看出台站的垂直方向波动比水平方向波动要大，并且具有明显的周期特性，类似于正弦函数，波峰多出现在每年的6、7月，波谷多出现在每年的12月与次年的1月，这说明环渤海区域的GPS台站周年运动较为规则，并且季节性是影响该区域台站高程方向运动的主要因素。据以前学者研究，台站高程方向波动的周期性常表现出周年性和半周年性及季节性，有些台站会表现出2年的周期性。这种波动特性是由地球周期性的膨胀和收缩，还是由周围地块间的相互作用，还是由该区域复杂环境和气候造成，需要进一步的研究。

台站的水平方向同样在做微弱的运动，虽然海洋潮汐负荷与海洋非潮汐负荷、大气潮汐负荷与大气非潮汐负荷及固体潮负荷等对测站水平方向的运动有影响，但影响较小，因此台站坐标的水平方向的变化仍主要以板块运动为主。从时序图中可以看出台站有向北和西的微弱运动趋势，运动速度分别为1.10和3.15 mm/a，从而可得在2010—2014年，该地区板块有往西北向的微弱运动，但运动量级较小，为3.34 mm/a，并且该区域内所细分的次级块体间相对运动不是很活跃。总而言之，位于该区域的GPS台站在水平方向比垂直方向稳定。

2.2 台站噪声特性分析

对台站HETS和台站LNYK做lomb周期图法计算并画出频谱和解算天数的对应关系，如图2所示。功率谱密度的计算结果一般会随所使用的方法和计算中选定的参数不同而有差别，但总体上都可以反映台站的频谱关系。

图3是使用Lomb周期图法经计算后绘出的频谱图。N、E是水平面在南北、东西两个方向的分量，U是高程方向。从图中可看出N、E两个方向的坐标序列表现出了微弱的周期特性，而U方向的坐标序列则表现出了明显的周期特性。正如上所述，台站水平方向和高程方向的时间序列周期都存在一年和半年两种形式，除此之外高程方向还呈现出了季节特性，这说明台站垂直方向比水平方向容易受到季节性影响，这与张鹏等用小波分析法对IGS站坐标序列频谱分析所得出的结论一致[6]。

测站在水平方向的周期性弱于测站高程方向的周期性，是因为板块的运动基本可分解为两个方向上的运动：一是在水平方向的运动，二是在垂直方向的运动。板块在水平方向的运动主要表现为整体的移动，除此之外由于受到区域内次级板块的相互作用以及其他因素的影响还夹杂着微弱的周期振荡运动。板块在垂直方向主要表现周期性的上升或者下沉，这种运动与重力激发、热效应、水文动力学等有关，由于该区域测站位于渤海附近，容易受到海洋潮汐负荷的影响，因此，该区域台站的周年和半周年运动特性与海洋潮汐密切相关。

图2 HETS频谱分析图

图3 LNYK频谱分析图

从图3中各测站振幅分析可得，台站HETS水平方向分量在2010—2013年间振幅几乎呈现正弦形式的递增和递减，并且周期约为1年，在2013年后的时间段内台站HETS在南北向上一直在做急剧的向北运动，而东西向上在做向西的运动，则反映出该台站所在地块可能受到了某些外力作用发生了明显的运动趋势变化。对台站LNYK在水平方向分析可得，2010—2011年时间段内，该站所在的地块在做不明显的周期运动，且整体表现出向西南方向的运动趋势，但在2011年后，该地区板块呈现了明显的东北向运动。由此可得两个测站虽都位于环渤海地区，但位于两个不同的次级块体上，充分反映了环渤海区域地质、地貌复杂，断裂带较多，次级块体的运动差异较大的特性。据不完全统计，仅在北京、天津、辽宁、河北、山东等地出现了多达十几处断裂。

利用Lomb周期图法计算出功率谱后，可通过式(5)对常数P0和噪声类型参数α采取线性拟合，见表1。表1中列出了HETS和LNYK这2个测站拟合后的值。由表1可得出α的值在-1和1之间，则为平稳过程，且α值越大，所代表的随机过程越平滑。由上述可知环渤海地区的台站主要受分形白噪声、白噪声、闪烁噪声的影响，而随机漫步噪声(α=-2)的成分可能性不大，这与黄立人对中国地壳运动观测网络研究得出的结论基本一致[13]。但是α值在0附近波动，因此可从三者的权重分析得出台站受到的白噪声影响较大，而分形白噪声和闪烁噪声所占的成分相对较小。因此可得，台站主要受到白噪声的影响。

3 结 论

文中分析了环渤海地区GPS台站速度场及各台站的坐标时间序列信息，得出以下结论：

表1 坐标时间序列拟合结果

(1) 该区域台站高程方向呈现明显的周期振荡特性，周期约为1年，振幅最大可达2 cm，波峰和波谷分别出现在每年的6、7月和12月及下一年的1月；从U方向可以看出，在所观测的时间段内，该区域台站高程方向的波动较为规则，即此处的板块相对而言较为稳定。

(2) 该处的地壳有向西北方向运动的趋势，但在2010—2014年间运动量级较小，运动值仅在3 mm左右。因此可判定这段时间内该区域块体较为稳定。

(3) 使用Lomb周期图法对环渤海地区的台站坐标时间序列进行了分析，得出台站水平方向和高程方向都表现出了一年、半年的周期性，但高程方向的周期性相对更加明显，并且高程方向同时表现出了季节性变化，这是水平方向不具有的，进一步证实了结论(1)。

(4) 通过对α和P0线性拟合，得出α主要在0附近波动，并对影响台站的噪声做了类型分析，结果表明台站主要受到白噪声、分形白噪声、闪烁噪声的影响，在这3种噪声中，最主要的是白噪声的影响。但此处所有的台站几乎不会受到随机漫步噪声的影响。总的来说，在求高精度测站坐标和速度场信息时，为了消除季节性对测站造成的影响，需要选择相同季节进行GPS重复观测和解算；考虑影响台站高程方向的因素表现出周年性和半周年性，因此在对环渤海地区高精度地壳运动监测时，最佳方案是通过长期跟踪观测，从而可以消除周期性误差。

参考文献：

[1] 姜卫平，夏传义，李昭，等.环境负载对区域GPS基准站时间序列的影响分析[J].测绘学报，2014,43(12)：1217-1223.

[2] 陈俊勇，张鹏，武军郦，等.关于在中国构建全球导航卫星国家级连续运行站系统的思考[J].测绘学报，2007,36(4)：366-369.

[3] 符养.中国大陆现今地壳形变与GPS坐标时间序列分析[D].上海：中国科学院上海天文台，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基准站坐标序列特征分析[J].地球物理学报，2008,51(5)：1372-1384.

[6] 范朋飞.高精度GPS站点坐标时间序列分析与应用[D].西安：长安大学，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跟踪站数据处理与时间序列特征分析[J].武汉大学学报，2007,32(3)：251-254.

[9] 朱文耀，符养，李彦.GPS高程导出的全球高程振荡运动及季节变化[J].中国科学，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].武汉：中国地震局地质研究所，2011.

[12] 赵斌，聂兆生，黄勇.大规模GPS揭示的华北地区现今垂直运动[J].大地测量与地球动力学，2014，34(5)：35-39.

[13] 黄立人.GPS基准站坐标分量时间序列的噪声特性分析[J].大地测量与地球动力学，2006,26(2)：31-38.