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

1、1实时高频 GPS 在地震学中的应用研究摘要：以 TrackRT 软件为例介绍了实时高频 GPS 技术的数据处理方法，通过与测震（强震）的比较，认为该技术可获取可靠的长周期位移，对测震（强震）资料进行有效弥补。并系统地阐述实时高频 GPS 技术在地震学中的应用方法，表明在我国开展实时高频 GPS 方面的研究具有十分重要的实际意义。 关键词：实时高频 GPS；强震；地震预警；地表位移 中图分类号：P315727 文献标识码：A 文章编号：1000-0666（2013）03-0330-07 0 引言 GPS 技术作为地学研究的一个有效途径，具有全天候、高精度、自动化、高效益等显著特点，在地壳运动监

2、测和地球动力学等领域发挥着越来越重要的作用（Wang et al，2001；Gan et al，2007；Zhang，Gan，2008） 。目前 GPS 通常利用日均值来分析和研究长周期的地壳形变信息，可为中长期地震预测提供可靠的依据，但它不能捕捉内地壳瞬时形变信息。 随着接收机技术和存储能力的提高，我们可以获取高频（1 Hz）和超高频（2050 Hz）GPS 观测数据，而单历元 GPS 处理技术的逐渐成熟，也使高频 GPS 有能力观测到周期大于 1 s 的真实地表位移。目前国内外学者已经利用高频 GPS 观测网成功的获取了几次大地震所引起的瞬时地2表运动状态，如 2003 年加州 San S

3、imeon 65 级地震（Ji et al，2004） ，2004 年 SumatraAndaman 93 级巨震（Ohta et al，2006） ，2008 年汶川80 级地震（Yin et al，2010）和 2011 年日本 90 级地震（Grapenthin，Freymueller，2011）等，这为我们更好地理解地震破裂过程以及地震波在地表的传播方式开辟了新的途径。 通讯技术的飞速发展使高频 GPS 数据的实时传输成为了可能，可获取实时地表位移状态，弥补现有技术手段的不足。当前实时高频 GPS 技术成为国际研究的新热点（Allen，Ziv，2011） 。本文介绍了实时高频GPS 技

4、术及其数据处理方法，并与传统地震（强震）仪进行了对比分析，论述其如何应用于地震预警等实时地震学领域。 1 实时高频 GPS 及处理方法 GPS 全面建成至今只有不到 20 年的时间，但却广泛地应用于导航、大地测量、摄影测量、地震学等多个领域。随着技术和方法的不断发展，其定位精度和采样频率有了明显提高，目前 GPS 定位精度可达亚毫米级，相对定位精度可达 10-9。而数据采样率从 30 s 提升到了 0021 s，甚至 001 s，极大拓展了观测的时间尺度。实时高频 GPS 技术（简称 RH GPS）是目前 GPS 发展的最新阶段，主要是通过接收机和 4 颗以上的卫星来实时获取高频数据（采样率高

5、于 1 s） ，并通过通讯网络和处理中心实时解算出地面点亚厘米级的位移。其数据处理方法主要有差分模式和非差模式两种，差分模式必须要选择一个参考站点，然后对其他站点进行相对定位，非差模式主要依赖于高精度的卫星轨道信息。两种方法在后处理情况下精度相当，但在高精度的实时应用中，非差模式由于无法及3时获取可靠的轨道信息而受限。 目前国际上的 GAMIT/GLOBK、RTD 和 GIPSY 软件均可做实时高频 GPS数据处理，前两种软件为差分模式定位，GIPSY 为非差模式定位。本文主要详细阐述 GAMIT/GLOBK 的实时处理方法，该软件中的 TrackRT 模块是由麻省理工学院的 Thomas H

6、erring 教授于 2010 年开发出来的，主要用来进行实时 GPS 数据处理。该模块可直接得出站点三维位移时间序列和对流层延迟，其水平位移精度优于 1 cm，垂直方向精度优于 2 cm。目前该软件已应用于 UNAVCO，USGS，PBO 等 GPS 观测网的实时解算和发布。虽然实时高频 GPS 数据的解算方法与高频 GPS 类似，但在其软件安装、数据获取、卫星轨道的选择等方面有很大差异。 1 地震研究 136 卷第 3 期殷海涛等：实时高频 GPS 在地震学中的应用研究 11 实时数据获取 在安装、运行 TrackRT 模块前，必须要安装 BKG NTRIP Client（BNC）和 QT

7、 软件的 libraries 和 include 文件。BNC 软件主要是利用互联网协议和客户端来传输 RTCM 数据，可通过互联网端口为trackRT 模块提供 GPS 原始数据。BNC 软件除了可以实时传输 GPS 原始数据外，还具有将原始数据转换为 RINEX 格式数据、存储星历文件、发布数据、精密单点定位等功能，可通过http：/igsbkgbundde/ntrip/download 下载。 12 卫星星历文件 实时 GPS 数据处理时所需的卫星轨道文件也必须满足实时的要求。目前国际 GPS 服务机构（IGS）提供了 3 种卫星星历：最终星历（IGF） 、4快速星历（IGR）和预报星历

8、（IGU） 。对于 IGS 最终星历来说，精度虽然很高，但时间延迟 1218 d，这限制了实时或准实时用户的使用。IGS提供的快速精密星历（时间延迟 1741 h） ，也不能满足实时的要求。为了实现实时应用，IGS 机构于 2000 年 3 月 5 日（GPS 周 1 052 周）提供了超快星历。超快星历和超快卫星钟差的更新率为 6 h，在每天UTC（Coordinate Universal Time）3：00、9：00、15：00 和 21：00 各发布一次。超快星历的历元间隔为 15 min，共包括 48 h 的轨道信息，前24 h 是基于跟踪站的观测值计算得到，后 24 h 是外推预报得

9、到的。实时GPS 数据处理只能利用后 24 h 外推结果来进行计算。除了 IGS 发布的精密星历外，观测数据本身也可实时获取导航星历，但其精度较低，不建议使用。表 1 对现有 GPS 卫星星历文件进行了比较，发现只有预报星历中外推的 24 h 结果在精度和时间延迟上满足实时高频 GPS 数据处理的要求。 13 实时处理方法 实时处理 GPS 数据时，TrackRT 模块中的部分命令与 Track 模块相同，需要在命令文件中设置一表 1GPS 星历文件对比 Tab1Comparison of GPS satellite ephemeris files 星历名称 1 轨道误差/cm1 卫星钟差 1

10、 时间延迟 1 更新时间导航星历1100 cm15 ns RMS 25 ns SDev1 实时 1 实时预报星历 （外推的 24 h）15 cm13 ns RMS 15 ns SDev1 实时 1UTC 3：00，9：00， 515：00，21：00 预报星历 （观测的 24 h）13 cm1150 ps RMS 25 ps SDev139 h1UTC 3：00，9：00， 15：00，21：00 快速星历 125 cm175 ps RMS 25 ps SDev11741 h1UTC 17：00 最终星历 125 cm175 ps RMS 20 ps SDev11218 d1 每周四 个参考站

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

12、甚宽相（EX-WL）和消电离层浮点估计（LC）来计算整周模糊度。设 L1 和 L2 的模糊度整周数为N1，N2，则： （1）MW-WL 是基于相位和距离数据来估计 N1-N2； （2）EX-WL=N1-（f1/f2）*N2 是 L1 的整周数，但 L2 的周数变为了1283*N2； （3）LC=2546*N1-1984*N2。 6EX-WL 不受几何距离变化的影响，但是依赖于电离层延迟。对于短基线来说，EX-WL 应该接近于 0 才符合 N1 和 N2 的选择。当 N1 和 N2 正确的时，LC 残差也应该接近于 0。 在定位精度方面，目前实时高频 GPS 的定位精度在亚厘米级，这表示它可以为

13、大地震的震级估计提供约束。在解算时间方面，实时高频 GPS的数据处理在理论上可达到实时解算，但在实际操作过程中，其数据采集、传输和处理受到通讯设备、处理机配置等硬件因素的影响，时间延迟约为 25 s，可为地震预警等实时地震学应用提供帮助。 2 与地震（强震）仪的比较 震时地表位移是估计地震破裂过程和地震震级的一项重要参量，传统的方法是对地震仪或强震仪所测速度（加速度）进行积分得到（殷海涛等，2009） 。虽然地震仪（强震仪）和高频 GPS 都可以获得震时的地表位移信息，但两种技术在观测方法上有着本质的区别：惯性地震仪的初始状态是平静的（不工作的） ，当地面震动时则产生测量值，震动结束后恢复平静

14、。数字地震仪提供的是一系列速度或加速度观测值，其值是在一个惯性框架下产生的，然后再经过积分（二次积分）得到站点位移。相对而言，GPS 记录的是由 GPS 卫星来确定的 GPS 天线的位置，其位置时间序列是在地球参考框架下计算出来的。表 2 列举了两种技术在性能上的差异（殷海涛等，2012） 。表 2 高频 GPS 与地震仪（强震仪）性能对比 Tab2Performance comparison of high rate GPS and seismometer（strong motion seismography） 性能 1 地震（强震）仪 1 高频 GPS 参考框架 1 惯性坐标系 1 地球参

15、7考框架观测量 1 测量速度或加速度， 通过积分可得到位移 1 观测卫星与接收机间的距离， 直接获得地表点的位移振幅 1 会出现饱和与限幅的情况 1 振幅不受限适用范围 1 对微震，远震都非常敏感 1 仅对大震和近震敏感 通过分析对比高频 GPS 与地震（强震）仪的工作原理及性能，可以发现高频 GPS 技术可以为传统地震学技术进行有效补充： （1）位移时间序列是 GPS 的直接观测量，但是对于地震（强震）仪来说，其观测量为速度（加速度） ，所以必须要对其进行积分处理才能得到位移，而积分处理经常伴随着传感器旋转、倾斜、滞后造成的漂移和积分过程中出现的不确定性，很容易产生放大噪声和扭曲真实信号的情

16、况（Boore et al，2002） ； （2）当大地震发生时，地震仪能够产生振幅饱和现象，为了不记录到满幅的速度和加速度，采取限幅的方法，而 GPS 在振幅方面不会产生饱和，它没有仪器响应来限制接收机的观测能力，相对而言，位移越大其定位相对精度会越高，所以 GPS 技术更适用于大震造成的地表形变； （3）地震仪只有在地面震动后才能产生观测值，但很多地震事件如：断层滑移，岩石粘性变化，流性变化和震后形变等，虽然产生了地壳形变，但是并未产生地震波。在火山活动和非构造性变形过程等，由于形变的时间尺度较大而未必产生地震波。利用高频 GPS 监测这些事件，可以较好地获得低频形变信息，进而了解慢形变过

17、程与机理； （4）目前的地震（强震）仪在大震中很难提取出长周期位移。而Allen 和 Ziv（2011） ，Bock 等（2011）认为：实时 GPS 定位得到的最稳8定的信息就是水平位移，可有效弥补传统地震技术的不足，为大地震的震级估计和地震破裂过程研究提供非常重要的信息。 虽然高频 GPS 在地震学应用上具有很多的优势，在大地震发生时可以记录到明显的位移等，但它也存在一些限制因素：（1）相对于地震仪而言，高频 GPS 的噪声基底较大，Ji 等（2004）分析认为，1 Hz GPS 数据计算的 EW 向，NS 向和垂直方向位置标准差可达到 3 mm，7 mm 和 11 mm。Langbein

18、 和 Bock（2004） ，Bock 等（2011）研究认为 1 Hz GPS 观测数据可以在 99%的置信水平，估计振幅超过 6 mm 的水平位移，因此目前高频 GPS 还不足以记录小震所造成的微弱地表形变。 （2）目前 GPS 接收机仅能获得 1 Hz 或者超高频（2050 Hz）的数据，但是对于更高频（约 200 Hz）的信号，它最终还是要受限于接收机的频带宽度。 因此，结合 GPS 和地震仪资料各自发挥其优势，可起到相互补充的作用，目前 GPS 观测技术与地震学的观测谱范围逐渐合并，逐步发展成了一门新兴的学科GPS 地震学（Nikolaidis et al，2001） ，将会把传统地

19、震学向前推进一大步。 3 实时地震学应用 目前地震预警系统（Earthquake Early Warning，简称 EEW）等实时地震学应用大都只基于地震仪观测网，利用 P 波检测来估计震级，这对快速估计中小地震震级是非常有效的。但目前有可靠的证据证明（Yamada et al，2008） ，在利用前几秒的 P 波信息估计震级时，当 M7地震发生时，震级估计会出现饱和。也就是传统的 P 波检测方法无法区分所探测到的地震是 68 级还是 78 级，虽然说这并不影响地震预警的发9布，但是以 2011 年 3 月 11 日日本的 90 级巨大地震为例，就是因为低估了震级，才导致了对地震和海啸没有做出

20、相应的防御措施，因此无论从科学研究还是从实际应用的角度来讲，都需要继续开展这方面的研究。 31 高频 GPS 与强震仪的结合 地震震级可以通过震时位移的振幅反演来进行估计，在获取震时地表位移的方法中，只有有效结合地震仪和高频 GPS 才能得到可靠的高精度强震动地表位移。Bock 等（2011）引入了连续 Kalman 滤波（卡尔曼滤波）和最优平滑滤波等方法来合并强震仪与高频 GPS 结果。 卡尔曼滤波器包括两个主要过程：预估与校正。预估过程主要利用时间更新方程建立对当前状态的先验估计，及时向前推算当前状态变量和误差协方差估计的值，以便为下一个时间状态构造先验估计值；校正过程负责反馈，利用测量更

21、新方程在预估过程的先验估计值及当前测量变量的基础上建立起对当前状态改进的后验估计。 我们假定每个站点的坐标方向（N，E，U）都是独立的一维运动，将每个方向的运动作为一阶线性差分等式，利用前向卡尔曼滤波来融合 GPS和强震仪数据并监测状态变量。在每个时间步长中，每个坐标方向的位移和速度可表示为初始值：X0=0； P0=I， （1）时间更新（先验估计）为P-k+1=AkPkATk+Qk， X-k+1=AkXk+BkUk， （2）测量值更新（后验估计）为 Pk+1=（P-k+1）-1+HTk+1R-1k+1Hk+1-1， Xk+1=X-k+1+Pk+1HTk+1R-1k+1（Zk+1-Hk+1X-k

22、+1） （3）其中，X为状态变量，k 为时间历元，P-为先验估计误差协方差矩阵，P 为后验估10计误差协方差矩阵；A、B 和 H 为状态变换矩阵，是状态变换过程中的调整系统，数值是从建立的系统数学模型中导出来的；Q 为过程激励噪声协方差矩阵，R 为观测噪声协方差矩阵；观测变量 Z 表示 GPS 位移，U 为系统控制输入，表示强震仪数据。 强震仪时间序列提供了系统输入，其采样率通常为 80250 Hz；而GPS 位移控制测量过程，其采样率通常为 15 Hz。这样式（3）在每一个时间步长内进行时间和测量值更新，如果强震仪采样率是 GPS 采样率的整数倍，则可以简化公式。如果采样频率自始至终保持不变

23、，其噪声特性也不改变，则式（3）中的各项参数 A，B，Q 和 R 均保持不变。此方法的特点是矩阵的维数较少，可以使数值计算更加简洁。此外，卡尔曼滤波只需要当前的数据，这样可以实现实时应用。 应用前向卡尔曼滤波可以解决预测问题，但如果已经存在了一部分数据，在计算中会用到所有存在的和未来的数据，在任意点的估计值可以进行优化。虽然平滑滤波通常应用于非实时处理，但可以通过约束平滑的时间间隔来实现近实时的处理。固定时间间隔的平滑滤波包括 3 步：（1）前向卡尔曼滤波；（2）后向卡尔曼滤波，将时间顺序转换到整个时间序列中。前两步提供的分散信息流在第三步中结合，并计算出最终状态值。结果如图 1 所示，通过滤波处理，两种技术手段可以进行有效结合，获取更精确的宽相位移，而且也减少了处理过程中的主观因素。 32 地震矩震级估计 快速的矩张量（Moment Tensor）估计为地震响应、海啸预警等提供了有价值的信息。目前所有的方法基本都使用地震仪资料，利用时间域