目前地震预报尤其是短临预报的水平还很低，远远不能满足公众日益增强的地震安全需求. 解决此问题的关键在于发现震前普遍存在的、 物理意义明确的短临前兆异常信号，而地震“前驱波”似乎具备这一性质. 对“前驱波”的研究已成为地震预报工作者关注的热点.

“前驱波”一词最早由金森博雄等(Kanamori，Cipar，1974)提出. 1974年，他在距震中近万公里的美国加州帕萨迪纳(Pasadena)地震台长周期地震仪记录中发现，1960年5月22日智利8.3级大地震前l5分钟有振荡周期为300—600 s的长周期波，称为“前驱长周期波”. 此后，一些地震工作者对“前驱波”展开了研究与讨论. 冯德益等研究的震前“长周期形变波”与“前驱波”类似，认为长周期观测仪器能够直接观测到这种形变波(冯德益等，1984). 杨又陵等(2003)认为，在2001年昆仑山口西M_S8. 1地震前存在“前驱波”. 赵根模等(2001)及陈德福(2006)对20世纪60年代以来出现的“前驱波”震例进行了不完全整理，所用震例虽有不同，但得出的“前驱波”时空分布特征却基本相似： ① 一般在震前7天至发震前出现； ② 中强震(M_S>6)震前记录较为明显； ③ 震级越大，有效观测距离越长. 陈德福(2006)还总结了“前驱波”应变变化的幅度，多在10^-8—10^-7，个别达到10^-9.

到目前为止，对“前驱波”还没有严格定义. 各种讨论“前驱波”的文献(冯德益等，1984； 杨又陵等，2003； 赵根模等，2001； 陈德福，2006； 张晃军等，2005； 车用太等，2002； 王武星等，2007； 王梅等，2002)给出了震前多种观测仪器记录的“前驱波”图象. 归纳起来，一般系指震前周期为几十秒到几十分钟的地面波动，其主要特征为持续时间较长，振幅较大. 相对于潮汐形变的时间尺度而言，这个过程很短； 与地震的过程相比，这种事件的特征时间尺度则很长(吴忠良，2001).

以往各种文献(冯德益等，1984； 杨又陵等，2003； 赵根模等，2001； 陈德福，2006； 张晃军等，2005； 车用太等，2002； 王武星等，2007； 王梅等，2002)讨论强震前驱波，曾举出过一些震例. 但是，鉴于前驱波的重要性，这种论证是不够的，缺乏系统性，无法说明这种现象是偶然的还是必然的. 例如，M_S7.5地震前是否都存在“前驱波”？ 解决此问题应该对地震目录中全部M_S7.5地震“前驱波”可能存在的时段做系统的分析.

车用太等研究者的结果(车用太等，2002； 高金哲等，2005； 王梅等，2002； 王庆良等，2005； 王武星等，2007； 许昭永等，2003； 张晁军等，2005)表明，倾斜仪(水平摆、 水管仪、 垂直摆)、 重力仪、 应变仪(伸缩仪、 钻孔应变仪)、 水位仪、 长周期地震仪、 海洋潮位仪、 脉动仪、 应力仪等多种仪器能够记录到“前驱波”. 高精度钻孔体应变仪是一种理想的长周期地壳形变观测仪器，已被大量用于美国正在开展的板块边界观测(邱泽华，石耀霖，2004； 邱泽华等，2007). 中国的钻孔应变观测台网有30多个观测点，其中山东省泰安台的观测资料质量最好，2001—2005年在同类观测资料评比中连续5年获得全国第一. 该台站的钻孔体应变仪于1998年9月启用，探头安装在70 m深的花岗岩钻孔中，具备以下特点： ① 观测精度达到10^-9量级； ② 数据记录间隔1 min； ③ 记录到的固体潮分钟值曲线光滑、 平稳； ④ 观测资料连续(王梅等，2004). 该台站有同步的气压、 水位、 气温等辅助观测，为消除所有这些因素对体应变观测的影响提供了条件(邱泽华等，2007).

我们利用山东省泰安台的钻孔体应变仪观测资料，深入分析了2001—2005年全球M_S≥7.5强远震前应变变化，用超限率分析方法对“前驱波”进行了客观、 系统地检验.

根据前人的描述，强震“前驱波”一般在震前7天内出现. 为适当扩大“前驱波”可能出现的时段，本文选取地震发生时刻前30天(43 200个分钟值)数据进行处理.

根据中国地震台网(CSN)地震目录，2001—2005年M_S≥7.5所有震例如表 1所列，总计28次地震，分布于17个月. 对所分析的数据需要说明以下几点： ① 所收集数据起始于2001年1月1日，发生于2001年1月10日、 13日的2次强震前只有10—13天数据，可能影响震前正常和异常情况的对比，故不进行分析； ② 由于当月数据缺失，对发生于2001年7月7日及2002年1月2日的M_S≥7.5地震不作处理； ③ 强震后一定时段可能存在余震，若将相距较近的强震分别处理，可能混淆余震与“前驱波”信息，因此将相距3天之内的强震看作1个事件处理； ④ 对2001年1月26日发生的强震，只取前26天的数据进行分析. 未作处理的震例见表 1备注.

表 1(Table 1)

表 1 2001—2005年Ms≥7.5强震目录Table 1 Parameters of great earthquakes(Ms≥7.5）during 2001—2005

表 1 2001—2005年Ms≥7.5强震目录 Table 1 Parameters of great earthquakes(Ms≥7.5）during 2001—2005

图 1为表 1中28次地震与泰安台的分布情况．5号强震发生在泰安台以西，震中距17748．4km，为震中距最大的处理震例；17号强震发生在泰安台以东，震中距2336.96km，为震中距最小的处理震例．由图 1可知，泰安台周围各个方向均有强震发生．

图 1(Fig.1)

图 1 MS≥7.5震中与泰安台站分布情况 Fig.1 Location of Taian station and epicenters of MS≥7.5 earthquakes

在钻孔体应变观测数据中，日月引力作用引起的固体潮汐和气压变化对构造应变的观测影响最大. 就本文采取的“前驱波”检验方法而言，这两种信号为主要干扰，可能掩盖“前驱波”信息，应将其影响消除.

为达到上述目的，对数据按以下步骤进行预处理：

1)预处理分钟值原始数据. 对缺数据利用插值方法进行拟合，消除由于标定、调仪器等原因引起的突跳变化.

2)为消除固体潮影响，对观测数据进行高通滤波. 我们采用8阶的零相位Butterworth高通滤波器. 该滤波器的归一化截止频率设定为0.015 6. 对应这个截止频率，滤掉了周期在128 min以上的信号，包括固体潮中绝大部分优势成分，得到的滤波结果周期成分为2—128 min，基本包括了“前驱波”出现的优势周期.

3)消除气压干扰. 对气压数据经过高通滤波后发现，气压对体应变数据影响较大，需要在高通滤波结果中将其影响消除. 采用了一元线性回归分析方法消除气压的干扰.

我们利用上述方法处理了表 1中18次地震前30天体应变与气压数据. 图 2给出了一个实例. 从图 2a中的体应变曲线可看出明显的固体潮变化，一方面反映该数据质量较好，另一方面也说明固体潮成分可能掩盖我们所关心的“前驱波”信号. 图 2b给出的气压曲线

图 2(Fig.2)

图 2 23号强震前30天数据预处理结果. (a) 体应变观测曲线； (b) 气压观测曲线；(c) 气压高通滤波后曲线； (d) 体应变高通滤波后曲线； (e) 消除气压干扰后的体应变高通滤波曲线 Fig.2 Pretreatment result of the data within 30 days before No.23 earthquake. (a) Volumetric strain observation; (b) atmospheric pressure change observation; (c) high-pass filtered atmospheric pressure change; (d) high-pass filtered volumetric strain observation; (e) high-pass filtered volumetric strain variation after removing the strain changes induced by atmospheric pressure effect

与图 2a虽然没有明显的相关性，但经过高通滤波处理后的两组数据图 2c和d表现出相似的变化趋势，说明在2—128 min周期段内气压对体应变数据干扰很大，必须予以滤除. 图 2e显示了滤除气压干扰后的体应变数据，与未消除干扰前信号图 2d相比，表明所做的处理是有效的.

小波变换首先由Grossmann和Morlet(1984)提出. Grossmann采用平移和伸缩不变性建立了其理论体系(Morlet et al，1982). 1985年法国数学家Meyer(1985)构造出一定衰减性的光滑小波. 1988年比利时数学家Daubechies(1988)证明了紧支撑正交标准小波基，提出了离散小波分析方法. 1989年Mallat(1989)提出了二进小波的快速算法，使小波变换得到广泛应用. 杨选辉等(2005)用小波方法对测震记录中的核爆炸信息进行了研究. 李琪等(2006)用小波方法提取了震磁效应. 我们率先用小波方法对强震“前驱波”进行了检验.

对于数字信号f(x)，应用离散小波变换DWT(discrete wavelet transform)作为不同频率的信息识别基础，即

式中，a为尺度因子(伸缩因子)，b为平移因子 (a，b∈R). 在计算中，采用a=2^k. 随着k的增加，信号从最高频向低频分解. 当k=0时，信号为采样频率； 当k=1时，将频率2等分； 依此类推(宋治平等，2003).

f(x)可以近似地表示为

式中，a^j_kf(x)与d^j_kf(x)分别为信号f(x)在分辨率为j时的近似部分与细节部分.

设f(x)采样后的信号{f_i}(i=0，1，…，N-1)的频带为[0，Ω]，则a⁰_k频带为[0，Ω]，a^l_k的频带为[0，Ω/2]，d^l_k的频率为[Ω/2，Ω]. 以此类推，a^j_k与d^j_k的频带分别为[0，Ω/2^j]与[Ω/2^j，Ω/2^j-1](图 3).

图 3(Fig.3)

图 3 小波分解的近似部分和细节部分与频率的关系 Fig.3 Approximate part and detailed part of wavelet decomposition result

根据式(2)以及小波分解的近似部分和细节部分与频率的关系(图 3)，对观测资料进行小波分析. 小波分解层数取6，即将2—128 min的信号分为6个不同频段，周期分别为2—4 min、 4—8 min、 8—16 min、 16—32 min、 32—64 min和64—128 min，进一步检验“前驱波”信息. 利用不同的小波函数分解高通滤波结果后发现，选用正则性较好的Daubechies 1小波(简写为db1)分析结果较理想(李杰等，2005)，能够凸显体应变信号细微变化信息.

我们对表 1中所处理的18次强震前30天数据作了小波分析，图 4是23号强震前30天的小波分析结果. 在未作小波处理的曲线中，类似“前驱波”的异常波形并不明显. 但是经过小波分解，细节部分(图 4d₁—d₄层标注“？”处)有类似“前驱波”变化波形的信号出现，波形振幅主要变化在±0.5×10^-9 之间，主要周期成分在2—4 min之间，持续时间1天左右. 信号在图 4d₂—d₄层逐渐减弱，到d₅层消失. 可以看出，体应变观测信号中确实存在类似“前驱波”的波形变化，而小波变换则可有效提取出这种信号，并给出其频率分布特征.

图 4(Fig.4)

图 4 23号强震前30天预处理数据db1小波分析结果. R图表示数据预处理后的信号，d1—d6图为小波分解后的第1层至第6层细节(高频)部分，a6图为小波分解后的第6层近似(低频)部分 Fig.4 Wavelet analysis result during 30 days before No.23 great earthquake. ‘R’ represents the signal after pretreatment, d1—d6 show detailed (high-frequency) part. ‘a6’ is No.6 scale of approximate part

根据车用太等研究者(车用太等，2002； 高金哲等，2005； 王梅等，2002； 王庆良等，2005； 王武星等，2007； 许昭永等，2003； 张晁军等，2005)给出的“前驱波”曲线，其形态与地震脉冲信号不同，是波幅变化明显超过背景噪声的持续时间较长的长周期形变波，与图 4中标“？”处变化类似. 针对“前驱波”这一主要特征，我们提出了“超限率分析”方法来定量检验此类波形.

经过高通滤波，并消除了气压影响后，除地震信号外，观测曲线整体变化相当平稳，可以用算术平均值反映信号值的平均分布情况，而用其标准差描述平时的波动水平. 我们称超过这种波动水平的数据为“超限点”，将单位时段内的超限点数定义为“超限率”. 如果震前存在“前驱波”，根据其波形特征，在相应时段内超出平时波动水平的超限点必然明显多于正常时段. 相比之下，虽然一些中小地震信号也超出平时的波动水平，但是因为其延续时间短，所以不产生很多超限点，超限率会低于“前驱波”.

根据前人的研究，“前驱波”一般在震前7天内存在. 为了不遗漏“前驱波”信号，我们选取了震前30天的数据对其进行统计检验. 以日为单位作超限率的时间序列，作为一级近似，我们用直线来拟合超限率时间序列. 只要震前15天任意时段内存在“前驱波”，回归直线斜率的统计结果就应该明显大于0.

根据这一思路，我们具体按以下步骤处理数据：

1)对小波分解的d₁—d₆层得出每一层数据值的算术平均值和标准差. 需要注意的是，在求这两个值的过程中必须将明显超过标准差范围的“坏点”剔除，否则其均值与标准差的值会偏大，可能将背景噪声值放大到“前驱波”幅度，导致无法区别. 本文将大于3倍标准差的数据值剔除，然后再求其均值与标准差.

2)计算“超限率”. 为不遗漏可能存在的前驱波信息，我们采用未剔除“坏点”的小波分解曲线与步骤1所求剔除“坏点”后的平均值和标准差来计算超限率，计算时段选为1天(1 440个分钟值)，即将一天的“超限点数”取一个值.

3)线性回归. 震前30天共取30个超限率值，对其进行直线拟合，用回归方法求出直线斜率.

4)统计分析. 对所有震例的超限率时间序列的回归直线斜率进行统计.

我们用上述“超限率”分析方法处理了18个震例的小波分解细节部分(图 4d₁—d₆). 图 5为23号强震前30天超限率分析结果. 图 5f₁—f₄中标注“11”处(对应时段为图 4d₁—d₄层标注“？”处)的“超限率”分别为500个、 672个、 624个和624个，标注“12”处的超限率值分别为464个、 532个、 728个和688个. 与该层中的其它值比较，两个值比较大，说明“？”处的超限率分布在11和12两天中，这与我们对该处异常波形的定性认识是一致的； 而图 5f₁中标注“3”处(对应图 4中d₁层第3天)的“超限率”达到586个，是该层超限率的最大值，说明超限率分析不仅能刻画出定性认识到的结果，而且可以定量描述出表现不明显的类似“前驱波”形态的细微变化. 由于这两处超限率高值都存

图 5(Fig.5)

图 5 23号强震前30天超限率分析结果. f1—f6分别为小波分解细节部分d1—d6层通过超限率分析得到的超限率时间序列( )及回归直线 Fig.5 The overrun point ratio analysis result before No.23 earthquake; f1—f6 are the overrun point ratio time series ( ) and the straight lines fitting them

在于震前15日之外，受其影响，超限率时间序列的回归直线斜率k从图 4d₁—d₄层分别为-1.63，-3.77，-1.84和-5.41，全部小于0. d₅和d₆层的超限率为-3.39和-4.40，可能也是由于震前15日外存在的超限率高值，使得这两层的斜率小于0.

需要说明的是，图 5中虽然也分布着一些看似没有规律的超限率高值，但我们的分析最终落足在对超限率时间序列回归直线斜率的统计结果上，其平均值能够反映出超限率曲线的整体变化趋势. 个别超限率高值也可能由偶然原因引起，不影响最终的统计结果，因此不再对其作具体分析.

利用超限率分析方法检验大震“前驱波”的统计结果见表 2. 由表 2可见，各层小波细节部分的超限率回归直线的斜率都很小，甚至比标准差小很多. 并且，对所有震例的斜率 所求算术平均值多数小于0. 对于各层小波细节部分的超限率回归直线，所有震例的斜率取正值和负值的个数大体相当，甚至负值个数稍多，这也表明斜率没有明显的取正值的趋势. 超限率回归直线斜率的平均值、 标准差和正负比这3项统计数据说明，这些强远震前15天内不能普遍检测到“前驱波”.

表 2(Table 2)

表 2 小波分解各层的超限率回归直线斜率统计表Table 2 Statistics of fitting line slopes of the overrun point ratio time series of all wavelet-decomposition detailed parts

表 2 小波分解各层的超限率回归直线斜率统计表 Table 2 Statistics of fitting line slopes of the overrun point ratio time series of all wavelet-decomposition detailed parts

表 2中的4号、 15号和26号强震k(d₁)—k(d₆)全部大于0； 21号和22号强震有5层的斜率大于0; 10号、 14号和20号强震有4层的斜率大于0. 共有8个序号强震与我们对“前驱波”的判定标准基本相符，这种现象不容忽视，需要深入研究.

需要特别说明的是，所有斜率值皆为正并且明显大于标准差的26号强震. 该震前6天左右有应变变化幅度较大、 持续时间较长的异常信号出现，其小波和超限率分析结果见图 6. 经查中国地震台网(CSN)地震目录，标注“？”处时段没有6级以上地震发生，应该不是地震的影响，而且对观测曲线也进行了气压消除，不是气压的影响. 经查当月的原始记录说明，原来是整理仪器电缆的人为干扰所致.

图 6(Fig.6)

图 6 (a) 26号强震前30天db1小波分析结果； (b) 26号强震前30天超限率分析结果； 标注含义同图4、图5 Fig.6 (a) Wavelet analysis result during 30 days before No.26 earthquake. (b) overrun point ratio analysis result during 30 days before No.26 earthquake; tagging is the same with Fig.4a and Fig.5

本文利用山东省泰安台体应变观测资料，采取多种数学方法处理了2001—2005年M_S≥7.5强震前体应变变化，客观、 系统地对强震“前驱波”进行了检验.

我们从震前原始观测曲线没有看出明显异常，因此认为“前驱波”如果存在，其变化可能很微弱，容易被信号中的其它优势成分所掩盖. 通过高通滤波及线性回归等数学方法对数据进行预处理，消除了固体潮汐、 气压等可能影响到“前驱波”检验的各种干扰因素.

用小波分解方法对预处理后数据进行了分析．由于小波变换可给出不同层数（不同频段）下信号分解的近似和细节部分，使我们进一步了解数据在不同频段内的变化．根据细节部分不同频段的变化特点，分析数据时可以正常背景为依据，提取信号中蕴含的微弱异常变化．

针对小波变换提取的振幅较大、持续时间较长的类似“前驱波”的形态特征，我们提出了“超限率分析”方法对这种现象进行检验．统计分析所有震例，若强震前15日内存在类似“前驱波”的异常变化，则超限率在其存在时段内必然提高，回归直线斜率的平均值明显为正，数值应该超过其标准差．表2中最终的统计结果显示，超限率时间序列的回归直线斜率的平均值接近为0，数值远小于标准差，并且对于所有震例的小波细节部分，其斜率正负个数之比接近为1．所有这些都显示了平稳随机过程的特点．这3项统计数据说明，在这些M_s≥7.5强远震前15日内不能普遍检测到“前驱波”．表2中虽然有8个序号超限率的回归直线斜率基本为正，但是否为前驱波影响所致，仍需深入地进行个案分析．

我们提出超限率分析方法，检验了强远震“前驱波”，给出了一个客观的统计结果，说明震前15天内不能普遍检测到“前驱波”，但这也不排除某些地震可能有“前驱波”．并且本文的研究结果与选取的震例均与强远震有关．“前驱波”产生后有可能迅速衰减，从而使距震中较远的钻孔体应变仪不能记录到．我们将对M_s＜7.5的近震展开进一步研究