地震波初至波与scholte波的时差提取研究

沈红薇，李 环

(沈阳理工大学 信息科学与工程学院，沈阳 110159)

地震波初至波与scholte波的时差提取研究

沈红薇，李 环

(沈阳理工大学 信息科学与工程学院，沈阳 110159)

通过对初至波与scholte波到达时刻提取的方法研究，提出计盒分形维数算法和修正后的平滑伪魏格纳分布时频分析算法分别对初至波和scholte波到达时刻进行提取，两者时差主要应用于水下舰船目标探测研究。在此背景下，通过模拟浅海海底环境进行地震波数据采集，并拾取初至波与scholte波，获取两者时差，从而验证了所提算法的有效性。

初至波；scholte波；计盒分形维数法；修正的平滑伪魏格纳分布

随着现代舰船防护和隐身能力的不断增强，水雷及声呐等传统探测技术对舰船目标识别与探测难度日益加大，迫切需要某种新型的水下舰船探测及定位技术来适应现代科技的发展。本文利用舰船地震波测距这种目前在水下探测技术方面鲜少用到的方法进行目标识别与定位。

舰船在水下行进时会产生巨大的振动，进而产生噪声信号，这种振动波会以水体为媒介传播到海底，引起海床振动，并在海底地层表面以弹性波的形式向四周传播。这种以行进的舰船为震源所产生的弹性波和普通的地震波相同，被称为舰船海底地震波。地震波又分为体波和面波两类，体波又分为横波和纵波[1]。在海洋地震测量中，均匀地层中地震观测点最先接收到的由震源直接传过来的地震波称为初至波[2]。而沿着弹性介质表面或两个不同弹性介质表面上传播的地震波称为面波，它是介质中横波与纵波相互叠加与干涉所产生的波。本文选用初至波与面波的时差进行测距进而实现水下舰艇目标定位。初至波相对于横波和面波起跳时间最早，到达时间最快，频率最高，幅值最低。正因为初至波具有这些特性，因而大部分初至波拾取算法中均利用这些特点来进行研究。面波主要能量分布在10～80Hz频带内。主要在地表传播，其能量约占地震波总能量的百分之七十左右，波速约为1～3km/s，相比于体波波速较低。面波根据其传播时媒介特质的不同而有不同的面波分类，当仅在固体表面上传播时称为瑞利波，在固体与固体界面上传播时称为斯通利波，而沿着液体与固体界面上传播时称为scholte波[3]。由于本文研究所得结论旨在用在海底舰船目标探测及定位上，故重点研究scholte波的特性及拾取。scholte波沿着海底表平面传播，离开界面时波在这两种介质中逐渐消失，幅度则按指数规律进行衰减。质点在深度-距离平面内做椭圆运动，没有低频截止频率。且具有到达时间最晚、频率低、幅值居中、能量最大等特点。scholte波所在频带内40～70Hz的能量占主要地位，故scholte波的提取通常采用时频分析的方法，观察其能量及频率分布来提取其到达时刻，进而测量其距离可实现水下舰船目标定位。

1 计盒维数法提取初至波原理

1.1 分形维数算法

本文采用计盒分形维数算法来提取初至波的到达时刻。分形维数由法国数学家Mandelbrot于1973年率先引入到信号处理方面进行研究，它研究极不规律、极不连续的几何问题。它不同于以往人们在欧式空间所研究问题的习惯，而是用极简单的方法描述复杂现象，用规律的方法来处理不规律的问题，它为一些非均匀型、突变型、差异性和间断型的研究提供了一种非常有效的方法[4]。分形维数的计算方法有许多，方法不同则得到的维数值的含义也不尽相同，在关于地理问题的研究中，计盒维数法是相对用的较多的一种分形维数算法。取边长为r的多个盒子，用其将分形曲线覆盖起来。覆盖全部曲线所得的盒子数记为N，由此计算其分形维数及其变化趋势。通过几中不同的分形维数计算方法的对比，不难发现计盒维数法更适用于地震道信号的提取与研究。

1.2 地震数据的分形维数计算

当初至波到达以后，地震记录的信号变成了噪声信号与初至波信号的叠加，其能量更强，振幅也较大。通过对比初至波到达前、后计盒分形维数的变化，可以较准确地描绘出其对应时间序列的能量变化趋势，分形维数的第一个突变点即指出了地震波初至到达相对应的时间位置。

定义n个边长为r的盒子覆盖地震道曲线，即其测量长度为L=n×r，改变盒子边长r，则测量的分形曲线长度L也随之变化，发现盒子边长r与分形曲线长度L存在一定的指数关系，即：

L=C×r1-D

(1)

式中,C是比例常数，D为所要求得的分形维数。

在实际算法运用中即在某个需要处理的地震道上确定长度为r的一个工作窗，将工作窗逐个沿地震道向前移动直至覆盖整个地震道曲线，计算每个工作窗内地震道采样序列的分形维数，得到一条分形维数变化轨迹趋势图，此轨迹第一个突变点即为初至波到达的时间点。改变窗口长度r，不同的窗口长度会对应得到不同的分形维数，r越小则计算精度越高，但会有计算时间上的限制。将窗口长度r和分形曲线长度L绘制在同一双对数坐标平面上，代入式(1)，得到分形维公式如下：

logL=(1-D)logr+logC

(2)

拟合各点曲线的斜率和它所分析的分形维数，不难发现两者间存在一定的关系，如图1所示。

图1 窗口长度与分形曲线对数关系

即斜率k=logL/logr，故分形维数与此拟合直线的斜率关系为D=1-k，因此分形维数也可以利用分形曲线长度L与窗口长度r的对数比来求得。这和利用最小二乘法对式(2)来估计分形维数D的原理相同。由此分形维数理论可知，利用计盒维数算法提取初至波可行。系统会自动计算出分形维数突变点所在位置，再根据地震道数据的采样序列查找出突变点所对应的初至波到达时刻。

2 修正的平滑伪魏格纳分布提取scholte波原理

时频分析的主要研究对象是非平稳信号或时变信号，它能清晰地描绘出信号的时变频谱特性及能量分布，在信号处理中有着非常重要的作用。典型的时频分析方法有S变换、小波变换、短时傅里叶变换、魏格纳分布、平滑魏格纳分布等。

时频分析是利用时间和频率的联合分布函数来表示某一非平稳信号，并对此非平稳信号进行分析和处理的一个过程[5]。地震波信号是一种非平稳信号，因此利用时频分析能更准确地将时间、频率和能量清晰地表现出来，以此最终精确地计算出面波scholte波的到达时刻。

本文采用修正的魏格纳分布来对截取后的滤波信号进行时频分析，观察其频率、时间及能量分布的情况来确定scholte波到达时刻。

魏格纳分布是由魏格纳在二十世纪初提出的，起初被应用到量子力学方面，此后由Ville将它引进到信号分析处理方面进行算法研究。魏格纳分布是完全非局部的，能更好地适用于非平稳信号的时频分析。

(3)

式中:W(t,ω)是信号s(t)的魏格纳分布;t为时移;τ为积分变量。

根据希尔伯特变换对s(t)构造解析信号z(t)形式。希尔伯特变换公式如下：

z(t)=s(t)+jH[s(t)]

(4)

式中H[s(t)]定义为

(5)

而在频域上，由实信号s(t)构造的解析信号z(t)魏格纳分布为

(6)

由于魏格纳分布具有良好的时频聚焦性，有时魏格纳分布在时间和频率上把某些值置于两个信号的中间，有时这些值又处于时频平面与所希望得到的信号争夺位置，这样就产生了交叉项，是指不同信号分量之间的相互交叉作用所产生的“虚假信号”[6]。交叉项的幅值可以达到其信号本值的两倍，该现象会严重影响到信号的时频分析特征[7]。所以，为抑制这种“虚假信号”人们提出多种解决办法，如平滑伪魏格纳分布、choi-william分布等。目前应用较为广泛的是平滑伪魏格纳分布。

信号s(t)的平滑伪魏格纳分布为

(7)

式中g(u-t)、h(τ)均为窗函数，用来抑制交叉项，其中当且仅当g(0)=h(0)=1时，g(u)和h(τ)为两个实偶窗函数，魏格纳分布经过加窗处理后即为平滑伪魏格纳分布。

若将平滑伪魏格纳分布适当重新排列，则会进一步改善其时频分析的性能，这种特殊形式的平滑伪魏格纳分布被称作修正的平滑伪魏格纳分布。其分布公式如下所示：

(8)

3 初至波与scholte波的时差提取实验与分析

3.1 采集系统的组成

由于条件限制，无法利用真正的舰船来拾取其舰船地震波进行试验，故模拟了一个浅海海底环境来进行测试。

实验场地：某水声研究所院内，选择外界干扰较少时进行。

实验器材：7m深水池，方形铁板(约70cm×80cm)，18.6kg大铁锤(用于敲击铁板作为震源)。

实验过程：将方形铁板浮在水池上，利用大铁锤从高于水面5m左右的空中自由落体砸在铁板上，以此作为震源放在右侧声源口上。传感器置于地面连接采集系统进行数据采集。

采样频率：2048Hz。

其他硬件实验器材：CA-YD-189型传感器、pxi-3342型采集卡、泛华PXI-9106型机器、YE-3821型信号调理器。

所用软件：Lzbview2013、Matlab2012b。

3.2 数据采集和数据处理

图2为采用泛华 PXI-9106型号机器和Labview2013软件编译的数据采集系统前面板。

图2 数据采集程序前面板

测量幅值范围根据采集卡固有的电压范围，设置为最小-5V最大+5V,采样频率为2048Hz。利用上述采集系统采集到的原始地震数据波形如图3所示。

图3 采集到的原始数据波形

图4为截取其中部分有用信号所得数据波形。

图4 截取后的数据波形

3.3 初至波到达时刻提取

利用计盒分形维数算法对截取后的部分原始地震波数据进行初至波提取。在提取过程中改变盒子边长即窗口长度r分别计算其分形维数。并通过其变化趋势来确定最适合的窗口长度，进而更为准确地提取初至波。通过实际程序中r的更改对比，当r=20时，既是采样频率的整数倍，且能在计算机计算精度较高的情况下计算时间较短，因此选则窗口长度r为20进行初至波到达时刻拾取。

图5为对截取后的数据进行无量纲化处理后，利用计盒维数法得到的窗口长度为20时地震信号分形维数趋势图。

由图5可见，图中程序自动标出的第一个突变点即为初至波到达时刻所对应的分形维数离散点。根据程序计算结果如图6所示，按其窗口长度为20计算其分形维数突变点在第15点。其对应的实际地震信号的离散点为15+1再乘以窗口长度20即320点，对应实际信号数据如图7所示，可得其初至波到达时刻为17.055664s。

3.4 scholte波到达时刻提取

由于scholte波主要能量分布在10～80Hz之间，故对图4的截取后的信号波进行滤波去噪，采用阶数为3阶的巴特沃兹带通滤波器进行滤波，通带范围为10～80Hz，由此截取其频带在10～80Hz之间的频率信号。滤波后的截取信号如图8所示。

图5 地震信号分形维数趋势图

图6 分形维数突变处对应的离散点

图7 初至波到达时刻数据

图8 滤波后信号

再将图8所示滤波后的有用信号利用修正后的平滑伪魏格纳分布进行时频分析，得到图9所示的时-频能量分布图。可以准确地看出10～80Hz频带内信号的时间、频率及能量密度的分布情况。颜色越深则表示能量越大。由于scholte波同其他干扰波相比能量最大，频率居中，且已经10～80Hz外的信号滤掉，因此将其时频分析能量分布图转换视角并进行部分放大如图10后，可明显观测到在6.72s时刻scholte波到达，成功接收到scholte波。

图9 时频分析能量分布图

图10 放大后的能量分布图

3.5 时差提取

由以上结论可知，对于所采集的实际地震信号进行数据处理后，提取到其初至波到达时刻为17.056s，其scholte波到达时刻为17.42s，两者时差为0.364s。在此初至波与scholte波时差提取成功的基础上，可以设置多个传感器在不同方位进行数据采集，利用时延定位法计算每个传感器对该数据的接收时间时差和已知的波速得到目标到传感器的距离进而对目标进行定位。

4 结论

(1)利用计盒分形维数法能较准确地提取初至波到达时刻，且效率明显提高。

(2)利用修正后的平滑伪魏格纳分布进行scholte波的提取，能很好地抑制交叉项，且很好地表现了能量的分布情况。结合scholet波能量最强频率居中等特点进行时频分析，这样可以非常清晰地提取scholte波到达时刻，进而得到初至波与scholte波的到达时差，以便实现水下目标定位。

[1]李响,颜冰,周穗华.基于直达波与地声界面波时延的目标定位方法[J].探测与控制学报,2010,32(6):50-53.

[2]曾富英,李敏锋,申维.地震波初至波拾取的分形研究[J].现代地质,2002,16(2):209-213.

[3]张海刚.浅海甚低频声传播建模与规律研究[D].哈尔滨:哈尔滨工程大学,2010.

[4]钱光萍.基于分形维地震道初至拾取方法研究[D].成都:成都理工大学,2001.

[5]马见青.多分量地震自适应极化滤波方法[D].西安:长安大学,2012.

[6]田琳,陈颖频,梁华兰.平滑伪Wigner-Ville分布在地震信号处理中的应用[J].新疆师范大学学报:自然科学版,2013,32(3):1-4.

[7]金银燕,于凤芹.Gabor变换与双线性时频分布的时频结构[J].计算机工程与应用,2011,47(25):146-148.

(责任编辑：马金发)

Research on Time Difference Extraction of Seismic Preliminary Wave and Scholte Wave

SHEN Hongwei,LI Huan

(Shenyang Ligong University,Shenyang 110159,China)

By researching the method of preliminary wave and scholte wave of arrival time extraction,a box counting fractal dimension algorithm is proposed to extract the time of preliminary wave arrival and give a correction to extract the time of scholte wave arrival after the smoothed pseudo Wigner Ville distribution frequency analysis method.The time difference between the two waves is mainly used in underwater ship target survey research.In this context,the shallow seabed environment model system is taken to extract preliminary wave and scholte wave with time difference for seismic wave data collection，which verifies the algorithm validity.

preliminary wave；scholte wave;box counting fractal dimension;modified smoothed pseudo Wigner-Ville transform

2016-01-08

沈红薇(1990—)，女，硕士研究生；通讯作者：李环(1964—)，女，教授，研究方向：扩频通信技术及应用等。

1003-1251(2016)06-0082-06

TP391

A