传统的地震仪在地震监测和研究中发挥着重要作用^[1-3]。但在大地震发生时，震源附近的地震仪可能会受到地震波冲击导致振幅满格而难以完整地记录地震波信号。同时，大地震造成的近场地表永久性变形及地面倾斜会导致地震仪记录的地震波信号产生扭曲、变形。而高频GPS观测资料具有能够获取宽频带信号、可以直接获取高精度水平位移信息、不会出现限幅现象以及对地面倾斜不敏感等优势。因此，充分利用高频GPS观测资料中记录的地震信息已被越来越多的学者所重视，GPS地震学^[4]随之应运而生。

文献[5]利用1 Hz的GPS观测数据获取了FAIR测站记录到的2002年的Denali地震波形，并与地震仪积分结果比较，两者具有较好的一致性。此后，越来越多的学者将高频GPS技术广泛应用于其他地震的研究，并取得了大量的研究成果^[6-8]。由GPS获得的同震静态永久变形信息(static offset)可用于约束震源参数如震级、震源深度等的反演，地震波信息可作为地震仪、强震仪等观测数据的重要补充，用于反演地震断层破裂过程。但是，震后的高频GPS观测资料包含了由地震引起的静态永久变形和地震波的详细信息。当需要永久性变形信息反演震源参数时，地震波信息会影响准确的获得永久性变形，反之，当需要地震波信息反演断层破裂过程时，永久性变形信息也会对地震波信息的准确性产生影响。因此，实现两者的有效分离是准确获取两种信息进而反演地震参数的关键因素。以往为了从GPS观测资料中获取同震静态永久变形，通常采用震前和震后一段时间的资料进行计算，并需要剔除地震发生或以后一段时间的数据以避免地震波位移影响。文献[9]利用卡尔曼滤波方法估计了2003年日本Tokachi-Oki地震主震和震后静态永久性变形效应，但是在应用该方法时，为了避免动态位移对滤波结果的影响，在估计中必须剔除地震发生时的观测数据，这样不仅增加了数据处理的工作量，而且也会丢失地震波信息。文献[10]利用小波分解重构方法提取高频GPS数据静态永久性变形信息，但是在利用小波分层技术时，无法完整保存地震波到时及振幅等信息，而到时和振幅信息对于准确确定震中和震级等参数具有重要意义。

为充分发挥高频GPS资料的优势，本文旨在研究在快速获取高频GPS数据静态永久变形信息的同时，如何有效分离并完整保存高频GPS资料中记录的地震波信息。本文选用平滑先验信息方法(SPM)^[11-14]来实现地震引起的高频GPS观测结果中静态形变和地震波信息的分离。该方法是一种高通滤波器，对具有周期特性的高频信号最为敏感，能够有效地识别具有周期特性的高频周期信息，从而实现周期信号与非周期信号的有效分离，并且不需要剔除数据。由于地震波具有一定周期或准周期性而表现为高频信息，永久性变形则不具有周期性，正是针对高频GPS资料中地震波与永久性变形的这一特点，本文将SPM方法引入高频GPS数据后处理中，分别获取了静态永久性变形信息和地震波信息。

1 数据处理与方法 1.1 数据及地震位置

本文采用的1 Hz的GPS观测数据由美国地学研究组织UNAVCO网站提供(http://www.unavco.org/data/data.html)。本文选取了近几年的3次地震，地震的详细信息见表 1。本文依照两个原则，每个地震对应选取了两个观测站：①测站到震中的距离；②在地震发生前后，观测站的数据记录是否完整。测站分布与震中位置分布情况如图 1所示。

1.2 高频GPS数据处理

本文采用Bernese 5.2软件进行动态精密单点定位，卫星轨道及钟差采用欧洲定轨中心(CODE)提供的15 min间隔的精密星历与5 s间隔的精密卫星钟差^[15]，并且采用IERS推荐的最新的天线相位中心模型、大气潮汐及海潮模型，卫星截止高度角设置为3°。同时，估计了载波相位模糊度、卫星钟差及接收机钟差。电离层误差通过采用双频电离层无关组合消除，天顶对流层延迟采用WET_GMF模型估计。计算得到的动态位移解仍然受多路径效应等周期性噪声的影响，本文采用改进的恒星日滤波(modified sidereal filtering, MSF)方法^[16]进行处理。对P496和P744测站选用震前2 d(P496和P744测站震前只有2 d的数据)，QSEC、SAJU、MIZU和USUD测站选用震前4 d相同恒星日段的时间序列进行加权平均，得到背景噪声。图 2以MIZU测站N方向为例，对比了分别采用震前2 d和震前4 d的观测数据计算得到的背景噪声。其中绿色曲线分别对应2011-03-07(年积日066)到2011-03-10(年积日069)动态位移时间序列，红色曲线为滑动平均后的结果，消除了高频噪声。蓝色曲线即为所得的背景噪声。两组背景噪声的相关系数达到了90.4%，从图中也可以看出，在震前4 d该台站周围环境稳定，没有受到明显的噪声事件干扰。从图 2结果来看，由于背景噪声只与观测台站质量和周围环境的噪声影响有关，在没有突然噪声的干扰下，选用震前2~4 d的数据处理的结果是基本相同的。从地震当天的动态位移时间序列中扣除背景噪声，扣除后的坐标时间序列基本上不再含有多路径效应等周期性噪声的影响。

1.3 信号分离

1.3.1 SPM方法

对于任意时间序列z，均可表示为z=[z₁ z₂ … z_n]^T。其中，n为采样点个数。z可看作近似稳态的周期分量和非周期趋势项之和

式中，z_stat表示近似稳态的周期分量(即地震波)；z_nonstat表示非周期趋势分量(即静态永久变形)，可表示为

式中，H∈R^n×m为观测矩阵；θ∈R^m为回归参数；m是回归参数的个数；ε为观测误差。求解回归参数的估计值时，常用最小二乘准则来估计，但由于静态永久变形信号为非稳态信号，因此，本文引入了正则化估计方法

式中，λ是正则化参数; D_d为第d次导数操作的离散逼近。显然，式(3)是在普通最小二乘方法的基础上增加了一个条件，属于条件极值问题。极值条件的目的是使式(3)的解在满足‖Hθ-z‖取得极小值的同时，使得‖D_d(Hθ)‖取值也尽可能小。这样就可以在推估线性项Hθ的过程中，加入一些已知的先验信息。式(3)的解可表示为

式中，为非周期分量的估计值。观测矩阵H的选择可以包含时间序列z的一些已知信息，例如使用一些高斯型函数或三角函数，这里选取单位阵I作为观测矩阵，即H=I∈R^n×n，详细参见文献[13]。式(3)中正则化部分可以理解为求解正则化矩阵D_d零空间里的解，二阶差分矩阵的零空间包含了所有的一阶曲线，因此D₂是用来估计非周期项的一个不错的选择^{[11, 13]}。二阶差分矩阵D₂∈R^(n-2)×n形式如下

有了这些特殊选择后，该方法称为平滑先验信息方法(smoothness priors method, SPM)^[11-14]。这样，周期信号可以表示为

1.3.2 频率响应及平滑参数的选取

式(7)可以写为，其中L=I-(I+λ²D₂^TD₂)^-1，相当于一个时变高通滤波器，对于每一个离散时间点，L的频率响应(以其行的傅氏变换形式获得)如图 3(a)所示，本文给出了前半部分的频率响应，后半部分是完全相同的。可以看到，滤波器基本为一常数，仅对起始和结束的信号做了不同处理。其滤波效果仅对数据的起始和末尾有扭曲现象，最大程度的避免截断效应的影响。对于信号z的第一个和最后一个元素，滤波效果是衰减的，因此避免了数据终点的扭曲。平滑参数λ对于频率响应的影响如图 3(b)所示，截止频率随λ增大而减小，图中平滑参数λ分别为1、2、4、10、20、50、300，相应的3 dB截止频率分别是0.189、0.132、0.093、0.059、0.041、0.025 Hz和0.011 Hz(对应1 Hz采样频率)。除了平滑参数λ，频率响应当然还与信号z的采样率有关。因此在给定信号z的情况下，其频率响应的调整仅依赖于一个平滑参数λ，受人为因素影响较小，使用起来非常简单。

对于地震波来说，体波的频率范围约为0.02~100 Hz^[17]，一般地震波首波都是高频信号。为了进一步验证，笔者对P744测站E-W方向分量作了功率谱分析，结果如图 4所示。从图中可以看出，高频地震波的能量主要集中在0.05~0.40 Hz频段上，在一般理论地震波中的体波频率范围之内，因此本文假设频率大于0.04 Hz的信号是所需要保留的高频地震波信号，则所设计滤波器的截止频率为0.04 Hz。对应频率响应关系，本文选取λ=30。在一般的正则化方法中，正则化参数的选择是难点，而本文根据实际物理问题本身的特点(地震波传播频率固定，与震级大小、震中距远近无关)选取了截止频率，进而确定正则化参数。这也解决了应用正则化方法时平滑因子的选取往往会受到人为因素影响的问题。根据地震波的频率范围，同时结合FFT结果，笔者认为在以后使用SPM方法处理1 Hz的GPS数据时，可选取30作为滤波参数的经验值。

2 结果

图 5是3次地震前后相应测站高频GPS动态位移时间序列。从图中结果可以明显看出，高频GPS动态位移记录了地震导致的地壳永久变形和地震波，但两种信号是混叠在一起的。在利用GPS资料反演地震参数或地震破裂过程时，通常需要永久性变形信息或地震波信息作为约束条件。如果参数不够准确，势必会导致反演结果存在一定的偏差，因此，在实际应用中，需要将永久性变形信息和高频地震波信息精确分离。

图 6至图 8分别给出了3次地震前后相应测站高频GPS动态位移的SPM信号分离结果。从图中可以看出，地震引起的永久变形和地震波得到了有效分离。将SPM方法获取的静态永久变形与用地震前后静态天解计算的结果进行了比较，结果见表 2和表 3。从表中的结果可以看出，SPM分离结果与GPS静态天解结果有一些差异。产生该差异的原因有二。一是在地震发生后，震中距较近的GPS站产生较大的静态永久变形，震后一段时间由于断层余滑的影响，在静态天解的静态永久变形结果中包含了主震与震后余滑的共同效应。因此SPM结果与GPS静态天解结果具有差异。如2010-04-04 EI Mayor-Cucapah地震P496和P744测站南北分量，2012-09-05哥斯达黎加地震QSEC和SAJU测站东西分量及南北分量。二是在主震发生后，强余震也可能会进一步引起震中距较近的GPS站产生较大的静态永久变形。因此，在静态天解结果中包含了主、余震引起的静态永久变形的共同效应，使得SPM结果与GPS静态天解结果具有差异。例如2011-03-11东日本大地震在M_w 9.1级主震发生后约30 min又发生了M_w 7.9级强余震，此后约10 min，M_w 9.1震中东部外海又发生了一次M_w 7.7级强余震。由图 8和表 3可以看出，M_w 7.9级强余震发生后，分别引起MIZU和USUD有2.3 cm和2.0 cm的东西向永久性变形。由于M_w 7.7级强余震距离陆地较远，两个测站都没有观测到明显的静态永久变形。

从图 6-图 8中可以看到，利用SPM方法不仅能够准确获得主、余震引起的静态永久变形，同时，得到的高频地震波信号(图 6(b)、6(d)，图 7(b)、7(d)，图 8(b)、8(d))比未分离前(图 5(a)-(f))更加清晰。

如图 8(b)和图 8(d)给出的2011-03-11东日本大地震经过SPM信号分离得到的MIZU和USUD测站高频地震波信号。其中黑色虚线对应的时间为M_w 9.1级主震P波到时，紫色虚线对应的时间为M_w 7.9级强余震P波到时。从图中可知，MIZU测站在M_w 9.1级主震发生后约40 s记录到了地震波信号，在M_w 7.9和M_w 7.7级余震发生后并没有记录到地震波信号；USUD测站在M_w 9.1级主震发生后约105 s记录到了地震波信号，在M_w 7.9级余震发生后约50 s记录到了地震波信号，而在M_w 7.7级余震发生后没有记录到地震波信号。

从以上的分析结果可以得出，利用SPM方法，在完整、清楚地获取地震波信息的同时，还可以清晰地分离主震引起的静态永久变形、余震引起的静态永久变形和震后断层余滑量。

3 结论

(1) 利用GPS静态天解解算地震引起的静态永久变形，不仅需要几天时间才可获得稳定解，而且解算结果包含了主震、余震及震后断层余滑位移，仅靠GPS静态天解结果是难以区分的。而高频GPS数据可以清晰地记录地震引起的静态永久变形和高频地震波信息，采用SPM方法，无须剔除地震时段的数据，因而在完整保留地震波信号的同时，还可以迅速提取静态永久变形信息。分离得到的静态永久变形信息可用于震源参数如震级、震源深度等的反演，而高频地震波信息可以结合地震仪、强震仪等观测资料反演断层破裂过程，为震后快速反应及灾后评估提供保障。

(2) 相比于小波等其他信号分离的方法，SPM滤波方法最大优点是其频率响应的调整仅依赖于一个平滑参数λ，使用起来非常方便简单。当利用SPM方法分离高频GPS动态位移中的静态永久变形和高频地震波信息时，不必对其做频谱分析，λ=30可作为经验值直接使用。

(3) 本文选取的3次地震均为大地震，从分离效果来看，能够实现快速、准确地分离，与前人结果也比较一致。对于小地震引起的静态永久变形和地震波信息，只要高频地震波信号能量在噪声水平以上，也可采用SPM方法快速分离。