﻿ 基于HHT时频分析的水声信号特性提取与建模
 舰船科学技术  2023, Vol. 45 Issue (23): 122-126    DOI: 10.3404/j.issn.1672-7649.2023.23.021 PDF

Underwater acoustic signal characteristic extraction and modeling based on HHT time-frequency analysis
ZHU Yong-yong, LI Zong-ji, WANG Shi-zhe
School of Weaponry Engineering, Naval University of Engineering, Wuhan 430033, Hubei, China
Abstract: In order to complete the modeling of the measured underwater impact signal, it is necessary to fully extract the acoustic characteristics of the signal. In this paper, the measured signal of the pool is taken as the research object, and different analysis methods are used to process the signal. Finally, it is found that the Hilbert Huang transform (HHT) time-frequency analysis method has obvious advantages, and very clear time and frequency distribution information can be seen on the HHT time-frequency spectrum. On this basis, the multi - segment linear frequency modulation (LFM) signal is used to simulate the measured underwater impact signal. Comparing the HHT spectra of the measured signal and the modeled signal, it can be concluded that the spectrum trend of the modeled signal is the same as the measured signal, and has the time-frequency characteristics very similar to the measured underwater acoustic signal, which verifies the advantages of HHT time-frequency analysis and the modeling effect.
Key words: signal time-frequency analysis     short time Fourier transform     Hilbert Huang transform     linear frequency modulation signal     feature extraction
0 引　言

1 希尔伯特黄变换时频分析

1.1 EMD原理

1）求出信号 $x(t)$ 的所有局部极小值点与局部极大值点，利用三次样条插值对求得的所有极值点进行拟合，得出极大值与极小值的包络曲线 ${x_{\max }}(t)$ ${x_{{\rm{min}}}}(t)$

2）对极大值、极小值包络曲线求平均，记为 ${a_{1,1}}(t)$

 ${a_{1,1}}(t) = \frac{{{x_{\max }}(t) + {x_{\min }}(t)}}{2} 。$ (1)

3）将原始信号与 ${a_{1,1}}(t)$ 作差，得到 ${h_{1,1}}(t)$

 ${h_{1,1}}(t) = x(t) - {a_{1,1}}(t)。$ (2)

 ${h_{1,k}}(t) = {h_{1,k - 1}}(t) - {a_{1,k}}(t)。$ (3)

 ${c_1}(t) = {h_{1,k}}(t) 。$ (4)

4）将第1个IMF分量 ${c_1}(t)$ 在原始信号中去除，剩余序列为 ${r_1}(t)$

 ${r_1}(t) = x(t) - {c_1}(t) 。$ (5)

5）将剩余序列看做一个新的原始信号，进行步骤1～步骤4，重复 $n$ 次，可得出第2个IMF分量以及其他IMF分量。

 $\begin{gathered} {r_1} - {c_2} = {r_2} ，\\ \cdots \\ {r_{n - 1}} - {c_n} = {r_n}。\\ \end{gathered}$ (6)

6）直到剩余序列 ${r_n}(t)$ 满足一定条件，比如小于一定的值或成为单调函数，停止分解过程。此时，原始信号被分解为 $n$ 个IMF分量与一个剩余分量 ${r_n}(t)$

 $x(t) = \sum\limits_{i = 1}^n {{c_i}(t) + {r_n}(t)} 。$ (7)

1.2 希尔伯特谱分析

EMD分解方法是基于信号局部特征的时间尺度，将信号“自适应”地分解为有限数目个IMF分量，给信号“瞬时频率”的概念赋予了实际意义。希尔伯特谱分析是希尔伯特黄变换理论中的第2步骤，是在EMD分解基础上进行的。根据式(7)，得知EMD方法将原始信号分解为一系列IMF分量与一个剩余分量，对每个IMF分量 ${c_i}(t)$ 进行希尔伯特变换处理，可得出：

 $\bar {{c_i}} (t) = \frac{1}{\text{π} }\int_{ - \infty }^\infty {\frac{{{c_i}(t)}}{{t - \tau }}} {\rm{d}}\tau 。$ (8)

 ${z_i}(t) = {c_i}(t) + j\bar {{c_i}} (t) = {a_i}(t){e^{j{\phi _i}(t)}} 。$ (9)

 $\left\{ \begin{gathered} {a_i}(t) = \sqrt {{c_i}^2(t) + {{\bar {{c_i}} }^2}(t)} ，\\ {\phi _i}(t) = \arctan \frac{{{c_i}(t)}}{{\bar {{c_i}} (t)}}。\\ \end{gathered} \right.$ (10)

 ${f_i}(t) = \frac{1}{{2\text{π} }}{\omega _i}(t) = \frac{1}{{2\text{π} }} \times \frac{{{\rm{d}}{\phi _i}(t)}}{{{\rm{d}}t}} 。$ (11)

 $x(t) = {{\rm{Re}}} \sum\limits_{i = 1}^n {{a_i}(t){e^{j\int {{\omega _i}(t){\rm{d}}t} }}} 。$ (12)

 $x(t) = {{\rm{Re}}} \sum\limits_{i = 1}^n {{a_i}(t){e^{j{\omega _i}t}}} 。$ (13)

 $H(\omega ,t) = {{\rm{Re}}} \sum\limits_{i = 1}^n {{a_i}(t){e^{j\int {{\omega _i}(t){\rm{d}}t} }}} 。$ (14)

$H(\omega ,t)$ 进行从时间0～T上的积分，得出希尔伯特边际谱（HHT边际谱）：

 $h(\omega ) = \int_0^T {H(\omega ,t)} {\rm{d}}t 。$ (15)
2 实验研究与分析

 图 1 实测信号时域图与功率谱图 Fig. 1 Time domain diagram and power spectrum diagram of measured signal

 图 2 短时傅里叶变换时频图 Fig. 2 Short time Fourier transform time-frequency diagram

 图 3 滤波处理的短时傅里叶变换时频图 Fig. 3 Short time Fourier transform time-frequency diagram of filtering processing

 图 4 HHT时频分析图 Fig. 4 HHT time-frequency analysis diagram

 图 5 滤波与阈值处理后信号的HHT时频图 Fig. 5 HHT time-frequency diagram of filtered and threshold processed signals

 图 6 滤波与阈值处理后信号的放大HHT时频图 Fig. 6 Amplified HHT time-frequency diagram of filtered and threshold processed signals

 $x(t) = A\cos \left(2\text{π} *\left({f_0}t + \frac{\mu }{2}{t^2}\right) + \varphi \right)。$ (16)

 ${s_1} = \left\{ \begin{gathered} 0,t \in [0,1{\rm{ms}}) \\ A = 0.05,{f_0} = 4\;000\;{\rm{Hz}},\mu = 4.8 \times {10^5},\varphi = {90^ \circ },\\ t \in [1\;{\rm{ms}},3.5\;{\rm{ms}}) \\ A = 0.06,{f_0} = 5\;200\;{\rm{Hz}},\mu = \frac{2}{7} \times {10^6},\varphi = {72^ \circ },\\t \in [3.5\;{\rm{ms}},7\;{\rm{ms}}) \\ A = 0.05,{f_0} = 6\;200\;{\rm{Hz}},\mu = 4 \times {10^5},\varphi = - {18^ \circ },\\ t \in [7\;{\rm{ms}},8.5\;{\rm{ms}}) \\ A = 0.025,{f_0} = 6\;800\;{\rm{Hz}},\mu = - 1.12 \times {10^6},\varphi =\\ {162^ \circ },t \in [8.5\;{\rm{ms}},11\;{\rm{ms}}] \\ 0,t \in (11\;{\rm{ms}},20\;{\rm{ms}}]。\\ \end{gathered} \right.$ (17)

 图 7 信号模拟函数的时域图与HHT时频图 Fig. 7 Time domain diagram and HHT time-frequency diagram of signal analog function

3 结　语

 [1] 陈博涛. EMD分解联合时频分析在阵列声波信号中的应用[D]. 长春: 吉林大学, 2010. [2] 肖祺. 低信噪比信号调制识别中的时频分析与应用[D]. 杭州: 杭州电子科技大学, 2021. [3] 黄会婷. 基于时频分析的超声信号处理方法研究[D]. 荆州: 长江大学, 2020. [4] 曹晓勇. 基于时频分析的砼脱空声频信号特征研究[D]. 西安: 长安大学, 2014. [5] 马妍. 基于时频分析的地震信号瞬时参数提取方法研究[D]. 大庆: 东北石油大学, 2013. [6] BABA T. Time-frequency analysis using short time Fourier transform[J]. The Open Acoustics Journal, 2012, 5(1): 245–253. [7] DJEBBARI A, REGUIG F B. Short-time Fourier transform analysis of the phonocardiogram signal[C]. Proceedings of ICECS, IEEE, 2000, 2: 844–847. [8] ZHONG J, HUANG Y. Time-frequency representation based on an adaptive short-time Fourier transform[J]. IEEE Transactions on Signal Processing, 2010, 58(10): 5118-5128. DOI:10.1109/TSP.2010.2053028 [9] MATEO C, TALAVERA J A. Short-time Fourier transform with the window size fixed in the frequency domain[J]. Digital Signal Processing, 2018, 77: 13-21. DOI:10.1016/j.dsp.2017.11.003 [10] PORTNOFF M. Time-frequency representation of digital signals and systems based on short-time Fourier analysis[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1980, 28(1): 55-69. DOI:10.1109/TASSP.1980.1163359 [11] 王红萍. 非平稳信号时频分析方法性能研究[D]. 南京: 南京航空航天大学, 2008. [12] 张旗. 非平稳信号的改进时频分析方法研究[D]. 西安: 西安电子科技大学, 2020. [13] AALAM M K,. Shubhanga k n. emd based detrending of non-linear and non-stationary power system signals[C]//Proceedings of INDICON, IEEE, 2021: 1–6. [14] HUANG N E, CHERN C C, HUAng K, et al. A new spectral representation of earthquake data: Hilbert spectral analysis of station TCU129, Chi-Chi, Taiwan, 21 September 1999[J]. Bulletin of the Seismological Society of America, 2001, 91(5): 1310-1338. [15] HUANG N E, SHEN Z, LONG S R. A new view of nonlinear water waves: the Hilbert spectrum[J]. Annual Review of Fluid Mechanics, 1999, 31: 417-457. DOI:10.1146/annurev.fluid.31.1.417 [16] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society of London. Series A:mathematical, physical and engineering sciences, 1998, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193 [17] 王黎黎. 基于希尔伯特—黄变换的时频分析算法研究[D]. 西安: 西安电子科技大学, 2009. [18] 宋倩倩. 基于Hilbert-Huang变换的语音信号时频分析[D]. 无锡: 江南大学, 2009. [19] 孙涛, 刘晶璟, 孔凡, 等. 小波变换和希尔伯特—黄变换在时频分析中的应用[J]. 中国水运(理论版), 2006(11): 111-113. [20] 杨培杰, 印兴耀, 张广智. 希尔伯特-黄变换地震信号时频分析与属性提取[J]. 地球物理学进展, 2007(5): 1585-1590. DOI:10.3969/j.issn.1004-2903.2007.05.037 [21] HUANG N E. Review of empirical mode decomposition[J]. Proceedings of SPIE - The International Society for Optical Engineering, 2001, 4391: 71-80. [22] HUANG N E. New method for nonlinear and nonstationary time series analysis: empirical mode decomposition and Hilbert spectral analysis[J]. Proceedings of SPIE - The International Society for Optical Engineering, 2000, 4056: 197-209.