被淹没地震信号的小波熵检测与自动识别方法

杨建平，帅晓勇，陶黄林



被淹没地震信号的小波熵检测与自动识别方法

*杨建平1，帅晓勇1，陶黄林2

（1. 井冈山大学电子与信息工程学院，江西，吉安 343009；2. 井冈山大学数理学院，江西，吉安 343009）

为探测大震前的微震，保护大型煤矿、油田和矿山等重要设施，急需地震信号的实时处理、自动识别和提取地震初至点等地震数据处理技术。采用了小波变换和信息熵理论相结合的一种具有多分辨率的复杂度参数——小波熵，该参数能够从被淹没环境中清晰地显示出勘探数据中地震波到达所带来的变化。结合实测数据进行了仿真，并对比了单一的小波变换、数字带通滤波器的监测效果，结果表明小波熵参数能够更好地自动识别微震初至点。

小波变换；信息熵；FIR带通滤波器；地震波；初至点

0 引言

通过地震观测仪记录、采集地震波并设法从波形中获取地震的关键参数是地震勘探中一项最基本、也是最重要的内容之一。分析和处理地震勘探数据，能够及时检测和识别大地震前的微震以及制定对策、及时保护大型煤矿、油田和矿山等重要设施带来帮助。然而，通常在地震勘探数据的采集过程中，由于受到探测区内以及周边各种机械设备的影响，采集数据中不可避免地包含着各种噪声的干扰，有时干扰的强度甚至会很大（地震信号完全被淹没），因此如何从地震勘探数据中自动检测出地震波初至点、发震时间段等参数将非常有利于地震勘探的研究，减少地震灾害带来的影响。

通常，所勘探的地震波数据为非平稳时间信号，变化强度、变化频率都具有典型的非平稳信号特征，因而采用传统的、基于傅里叶分析的、具有全局性意义的功率谱方法无法表征非平稳信号中最根本的特征——时频的局域化特性。目前多用时变谱来刻画时频的非平稳特性，如Iyama、Liu X以及Basu等人[1-6]分别利用在各尺度下能量无泄漏的正交小波基（或双正交小波基）对地震波作小波变换，得到地震波的小波时变谱来探讨具有非平稳特性的地震波，分析小波变换系数与地震波演变之间的关系，取得了一定的研究成果，这些成果也表明小波变换能够较好地处理地震波信号。另一方面，由于地震波数据是地震沿地球表层的复杂振动过程中的一个观测资料，是复杂地质系统中的复杂现象，目前各种复杂现象、复杂系统的研究已成为非线性科学的前沿课题，如各种研究中经常使用的关联维、最大李雅普诺夫指数[7]、Lempel-Ziv方法的算法复杂度[8]、近似熵复杂度[9]等一系列描述非线性特征的参数来定量评价复杂系统的动力学规律。本文欲用小波变换的时频特性，结合描述系统复杂度的熵特征，首先提出并使用小波熵[10-11](wavelet entropy)分析方法对地震波的发生时间进行提取，目前尚没有发现其他文献已采用该方法进行地震淹没信息的识别，该方法有利于计算机系统自动识别、实时处理等操作的执行。

1 地震勘探数据的小波分解思路

小波分析中对地震波信号的分解是基于小波变换的多分辨率分解算法[12]，算法采用正交小波基（或者是双正交小波基）将地震波信号分解为多个不同频带上的分信号之和，这个分解过程本质上是地震波信号的滤波过程，即将地震波信号重复通过一组高通滤波器和一组低通滤波器，逐个进行滤波以达到将地震波信号逐级分解，由高通滤波器输出地震波信号的高频细节分量，低通滤波器输出地震波信号的低频粗略分量。经过第一次滤波后得到的两个地震波分信号所占频带宽度相等，即分别占原地震波信号总频带的二分之一，分解后地震波信号的采样频率降低了一倍；后面的每一步分解都是对其上一步的低频段地震波信号重复使用一对高通和低通滤波器，进而得到下一低频层次上的两个分解分量（低频中的高频和低频中的低频分量）。

2 勘探信号的小波熵

为了将地震勘探信号中地震时间段内、外的小波系数的差异体现出来，而背景信号和地震波信号的各频带信息也不尽详细，为表征这种差异熵无疑是一个最好的参数。信息熵是在一定的状态下定位系统的一种信息测度，可以对勘探信号的未知程度进行度量，可以用来估计这类随机信号的复杂性。对于地震这个不确定性系统，若用一个取有限个值的随机变量表示其状态特征，取值为的概率为，且，则的某一结果所得到的信息用表示，的信息熵为：。

地震勘探信号小波熵具有鲜明的物理意义，是通过计算地震勘探信号的小波分解系数的信息熵所计算出来的一种复杂度，其值的大小与小波系数大小的分布特征有关，小波系数的大小又与地震勘探信号各频率成分有关，因而其小波熵值的大小与地震波各频率的分布有关，但与幅值的绝对大小无关，能反映地震勘探信号的不确定性（变化），小波熵大，则表明地震勘探信号中出现的变化多，地震勘探信号变化的速率快，即大的小波熵值对应地震勘探信号变化是无序而复杂的，其复杂度也大；反之，复杂度越小，则说明发生变化的速率慢，变化是规则的，周期性较强。

3 仿真研究

图1为中国数字地震台网的长春（ChangChun）地震信号波形图，采样间隔为0.02 s，采样频率为50 Hz，从图中可以看出地震信号淹没在振动很强、长周期的主干扰成分之中，另外还包含有随机噪声，很难识别地震波信号的初至点及其结束时刻。

图1 淹没在长周期干扰中的地震信号

3.1 采用数字滤波器的波形识别

长春地震信号主要受简单的单一长周期成分干扰，其所含随机噪声与地震波相似（全部频率成分的分布相对均衡），因而若能消除此长周期干扰，则应能更清晰地观察出地震波的初至点。

这里设计一个有限长脉冲响应数字滤波器——FIR带通滤波器，其阻带的通带边界频率（归一化频率）为wp=[0.032 0.2]，以滤除狭窄带宽的具有单一频率的周期成分（接近正弦信号的波形），获取消除长周期干扰成分的地震波信号，带通滤波器的幅度响应曲线和相位响应曲线如图2中的（a）、(b)所示。将选用的长春台地震勘探信号输入此带通滤波器以滤除长周期干扰成分，获得图3所示的输出波形（消除主干扰成分的地震勘探信号），观察图形可以粗略看出地震的初至点——图中60 s附近波形幅度明显变大的地方。但是，这种变化信息所反映的地震初至点、终止时刻，用人眼仅可粗略判别，若用自动控制系统进行自动识别则非常困难，其原因有两个：（1）如在滤波输出波形的29 s附近处，随机成分的振动也比较强，自动识别时很难设计出合适的算法将此种情况排除在外。（2）在60 s附近地震波形的随机波动性很强，借助计算机进行自动识别也将非常困难。

图2 带通滤波器的幅度和相频响应曲线：(a)幅度响应曲线 (b)相位响应曲线

图3 带通滤波器的输出结果

3.2 地震信号的小波分解与识别

3.2.1 地震信号的小波系数

所勘探的地震波信号实际上可看作由两部分叠加而成，一是背景信号（无地震波成分），另一是真正的地震波信号。对勘探信号进行小波分解可分两步来理解：背景信号的小波分解，其各频带成分体现了背景信号中各种频率的分布特征；地震波信号的小波分解包含了地震波的各种频率信息，因为地震波由到达到结束只经历了采集信号的部分时间段，也就是地震波信号在震前、震后的振幅趋于0。对背景信号与地震波信号的合成信号（勘探信号）进行小波分解，地震时间段上的小波分解系数（相比于地震时间段外）在各频带上变大。由于小波变换能突出信号的变化细节，这种变化比勘探信号中地震时间段内原信号的变化更明显，而地震时间段外的小波分解系数与背景信号的分解情况一致，因而可以利用地震时间段内和地震时间段外的小波系数的差异来更有效地提取地震时间段的信息。

3.2.2 基于小波变换的地震信息检测

对地震勘探信号进行5层小波分解，图4为低频重构信号及1～5层高频重构信号，低频信号主要由长周期干扰成分构成，各高频为噪声及地震波的分解成分，若对不同频带取阈值滤除噪声部分，如文献[6]那样可以达到很好的消噪效果。

图4 地震勘探信号的5层小波分解，横轴为时间，纵轴为幅度

小波变换将地震波分解成各频带的不同分量，根据小波变换的特性，地震波到达时间段上呈现更大的变化细节。从2～5层高频重构信号中均可观察出地震波时刻在60 s附近，但在每一层信号变化的时间点上都有差异。这样差异对于地震勘探仪来说无论采用哪一层高频信号都很难借助自动识别系统进行自动检测并给出初至点及地震时间段，因为各高频信号分量包含了地震波信号中的不同频率成分，频率变化信息所对应的时间位置也各不相同。若用任一高频段去分析地震波信号都将不够全面，得不到反映地震情况的准确信息，因而不能采用单一的小波变换方法进行较为准确的检测、识别及提取相关参数。

3.3 小波熵算法自动检测地震波的初至点

由于采集地震勘探信号时已将周围各种环境信息一并收集，这其中包含有周围的环境信号，想直接从中检测出地震的初至点尤其是当地震信号微弱时是非常困难的。当环境信息相对稳定时，随着地震波到达使得整个勘探信号的成分增多，小波变换的高频系数的幅度将增大，因而信号的小波熵复杂度变大；当地震波停止时，不含地震波的勘探信号小波熵复杂度将变小。因此，可以通过对勘探信号进行小波变换求取各分量系数（或重构分信号），再根据各分量系数并利用信息熵理论的计算公式求得小波熵值，该小波熵值的大小情况能够反映地震波的到达信息。若采用一个移动的窗口并依次计算各移动窗口内的小波熵，可获得一小波熵的变化曲线，该曲线上不同点的小波熵值大小能够反映不同窗口内的勘探信号是否含有地震波。于是，为了探明地震波初至点、地震波时间段，可设计算法对小波熵曲线上的不同点进行识别，以检测出地震波的初至点等参数。通过上述的分析，设计识别算法步骤如下：

步骤1 采用db6小波基（选取原则是正交性好，信号分解后的细节更明显）对地震波信号进行三层小波分解，信号分解为低频系数和一、二、三层高频系数，进行单支重构得a3、d3、d2、d1等四个分信号，即；

步骤2 滑动窗口法计算各窗口内的小波熵：用移动窗口逐步计算每一个窗口内的小波熵，滑动步长为1，即用一个宽度为的时间窗从四个分信号的第1点开始截取长为的一段，根据小波熵定义求出此时间窗的小波熵；

步骤3 将时间窗向右移动一个点，得到四个分信号在第二窗口的小波熵；如此下去，可得到整个地震勘探数据各个时间点上的小波熵，即获得勘探数据的小波熵曲线；

步骤4 计算小波熵曲线的差分获取波形变化率的正、负值，正值为上升沿，负值对应下降沿；以（最大值+平均值）/2为阈值搜索小波熵曲线中的阈值点，阈值点和上升沿重合时为地震波初至点，紧接着的阈值点和下降沿重合点为地震波终止点，并计算两阈值点间的时间宽度；继续按此方式搜索，寻找阈值点和上升沿（或下降沿），获取下一个地震初至点，直到所有勘探信号结束。

图5 地震勘探信号的小波熵曲线

图5为地震勘探信号的小波熵曲线，其波形中小波熵“山峰”清晰可见，这段时间内的小波熵大，包含有地震波到达所致的变化信息，即这个过程对应地震波时间段。根据设定的阈值和上升沿（或下降沿）可检测出地震到达时间点即初至点为56.82 s，终止点为68 s，地震时间为11.18 s。

4 讨论与工作展望

本文虽然仅以淹没在长周期为主干扰成分的地震勘探信号为例，但文中的多分辨率小波熵算法能够广泛地适用于检测地震波初至点的各种不同干扰情形，这个结论可以从小波熵处理信号的理论特点进行推断：只要地震勘探信号中的各种干扰成分在地震波到达前后性质基本相同，则由于地震波到达时刻信号的成分增多、复杂度变大，此时的小波熵必定变大；而当地震波消失时，信号的成分减少、复杂度变小，因而地震波时间段的复杂度大，其小波熵值也大，能够从中识别和自动提取地震波的初至点。该算法在实时监控、自动识别中能够准确地判断微震初至点，但程序在计算中存在计算量较大的特点，今后的研究将致力于减少运算量。

[1] Iyama J, Kuwamura H. Application of wavelets to analysis and simulation of earthquake motions [J]. Earthquake Engineering and Structural Dynamics, 1999, 28(3)：255-272.

[2] Liu X. Time-Arrival Location of Seismic P-Wave based on Wavelet Transform Modulus Maxima [J]. Journal of Multimedia, 2013, 8(1): 32-39.

[3] Basu B, Gupta V K．Seismic response of SDOF systems by wavelet modeling of non-stationary processes [J]. Journal of Engineering Mechanics, 1998, 124（10）： 1142-1150.

[4] Ao M S, Hu Y J, Zhao B, et al. Extraction of High-Rate GPS Seismic Wave Signals with Wavelet Packets Decomposition [J]. Earth Science/Diqiu Kexue, 2012, 37(5).

[5] Mukherjee S, Gupta V K. Wavelet-based characterization of design ground motions [J]. Earthquake Engineering and Structural Dynamics, 2002,31(5):1173-1190.

[6] En D, Wei N N, Wei H H, et al. Research of the Seismic Wave Data Processing Based on the Wavelet Transform[J]. Applied Mechanics and Materials, 2013, 273: 184-187.

[7] Krystal AD，Zaidman C，Greenside HS，et al．The largest Lyapunov exponent of the EEG during ECT seizures as a measure of ECT seizure adequacy[J]. Electroencephalography and Clinical Neurophysiology，1997,103 (6):599-606．

[8] Shaw FZ，Chen RF，Tsao HW，et al．Algorithmic complexity as an index of cortical function in awake and pentobarbital-anesthetized rats[J].Journal of neuroscience methods，1999，93(2)：101-110．

[9] Pincus S．Approximate entropy as an irregularity measure for financial data [J]．Econometric Reviews, 2008, 27(4)：329-362.

[10] 杨建平. 地磁场变化的小波熵复杂度分析方法[J]. 地球物理学进展, 2010, 25(5): 1605-1611.

[11] 杨建平,朱平．分析Lorenz系统动力学特征的新方法[J]. 计算机工程与应用,2012,48(23):230-233．

[12] Mallat SG．A theory for multi-resolution signal decomposition: the wavelet representation [J]. IEEE Transaction on Pattern Analysis and Machine Intelligence, 1989,11:674-693.

METHOD OF DETECTION BY WAVELET ENTROPY AND IDENTIFICATION AUTOMATICALLY FOR SUBMERGED SEISMIC SIGNAL

*YANG Jian-ping1, SHUAI Xiao-yong1, TAO Huang-lin2

(1. School of Electronics and Information Engineering, Jinggangshan University, Ji’an, Jiangxi 343009, China 2. School of Mathematics and Physics, Jinggangshan University, Ji’an, Jiangxi 343009, China)

In order to detect the micro-seismic before large earthquake, protect the important facilities, such asthe large coal, oil, mine and so on. It’s an urgent need for seismic data processing technique, such as real-time process, recognize automatically and extract the submerged seismic onset point. A multi-resolution complexity parameter was acquired based on the wavelet transform and the theory of information entropy, the parameter can clearly shows the change in the exploration data from the arrivals of seismic waves. A simulation was done with the exploration data, Comparison of the monitoring effect of wavelet transform or digital band-pass filter, the results show that the parameter can be very good at the micro seismic onset point for automatic identification.

wavelet transform；information entropy；FIR band-pass filter；seismic wave；onset point

1674-8085(2015)04-0043-06

P318

A

10.3969/j.issn.1674-8085.2015.04.008

2015-03-30；修改日期：2015-06-22

国家自然科学基金项目（31260238）；江西省自然科学基金项目（20151BAB207063）

*杨建平(1970-)，男，江西吉安人，副教授，主要从事信号与信息处理研究(Email：yangjp9273@163.com);

帅晓勇(1977-)，男，江西吉安人，讲师，主要从事通信信息处理(Email：jgsxythl@126.com);

陶黄林(1981-)，女，湖北黄冈人，讲师，主要从事智能ERP、数据挖掘研究(Email：jgsxythl@126.com).