0 引　言

为了对水下实测信号进行建模，充分提取信号声学特性十分必要，这对水下信号处理、水下攻防、水下探测与识别等研究领域具有重要意义。时间和频率是描述信号的2个重要物理量，时频特征是声信号的重要特性，因此从时域、频域2个角度进行信号处理是提取信号声学特性研究中较常见的方法。

时域分析是在时间的维度对信号进行分析与处理，但实际中采集到的声信号一般较为复杂，存在信息相互叠加的现象，仅仅时域分析很难从中提取到有效信息。傅里叶变换解决了这个问题，使很多在时域中难以看到的信号特征在频域下很好的体现出来^[1-5]。但傅里叶变换存在一定局限性：它是一种整体的、全局的积分变换，要么完全在时域，要么完全在频域，无法得出信号的时频局部特性。也就是说，通过单独的时域分析与频域分析，并不能知道信号在某一特定时间点的对应频谱，想要完成水声信号的建模，仅靠时域图与频域图得出的声特性信息还远远不够，需更清楚的了解水声信号的时频特性，这就要对信号进行时频分析。

短时傅里叶变换^[5-12]时频分析理论是由Gabor在20世纪40年代提出的，是在傅里叶变换基础上的推广改进，其原理简单、功能强大，在今天仍然是信号分析与处理的有力工具。但针对某一特定信号进行短时傅里叶变换时频分析时，选择何种宽度、何种类型的时间窗使信号分析效果更好，这是不能确定的，也就是说，时间分辨率与频率分辨率之间是相互矛盾的。为解决这个问题，Norden. Huang提出了一种全新的信号处理时频分析方法——希尔伯特黄变换方法（HHT），它由EMD分解和Hilbert变换两部分构成，非常适用于处理水下冲击信号等这类非平稳信号^[13-16]。

本文采用希尔伯特黄变换（HHT）时频分析方法，首先将信号进行HHT时频分析，得到信号声学特性；然后进行数学建模，构造函数表达式表示原始声信号。最后，比对模拟信号与原始信号的HHT分析结果，验证建模效果。

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

希尔伯特黄变换方法（HHT）是基于经验模态分解与希尔伯特谱技术，将复杂水声信号分解为有限数目个固有模态函数，赋予了信号瞬时频率合理的意义^[17-20]。该方法分析信号的基本过程由两部分组成：首先，信号经过经验模态分解（EMD）后，自适应的按序排列为一系列有限数目的固有模态函数（IMF）；然后，将固有模态函数经过希尔伯特变换得出时频谱图。

1.1 EMD原理

在希尔伯特黄变换分析信号的2个步骤中，经验模态分解（EMD）是其中的核心部分，是将待分析信号分解为一系列固有模态函数（IMF）。EMD分解方法体现了“筛分”思想，基于上述假设，对信号 $x(t)$ 进行分解的步骤如下：

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,1}}(t)$ 是否满足固有模态函数的2个限制条件：在整个时域内，信号过零点的数目与极值点的数目必须相等或最多相差1；在信号的任一数据点，极大值包络线和极小值包络线的均值为0。如果不满足这2个条件，需重复进行步骤1～步骤3。在这个重复过程中，会得到新的差值序列 ${h_{1,k}}(t)$ ， $k$ 为求取第1阶固有模态函数需要的迭代次数。

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

(3)

直至 ${h_{1,k}}(t)$ 满足上述2个限制条件为止，则 ${h_{1,k}}(t)$ 为第1阶固有模态函数。

$ {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)

分析EMD分解的整个过程，可看出此分解是一种具有经验性、自适应性与滤波性的过程^[21]。也就是说，整个分解过程可看做是一种自适应的“筛选”过程^[22]，且得到的一系列IMF分量是按照频率成分从高到低排列的。经过EMD分解得出的IMF分量非常适合进行希尔伯特变换，并且通过IMF的2个限定条件，可有效保证EMD分解得到的每一阶IMF分量经过希尔伯特变换求得的瞬时频率具有合理的物理意义。

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)

结合原始IMF分量，可构造出解析信号表达式：

$ {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)

根据式(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)

式中， ${{\rm{Re}}} $ 的含义为取实部。剩余分量 ${r_n}(t)$ 一般是单调函数或者是一个常量，因此这种表示形式忽略了剩余分量的影响。该表达形式用傅里叶级数可表示为：

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

(13)

式中， ${a_i}$ 、 ${\omega _i}$ 均为常量。分析式(12)，发现信号的幅值是以时间 $t$ 与瞬时频率 ${\omega _i}$ 为自变量的一个函数，这种信号幅值的时间-频率分布叫做希尔伯特谱（HHT谱）：

$ 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所示。

分析图1，发现单独对信号的时域与频域进行分析时，能抓取的信号特征很少。观察时域图，得出该信号为冲击信号，是非平稳信号，重物落水瞬间的信号能量较高。观察频域图，发现该信号是一个频带较宽的宽带信号，信号各频带对应的能量不同。对发现在2～7 kHz频带内能量较强，但找不出更多信息和更好规律来模拟该信号，

对信号进行短时傅里叶变换时频分析，结果如图2所示。分析时频图，可得知该信号的大部分能量集中在100 Hz～10 kHz范围内，进行滤波处理，去除频率小于100 Hz和大于10 kHz的部分，得出的时频图如图3所示。

分析可知：在4 ms处有持续1.5 ms左右的从3～7 kHz的能量强区。从时频图可看到信号时频平面上的整体图，但仅可拿到较少的信号时间-频率的分布信息，无法给出时间-频率的函数关系，因此想要建立数学模型，拿到的这些信息还远远不够。

对100 Hz～10 kHz滤波处理后的实测信号进行HHT时频分析，结果如图4所示。

分析HHT时频谱图可发现：对于该信号某一频率成分出现的时刻，希尔伯特黄变换给出了精确结果。也就是说，可以给出瞬时频率与能量的分布情况，即非常清晰的时间、频率分布信息。在这些时频信息的基础上，可进行信号的建模研究。

对时频图设置一个阈值，去除能量较低的部分，这样使得时频图中时间与频率的关系更加明显。设置阈值限制后，选取其中能量强的部分，100 Hz～10 kHz滤波与阈值处理后该信号的HHT时频图如图5所示。观察图5可得：时间与频率基本满足线性关系，可将频率看做是一个自变量为时间的分段函数，因此该信号等同于一系列线性调频信号的“加和”。

将HHT谱图放大并标示坐标，如图6所示。

该信号的时间-频率分段函数可看做由4部分线性成分构成，因此可使用线性调频（LFM）信号模拟该信号。由图可知，该信号冲击时间大概持续10 ms，因此构造如下函数表达式：在0～0.02 s时间段内，包含10 ms的线性调频信号，其他时间段进行补零处理。构造信号的采样率与原始实测信号采样率保持一致，均为128 kHz。LFM信号的实数形式表达式如下式：

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

(16)

式中： $A$ 为信号幅度； ${f_0}$ 和 $\mu $ 分别为信号初始频率和调频率； $t$ 为时间； $\varphi $ 为初始相位。在进行建模时，需保证建模函数与原始实测信号的HHT分析频谱走势相符，能量大致相同。构造函数为分段函数，各段的相关参数值如下式：

$ {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谱图可知，模拟函数较好体现了时频变化规律，建模信号的频谱走势与实测信号相同，具有与实测信号非常相似的时频特性。在各段信号接头处出现异常现象，是因为接头处并未完全光滑，频率发生了瞬间跳变与振荡。

3 结　语

本文以实测水下冲击信号为研究对象，从时域与频域2个角度出发，对信号进行分析处理以提取声学特性。最终发现，希尔伯特黄变换时频分析方法可得到非常清晰的时间、频率分布信息，可以完成实测信号的数学建模。比较实测信号与模拟信号的HHT时频分析结果可知：用分段线性调频信号来模拟实测信号，可以得到相同的频谱走势与相似的时频特性，验证了HHT时频分析的优势与建模的可靠性。