实时高频GPS在地震学中的应用研究*

殷海涛，刘希强，甘卫军

0 引言

GPS技术作为地学研究的一个有效途径，具有全天候、高精度、自动化、高效益等显著特点，在地壳运动监测和地球动力学等领域发挥着越来越重要的作用 (Wang et al，2001;Gan et al，2007;Zhang，Gan，2008)。目前GPS通常利用日均值来分析和研究长周期的地壳形变信息，可为中长期地震预测提供可靠的依据，但它不能捕捉内地壳瞬时形变信息。

随着接收机技术和存储能力的提高，我们可以获取高频 (1 Hz)和超高频 (20～50 Hz)GPS观测数据，而单历元GPS处理技术的逐渐成熟，也使高频GPS有能力观测到周期大于1 s的真实地表位移。目前国内外学者已经利用高频GPS观测网成功的获取了几次大地震所引起的瞬时地表运动状态，如2003年加州San Simeon 6.5级地震 (Ji et al，2004)，2004年 Sumatra-Andaman 9.3级巨震 (Ohta et al，2006)，2008年汶川8.0级地震(Yin et al，2010)和 2011年日本 9.0级地震(Grapenthin，Freymueller，2011)等，这为我们更好地理解地震破裂过程以及地震波在地表的传播方式开辟了新的途径。

通讯技术的飞速发展使高频GPS数据的实时传输成为了可能，可获取实时地表位移状态，弥补现有技术手段的不足。当前实时高频GPS技术成为国际研究的新热点 (Allen，Ziv，2011)。本文介绍了实时高频GPS技术及其数据处理方法，并与传统地震 (强震)仪进行了对比分析，论述其如何应用于地震预警等实时地震学领域。

1 实时高频GPS及处理方法

GPS全面建成至今只有不到20年的时间，但却广泛地应用于导航、大地测量、摄影测量、地震学等多个领域。随着技术和方法的不断发展，其定位精度和采样频率有了明显提高，目前GPS定位精度可达亚毫米级，相对定位精度可达10－9。而数据采样率从30 s提升到了0.02～1 s，甚至0.01 s，极大拓展了观测的时间尺度。实时高频GPS技术 (简称RH GPS)是目前GPS发展的最新阶段，主要是通过接收机和4颗以上的卫星来实时获取高频数据 (采样率高于1 s)，并通过通讯网络和处理中心实时解算出地面点亚厘米级的位移。其数据处理方法主要有差分模式和非差模式两种，差分模式必须要选择一个参考站点，然后对其他站点进行相对定位，非差模式主要依赖于高精度的卫星轨道信息。两种方法在后处理情况下精度相当，但在高精度的实时应用中，非差模式由于无法及时获取可靠的轨道信息而受限。

目前国际上的GAMIT/GLOBK、RTD和GIPSY软件均可做实时高频GPS数据处理，前两种软件为差分模式定位，GIPSY为非差模式定位。本文主要详细阐述GAMIT/GLOBK的实时处理方法，该软件中的TrackRT模块是由麻省理工学院的Thom-as Herring教授于2010年开发出来的，主要用来进行实时GPS数据处理。该模块可直接得出站点三维位移时间序列和对流层延迟，其水平位移精度优于1 cm，垂直方向精度优于2 cm。目前该软件已应用于UNAVCO，USGS，PBO等GPS观测网的实时解算和发布。虽然实时高频GPS数据的解算方法与高频GPS类似，但在其软件安装、数据获取、卫星轨道的选择等方面有很大差异。

1.1 实时数据获取

在安装、运行TrackRT模块前，必须要安装BKG NTRIP Client(BNC)和QT软件的libraries和include文件。BNC软件主要是利用互联网协议和客户端来传输RTCM数据，可通过互联网端口为trackRT模块提供GPS原始数据。BNC软件除了可以实时传输GPS原始数据外，还具有将原始数据转换为RINEX格式数据、存储星历文件、发布数据、精密单点定位等功能，可通过 http://igs.bkg.bund.de/ntrip/download下载。

1.2 卫星星历文件

实时GPS数据处理时所需的卫星轨道文件也必须满足实时的要求。目前国际GPS服务机构(IGS)提供了3种卫星星历:最终星历 (IGF)、快速星历 (IGR)和预报星历 (IGU)。对于IGS最终星历来说，精度虽然很高，但时间延迟12～18 d，这限制了实时或准实时用户的使用。IGS提供的快速精密星历 (时间延迟17～41 h)，也不能满足实时的要求。为了实现实时应用，IGS机构于2000年3月5日 (GPS周1 052周)提供了超快星历。超快星历和超快卫星钟差的更新率为6 h，在每天 UTC(CoordinateUniversalTime)3:00、9:00、15:00和21:00各发布一次。超快星历的历元间隔为15 min，共包括48 h的轨道信息，前24 h是基于跟踪站的观测值计算得到，后24 h是外推预报得到的。实时GPS数据处理只能利用后24 h外推结果来进行计算。除了IGS发布的精密星历外，观测数据本身也可实时获取导航星历，但其精度较低，不建议使用。表1对现有GPS卫星星历文件进行了比较，发现只有预报星历中外推的24 h结果在精度和时间延迟上满足实时高频GPS数据处理的要求。

1.3 实时处理方法

实时处理GPS数据时，TrackRT模块中的部分命令与Track模块相同，需要在命令文件中设置一个参考站，然后进行差分处理，但由于实时数据流中并不包括站点的位置坐标，所以必须要预先给定站点坐标值。在命令文件中还需利用SITE_STATS和ATM_STATS命令对站点的坐标噪声和大气延迟进行先验约束，确保其解算精度。参考站的选择方法可参考殷海涛等 (2012)，为了使参考站“固定不动”，通常将参考站的这两组值设置为0。在命令文件中还需设定站点的天线和接收机型号，以及更新数据码偏心文件 (DCB)，这样不仅可以提高解算精度，还可减少计算量，提高计算效率。

表1 GPS星历文件对比Tab.1 Comparison of GPS satellite ephemeris files

实时GPS数据处理的核心问题也是模糊度的快速解算，TrackRT模块使用Melbourne-Wubbena宽相 (MW－WL)，甚宽相 (EX－WL)和消电离层浮点估计 (LC)来计算整周模糊度。设L1和L2的模糊度整周数为N1，N2，则:

(1)MW－WL是基于相位和距离数据来估计N1－N2;

(2)EX－WL=N1－ (f1/f2)*N2是L1的整周数，但L2的周数变为了1.283*N2;

(3)LC=2.546*N1－1.984*N2。

EX－WL不受几何距离变化的影响，但是依赖于电离层延迟。对于短基线来说，EX－WL应该接近于0才符合N1和N2的选择。当N1和N2正确的时，LC残差也应该接近于0。

在定位精度方面，目前实时高频GPS的定位精度在亚厘米级，这表示它可以为大地震的震级估计提供约束。在解算时间方面，实时高频GPS的数据处理在理论上可达到实时解算，但在实际操作过程中，其数据采集、传输和处理受到通讯设备、处理机配置等硬件因素的影响，时间延迟约为2～5 s，可为地震预警等实时地震学应用提供帮助。

2 与地震 (强震)仪的比较

震时地表位移是估计地震破裂过程和地震震级的一项重要参量，传统的方法是对地震仪或强震仪所测速度 (加速度)进行积分得到 (殷海涛等，2009)。虽然地震仪 (强震仪)和高频GPS都可以获得震时的地表位移信息，但两种技术在观测方法上有着本质的区别:惯性地震仪的初始状态是平静的 (不工作的)，当地面震动时则产生测量值，震动结束后恢复平静。数字地震仪提供的是一系列速度或加速度观测值，其值是在一个惯性框架下产生的，然后再经过积分 (二次积分)得到站点位移。相对而言，GPS记录的是由GPS卫星来确定的GPS天线的位置，其位置时间序列是在地球参考框架下计算出来的。表2列举了两种技术在性能上的差异 (殷海涛等，2012)。

表2 高频GPS与地震仪 (强震仪)性能对比Tab.2 Performance comparison of high rate GPS and seismometer(strong motion seismography)

通过分析对比高频GPS与地震 (强震)仪的工作原理及性能，可以发现高频GPS技术可以为传统地震学技术进行有效补充:

(1)位移时间序列是GPS的直接观测量，但是对于地震 (强震)仪来说，其观测量为速度(加速度)，所以必须要对其进行积分处理才能得到位移，而积分处理经常伴随着传感器旋转、倾斜、滞后造成的漂移和积分过程中出现的不确定性，很容易产生放大噪声和扭曲真实信号的情况(Boore et al，2002);

(2)当大地震发生时，地震仪能够产生振幅饱和现象，为了不记录到满幅的速度和加速度，采取限幅的方法，而GPS在振幅方面不会产生饱和，它没有仪器响应来限制接收机的观测能力，相对而言，位移越大其定位相对精度会越高，所以GPS技术更适用于大震造成的地表形变;

(3)地震仪只有在地面震动后才能产生观测值，但很多地震事件如:断层滑移，岩石粘性变化，流性变化和震后形变等，虽然产生了地壳形变，但是并未产生地震波。在火山活动和非构造性变形过程等，由于形变的时间尺度较大而未必产生地震波。利用高频GPS监测这些事件，可以较好地获得低频形变信息，进而了解慢形变过程与机理;

(4)目前的地震 (强震)仪在大震中很难提取出长周期位移。而Allen和Ziv(2011)，Bock等(2011)认为:实时GPS定位得到的最稳定的信息就是水平位移，可有效弥补传统地震技术的不足，为大地震的震级估计和地震破裂过程研究提供非常重要的信息。

虽然高频GPS在地震学应用上具有很多的优势，在大地震发生时可以记录到明显的位移等，但它也存在一些限制因素:(1)相对于地震仪而言，高频GPS的噪声基底较大，Ji等 (2004)分析认为，1～Hz GPS数据计算的EW向，NS向和垂直方向位置标准差可达到3 mm，7 mm和11 mm。Langbein和 Bock(2004)，Bock等 (2011)研究认为1 Hz GPS观测数据可以在99%的置信水平，估计振幅超过6 mm的水平位移，因此目前高频GPS还不足以记录小震所造成的微弱地表形变。(2)目前GPS接收机仅能获得1 Hz或者超高频(20～50 Hz)的数据，但是对于更高频 (约200 Hz)的信号，它最终还是要受限于接收机的频带宽度。

因此，结合GPS和地震仪资料各自发挥其优势，可起到相互补充的作用，目前GPS观测技术与地震学的观测谱范围逐渐合并，逐步发展成了一门新兴的学科——GPS地震学 (Nikolaidis et al，2001)，将会把传统地震学向前推进一大步。

3 实时地震学应用

目前地震预警系统 (Earthquake Early Warning，简称EEW)等实时地震学应用大都只基于地震仪观测网，利用P波检测来估计震级，这对快速估计中小地震震级是非常有效的。但目前有可靠的证据证明 (Yamada et al，2008)，在利用前几秒的P波信息估计震级时，当M＞7地震发生时，震级估计会出现饱和。也就是传统的P波检测方法无法区分所探测到的地震是6.8级还是7.8级，虽然说这并不影响地震预警的发布，但是以2011年3月11日日本的9.0级巨大地震为例，就是因为低估了震级，才导致了对地震和海啸没有做出相应的防御措施，因此无论从科学研究还是从实际应用的角度来讲，都需要继续开展这方面的研究。

3.1 高频GPS与强震仪的结合

地震震级可以通过震时位移的振幅反演来进行估计，在获取震时地表位移的方法中，只有有效结合地震仪和高频GPS才能得到可靠的高精度强震动地表位移。Bock等 (2011)引入了连续Kalman滤波 (卡尔曼滤波)和最优平滑滤波等方法来合并强震仪与高频GPS结果。

卡尔曼滤波器包括两个主要过程:预估与校正。预估过程主要利用时间更新方程建立对当前状态的先验估计，及时向前推算当前状态变量和误差协方差估计的值，以便为下一个时间状态构造先验估计值;校正过程负责反馈，利用测量更新方程在预估过程的先验估计值及当前测量变量的基础上建立起对当前状态改进的后验估计。

我们假定每个站点的坐标方向 (N，E，U)都是独立的一维运动，将每个方向的运动作为一阶线性差分等式，利用前向卡尔曼滤波来融合GPS和强震仪数据并监测状态变量。在每个时间步长中，每个坐标方向的位移和速度可表示为

时间更新 (先验估计)为

测量值更新 (后验估计)为

其中，X为状态变量，k为时间历元，P－为先验估计误差协方差矩阵，P为后验估计误差协方差矩阵;A、B和H为状态变换矩阵，是状态变换过程中的调整系统，数值是从建立的系统数学模型中导出来的;Q为过程激励噪声协方差矩阵，R为观测噪声协方差矩阵;观测变量Z表示GPS位移，U为系统控制输入，表示强震仪数据。

强震仪时间序列提供了系统输入，其采样率通常为80～250 Hz;而GPS位移控制测量过程，其采样率通常为1～5 Hz。这样式(3)在每一个时间步长内进行时间和测量值更新，如果强震仪采样率是GPS采样率的整数倍，则可以简化公式。如果采样频率自始至终保持不变，其噪声特性也不改变，则式(3)中的各项参数A，B，Q和R均保持不变。此方法的特点是矩阵的维数较少，可以使数值计算更加简洁。此外，卡尔曼滤波只需要当前的数据，这样可以实现实时应用。

应用前向卡尔曼滤波可以解决预测问题，但如果已经存在了一部分数据，在计算中会用到所有存在的和未来的数据，在任意点的估计值可以进行优化。虽然平滑滤波通常应用于非实时处理，但可以通过约束平滑的时间间隔来实现近实时的处理。固定时间间隔的平滑滤波包括3步:(1)前向卡尔曼滤波;(2)后向卡尔曼滤波，将时间顺序转换到整个时间序列中。前两步提供的分散信息流在第三步中结合，并计算出最终状态值。结果如图1所示，通过滤波处理，两种技术手段可以进行有效结合，获取更精确的宽相位移，而且也减少了处理过程中的主观因素。

3.2 地震矩震级估计

快速的矩张量 (Moment Tensor)估计为地震响应、海啸预警等提供了有价值的信息。目前所有的方法基本都使用地震仪资料，利用时间域的波形匹配反演方法来实时估计 (Dreger，2003)，这对中小地震来说非常有效，但对较大地震来说，仍然是研究的热点。

Kanamori和Rivera(2008)发展了W相位反演方法，这使大地震CMT(Centroid moment tensor)的计算前进了一大步，目前该方法正应用于USGS等机构。虽然W相位反演方法非常稳定，但也需要长周期位移资料，而目前的地震 (强震)仪很难在大震中提取出长周期位移，限制了其进一步发展。Diego等 (2012)提出了一种利用实时高频GPS反演地震矩张量的方法fastCMT，该反演方法主要利用地表的同震位移来反演深部的震源参数，主要由以下4部分组成:

(1)在指定区域构建反演节点，其节点格网的构建方法有两种:一是需要预先设定断层面模型，这就要求对该地区的地质条件有很深入的了解;二是不需要已知断层模型，直接根据格网搜索法来进行估计，这对断层复杂的地区非常有帮助，但是必须要有足够密度的GPS观测网才能实现;

图1 卡尔曼滤波前后位移时间序列①汶川地震对XYAN站点观测数据.(a)南北向分量;(b)东西向分量;(c)高程分量Fig.1 Displacement time series before and after Kalman filter(a)NS-component;(b)EW-component;(c)evaluation component

(2)根据地壳速度模型计算格林函数，所用程序为 EDGRN/EDCMP。从Hankel变换得到的部分差分运动等式的结果开始，然后应用Thomson－Haskell传播矩阵将地面和深部的变形进行关联，其简化公式为

结合公式 (4)和 (5)，则得到

(3)高频GPS数据的获取及处理，首先利用移动平均法对获取的高频GPS时间序列进行滤波处理，得到永久位移，然后再对其结果进行加权处理。在反演过程中，就使用了两个权重，一是观测噪声权重，用来提高观测值的可靠性。二是震中距权重，为了避免近场最大同震形变主导整个反演过程。其观测噪声加权公式为

位移分量震前60 s的观测噪声标准差。由于实时GPS位移噪声在几分钟尺度内基本保持稳定，而且Bock等 (2011)利用震动台试验表明震动前和震动过程中的噪声水平并没有明显的改变，所以我们可以设定观测噪声和权重矩阵在整个反演过程中保持不变。

为了避免地面运动衰减产生的影响，我们利用站点震中距作为另一个权重。Dreger(2003)给出

(4)地震矩震级估计，在实时应用过程中，每获取一个历元的观测值就反演一次，产生新的MT值

然后，可以利用 Hanks和 Kanamori(1979)提出来的关系式计算矩震级

其中，M0单位为dyn·cm。

通过上述计算可以实时估计出地震矩震级，并可通过格网搜索法估计地震矩心的位置，利用方差减少法计算反演差值，将方差减少值VR最大的节点视为矩心，公式为

式中，d为GPS站点到节点的距离。

fastCMT方法的特点是可以在没有任何断层模型的情况下实时估计大震矩震级和矩心位置，但是前提条件是近震区装备有足够密度的高频GPS站点。Diego等 (2012)通过2003年日本Tokachioki MW8.3地震和2010年EI Mayor-Cucapah MW7.2地震的实际数据分析认为:在地震发生2～3 min，该方法便可计算出精确的地震矩震级，可为地震快速响应和海啸预警等提供有力支持。

4 讨论与结论

实时高频GPS是目前GPS技术发展的最新阶段，大大拓展了其观测信息频域和应用范围，也将使现有的GPS观测网发挥更大的作用。但由于其研究刚刚起步，在数据处理和地震应用过程中仍存在很多问题亟待解决:

(1)模糊度的快速解算仍然是实时GPS数据处理的核心问题，也是一直以来国际研究的热点问题。除此以外，接收机进行观测时易受外部环境因素干扰，产生突跳或卫星失锁的情况，如何消除突跳和快速收敛也是需要解决的技术难点;

(2)对于差分模式定位来讲，如何选择参考站是主要问题，而非差模式的问题在于如何获取高精度的卫星轨道信息。对于地震预警系统来说，两者各有利弊，差分模式可应用于区域预警，而非差模式更适用于现地预警;

(3)如何消除实时GPS观测中的噪声影响，得到更可靠的位移信息，是目前所要面临的主要问题;

(4)如何更好的结合地 (强)震仪与高频GPS资料，获取稳定性高敏感性高的震时地表位移是下一步应用研究的重点;

(5)在获取地表位移后，如何快速提取地表永久位移是估计地震震级的关键因素，也是以后研究的主要问题。

在地震学应用中，与传统地震 (强震)仪相比，RH GPS可直接测量震时地表位移，增强震级估计 (特别是大地震)的可靠性，避免了传感器出现的旋转、倾斜、滞后造成的漂移和积分过程中出现的不确定因素。而且它可以不受振幅限制，有效的弥补了地震 (强震)仪的不足，但由于其定位精度要低于地震 (强震)仪，对微小地震不敏感，所以结合两种手段才能达到最佳的监测效果。

利用实时高频GPS技术能快速获取震时地表位移，可以计算出精确的地震矩震级可以为地震破裂过程研究、地震预警、海啸预警等地震研究和应用提供一项新的手段，因此在我国及时开展此领域的研究工作是十分必要和迫切的。

殷海涛，甘卫军，肖根如，等.2009.利用高频GPS技术进行强震地面运动监测的研究进展［J］．地球物理学进展，24(6):2 012－2 019.doi:10.3969/j.issn.1004－2903.2009.06.011.

殷海涛，肖根如，张磊，等.2012.TRACK高频GPS定位中震时参考站的选取方法［J］．大地测量与地球动力学，32(4):15－19.

Allen R M，Ziv A.2011.Application of real-time GPS to earthquake early warning［J］．GeophysResLett，38(16):1 － 6，doi:10.1029/2011GL047947.

Bock Y，Diego M，Crowell B W.2011.Real-time strong-motion broadband displacements from collocated GPS and Accelerometers［J］．BSSA，101(6):2 904 －2 925，doi:10.1785/0120110007.

Boore D M，Stephens C D，Joyner W B，et al.2002.Comments on baseline correction of digital strong-motion data:examples from the 1999 Hector Mine，California earthquake［J］．BSSA，92(4):1543 － 1560，doi:10.17850120000926.

Diego M，Bock Y，Crowell B W.2012.Real-time centroid moment tensor determination for large earthquakes from local and regional displacement records［J］．Geophys J Int，188(2):703 － 718，doi:10.1111/j1365－246X.2011.05297.x.

Dreger S D.2003.TDMT_INV:Time domain seismic moment tensor inversion［M］//William H K L，Hiroo K，Paul J，et al.International handbook of earthquake and engineering seismology.American:Academic Press，81B:1 627 －1 631.

Gan W J，Zhang P Z，Zheng K S，et al.2007.Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements［J］．JGR，112(B8):1 －14，doi:10.1029/2005JB004120.

Grapenthin R，Freymueller J T.2011.The dynamics of a seismic wave field:Animation and analysis of kinematic GPS data recorded during the 2011 Tohoku-oki earthquake，Japan［J］．Geophys Res Lett，38(18):1 －5，doi:10.1029/2011GL048405.

Hanks T C，Kanamori H.1979.A moment magnitude scale［J］．JGR，84(B5):2 348－2 350.

Ji C，Larson M，Tan Y，et al.2004.Slip history of the 2003 San Simon earthquake constrained by combining 1 － Hz GPS，strong motion，and teleseismic data［J］．Geophys Res Lett，31(17):1 － 4，doi:10.1029/2004GL020448.

Kanamori H，Rivera L.2008.Source inversion of W phase:speeding up seismic tsunami warning［J］．Geophys J Int，175(1):222 － 238.

Langbein J，Bock Y.2004.High-rate real-time GPS network at Parkfield:Utility for detecting fault slip and seismic displacements.Geophys Res Lett，31(15):1 － 4，doi:10.1029/2003GL019408.

Nikolaidis R，Bock Y，Jonge P J，et al.2001.Seismic wave observation with the Global Positioning System［J］．JGR，106(B10):218 097 －21 916.

Ohta Y，Meiano I，Sagiya T，et al.2006.Large surface wave of the 2004 Sumatra-Andaman earthquake captured by the very long baseline kinematic analysis of 1 － Hz GPS data［J］．Earth Planets Space，58(2):153－157.

Wang Q，Zhang P Z，Jeffrey T F，et al.2001.Present-day crustal deformation in continental China constrained by Global Positioning System measurements［J］．Science，294(5542):574 － 577.

Yamada M，Heaton T，Beck J.2008.Real-time estimation of fault rupture extent using envelopes of acceleration［J］．BSSA，98(2):607－619.

Yin H T，Zhang P Z，Gan W J，et al.2010.Near-field surface movement during the Wenchuan MS8.0 earthquake measured by high-rate GPS［J］．Chinese Sci Bull，55(23):2529 － 2534，doi:10.1007/s11434－010－4026－2.

Zhang P Z，Gan W J.2008.Combined model of rigid-block motion with continuous deformation:Patterns of present-day deformation in Continental China［J］．Geological Society of America Special Paper，444:59 －71，doi:10.1130/2008.2444(04).