基于HHT的舰船地震波信号实时检测算法

孙庆鹏，周海文

（1.中国人民解放军92213部队，广东 湛江 524064；2.海军工程大学 电器工程学院，湖北 武汉 430033）

0 引言

随着反水雷技术的不断进步，水雷引信中常用的声、磁、水压等物理场信号不断被抑制，水雷引信对目标舰船的检测识别变得越来越困难。舰船海底地震波是近些年发现的新物理场，其是舰船在浅海航行时引起的振动、噪声和水体扰动以压力波的形式经流体介质传播到海底，进而引起海底介质的振动，在海底界面和海底以弹性波的形式向周围传播[1]。舰船地震波信号的成分集中于低频段，其具有衰减小、传播距离远、震动幅度大等特点，可以作为水雷引信的补充[2]。

目前，俄、美两国已经在水雷引信系统中使用舰船地震波场联合其他物理场对目标进行检测和识别[3-4]，但出于保密，其相关技术并未报道。国内的研究主要集中于解释舰船海底地震波形成机理和传播规律上[5-8]，着重分析了舰船地震波波动成分，并研究了其在不同海底介质中的传播特性。对于舰船地震波信号的应用，国内只初步探讨了其可用于舰船目标检测和识别的可能性，例如：文献[9]初步探讨了舰船地震波信号用于舰船目标检测的可行性；文献[10]指出舰船地震波信号能够进行远距离的检测；沈阳理工大学将舰船等效为固定的点震源，利用舰船地震波中的Scholte波特点对舰船目标进行探测和识别[11-12]。

舰船地震波信号是一种非平稳信号，文献[13]基于短时傅里叶变换提出一种实时的舰船检测算法，但短时傅里叶变换的时域和频域分辨率固定，而且其所用检测阈值是固定的，检测效果受环境影响较大；文献[14]中首次运用了小波变换对舰船地震波信号进行分析，但只探讨了小波变化在特征提取中的优势，并未研究具体的舰船地震波信号检测算法，并且小波变换也是以傅里叶变换为基础，其小波基的选择也会影响信号分析的结果。HHT具有简单高效、自适应性强、时频分辨率高等优点，更适用于处理非线性、非平稳信号[15]。本文将利用HHT方法对实测舰船地震波信号进行分析，然后提出基于 Hilbert能量谱的舰船地震波信号实时检测算法。

1 舰船地震波信号HHT分析

1.1 HHT信号分析原理

HHT分析是由经验模态分解（EMD）和Hilbert变换组成。EMD方法能够自适应地将复杂信号分解成一系列具有不同时间尺度的IMF分量和1个残余分量，从高频到低频的顺序依次分解，如下式[16]：

式中：IMFi为EMD分解得到的一些本征模态分量；rn(t)为残余分量。

对EMD分解得到的IMFi分量进行Hilbert变换，可以在不造成信息缺失的前提下，得到一个复IMF分量信号，即

式中：ai(t)为复IMF分量的瞬时幅值；θi(t)为瞬时相位。根据瞬时频率定义可得瞬时频率：

因此，可以得到希尔伯特时频谱，也称Hilbert谱，记为

如果将 H(w,t)对时间积分，就能得到 Hilbert边际谱，即

边际谱提供了对每个频率的幅值测量，表达了整个时间段内幅值的累积。

将式（4）幅值取平方，对其频率积分，可得Hilbert瞬时能量为

Hilbert瞬时能量提供了信号能量随时间变化的情况。将幅值的平方对时间积分，便可得Hilbert能量谱：

Hilbert能量谱提供了对每个频率的能量大小，表示的是每段频率再整个时间长度内所累计的能量。

1.2 实测舰船地震信号HHT分析

根据上述HHT分析原理，现对一组实测的货船地震波信号进行HHT分析。该数据是在汉江某处由地震波采集系统（采集加速度矢量）测量获得，试验时采样频率为 460 Hz，周围无明显的工业噪声污染，测量目标为江上正在行进的货船，航速大约为10 kn，沿水流方向直线行驶。图1为当时实测的试验场景。目标A/B船的吨位大约为2 000 t，航速分别为10 kn、13 kn，采集系统放在沙地中，与目标的正横距离大约为50 m。

图1 试验场景图Fig.1 Experimental environment

文献[10]中对舰船地震波信号的时频进行分析，指出舰船地震波信号在垂直方向上的检测距离最远，频域线谱衰减时间长。图2为实测的货船垂直轴方向的地震波信号，从图中可知当舰船经过测量系统时，加速度有明显的幅度变化。同时，考虑到数据较多，在这截取垂直轴 az方向的一部分加速度数据进行HHT分析，如图2（a）中的红色框中的一段数据，进行HHT分析。

图2 实测货船地震波数据Fig.2 Measured seismic wave data of cargo ship

图3为信号经EMD分解后得到的13个IMF分量和1个残余项（Res）。从图3中可以看出，随着EMD分解的进行，本征固有模态函数IMF的频率从高到低，并且幅值越来越小。最后1项为分解残余项，该部分表征信号总体变化趋势，可以通过消除该项进行零点校准和温度校正。

图3 货船地震波信号EMD分解Fig.3 EMD decomposition of cargo ship seismic wave signal

图4为信号的HHT信号时频图，该图很好地表达了信号幅值随时间和频率的变化过程，用HHT时频图可以很好地完成对突变信号的检测，如图4中所示，突变信号发生在第20 s，突变信号的频率集中在10～30 Hz之间。

图4 信号HHT时频图Fig.4 Time-frequency diagram of HHT signal

图5为信号利用周期图法得到的信号功率谱。图 6-8分别为根据公式（5）-（7）得到的地震波信号的边际谱、瞬时能量谱及 Hilbert能量谱。对比图5和图6可知，信号的功率谱往往会低估低频部分的幅值；瞬时能量谱表达了信号能量随时间变化情况，从图7中可以明显看出在20 s左右信号能量发生变化，这与信号（图2（b））动态变化相一致；Hilbert能量谱反映信号每个频段在整个时间长度内所累积的总能量，从图 8中可知，信号能量集中在10～30 Hz之间，在100～200 Hz频段内存在一些线谱，这是由于信号中含有目标水中传播的直达噪声信号，对比图5和图8，可以看出 Hilbert能量谱能够更好地反映信号高频段的特征。

图5 信号功率谱Fig.5 Signal power spectrum

图6 信号边际谱Fig.6 Signal marginal spectrum

图7 信号瞬时能量谱Fig.7 Signal instantaneous energy spectrum

图8 信号Hilbert能量谱Fig.8 Signal Hilbert energy spectrum

从功率谱、边际谱以及信号的 Hilbert能量谱都能清楚地看出，信号的能量集中在10～30 Hz之间，该明显的频域特征能够用于舰船目标的地震波信号检测。

2 基于Hilbert能量谱的实时检测算法

2.1 构造检测特征量

由前面的实测舰船地震波面HHT分析可知，舰船地震波信号的能量集中于10～30 Hz之间，因此，可将信号的10～30 Hz频段内的能量和作为检测特征量，对舰船目标进行检测。检测特征量计算公式如下：

式中：Ei(a)为 i时刻的 Hilbert能量谱；w1和 w2为选取的频段的下限截止频率和上限截止频率，这里可取w1为10 Hz，w2为30 Hz。

2.2 实时检测算法

为实现对舰船地震波信号的实时检测，基于上述构造的特征量，利用滑动窗进行数据截取和检测特征量更新。检测时，设置检测阈值，当特征量连续多次大于检测阈值时，判断有信号出现。算法的原理框如图9所示。

图9 地震波信号实时检测算法Fig.9 Detection algorithm for seismic wave signal

算法的具体步骤如下：

1）信号预处理，利用带通滤波器对信号进行滤波，带通滤波器的截止频率为10 Hz和30 Hz。

2）Hilbert能量谱计算。利用滑动窗W截取一段数据，对该段数据求解Hilbert能量谱：

式中，H(w,t)为信号的HHT时频谱。

3）特征提取。提取w1～w2频段内的能量，按照式（8）计算i时刻特征量。

4）自适应阈值确定。从初始时刻开始，选取本时刻 i之前的K个特征量平均值作为基阈值U0(i)，后期可根据实际情况利用阈值调整因子μ（＞1）进行调整，即

5）目标检测。从初始时刻开始，若从 i时刻起连续 3次特征量大于阈值，即当iT＞Ui-1、Ti+1＞Ui、且Ti+2＞Ui+1时，则判定存在目标信号。

3 算例

3.1 实测算例

根据上述的地震波信号实时检测算法，对2种货船实测的地震波数据进行试验，该数据是在前面的同样试验条件下获取的。

检测试验时按照算法步骤进行，其中数据滑动窗的窗长N为1 024、步长L为460，信号的采样率Fs为460 Hz，即每秒更新1次特征量，阈值调整窗口K为30，阈值调整参数μ为1.5。

检测的结果如图10-11所示，图10为A货船的检测结果，图11为B货船检测结果。从两张图中可以看出，当目标接近测量系统时，特征量均发生较为明显的变化，两张图中的红色实心“□”为信号检测时刻点，竖直虚线为检测点的对应时刻，A货船中心通过是大约在164 s左右通过，B货船中心大约在177 s左右通过。从图10-11可知，A货船在 153 s时被检测到，B货船在 166 s时被检测到，目标被检测到均在目标中心通过之前完成，且A货船在经过前12 s完成，B货船在前11 s完成，试验时货船的速度大约为10 kn，换算获得检测的距离大约在 60 m，该作用距离基本能够满足水雷的作战需求。

图10 A货船检测结果Fig.10 Detection results of cargo ship A

图11 B货船检测结果Fig.11 Detection results of cargo ship B

从上述结果可知，本文所提出的检测算法能够有效检测地震波信号，且其检测时刻点都在目标中心经过正横之前，具有较好的实时性和一定的检测距离。

3.2 仿真算例

3.1节中所用的实测数据信噪比较高，所以检测效果较为明显，而现实应用中信噪比往往无法控制，因此还需验证算法在低信噪比条件下的检测性能。由于实测时没有获得信噪比低的数据，现截取A货船Z方向上的一段实测数据叠加背景噪声代替低信噪比数据。

第1节中截取的Z方向上的数据长度为21 s，目标中心通过时刻大约在 10 s处。现假设环境背景噪声的长度为180 s，将目标信号每隔10 s叠加到其中，可叠加6段目标信号。环境背景可用高斯白噪声代替，其大小根据设定的信噪比SNR来计算，信噪比SNR的计算公式为10log（信号时域峰-峰值/环境背景噪声峰-峰值）。

仿真时，将信噪比设置为-6 dB，如图12为叠加前后的地震波信号，从图中可以看出，目标信号完全被淹没在背景噪声中。下面对仿真信号进行试验（检测参数按照3.1节中设置）。

图12 仿真信号Fig.12 Simulated signal

利用本文提出的检测算法进行检测，检测结果如图13所示。图中红色实体“□”为检测点，从图中可以看出，在目标出现处特征量都有明显的变化，且算法能够检测到目标信号。

图13 仿真信号检测结果Fig.13 Detection results of simulated signal

表1中列出了检测点和目标中心通过时刻。从表1中可知，在目标中心通过之前就能够完成对目标的检测，但是检测时刻都在目标经过前2 s左右，换算为距离大约是 10 m左右，即检测距离约为10 m。对比 3.1节中的检测结果，可以得出 SNR的变化会影响本算法的检测距离，SNR越小检测距离越近。

表1 仿真信号检测时间Table 1 Detection time of simulated signal

为计算检测算法在低信噪比情况下的检测概率，现按表2中信噪比生成仿真数据，各信噪比数据段中存在100段目标信号，检测结果如表2所示。从表中可知，在信噪比为-20 dB时，信号的检测概率大约为86%，虚警概率大约为6%，若信噪比大于-10 dB，检测概率为100%。因此，由仿真结果可以判断，本文提出的检测地震波信号检测算法在较低的信噪比环境中仍具有较高的检测概率。

表2 不同信噪比下的检测结果Table 2 Detection results of different SNR

4 结束语

本文针对舰船地震波信号的非平稳特点，利用非平稳处理方法HHT对信号进行分析，得到了舰船地震波信号的 Hilbert能量谱，并基于信号的Hilbert能量谱特点，提出基于滑动小波能量谱的舰船地震波检测算法，重点讨论了算法的检测效果及信噪比对检测算法性能的影响。基于实测和仿真数据试验表明：本文所提出的地震波信号实时检测算法能够有效地实时对目标信号进行检测；在正常环境中检测距离大约为60 m，在SNR为-6 dB时检测距离大约在10 m左右，信噪比越小检测距离越近。从不同信噪比仿真信号检测的结果中可以得到，在信噪比为-15 dB时，检测概率达到96%。

本文在验证检测算法性能和计算检测概率时，利用的是在汉江岸边所测数据，目标为普通的货船，对于应用于水雷中的地震波引信应该考虑接受的是舰船海底地震波信号。因此，为进一步验证算法，接下来将进行海底舰船目标检测试验。