引言

小波变换作为一种信号的时间-尺度分析方法,不仅具有多分辨率分析的特点,而且在时域和频域都具有表征信号局部特征的能力,很适合于探测正常信号中的局部奇异性和瞬态反常现象,并且还可以展示这些反常现象的成分,因此,被业界誉为分析信号的“显微镜”。20世纪90年代以来,小波变换理论在我国地球物理学领域得到了广泛的研究和应用。例如:宋治平等(2003)和刘强等(2007)通过对不同类别地震前兆资料(包括重力、形变、井水位、连续GPS观测等)的小波分析研究表明,小波分析方法对形变数字化资料中的干扰识别与消除、不同频率的潮汐波的分离、趋势异常与短期异常的提取等都具有较强的处理功能,是形变数字化资料分析处理的一种有效方法;邱泽华等(2010)对姑咱台形变资料小波分解的研究表明,在汶川地震发生前,较长周期的异常成分出现的比较晚,这是破裂尺度在增大的表现;池顺良等(2013)对芦山M7.0级地震的研究也表明,距离地震震中70km的姑咱地震台形变资料出现了明显的应变异常,这些异常信号与芦山地震在时间和空间上相关;蒲小武等(2014)通过对武都台应变在汶川M8.0级地震前不同阶段特征的分析后发现,应变异常变化均出现在震前2—3年,持续时间约1年;同时,小波变换理论在地震探油或探矿方面也得到了广泛的应用(孟昭波等,1995)。

钻孔形变观测是研究地壳运动的主要观测手段,同时其在地震前兆观测中也占有重要的地位。本文在南京市数字化钻孔形变观测资料的基础上,将小波变换方法应用于数字化形变资料的处理中,对本区域内的地震前兆异常信息进行了分析和探讨。

1 小波分析理论基础
1.1 小波变换

小波变换概念是法国地球物理学家Grossmann等(1984)在分析处理地球物理勘探资料时提出的,它包括离散小波变换和连续小波变换两种。其中,连续小波变换CWTa,b(Continue Wavelet Transform)的定义为:

(1)
CWTa,b= R f (tΨa,b(t)dt= 1
a
R f(t)Ψ | t-b
a
|dt

式中,Ψ(t)为小波;a为尺度因子,b为平移参数,且参数ab都是连续变化的值,满足a,b∈Ra≠0。

在实际应用中,考虑到信号ft)是离散序列,所以ab也需要离散化。对尺度因子a和平移参数b分别按下式进行离散采样:

a=a0m a0m>0,m∈z
b=nb0a0m b∈Rn∈z

这样,小波Ψa,b可变为:

Ψa,b(t)=a0-m/2(a0-mtnb0)

而任意函数ft)的离散小波可变换为:

DWTa,b=R f (tΨm,n(t)dt
(2)

通常,在对实际资料进行分析时采用a=2k。随着k的增加,信号从最高频向低频分解。当k=0时,信号为采样频率;当k=1时,可将频率二等分。此后依此类推。

1.2 高、低频信息的识别

从频谱分析角度看,正交小波变换DWTa,b是把信号分解到一系列相互独立的频带上,其分辨率j反映了频带的位置和带宽。而数字信号ft)可以表示为(郑治真,2001):

f (t) ≅ Aj f (x) = akj f (x) + dkj f (x)
(3)

式中,akj f (x)dkj f (x)分别为信号ft)在分辨率为j时的近似部分和分解部分,而系数akj f (x)dkj f (x)分别为离散和近似于离散的细节。

实际上,小波变换正交展开系数akj f (x)dkj f (x)是窗口小波变换。当j固定时,随着k的变化,系数akj f (x)dkj f (x)都占满了整个时间域,而且没有重叠。但随着j的变化,dkj f (x)将占满整个频率域,同时也没有重叠。

1.3 小波基的选取

小波基的选取应从一般原则和具体分析对象两个方面考虑。对地震信号进行小波分析时,选择小波基函数的一般原则是:选取的小波基函数具有紧支撑性、对称性和正则性(光滑性)。同时,选择与地震子波形状相近的小波基函数,或以模拟地震子波作为小波基函数。

函数按正交小波展开的分解算法和回复算法,统称为Mallat算法(徐长发等,2004)。Mallat算法为应用小波对信号分解重构提供了非常便捷的手段。而对于正交小波,由于期望它也同样具备有限支集、光滑和超强的时域以及频域的局部化能力等特性,以便使Mallat算法更加快捷,使其在信号分析处理中发挥更突出的作用。正是基于这样的期望(或要求),Daubechies(dbN)小波被构造了出来,该小波是Daubechies(1988)从两尺度方程系数{hi}出发,设计出来的离散正交小波。笔者尝试使用Daubechies小波进行信号处理。

2 南京数字化钻孔形变台站概括

南京地处华北地震区长江上游的黄海地震带内,属中强地震过渡带。在大地构造上,南京位于下扬子构造区宁镇隆起的东段,该地区地质构造复杂,断层比较发育,如图1所示。“十五”期间,南京市地震前兆台网进行了数字化技术改造,2008年底完成了六合台体应变、江宁台和高淳台四分量钻孔应变的数字化技术改造并投入使用,在地震预报方面取得了一定的实效(杨建军等,2008)。笔者选取了数据连续可靠和固体潮清晰的六合台、江宁台、高淳台(见图1和表1)的形变整点值资料,利用db4小波分解(个别缺数作了线性插值处理)并结合2006—2012年南京市周边300km范围内发生的地震事件,探讨和分析体应变、四分量钻孔应变在震前的短临异常信息。

图 1 南京市钻孔形变台站分布及地质构造图 Fig. 1Digital data of deformation stations and geological structure of Nanjing
表1 南京数字化形变台站 Table 1 Digital data of deformation stations of Nanjing
3 基于小波分析的震前地形变异常特征
3.1 体应变小波分析

体应变观测除了可清晰地记录到固体潮汐、地震波、震前异常等信息外,同时也记录到了大气压力、地下水位、降雨、抽水等对地应变场的干扰。因此,对干扰信息的甄别是进行数据分析和准备识别地震异常的前提和基础。而基于小波变换方法可以对形变资料中不同频率的潮汐波信息进行有效的识别,并对这些信息进行有效的分离。

笔者选取六合台2006—2007年形变整点值的观测资料进行分析,所采用的观测数据连续率超过99%(个别缺数作了线性插值处理),M2波潮汐因子相对误差小于0.04,并可记录到清晰的固体潮和远大震的地震波。表2是发生在六合台附近的2次地震事件,其中,安徽定远ML4.7级地震前,六合台体应变出现了明显的趋势性张-压性变换,从7月2日开始,体应变变化量达1088×10-8,日均变幅为43×10-8;响水ML4.0级地震前,在4月9日六合台体应变观测值也出现了明显的下降趋势,同时也呈现张-压性,体应变变化量为1021×10-8,日均变幅达38×10-8,这明显打破了六合台体应变年正常变化的规律。在这2次地震发生前,笔者均填写了临震预报卡。

表2 2次地震发生前六合台体应变异常特征 Table 2 Anomalies characteristics of body strain at Liuhe Seismostation before the two earthquakes

图2是六合台体应变2006—2007年形变整点观测值采用db4小波分解的结果。从图中可以清晰地看到,小波分解细节部分1阶是频率最高的部分,包含了突跳信息等;2—4部分出现了明显的潮汐现象;5—6阶出现了除相应的频率的随机噪声外,还有缺数引起的干扰;7—9阶是频率最低的部分,可能包含有地震前兆信息。通过小波分解得到的不同频率细节部分表明,在2次地震发生前,六合台体应变在7—9低频部分出现了明显的异常情况,这些异常信号的周期在24—30天左右,持续时间在1—2个月之间,之后就发生了地震。对比发生的2次地震可发现,安徽定远ML4.7级地震由于震级相对较大,且距离六合台较近,因而在7—9阶出现的异常信号明显。

图 2 南京六合台体应变2006—2007年形变整点观测值小波分解结果细节部分 Fig. 2Detail parts of wavelet decomposition of body strain at Liuhe Seismostation in 2006—2007
3.2 多分量钻孔小波分析

多分量钻孔应变仪是直接测量地层应变变化的新型地球物理观测仪器,由于分量钻孔应变仪观测数据比观测对象(平面应变3分量)多一路观测数据,构成联系4路观测量间的自检条件:1路+3路=2路+4路+任意常数,因此,观测数据的可靠性能得到很大的保证(邱泽华等,2005)。

笔者对江宁台2008—2009年形变整点观测值资料采用db4小波细节部分进行9阶分解,如图3和表3所示。从图中可以看到,2008年7月6日,句容与南京交界附近(=41.7km)发生了ML3.6级地震,在江宁台的多分量应变整点观测值小波分解第9层细节部分,在震前1个月左右出现了明显的异常信号。同时通过分析可发现,这种异常信号在1分向NS道、3分向NE道和4分向NW道更加明显,并在震后20天左右,这种异常信号逐渐消失。同样,笔者对2009年11月13日发生在江宁、句容交界的ML3.4级地震进行了小波分解,选用了2009年形变整点观测值进行处理分析后发现,从2009年10月15日开始,江宁台多分量应变1分向NS道、2分向NE道和4分向NW道监测到了明显的短临异常,持续时间在30天左右。


(a)一分向NS道;(b)二分向EW道;(c)三分向NE道;(d)四分向NW道
图 3 2008—2009年江宁台分量应变整点观测值小波分解第9层细节部分 Fig. 3Detail parts of wavelet decomposition at the 9th level of multi-component borehole at Jiangning seismostation in 2008—2009.
表3 2次地震前江宁台多分量应变异常特征 Table 3 Anomalies characteristics of multi-component borehole at Jiangning Seismostation before the two earthquakes

高淳台四分量钻孔应变位于县城东20km、茅山断裂带南端游子山脚下,在大地构造单元上属于南京凹陷的边缘地带,为茅西断裂带所控制。自2009年3月完成四分量钻孔应变仪的安装和调试以来,观测资料连续率和内精度都符合技术要求。

2012年7月20日20时11分,在江苏高邮、宝应交界处(北纬33.0°,东经119.6°)发生了M4.9级地震,震源深度5km。这是江苏省20多年来在陆地发生的最大一次地震,震中距高淳台240km。利用小波分析法并结合Daubechies正交小波,笔者对高淳台2012年形变整点观测值进行了分解,图4是高淳台分量应变整点值小波分解第9层细节部分,从图中可以看到,在高邮、宝应M4.9地震发生前一个多月,1分向NS道和2分向EW道出现同步的异常,3分向NE道和4分向NW道的异常趋势一致,这些异常信号不仅频段相同,而且出现的时间也基本一致;不同的是,距高淳台60km的江宁台分量应变在这次地震前4个测向出现的异常特征保持一致,异常时间周期则基本同步(图5)。这些变化可能反映了孕震过程中局部应力的释放并与时间尺度上的微破裂有关。


(a)一分向NS道;(b)二分向EW道;(c)三分向NE道;(d)四分向NW道
图 4 2012年高淳台四分量应变整点观测值小波分解第7—9层细节部分 Fig. 4Detail parts of wavelet decomposition at the 7-9 level of multi-component borehole at Gaochun seismostations in 2012

(a)一分向NS道;(b)二分向EW道;(c)三分向NE道;(d)四分向NW道
图 5 2012年江宁台分量应变分量应变整点值小波分解第9层细节部分 Fig. 5Detail parts of wavelet decomposition at 9th level of multi-component borehole at Jiangning seismostation in 2012
4 认识与讨论

(1)在小波分析理论的基础上,利用Daubechies小波对南京地区数字化形变整点观测值数据进行了信号处理分析,结果表明db4小波对于本区域地震异常信号的判别和分析具有一定的参考价值。

(2)对不同地震同一定点应变的响应研究表明,六合台体应变异常信号出现在小波分解7—9阶,异常信号的周期在24—30天左右,持续时间在1—2个月之间,之后就发生了地震;安徽定远MS4.7级地震异常信号较江苏响水MS4.0级地震明显,这表明震中距越小,震级越大,异常信号越明显。江宁台分量应变的整点值第9阶小波分解表明,在2次地震发生前,4个分向中至少有3个分向的分量出现了明显的异常情况,这种异常的周期通常在20—30天左右。这一异常信号出现的原因可能包含了两种情况:一是从异常起始时间与震中距的对应关系看,异常是由震源周围向外围扩散的过程;二是从异常信号出现的明显程度和时间不同看,这反映了震源机制的不同,所以在各频段呈现出的异常信息也不相同。

(3)对同一地震(高邮宝应M4.9)不同台站的对比研究表明,不同台站对同一地震的异常信息响应的频率和时间基本一致,但异常幅度并不相同;距离震中近的江宁台在震前4个分量出现的异常幅度基本保持一致,而距震中相对较远的高淳台则出现1、2分向和3、4分向一致的现象;但相同的是4个分向在震前的异常周期均为20—30天。

综上所述,应用Daubechies小波法提取地震异常具有明显的优势。对南京地区2006年以来震中距在300km范围内发生的几次地震的分析表明,经重构后的小波信号,在细节部分可以看到干扰被处理的很好,可得到非常理想的曲线;同时,随着分离层数的增加,高频信息所包含的干扰均能够被剔除掉,使得固体潮信号随时间尺度的变化趋势更加明显和直观。经Daubechies小波分解后的体应变、四分量应变仪等数字化形变信号,在震前均呈现出了中期、短期和临震异常特征,这对未来周边地区地震三要素的中期、短期和临震异常判定都具有一定的应用价值。

参考文献
[1]池顺良,刘琦,池毅等,2013.2013年庐山MS7.0地震前的震前及临震应变异常. 地震学报 35(3):296—303.[本文引用:1次]
[2]刘强,宋治平,2007.基于小波分析提取的云南强震数字化形变异常特征. 中国地震 23(3):310—318.[本文引用:1次]
[3]孟昭波,杨丽华,1995.地震资料的小波压缩. 石油地球物理勘探 30(增刊2):70—75.[本文引用:1次]
[4]蒲小武,陈征,高原等,2014.武都分量钻孔应变在汶川MS8.0级地震前的异常变化分析. 震灾防御技术 9(1):133—141.[本文引用:1次]
[5]邱泽华,石耀霖,欧阳祖熙,2005.四分量钻孔观测的实地相对标定. 大地测量与地球动力学 25:118—122.[本文引用:1次]
[6]邱泽华,张宝红,池顺良等,2010.汶川地震前姑咱台观测的异常应变变化. 中国科学:地球科学 40(8):1031—1039.[本文引用:1次]
[7]宋治平,武安绪,王梅等,2003.小波分析方法在形变数字化资料处理中的应用. 大地测量与地球动力学 23(4):21—27.[本文引用:1次]
[8]徐长发,李国宽,2004.实用小波方法. 武汉:华中科技大学出版社 1—191.[本文引用:1次]
[9]杨建军,冯志生,李飞等,2008.根据六合体应变异常所做两次地震短临预测意见的反思. 国际地震动态 350:12—18.[本文引用:1次]
[10]郑治真,沈萍,杨选辉等,2001.小波变换及其MATLAB工具的应用. 北京:地震出版社 110—115.[本文引用:1次]
[11]Daubechies I. ,1988.Orthonormal bases of compactly supported wavelets.Communication on Pure and Applied Math.41: 990—996.[本文引用:1次]
[12]Grossmann A., Morlet J., 1984.Decomposition of hardy function into square integrable wavelets of contant shape.Siam.J.Math.Anal.15: 723—726.[本文引用:1次]