2012, Vol. 55
地震发生后，对其做出快速､准确的评估，对于地震应急､震后减灾和救灾､震后趋势预测及强震引起的海啸预警等至关重要.基于短时速报的地表同震位移为研究震源参数､潜在的断层及地壳形变等事件提供了宝贵的信息.利用地震台网布设的地震仪能够在地震发生的数分钟内较为准确地确定震源及震级，然而对于一些破坏力极大的强震，地震仪容易产生振幅饱和､速度或加速度积分经常出错､噪声放大和扭曲真实信号等现象[1].此外，基于地震仪估计得到的震级较为保守，容易低估地震的大小及其破坏程度.例如2004年12月26日的印度洋海啸引发的苏门答腊地震，利用地震波短时速报的震级(M8.0~8.5)远小于其真实大小(M9.2~9.3)，从而严重影响了快速､准确的地震预警，造成了巨大的人员伤亡[2].近年来，随着GPS 观测精度及处理方法的不断改进与提高，尤其是高频(1 Hz)和超高频(20~50Hz)GPS 技术的出现，GPS 逐渐成为地震仪的重要补充，为地震学研究提供了一种新的途径.
利用高频GPS 进行地表位移监测及地震分析最早于1994年由日本Kyoto大学灾害预防研究所(DPRI)提出，其利用1 Hz的GPS 观测数据，采用事后相对定位方式获得了平面1~2cm 的定位精度，并得出利用GPS 能够有效监测近震､大幅的地表运动等结论[3].澳大利亚新南威尔士大学(UNSW)的Ge等将GPS 天线及地震仪安置在震荡的卡车顶部，利用10~20 Hz的GPS观测数据，采用RTK 方式测试了高频GPS用于监测表面位移的可行性，结果表明由GPS监测的位移与加速度计积分得到的位移两者吻合较好[4-5].此后，Larson等人基于高频GPS 对2002 年DenaliM7.9 级地震，2003年Tokachi-okiM8.1级地震，2003年加州SanSimeon M6.5 级地震，2004 年Sumatra-AndamanM9.3级巨震，2008年汶川M8.0地震进行了研究，获取了地震近场和远场的地表运动时序[6-13].在地震/海啸预警方面，Blewitt 等利用IGS 跟踪站网15min左右､30s 间隔的GPS 观测数据，研究了2004年Sumatra地震的震级､机制及其空间破裂尺度，得出其灵敏度(准确度)受卫星轨道精度及采样率影响等结论[2].然而，基于高频和超高频的GPS地震波震相识别及地震(海啸)预警研究仍然较少.
因此，本文首先利用2011年3 月11 日日本仙台港大地震当天若干测站1Hz的GPS观测数据分析了各站的同震位移，并在此基础上利用测站速度信息，采用S变换提取地震波P 波初至时刻(震相识别)，为强震预警提供支持和依据.
2 高频GPS 地震波获取及同震位移监测选取此次日本大地震近震及远震区的6个IGS跟踪站(MIZU､USUD､DAEJ､SHAO､PETS､GUAM)1Hz的GPS 观测数据，同时利用距东京15km 左右的JA01站1Hz的GPS观测数据，采用武汉大学研制的TriP 软件进行动态精密单点定位解算[14]，卫星轨道及钟差采用COD 中心提供的15min间隔的精密星历与5s间隔的精密卫星钟差，卫星截止高度角设置为5°.各站分布及震中位置(星形)如图 1所示.MIZU 站距震中最近，大约为140km;GUAM 站距震中最远，达到2700余千米.
|
图 1 震中及高频GPS测站分布示意图 Fig. 1 Location of the HR-GPS stations and the epicenter of the M9.0 Tohoku earthquake |
基于高频GPS可以直接获得地表的真实位移，捕获强震时复杂的地表运动时序，为研究断裂扩展机理提供可视化帮助.图 2 自上而下依次为所有测站北､东､垂直方向的地表同震位移时间序列，横轴代表UTC 时间，纵轴代表相应的位移量，为便于描述地震波的传播过程，将所有测站依震中距大小排列显示其震前､震时及震后的位移变化.从图中可以看出地震发生后，距离震中最近的MIZU 站最先感知，随后地震波传至其他各站.各站各分量在地震到时均有不同程度的震动，其震动大小主要取决于震中距及地震带的破裂方向.尤其以震中距较小的MIZU 站(~140 km)､JA01 站(~390 km)及USUD 站(~450km)最为明显，对于数千千米之外的DAEJ站､SHAO站､PETS站和GUAM 站利用高频GPS同样能够监测到其地震波信号.此外，由于破坏力极强的剪切波(S波)主要造成水平向的挤压和扭动，因此，水平向位移较竖向位移变化更加显著.
|
图 2 GPS同震位移时间序列(北/东/垂向) Fig. 2 Time series of coseismic displacement using HR-GPS (North/East/Up) |
尽管远震区(> 600km)的高频GPS能够监测到地震波信号，但其破坏力较小，同震位移容易被观测噪声所淹没，因此本文仅给出近震区测站的震时水平向运动轨迹(图 3所示)及其震后位移变化如表 1所示.在时域内，图 2 中已经详细反映了测站在N､E､U 方向的运动过程，因此，图 3 侧重于从空间域(水平向)描述测站的点位运动过程.
|
图 3 近震测站震时水平向运动轨迹(UTC 05 :46 :00—05 :56 :00) Fig. 3 Horizontal trajectory of near-field stations (during UTC 05 :46 : 00—05 : 56 : 00) |
|
表 1 近震区高频GPS测站同震位移 Table 1 Coseismic permanent displacements of near-field stations |
图 3中坐标原点为各站震前初始位置(Start)，地震波(尤其是横向剪切波)到达时测站开始朝东南方向发生剧烈的抖动，其中MIZU 站的最大位移达到东向3m､南向2m左右;随着地震波的衰减，测站位置发生回弹至最终点位(End).相对于MIZU站，JA01､USUD 站的位移量级较小，将其运动轨迹放大(Enlarge)进行图示，得到图 3中左半部分.表 1说明强震已经造成这些测站产生了永久性位移.其中MIZU 站向南偏移1.133m, 向东偏移2.130m, 垂直方向也向下沉降0.1m 左右.其他测站也呈向东漂移和向下沉降的趋势，这在一定程度上也反映强震造成日本所在的板块整体朝东运动，并略有沉降.
3 地震波震相识别方法地震预警的基本思想是利用地震P 波传播速度快于S 波与面波的特点，通过分析P 波初到时刻､振幅及频率等信息，估算随后而至的破坏性震动的大小及范围，进而确定预警范围与级别.尽管预警时间非常短暂(数秒至数十秒)，但却具有极其重要的社会效益.因此基于高频GPS的强震预警及震源参数反演首先需要解决的是初至波的拾取，即P波的震相识别.常用的震相识别方法有:人工经验法､STA/LTA 法､傅里叶变换法､小波分析法､AIC 方法､瞬时频率法等.
S变换(S-Transform)是由美国地球物理学家Stockwell等在1996年提出的一种无损可逆的线性时频分析方法，它综合了短时傅里叶变换和小波变换的优点，具有较高的时频分辨率，在信号处理领域得到了广泛应用[15].由于地震波初至前一般表现为随机噪声(背景噪声)，初至后则是地震信号与背景噪声的混合.因此，地震波初动时将形成振幅能量的突变.当信噪比较高时，在时域上很容易识别;当信噪比较低时，初至前后能量差异小，在时域难以识别，但在时频域内，由于地震波初至信号与背景噪声占用不同的时频区域，使得初至识别仍然可行.
连续时间序列函数h(t)的S变换公式如下:
|
(1) |
其中，t表示时间，f表示频率，τ 为高斯窗函数的平滑因子.对于1维的函数h(t)，变量τ 和其中一个频率fi可确定一个S(τ，fi)，我们称S(τ，fi)所确定的值为一道噪声.因此，当频率取不同值时，S 变换可以获得若干道噪声向量.
本文首先利用MIZU 和JA01 两个GPS 测站的速度信息，采用S 变换提取P 波初至时刻，并结合邻近地震仪的数据(K-NET 网，IWT011 站，距MIZU 站约2km)､美国地质调查局(USGS)提供的参考数据进行对比分析，最后利用拾取的地震波初至时刻与震中距反演出P波的传播速度.
4 震相识别结果及分析 4.1 背景噪声时频分析以MIZU 站为例，选取其震前900 个历元(15min)的位移序列进行一阶差分消除长周期低频噪声，得到测站速度波形，采用S变换将其速度波形转换至时频域进行频谱分析，估计背景噪声的频率范围及其能量大小，得到该站速度背景噪声的S变换谱，如图 4所示，其中横轴为历元数，纵轴为归一化频率.
|
图 4 MIZU站震前速度波形S变换谱 Fig. 4 S-transform spectrum of station MIZU waveform before earthquake |
图 4显示，震前速度波形表现为高频噪声，其频率主要分布在0.2~0.5 Hz之间，且竖向背景噪声的能量较水平向明显偏大，这也与GPS测量垂直方向较水平方向的观测噪声较大相符.北､东､垂直方向基底噪声的平均能量分别为0.001､0.0015 和0.002左右，这些数值将作为后续地震波P 波初至拾取的能量阈值.
4.2 GPS地震波P波初至拾取基于S变换去噪后的MIZU､JA01测站速度时频谱如图 5､6所示，图中从左向右､自上而下依次代表北､东､垂直三个方向对应的频谱图和速度时序图，横轴为UTC 时，纵轴为归一化频率f或速度V值.
|
图 5 MIZU站震时频谱图 Fig. 5 S-transform spectrum of station MIZU |
|
图 6 JA01站震时频谱图 Fig. 6 S-transform spectrum of station JA01 |
由图 5､6所知，地震波初至时，N､E､U 三个方向在时频域内均发生了较大的能量突变，且随着后续S波､面波的相继到达，能量逐渐增强至最大，并随着地震波的衰弱逐渐消退.基于S变换的P 波拾取结果如表 2所示，受GPS采样率(1Hz)的限制，P波拾取精度仅能精确到秒级.由N､E､U 三个方向拾取的P 波初至时刻不尽相同，差异达到数秒，且水平向较垂直方向的拾取时刻偏早.
|
|
表 2 各测站不同拾取类别对应的P波初到时刻 Table 2 Arrival time of P-wave of station MIZU，JA01 and IWT011 using S-transform and the IASP91 earth model (USGS method) |
为了验证GPS地震波P波拾取的准确性，采用日本K-NET 网络中与MIZU 站相近的IWT011站100Hz的地震仪数据对初到时间进行P波拾取.类似地，采用S变换获得了该站的震前､震时频谱图，如图 7所示.
|
图 7 IWT011站震时频谱图(A为加速度) Fig. 7 S-transform Spectrum of Seismograph IWT011 |
图 7 中从左向右､自上而下依次为N､E､U三个方向对应的时频谱图和加速度时序图.双向竖线表示P波初到时刻(具体时刻如表 2 所示，精确到0.01s)，且P波的归一化频率在各方向上具有较好的一致性(0.02~0.03Hz).
4.4 拾取结果验证及误差分析为验证本文所提方法的正确性及有效性，结合USGS网站提供的P 波初到时刻作为参考值，对比分析各站的拾取精度，结果如表 2所示.比较表 2中各站S 变换拾取的P 波初到时刻与USGS 提供的参考值发现，几乎所有测站拾取的结果均比USGS提供的结果晚5~8s 左右，对于高频GPS 测站(MIZU､JA01)，其水平向拾取精度明显优于高程方向.比较MIZU(高频GPS)与IWT011(地震仪)这两个相邻的测站，发现其地震波平均到时也相差2~3s.究其原因，这是因为:其一，受GPS本身观测精度的限制，基于GPS技术获取的地震波信噪比不及地震仪数据，导致使用S 变换不容易区分背景噪声与微弱的地震信号(尤其是垂直方向);而作为区分背景噪声与地震信号的能量阈值选取极其关键，较小的阈值差异可能引起较大的拾取误差(达到数秒);且不同方向对地震波初至的敏感程度不同，对其简单的求取平均到时也会引起拾取结果产生一定的偏差.其二，尽管已将MIZU 与IWT011 的时间归化至UTC 时，但时钟本身的误差导致两者并非严格同步，因此基于GPS和地震仪的拾取结果也会有所差异.此外，由于USGS 提供的参考值是基于IASP91地球模型计算得到的理论值，它并未考虑椭球改正及严密的介质传播速度等信息，可能会与实际情况存在较大的差异[16].因此，本文采用S变换拾取的地震波初至与USGS 提供的参考结果也并非完全一致.但总体而言，基于GPS 与地震仪数据在震相识别时具有较好的吻合性.
基于上述拾取的P 波，根据各站的震中距反演得到MIZU､JA01 两站的P 波传播速度分别为5.06km/s和5.82km/s.根据P波在不同的介质对应的传播速度一般为5~8km/s, 因此，这也在一定程度上验证了本文基于S 变换拾取的震相为P波.
5 结论本文以日本仙台大地震为例，首先利用近震区和远震区若干测站的高频GPS 观测数据采用动态精密单点定位技术获取了高频GPS地震波，分析了强震时复杂的地表运动时序，获得了地表的同震位移，结果表明:地震到时各站均有不同程度的震动，其震动大小主要取决于震中距及地震带的破裂方向，即使对于数千千米之外的DAEJ站､SHAO 站､PETS站和GUAM 站利用高频GPS同样能够监测到其地震波信号;由于破坏力极强的剪切波(S 波)主要造成水平向的挤压和扭动，导致水平向位移较垂直向位移变化更加剧烈;强震造成近震区内的部分测站已产生永久性变形，日本所在的板块整体朝东运动，并略有下沉.
提出利用S变换进行高频GPS 地震波的震相识别，基于MIZU 和JA01两个GPS测站的速度信息，采用S 变换提取了P 波初至时刻，并与USGS提供的参考值进行对比分析，同时选择距MIZU2km左右的IWT011站的地震仪数据进行S变换，将其P波拾取结果与USGS 提供的参考值进行对比分析.结果表明:受GPS 测量精度及采样率的限制，同一测站N､E､U 三个方向拾取的P 波初至时刻差异达到数秒，且水平向较垂直方向的拾取时刻早;不同测站基于S变换拾取得到的结果与USGS提供的参考值差异达到5~8s, 相邻GPS测站与地震仪台站的地震波初至时刻吻合较好，但仍有2~3s的差异.由地震波初到时刻及震中距反演得到的P波传播速度为5~6km/s, 与实际较为符合.
今后，随着高频和超高频GPS 技术的发展，尤其是高程方向定位精度的提高对于改善高频GPS地震波信噪比､提高震相识别精度，进而为地震应急､预警提供更加科学的决策依据具有十分重要的意义，也是今后急需解决的问题.
致谢衷心感谢IGS､USGS､K-NET 网络以及DLR 的Thomas Dautermann为本文提供了丰富的数据及产品服务;在此一并表示感谢!
| [1] | Boore D M, Stephens C D, Joyner W B. Comments on baseline correction of digital strong-motion data: examples from the 1999 Hector Mine, California, earthquake. Bull. Seismol. Soc. Am. , 2002, 92(4): 1543-1560. DOI:10.1785/0120000926 |
| [2] | Blewitt G, Kreemer C, Hammond W C, et al. Rapid determination of earthquake magnitude using GPS for tsunami warning systems. Geophys. Res. Lett. , 2006, 33: L11309. DOI:10.1029/2006GL026145 |
| [3] | Hirahara K, Nakano T, Hoso Y, et al. An experiment for GPS strain seismometer. // Proc. Japanese Symposium on GPS. Tokyo, Japan, 1994: 67-75. |
| [4] | Ge L. GPS seismometer and its signal extraction. // Proc. 12th Int. Tech. Meeting of the Satellite Division of the U. S. Inst. of Navigation (ION GPS '99). Nashville, Tennessee, 1999. |
| [5] | Ge L L, Han S W, Rizos C, et al. GPS seismometers with up to 20 Hz sampling rate. Earth. Planets and Space , 2000, 52(10): 881-884. DOI:10.1186/BF03352300 |
| [6] | Larson K M, Bodin P, Gomberg J. Using 1-Hz GPS data to measure deformations caused by the Denali fault earthquake. Science , 2003, 300(5624): 1421-1424. DOI:10.1126/science.1084531 |
| [7] | Bock Y, Prawirodirdjo L, Melbourne T I. Detection of arbitrarily large dynamic ground motions with a dense high-rate GPS network. Geophys. Res. Lett. , 2004, 31: L06604. DOI:10.1029/2003GL019150 |
| [8] | Irwan M, Kimata F, Hirahara K, et al. Measuring ground deformations with 1-Hz GPS data: the 2003 Tokachi-oki earthquake (preliminary report). Earth, Planets and Space , 2004, 56(3): 389-393. DOI:10.1186/BF03353070 |
| [9] | Miyazaki S, Larson K M, Choi K, et al. Modeling the rupture process of the 2003 September 25 Tokachi-Oki (Hokkaido) earthquake using 1-Hz GPS data. Geophys. Res. Lett. , 2004, 31: L21603. DOI:10.1029/2004GL021457 |
| [10] | Ji C, Larson K M, Tan Y, et al. Slip history of the 2003 San Simeon earthquake constrained by combining 1-Hz GPS, strong motion, and teleseismic data. Geophys. Res. Lett. , 2004, 31: L17608. DOI:10.1029/2004GL020448 |
| [11] | Ohta Y, Meilano I, Sagiya T, et al. Large surface wave of the 2004 Sumatra-Andaman earthquake captured by the very long baseline kinematic analysis of 1-Hz GPS data. Earth, Planets and Space , 2006, 58(2): 153-157. DOI:10.1186/BF03353372 |
| [12] | Shi C, Lou Y, Zhang H, et al. Seismic deformation of the Mw8.0 Wenchuan earthquake from high-rate GPS observations. Adv. Space Res. , 2010, 46(2): 228-235. DOI:10.1016/j.asr.2010.03.006 |
| [13] | 殷海涛, 甘卫军, 肖根如, 等. 利用高频GPS技术进行强震地面运动监测的研究进展. 地球物理学进展 , 2009, 24(6): 2012–2019. Yin H T, Gan W J, Xiao G R, et al. Progress on monitoring strong earthquake ground motions using high-rate GPS. Progress in Geophysics (in Chinese) , 2009, 24(6): 2012-2019. |
| [14] | 张小红, 刘经南, ForsbergR. 基于精密单点定位技术的航空测量应用实践. 武汉大学学报 (信息科学版) , 2006, 31(1): 19–22. Zhang X H, Liu J N, Forsberg R. Application of precise point positioning in airborne survey. Geomatics and Information Science of Wuhan University (in Chinese) , 2006, 31(1): 19-22. |
| [15] | Stockwell R G, Mansinha L, Lowe R P. Localization of the complex spectrum: The S transform. IEEE Transactions on Signal Processing , 1996, 44(4): 998-1001. DOI:10.1109/78.492555 |
| [16] | Kennett B L N, Engdahl E R. Travel times for global earthquake location and phase identification. Geophys. J. Int. , 1991, 105(2): 429-465. DOI:10.1111/gji.1991.105.issue-2 |