海底初至波特性及识别算法研究

李 环，刘云龙，吴 强，邵雨新

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

海底初至波特性及识别算法研究

李 环，刘云龙，吴 强，邵雨新

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

研究运用小波加分形维数的复合算法对地震波中的初至波进行提取。当水中目标振动以后产生海底地震波，最先到达观测点的是初至波。通过对初至波到达时刻的研究，可以对水中的目标进行探测，并通过浅海海底环境模型系统进行地震波数据采集，仿真结果验证了本文所提出的算法的可行性。

初至波；小波；分形维数；

当前水下目标定位技术越来越得到各个国家的关注，已经成为海底目标探测领域不可或缺的一门技术。当海中目标在水下前进时会产生振动，振动信号传到海底时引起海底振动，并在海底地层表面以弹性波的形式向四周传播。这种振动目标所产生的弹性波与普通的地震波相同，称之为海底地震波。根据地震波在介质中的运动特征及传播过程中的特性，通常将其分成两类，既体波和面波。体波分为S波(横波)和P 波(纵波)，两种类型的波在介质中各自独立传播[1]。在测量地震波时，最先被传感器接收到的地震波被称为初至波[2]。

初至波提取是一项基础而又重要的工作，初至波提取的准确性对水下目标的探测与定位有极其重要的意义。随着技术的发展，覆盖次数不断提高，初至波提取工作量越来越大，也占据了相当比重的处理成本。初至波的特点是能量小、频率高、传播速度快，这些特性决定了初至波提取的难度很高，很难精准地提取出初至波到达时刻，需要找到一种办法，既能从包含噪声的复杂地震波信号中检测出初至波到达时信号特征的变化，又能准确地采集到初至波到达时刻。

本文采用分形维数加小波的复合算法，结合小波的定位优势及分形维数处理突变信号的特点，形成分形维数加小波复合算法。本方法避免了小波无法处理分形维信号以及分形维数无法定位信号的缺点，达到提高初至波提取精度的目的。

1 分形维数加小波复合算法提取初至波原理

自从80年代初期提出“分形”一词后，在短短数十年中其概念已发展应用于工程、科学的诸多领域中。这主要是因为分形模型并不假定被研究对象有良好的连续性和平滑性，因此可以用来描述某些建立在欧式几何基础上的不易解决的问题。分形的严格数学定义比较抽象，这里只用简单且形象易懂的方法粗略说明。

分形最突出的特点是在不同尺度下表现出的自相似性，如图1所示。

图1 自相似图形

图1表示在不同尺度下观察图形都具有自相似特性。对于自相似图形，可以运用分形维数的方法对其进行处理。

引申到信号处理领域，如果一个信号具有自相似的特性，也可以用分形理论来对信号进行处理。对于地震波信号而言，从震源产生地震波信号，经过介质的传播后，会产生线性位移和信号的叠加，所以传感器接收到的地震波信号依然具有自相似特性，如图2～4所示。

图2 一尺度地震波信号

图3 二尺度地震波信号

图2～图4显示出不同尺度的地震波信号，可以明显看出，在不同尺度观察地震波信号，地震波信号具有自相似特性。在实际观测当中，传感器所接收到的信号由噪声和地震波信号所组成，噪声的分形维数和地震波信号的分形维数不同，当地震波信号到达时，分形维数将发生变化。如图5所示。

由图5可以清楚地看到，在信号初始阶段是噪声，当真实地震波到达时波形与噪声完全不同，也就是地震波信号的分形维数与噪声的分形维数不同。同时，噪声和地震波信号都属于自相似信号，所以可运用分形理论对地震波信号的突变点进行分析，也就是初至波到达时刻分析。

图4 三尺度地震波信号

图5 传感器接收到的地震波信号

1.1 分形维数原理

N=Kη-D

(1)

式中：K是一个常数；D=1时物体是直线；D=2时物体是平面，以此类推，D=n时是多维空间，既通常的n维概念，由式(1)可见

(2)

对于地震波数据，可用式(1)来估计维数。有

lnN=lnK-Dln(η)

(3)

用实测方法取得一组(N,η)对应值，画在对数坐标纸上(图6)，用最小二乘法对所得各点拟合一条直线，则曲线斜率便反映分形维数。

图6 分形维数

根据分型理论，可以用分形维数来检测信号的变化情况，由此来确定信号的变化趋势。首先模拟出一个具有自相似特性的信号，其中信号开始阶段是信噪比为-25dB的高斯白噪声，有效信号为正弦波信号，如图7所示。

图7 自相似信号

设计一个窗口，窗口大小为10，让窗口每次移动1个窗口的长度，直到窗口覆盖整个信号。窗口内，取5个不同的η值，画出一条拟合曲线图。给出第一个窗口的分形维数计算结果，如图8所示。

计算完所有窗口的分形维数后，可以得到一条分形维数曲线，如图9所示。

图8 拟合曲线图

图9 分形维数曲线

由图9能明显看到，分形维数曲线在第20个数据窗发生了突变。由于信号成分比较单一，可以清楚地看出分形维数曲线突变的时刻。但在处理其它成分复杂的信号时，很难从分形维数曲线中准确地找到信号的突变点。

由于地震波信号具有统计意义下的自相似性，则满足式(1)的信号x(t)可以推导出以下一些尺度上的自相似性：

(1)均值的自相似性。对式(2)求总集均值，得

E[(x(t))]=a(-H)E[x(at)]

(2)自相关函数的自相似性。根据定义

Rx(t1,t2)=E[x(t1),x(t2)]

=a(-2H)E[x(at1)x(at2)]

=a-2HRx(at1,at2)

(4)

如果过程是平稳的，Rx(t1,t2)=Rx(τ)，τ=t2-t1,则有

Rx(τ)=a-2HRx(aτ)

(5)

这是自相关函数的自相似性。

(3)功率谱的自相似性。对式(5)做傅里叶变换得

=a2H+1sx(aω)

(6)

这是功率谱的自相似性。不难证明，满足上式的功率谱具有下述形式：

(7)

可见自相似过程的功率谱有如下特点：

(8)

称为分维谱特性，或1/f谱。

1.2 小波原理

小波变换在信号处理中最重要的应用就是小波的定位功能。小波变换可以将信号的突变时刻明显地表现出来。由此可以发现，当分形维数曲线非常复杂时，可以运用小波的定位功能来对分形维数曲线进行分析，从而找到分形维数突变时刻。小波变换是将信号分解成一系列小波函数的叠加，而这些小波函数都是由一个母小波函数经过平移与尺度伸缩而得来。为了逼近地震波这种非平稳且尖锐变化的信号，不规则的小波函数是最好的选择。

连续小波变换用下式来表示：

(9)

表示小波变换是信号f(t)与被缩放和平移的小波函数φ()之积在信号存在的整个期间里求和的结果。小波变换的结果是许多小波系数C，这些系数是缩放因子(scale)和平移(position)的函数。

小波基的选取对小波分析来说是至关重要的。根据分型理论的要求，需要的小波基是一组两两正交归一、平滑并且具有紧支撑性的小波基，这样得到的小波系数才是互不相关且均值为零的系数。根据这些要求，找到了Daubechies小波和Symlet小波。但是分形维数曲线属于不断变化的带有突变时刻的振动信号，由于Symlet小波本身也属于振动信号，所以对于振动信号的处理，Symlet小波更适合于Daubechies小波，同时用这两种小波对分形维数曲线进行小波变换，发现经过小波变换后，Symlet小波的小波系数大于Daubechies小波，证明Symlet小波和分形维数曲线更加相似，所以选择用Symlet小波来进行小波变换。

1.3 小波与分形维的复合算法

根据上面对分形维数和小波的理论分析，可以先用分形维数算法对信号进行维数分析，作出信号的分形维数曲线，然后再用小波变换的方法对分形维数曲线进行小波变换，由于小波具有定位的功能，可以准确提取出分形维数曲线的突变时间点，得到的突变点即为初至波到达时刻。

本文选用的分形维数算法是周长面积法，根据数据量的大小来设定一个合适的窗，窗的大小为K，通过计算每个窗内的曲线的N个周长和N个面积，并在对数坐标系中对周长和面积作线性回归，即可得到分形维数D。对于地震波曲线S，根据分形维理论有

(10)

式中：L表示S的Hausdorff长度；A表示所包围的欧式面积；D表示S的分形维数；α0是无量纲常数，称为形状因子。

(11)

对于给定的几何图形，α0和D是常数，所以在不同尺度δi(i=1,2,……,n)下，测量地震波曲线S的N个周长值L和N个面积值Ai，在对数坐标系下对数据ln()、ln()作线性回归即可得分形维数D。设计一个窗口，用窗口对整个信号进行覆盖，通过计算每个窗内的分形维数，可以绘制出分形维数曲线。

再对分形维数曲线进行一维离散小波变换，计算出小波变换的系数，对每一层的系数再进行信号重构，即可得到突变时刻。

2 初至波到达时刻提取

在算法的实际应用前，需要先用仿真信号来验证复合算法的可行性。通过Matlab软件，仿真出自相似的信号,采样率为50Hz，如图10所示。

图10 仿真自相似信号

运用复合算法对信号进行处理。首先对信号进行分形维数的计算，设置窗的大小为10，通过计算每个窗的分形维数可以得到一条分形维数曲线，如图11所示。

图11 分形维数曲线

再对分形维数曲线运用小波变换，如图12所示。

如图11、图12所示，正弦波信号经过复合算法后，会在小波变换中准确体现出分形维数突变的时刻，即为正弦波开始时刻。根据正弦波信号的处理结果证明，复合算法是可行的。

根据上述仿真结果，可实现对真实地震波信号的初至波到达时刻的提取。首先，需要采集一段真实的地震波数据，由于条件限制，无法利用真正的大型水下目标来产生地震波数据，所以建立一个模拟浅海海底地震波模型来得到一段地震波数据。将约70cm×80cm的方形铁板浮在约7m深的水池中，用18.6kg重的大铁锤从5m左右的空中自由落体砸在铁板上，以此作为震源来产生水下地震波，通过传感器接收，可以得到一段地震波信号，如图13所示。

图12 小波变换

图13 原始地震波信号

然后对这段地震波数据进行分形维数的计算，取窗口大小为10，得到分形维数曲线如图14所示。

图14 地震波信号分形维数曲线

再对分形维数曲线进行一层离散小波变换既一尺度离散小波变换，得到小波变换的细节系数，如图15所示。

图15 小波变换系数

通过Matlab工具箱，根据得到的小波细节系数进行小波重构，得到一尺度小波变换的重构信号，如图16所示。

图16 一尺度小波重构信号

由图16可以清楚地看到，第51个点为幅值最大的点，即为分形维数曲线变化最剧烈的窗，也就是从第510个数据点开始，分形维数发生了剧烈变化，既初至波到达时刻。对应的原始地震波数据如表1所示。

表1 原始地震波数据

由表1可以看到初至波到达时刻为17.04834s。

3 结论

根据实验结果，证明小波加分形维数的复合算法是完全可行的，解决了小波无法处理分数维信号和分形维数算法无法定位的缺陷，从理论上让二者结合形成复合算法。更深一步探索了小波理论和分形维数理论，同时通过实验研究，验证了复合算法的可行性，提高了初至波提取精度。

[1] 师杰.地震信息的失真现象分析[J] .山西地震,2010(3):45-48.

[2] 杨婷.采空区存在条件下的地震波场特征[D].桂林:桂林理工大学,2013.

[3] 杨阳.广义S变换时频分析的应用研究[D].哈尔滨:哈尔滨工程大学，2011.

[4] 张丽娜.数字信号处理的时频分析方法综述[J].信息技术,2013(6):26-28.

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

[6] 王佳.逆合成孔径雷达对海面舰船目标成像技术研究[D].南京:南京航空航天大学,2010.

[7] 王萌,林云,郭凯丰,等.基于S变换的信号时频特性分析[J].电子世界,2013(14):85-85.

(责任编辑：马金发)

ResearchontheCharacteristicsandIdentificationAlgorithmoftheFirstArrivalWaveinSeabed

LI Huan,LIU Yunlong,WU Qiang,SHAO Yuxin

(Shenyang Ligong University,Shenyang 110159,China)

A new algorithm is studied mainly to extract the first wavefrom seismic wave.When the seismic wave is generated by the vibration target,the first arrival of the observation point is the first arrival wave,and the accuracy of the first arrival wave is directly affected by extraction accuracy.Through studying first arrival wave extraction,a composite algorithm of wavelet and fractal dimension is proposed to extract the arrival time of the first arrival wave.Through the self-made shallow sea environment model system,the seismic wave data acquisition is carried out,and the arrival time is picked up by the composite algorithm,which verifies the proposed algorithm effectiveness.

preliminary wave;wavelet;fractal dimension;

TP391

A

2016-12-26

辽宁省科学技术计划项目(2015020028)

李环(1964—)，女，教授，研究方向：扩频通信技术及应用等。

1003-1251(2017)05-0027-07