地球物理学报  2012, Vol. 55 Issue (2): 538-546   PDF    
引用本文
邱泽华, 唐磊, 张宝红 等. 用小波-超限率分析提取宁陕台汶川地震体应变异常. 地球物理学报, 2012, 55(2): 538-546, doi: 10.6038/j.issn.0001-5733.2012.02.016.  
QIU Z H, TANG L, ZHANG B H, et al. Extracting anomaly of the Wenchuan Earthquake from the dilatometer recording at NSH by means of Wavelet-Overrun Rate Analysis. Chinese J. Geophys. (in Chinese), 2012, 55(2): 538-546, doi: 10.6038/j.issn.0001-5733.2012.02.016.  
用小波-超限率分析提取宁陕台汶川地震体应变异常
邱泽华 , 唐磊 , 张宝红 , 宋茉     
中国地震局地壳应力研究所, 北京 100085
收稿日期 2011-04-06, 2011-08-06 收修定稿
基金项目: 中国地震局地壳应力研究所基本科研业务专项项目(ZDJ2008-40)资助.
作者简介: 邱泽华,男,1959年生,研究员,博士,主要从事地震、钻孔应变观测研究.E-mail:qzhbh@163.com
摘要: 姑咱台观测的汶川地震应变异常变得越来越显著.同时,在龙门山断裂带的另一端,原来不好判别的宁陕台钻孔体应变观测异常也逐渐清晰起来.利用小波分解和改进的超限率分析方法,可以将细微的地震异常变化精细地提取出来.一般地,小波分解得到的某个带通信号,可以分为正常背景和异常"坏点"两部分.通过对这两部分界线划分的研究,找到了客观、惟一确定超限率分析的阈值的方法.分析结果表明:沿整个汶川地震发震构造带,应变异常影响范围比我们以前知道的大很多;特别是,这种异常变化明显受固体潮影响,这对进一步揭示其机制提供了新的重要线索.
关键词: 地震前兆      钻孔应变观测      汶川地震      异常提取      超限率分析     
Extracting anomaly of the Wenchuan Earthquake from the dilatometer recording at NSH by means of Wavelet-Overrun Rate Analysis
QIU Ze-Hua, TANG Lei, ZHANG Bao-Hong, SONG Mo     
Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China
Abstract: The observed strain anomaly at station GZ has become quite evident. At the other end of Longmenshan Fault, an initially unseen anomaly at station NSH is more and more obvious in the dilatometer recording as the time lapses. By means of wavelet decomposition and the help of a modified Overrun Rate Analysis (ORA), the tiny anomaly can be extracted from the raw data. Generally, a band-passed signal obtained with wavelet decomposition can be divided into two parts: normal "background" and abnormal "bad points". An approach to locate the objective as well as exact distinction between the two parts has been developed for the determinate threshold of ORA. According to the analysis, the extent of strain anomaly along the seismic fault of the Wenchuan Earthquake is much larger than previously recognised. The anomaly is characterized by small pulses in consistence with peaks of Earth Tide, providing a clue to search for its mechanism..
Key words: Earthquake precursor      Borehole strainmeter      Wenchuan Earthquake      Anomaly extraction      Overrun Rate Analysis     
1 引言

地震前兆是地震预报的基础.地震前兆研究的首要问题是如何识别和提取前兆信号.中国的地震预报已经有数十年的历史，积累了大量观测资料，建立了正确的识别地震前兆的思路[1-3].但是，“九五计划"实现数字化之前，中国地震局地震前兆监测台网的观测数据或者是采样率很低(每日数次)，或者是采用模拟记录，并且缺少辅助观测.因此，难以从中细致地提取异常信息和确切地判定地震前兆.近年来，随着监测台网“九五"的数字化和“十五"的网络化，情况已经大为改观.用什么方法能更好地提取异常信息，已经成为摆在研究人员面前的重要问题.

根据现有的观测结果，大地震前的应力-应变异常变化在震中及邻近地区可能非常显著.例如，1976年唐山地震和1985年乌恰地震的前兆异常，用精度不很高的仪器就可以记录到直观可见的剧烈变化[4-5].但是，受目前的观测台网分布和仪器寿命所限，这种机会非常少.在其他情形里，台站可能距离震中相当远，只能记录到非常微弱的异常信号.例如姑咱台记录的汶川地震的异常变化，在原始曲线上很难看出来.但是用高通滤波处理，就可以将这些短周期的异常信号提取出来[6-8].在此基础上，可以进行各种处理，分析这种信号的特征.

地震是一种“非频发"事件[9].必须珍惜每一个地震(特别是大地震)的观测资料，努力从有限的观测资料获取更多信息，将更微弱的信息提取出来.超限率分析就是以此为目标提出的[7-8].

我们进一步深入研究了超限率分析的阈值问题，提出了客观､惟一的确定正常背景和异常信号的分界的方法.将经过完善的超限率分析与小波分析结合运用，可以最大限度地突出异常信息.利用这种方法，我们又从宁陕台钻孔体应变观测资料中新发现了与姑咱台观测异常可比的汶川地震的相关异常信息.

2 高通滤波和小波分解

宁陕台位于北纬33.33°，东经108.32°.虽然该台距汶川地震主震震中约500km, 距地震余震区也有300km, 但是正好位于汶川地震所在构造带的东北端.宁陕体应变观测点是“十五"建立的.仪器观测分辨率优于10-9，记录到清晰的固体潮.采样率为每分钟一次.该观测点开始产出数据的时间是2007年5月，2008 年5 月12 日汶川地震前有一年左右的观测记录.

在汶川地震后不久我们进行的分析中，由于资料记录比较短，没有看出宁陕台观测资料中存在明显异常.在原始曲线中，尽管同震阶跃十分清晰，但是看不出任何与该地震相关的异常变化.在高通滤波后的曲线中，异常变化也不明显.

至今，宁陕台已经有接近三年的观测资料(图 1中a1).观测系统漂移相当大，经过多次调零，中间还有停电等故障造成的数据中断(图中超幅直线)，这里给出的曲线是拼接而成的.尽管在原始观测曲线(图 1a1)中仍然看不出什么异常，但是在高通滤波信号(图 1b1)中，与汶川地震相关的变化已经隐约可见.这里需要说明的是，2007 年9 月出现的突变是调仪器格值(灵敏度)造成的.

图 1 宁陕体应变(al, bl)( 2007-05-2010-02)和姑咱面应变(a2，b2)(2006-12-2009-01)观测曲线 (al ，2)原始曲线(Sa) ； (bl, b2)髙通滤波信号. Fig. 1 Recordings of Ningshan volumetric strain (a1，b1: May 2007-Feb.2010) and Guza areal strain (a2，b2； Dec.2006-Jan.2009) (a1，a2) Original signal； (b1，b2) High-pass signal (cutoff period: 128 min).

相比之下，姑咱台的观测资料，在原始曲线(图 1a2)中看不出异常，然而在高通滤波(截止周期128min)信号中，汶川地震异常就非常清晰(图 1b2).

我们对姑咱台观测资料的研究实践表明，小波分解是一个更有效的选择[8].小波变换有离散小波变换和连续小波变换两种.我们采用的是由Daubechies[10]和Mallat[11-13]发展起来的离散小波变换.

在宁陕台体应变观测的高通滤波信号中，与汶川地震相关的异常变化表现得不十分明显，但是在小波分解信号中就比较清楚了(图 2).这里将应变信号分解到第9层，考察离散小波变换各层的细节系数，代表原始信号高频部分不同周期段的带通滤波信号.第一层细节部分的信号是2~4min的带通滤波信号，第二层是4~8 min 的带通信号，依此类推.

图 2 宁陕台体应变(Sa)及其小波分解细节部分水平灰线表示超限率分析的阈值，dl- d9分别表示第1到第9层.时间段:2007年5月-2010年2月. Fig. 2 Wavelet analysis of Ningshan volumetric strain (Dec.2006一Jan.2009)White solid lines indicate the thresholds for Overrun Rate Analysis.Sa: original signal； dl一d9 : detail coefficients.

我们看到，与姑咱台观测异常特征类似，宁陕台的很多层小波细节部分的信号在汶川地震前后有可识别的异常变化，波动幅度比其他时段大.用超限率分析可以进一步将这种信号突显出来.

3 超限率分析的阈值确定

超限率分析可以用于研究高通滤波或小波分解(带通)信号.在超限率分析中，关键是如何选取阈值.

在以往的研究中，我们简单地用标准差作为超限率分析的阈值.这样做的好处是，选取的阈值具有客观和惟一的性质.但是，观测资料通常存在一些“坏点".标准差的大小在很大程度上被这些“坏点"左右.这样选取的阈值，会丢失大量有用的信息，不能严格区别异常和正常.

要找到更好的阈值，就要去掉“坏点".这可以有多种方法.因为我们的超限率分析以日为时间单位，所以这里采用的是在日标准差的基础上去除“坏点"求得阈值的方法.

对于一个高通或小波分解信号，求出每日的标准差.数据“坏点"的存在必然使某些天的标准差变“坏"，成为标准差的“坏点".按照一定的百分比(这里记为pbp)，将“坏标准差"去掉，可以得到一个标准差的均值.以这个标准差的均值为阈值，我们可以得到一个数量超限率(每日的超限点数量，Ron)的时间序列，其均值记为mRon.显然，去掉的“坏点"越多，即百分比pbp越大，则阈值越小，超限点就越多，因而mRon越大.一般地，被这样去掉的数据“坏点"与超限点的关系是，前者是后者的一部分.

对于宁陕台体应变高通信号，数量超限率均值mRon与pbp的关系如图 3a所示.mRon随pbp的变化大体可分为一低一高两个阶段:第一个阶段，pbp很小，即去掉的“坏点"很少，mRon的水平也很低；第二个阶段，当去掉的“坏点"达到一定的百分比时，mRon陡然增高，然后趋于平稳.这种mRon 的曲线变化表明了这样的事实:对于我们的观测，“坏点"并不是相对的，而是有特殊的物理意义；“坏点"与“好点"在统计上有明确的界线.这为客观､惟一地严格区别“正常"和“异常"提供了一个依据.

图 3 宁陕台体应变高通信号的mRon-pbp关系 Fig. 3 Relation between mRon (mean of Number Overrun Rate) and pbp(percentage of bad-points) of high-pass signal of Ningshan volumetric strain

我们假定，地震前兆信号包含在“坏点"中.我们需要分析的“异常"，就是那些“坏点"部分，要从这些“坏点"部分中辨别是否存在地震前兆信号.

求出mRon对pbp的导数(图 3b)，我们看到有一个显著的最大值.与该最大值百分比对应的日标准差的均值就是超限率分析所需的阈值.这个阈值正好将“正常"背景和“异常"信号分开.

4 宁陕台的超限率前兆异常

在高通或小波分解信号中，我们把超出阈值范围的点称为超限点，把这个点的数值超出该范围的部分的绝对值称为超限强度，称单位时间内的超限点数为数量超限率(Ron)，称单位时间内所有超限点的绝对值超限强度之和为强度超限率(Ros)[7].

图 2中的水平灰线表示宁陕台观测体应变数据各层小波分解信号的超限率阈值.图 4 和图 5 分别给出了小波分解信号的强度超限率和数量超限率曲线.

图 4 宁陕台体应变小波分解细节信号的强度超限率(Ros)曲线 dl-d9分别表示第1到第9层.时间段:2007年5月-2010年2月. Fig. 4 Wavelet analysis of Ros (Strength Overrun Rate) of high-pass signal of Ningshan volumetric strain (Dec.2006一Jan.2009)Sa: original signal； dl- d9 : Ros detail coefficients.
图 5 宁陕台体应变小波分解细节信号的数量超限率(Ron)曲线 dl- d9分别表示第1到第9层，时间段:007年5月-2010年2月，纵坐标单位为每天的超限点个数，后同. Fig. 5 Wavelet analysis of Ron (Number Overrun Rate) of high-pa^s signal of Ningshan volumetric strain (Dec.2006-Jan.2009)Sa: the original signal； dl- d9 : Ros Detail coefficients.

在强度超限率曲线(图 4)上，原始观测曲线的超幅直线(仪器故障)依然存在，形成一定的干扰；但是，与汶川地震相关的异常已经变得相当明显.震前强度超限率有异常的增大；在地震发生的时刻，从原始曲线上就能看到显著的同震变化，小波分解信号短周期成分的超限强度也是非常大的；震后强度超限率逐渐减弱到通常水平.这些现象与姑咱台钻孔应变观测的异常都是相似的[8].

在数量超限率曲线(图 5)上，仪器故障断点的干扰变得不明显，这是由于这种断点实际上并不多，远远少于数据的“坏点".这就使与地震相关的异常显得更加清晰.由此可见用数量超限率分析地震相关异常的优点.

特别需要指出的是，我们看到宁陕台观测到的汶川地震异常变化与固体潮相关.

5 与附近台站观测的对比

目前，除宁陕外，陕西省还有乾陵､西安和安康三个体钻孔应变观测点.这些观测点开始正式观测的时间都是2007年.其中，除安康台是12月开始有数据外，其他两个台站都是5月就开始有数据了，2008年5月12日汶川地震前有一年左右的观测记录.至今，这些观测点积累了三年左右的观测资料.我们分析了所有这些台站的观测资料.尽管这些台站的体应变观测变化都存在因多次调仪器而出现的数据不连续等缺陷，但是小波分解的长期信号仍然具有重要的参考价值.

根据我们的分析结果，除宁陕台之外，乾陵､西安和安康三个台都没有观测到可识别的汶川地震异常变化(图 6~8).

图 6 安康台体应变小波细节信号的数量超限率(Ron)曲线 dl- d9分别表示第1到第9层，时间段:2007年8月-2010年10月. Fig. 6 Wavelet analysis of Ron (Number Overrun Rate) of high-pass signal of Ankang volumetric strain (Aug.2007-Oct.2010).Sa: original signal； dl- d9 : Ros detail coefficients.
图 7 乾陵台体应变小波细节信号的数量超限率(Ron)曲线 dl- d9分别表示第1到第9层.时间段:2007年5月-2010年10月. Fig. 7 Wavelet analysis of Ron (Number Overrun Rate) of high-pass signal of Qianling volumetric strain (May.2007-Oct.2010)Sa: original signal； dl- d9 : Ros detail coefficients.
图 8 西安台体应变小波细节信号的数量超限率(Ron)曲线 dl-d9分别表示第1到第9层，时间段:007年5月-2010年10月. Fig. 8 Wavelet analysis of Ron (Number Overrun Rate) of high-pass signal of Xian volumetric strain (May.2007-Oct.2010)Sa: original signal； dl- d9 : Ros detail coefficients.

对比汶川地震余震分布､龙门山断裂构造和这些台站的位置(图 9)，我们注意到，宁陕台实际上刚好位于发生了汶川地震的龙门山断裂带的北东端.这与姑咱台刚好位于该断裂的南西端形成一种对应.这样的两个台站都观测到相似的汶川地震前兆异常，可能不是偶然的.从断裂力学的观点看，断裂的端点正是应力集中乃至出现奇异的地方.根据这种看法，其他几个体应变观测点远离汶川地震断裂端点，没有观测到异常变化也是合理的.

图 9 汶川地震断裂与钻孔应变台站位置 * :分量应变台站；圆中菱形符号:体应变台站. Fig. 9 The seismic fault zone of Wenchuan earthquake and the observatory sites Marks * indicate tensor strainmeter sites and © dilatometers.
6 结语

新发现的宁陕台体应变仪观测到的汶川地震异常进一步肯定了姑咱台观测的异常的可靠性.伴随观测时间的延长，这种异常变得越来越清晰.宁陕台观测的异常在形态上与姑咱台观测的异常相似，在幅度上则小得多.这种幅度的差别与宁陕台和姑咱台在震中距上的差别一致.有限的震例表明，这种前兆异常主要发生在地震断裂带上.

当这种形态的地震异常十分显著时，例如1976年唐山地震陡河台和赵各庄台的资料和1985 年乌恰地震喀什台的观测资料，可以从原始观测曲线上直接看出[4-5].对于信号比较微弱的情况，不能直接从原始曲线上看出，例如汶川地震姑咱台的观测资料，可以用高通滤波或小波分解的方法将异常提取出来[7-8].对宁陕台观测资料的分析表明，当信号甚至在高通滤波和小波分解后也不易识别时，利用超限率方法可能提取异常.我们将高通滤波或小波分解信号分为正常背景和异常“坏点"两部分.通过对这两部分界线划分的研究，找到了客观､惟一确定超限率分析的阈值的方法.

超限率分析的作用一方面能最大限度地突出异常，压制正常背景，提取需要的信号.另一方面还能大大凝缩信息，减少运算的时间，提高分析效率.超限率分析还有一个优点，在对宁陕体应变观测数据的分析中表现得很清楚.这就是即使数据存在大量间断点，也仍然可以给出有价值的､令人信服的结果.

我们曾经提出，这种地震异常应变变化在发生机制上可能与岩石力学实验揭示的岩石破坏的声发射关系密切[8].一般的看法是，声发射的机制是微破裂的出现和扩展[14-19].宁陕台体应变仪观测的汶川地震前兆异常变化明显受固体潮影响，这对进一步揭示其机制提供了新的重要线索.

参考文献
[1] 陈运泰. 地震预测研究概况. 地震学刊 , 1993(1): 17–23. Chen Y T. Review of the earthquake prediction research. Journal of Seismology (in Chinese) , 1993(1): 17-23.
[2] 丁国瑜, 梅世蓉, 马宗晋. 地震预报方法. 北京: 地震出版社, 1981 : 174 -179.
[3] 梅世蓉, 冯德益, 张国民, 等. 中国地震预报概论. 北京: 地震出版社, 1993 . Mei S R, Feng D Y, Zhang G M, et al. Introduction to earthquake prediction in China (in Chinese). Beijing: Seismological Press, 1993 .
[4] 李茂玮. 乌恰7.4级地震的土层应力前兆. 内陆地震 , 1987, 1(1): 77–83. Li M W. Precursors of strain in soil before the M7.4 Wuqia earthquake. Inland Earthquake (in Chinese) , 1987, 1(1): 77-83.
[5] Qiu Z H, Zhang B H, Huang X N, et al. On the cause of ground stress tensile pulses observed before the 1976 Tangshan earthquake. Bull. Seism. Soc. Am. , 1988, 88(4): 989-994.
[6] 池顺良, 池毅, 邓涛, 等. 从5.12汶川地震前后分量应变仪观测到的应变异常看建设密集应变观测网络的必要性. 国际地震动态 , 2009(1): 1–13. Chi S L, Chi Y, Deng T, et al. The necessity of building national strain-observation network from the strain abnormality before Wenchuan earthquake. Recent Developments in World Seismology (in Chinese) , 2009(1): 1-13.
[7] 邱泽华, 周龙寿, 池顺良. 用超限率分析法研究汶川地震的前兆应变变化. 大地测量与地球动力学 , 2009, 29(4): 1–4. Qiu Z H, Zhou L S, Chi S L. Study on precursory strain changes of Wenchuan earthquake with ORA method. Journal of Geodesy and Geodynamics (in Chinese) , 2009, 29(4): 1-4.
[8] 邱泽华, 张宝红, 池顺良, 等. 汶川地震前姑咱台观测的异常应变变化. 中国科学 (地球科学) , 2010, 40(8): 1040–1047. Qiu Z H, Zhang B H, Chi S L, et al. Abnormal strain changes observed at Guza before the Wenchuan earthquake. Sci. China (Earth Sci.) (in Chinese) , 2010, 40(8): 1040-1047.
[9] 陈运泰. 地震预测: 回顾与展望. 中国科学D辑 (地球科学) , 2009, 39(12): 1633–1658. Chen Y T. Earthquake prediction: Retrospect and prospect. Sci. China Ser. D (Earth Sci.) (in Chinese) , 2009, 39(12): 1633-1658.
[10] Daubechies I. Orthonormal bases of compactly supported wavelets. Communication on Pure and Applied Math. , 1988, 41(7): 990-996.
[11] Mallat S. Multiresolution approximation and wavelets. Transactions of the American Mathematical Society , 1992, 315: 69-88.
[12] Mallat S G. A theory for multiresolution signal decomposition: the wavelet representation signal decomposition. IEEE Trans. on Pattern Analysis and Machine Intelligence , 1989, 11(7): 674-693. DOI:10.1109/34.192463
[13] Mallat S G. Multiresolution approximations and wavelet orthogonal bases of L2. Trans. of American Mathematical Society , 1989, 315(1): 69-87.
[14] Armstrong B H. Acoustic emission prior to rockbursts and earthquakes. Bull. Seismol. Soc. Amer. , 1969, 59(3): 1259-1279.
[15] Mogi K. Earthquake as acoustic emission-1980 IZU Peninsula Earthquake in particular. J. of Acoust. Emission. , 1982, 1(1): 37-44.
[16] Lockner D A, Byerlee J D, Kuksenko V, et al. Quasi-static fault growth and shear fracture energy in granite. Nature , 1991, 350(6313): 39-42. DOI:10.1038/350039a0
[17] Lockner D. Brittle fracture as an analog to earthquakes: Can acoustic emission be used to develop a viable prediction strategy. J. Acoust. Emission. , 1996, 14: S88-S101.
[18] 陈颙. 不同应力途径三轴压缩下岩石的声发射. 地震学报 , 1981(1): 41–48. Chen Y. Acoustic emission of rocks under triaxial compression along various stress paths. Acta Seismologica Sinica (in Chinese) , 1981(1): 41-48.
[19] 刘力强, 马胜利, 马瑾, 等. 岩石构造对声发射统计特征的影响. 地震地质 , 1999, 21(4): 377–386. Liu L Q, Ma S L, Ma J, et al. Effect of rock structure on statistic characteristics of acoustic emission. Seismology and Geology (in Chinese) , 1999, 21(4): 377-386.