0 引　言

助推过程中振动噪声一方面影响水下航行器的安全性，另一方面会影响航行器的水声设备正常工作。对助推过程中振动噪声特性的研究显得尤为重要，以振动特性为依据可以评估助推系统的振动性能^[1]、研究振动机理^[2]、研究振动控制方法^[3-4]。最后可将这些研究结果应用于系统减振降噪以及其它方面的优化设计。

信号处理方法是研究振动特性的常用方法，通常包括时域分析方法、频域分析方法，以及时频分析方法。时频分析方法不仅能够获取振动信号的频谱成分，而且能够获取某一频率成分出现的时间点。常见的时频分析方法有短时傅里叶分析、小波变换、HHT变换、VMD等。VMD方法最早由 Konstantin Dragomiretskiy 等^[5]在2014 年提出，大量实践表明该方法在处理非线性非平稳信号时能取得较好的结果。

李占龙等^[6]采用VMD方法对含噪声、冲击和间断信号这3类典型的非平稳信号进行分解，并结合Teage算子对分解得到的IMF信号进行解调得到瞬时的时频图。结果表明VMD方法抗混叠性能较好。杨宗林等^[7]利用GA-PSO算法对VMD算法中所需要预设的模态个数与惩罚因子进行了优化。柳絮等^[8]结合HHT和VMD对某海上平台进行振动分析。刘长福等^[9]结合VMD和FFT算法提取了颤振动特征。将VMD方法应用到助推过程并进一步分析振动源和振动传递等特性研究尚不充分。

本文基于变分模态分解方法结合希尔伯特黄变换^[10]分析系统助推时振动信号的时频分布情况，基于时频分布结果进一步分析振动源类别、振动传递特性。以所得到的振动特性为依据，结合主动^[11]或半主动方法^[12]抑制振动，在振动传递过程中对振动进行控制^[13-14]、对各低频^[15-16]或高频振动源提出针对性强的控制方法，综合评估系统的振动性能。

1 变分模态分解理论

变分模态分解方法是建立在维纳滤波、希尔伯特变换、解析信号、频率混合和外差解调基础上，在给定的变分模型下进行迭代从而将输入的实信号 $f$ 分解成调幅调频信号 ${u_k}$ 。通过处理观测到的信号 ${f_0}$ 恢复出原始的信号 $f$ 是一个典型的反问题，维纳滤波方法中采用Tikhonov准则^[17]处理这个病态问题。

实信号 $f\left( t \right)$ 的希尔伯特变换 $\mathop {f\left( t \right)}\limits^ \wedge $ 是将 $f\left( t \right)$ 通过一个特定的滤波器，在时域内 $ h\left( t \right) = {1 \mathord{\left/ {\vphantom {1 {\pi t}}} \right. } {\pi t}} $ ， $f\left( t \right)$ 的希尔伯特变换为 $f\left( t \right)$ 卷积 $ h\left( t \right) $ ，在频域内希尔伯特变换相当于一个乘算子。在实际信号处理中传感器测到的信号均为实信号，为了处理上的方便，常采用希尔伯特变换构造解析信号，解析信号在求取包络、相位、频谱搬移等方面比较方便。在变分模态分解方法中采用的混合器为乘法器，在傅里叶域中，将指数信号与解析信号相乘结果是简单的频移。

基于以上理论，对于每个 ${u_k}$ ，通过希尔伯特变换得到信号的单边频谱、通过混合指数信号估计中心频率 $ {\omega _k} $ 、通过平方范数 ${L^2}$ 估计带宽，引入二次惩罚项 $\alpha $ 和拉格朗日乘子 $\lambda $ 。最后所得到的变分模态分解模型如下式：

$ \begin{split}\mathop {\min }\limits_{\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\}} &\alpha \left\{ {\sum\limits_k {\left\| {{\partial _t}\left[ {\left( {\delta \left( t \right) + \frac{j}{{\pi t}}} \right) . {u_k}\left( t \right)} \right]{e^{ - j{\omega _k}t}}} \right\|_2^2} } \right\} +\\ &\left\| {f\left( t \right) - \sum\limits_k {{u_k}\left( t \right)} } \right\|_2^2 + \left\langle {\lambda \left( t \right),f\left( t \right) - \sum\limits_k {{u_k}\left( t \right)} } \right\rangle \end{split}，$

(1)

$ {\rm{s}}.{\rm{t}}.\sum\limits_k {{u_k} = f} 。$

(2)

式中， $ \delta \left( t \right) $ 表示Dirac函数， $ {u_k} $ 表示各本征模态函数。

对于 ${u_k}$ 的更新主要是基于维纳滤波器，更新的方式如下式：

$ \hat {u} _k^{n + 1}\left( \omega \right) = \frac{{{f\left( \omega \right) -\displaystyle {\sum}_{i \ne j} { {{\hat{u_i}}^{\left( \omega \right)}} + {\dfrac{{\lambda}}{2}} }} }}{{1 + 2\alpha {{\left( {\omega - {\omega _k}} \right)}^2}}} ，$

(3)

对于 $ {\omega _k} $ 的更新，其迭代收敛的准则如下式：

$ \omega _k^{n + 1} = \mathop {\arg \min }\limits_{{\omega _k}} \left\{ {\int_0^\infty {{{\left( {\omega - {\omega _k}} \right)}^2}} \left| {\mathop {{\hat{u_k}}\left( \omega \right)}\limits^ \wedge } \right|d\omega } \right\} ，$

(4)

$ {\omega _k} $ 的更新方式如下式：

$ \omega _k^{n + 1} = \frac{{\displaystyle\int_0^\infty {\omega {{\left| {{{\mathop u\limits^ \wedge }_k}\left( \omega \right)} \right|}^2}{\rm{d}}\omega } }}{{\displaystyle\int_0^\infty {{{\left| {{{\mathop u\limits^ \wedge }_k}\left( \omega \right)} \right|}^2}{\rm{d}}\omega } }} ，$

(5)

$\lambda $ 的更新方式如下式：

$ {\lambda ^{n + 1}} \leftarrow {\lambda ^n} + \tau \left( {f - \sum\limits_k {u_k^{n + 1}} } \right) ，$

(6)

迭代收敛的准则如下式： $\tau$ 为与拉格朗日乘子相关的一个参数，可以人为设置。

$ {{\sum\nolimits_k {\left\| {u_k^{n + 1} - u_k^n} \right\|} _2^2} \mathord{\left/ {\vphantom {{\sum\nolimits_k {\left\| {u_k^{n + 1} - u_k^n} \right\|} _2^2} {\left\| {u_k^n} \right\|}}} \right. } {\left\| {u_k^n} \right\|}}_2^2 < \varepsilon 。$

(7)

其中， $\varepsilon $ 为给定的误差值。

算法的主要流程为：初始化本征模态函数 ${u_k}$ 、中心频率 $ {\omega _k} $ ，拉格朗日乘子 $\lambda $ ，确定模态个数 $k$ 、惩罚因子 $\alpha $ ，然后按式(3)更新 ${u_k}$ ，按照式(5)更新 $ {\omega _k} $ ，按照式(6)更新 $\lambda $ 。最后，依据收敛准则(7)检查结果是否满足条件，若不满足输出条件则继续迭代直至迭代出满足该收敛条件的各个本征模态函数。算法的流程图如图1所示。

2 变分模态分解信号仿真

仿真信号由3个信号加上噪声组成： ${y_1} = 4 \times \sin (2 \times {\text{π}} \times 40 \times t)$ ， ${y_1}$ 表示初相位为 $0$ 、频率为40的正弦信号。 ${y_2} = 4 \times \sin (2 \times {\text{π}} \times 100 \times \left( {t + {t^2}} \right))$ ， ${y_2}$ 表示初相位为 $0$ 的频率随时间变化的信号， ${y_3} = 2 \times \sin (2 \times {\text{π}} \times 500 \times t) + 3 \times \sin (2 \times {\text{π}} \times 800 \times t)$ ， ${y_3}$ 由初相位为 $0$ 的频率分别为500和800的信号组成。 ${y_4} = 0.1 \times r\left( {size\left( {{y_1}} \right)} \right)$ ， ${y_4}$ 表示为噪声信号。合成信号 $ y = {y_1} + {y_2} + {y_3} + {y_4} $ ，其在时域内的图像如图2所示。

采用快速傅里叶变换方法对仿真信号 $y$ 进行变换，所得到的频域图形如图3所示。

采用VMD方法对仿真信号进行分解，由图3可得仿真时模态数量 $k = 4$ ，初始化 ${u_k}$ 和 ${\omega _k}$ 等参数，按式(3)、式(5)、式(6)对参数进行更新，检查迭代后的参数后代入式(7)最终得到的各模态函数如图4所示。

对分解后得到的各本征模态函数与原仿真信号做对比，imf₁和imf₂对应于分量信号 ${y_3}$ 。Imf₃与分量信号 ${y_2}$ 对应。Imf₄与分量信号 ${y_1}$ 对应。对分解后的信号进行合成与原仿真信号相比误差值图5所示，最大误差为 $3\% $ ，对应的时间点为0.00975。

基于VMD方法对本征模态函数做HHT变换得到的时频分布如图6所示。

3 实验测试

实验系统组成如图7所示。实验时在保证其他条件一样的情况下，助推系统动力源提供的动力大小分别设置为3 MPa，5 MPa，7 MPa。

实验获取了旋转机械1、减速器、旋转机械2、航行器器壁处振动信号。由于助推系统安装在航行器上，连接处的振动会经过航行器向外传递，实验时也测取了连接处的振动信号。

4 振动信号时频分布

旋转机械1振动信号时频分布如图8所示。

减速器振动信号时频分布如图9所示。

旋转机械2振动信号时频分布如图10所示。

助推系统与航行器连接处振动信号时频分布如图11所示。

5 时频分布结果分析

旋转机械1振动频率分布范围较广，主要集中在0～f₆ Hz之间，低频成分主要集中在f₁ Hz以下。频率较强的成分分布情况表1所示，从表中可以看出在所给出的3种工况下旋转机械1振动信号频率分布范围基本一致，较强的频率成分分布大体一致，相差较小。

减速器振动信号频率主要集中在f₂～f₃ Hz之间，振动频率相对集中，频率较强的成分分布如表2所示。工况1和工况3下减速器振动信号频率分布范围较为接近，较强的频率成分分布也较为一致。工况2下减速器振动频率分布范围较广，主要分布在f₂～f₄之间，较强的频率成分值比工况1和工况3高。

旋转机械2振动信号频率主要分布在f₂～f₄之间，频率成分较强的分布在f₃～f₄之间，工况2下对应的最强振动频率值最大，不存在f₁以下的低频信号。旋转机械2振动信号较强频率成分分布如表3所示。

助推系统与航行器连接处1振动信号较强频率成分分布如表4所示，较强的频率成分主要分布在f₄附近。从时频分布图中可以看出，工况压力与较强的频率值之间没有严格的对应关系，即助推时压力增大对应的最强的频率成分值不一定增大。

助推系统与航行器连接处2振动信号较强频率成分如表5所示。振动频率主要集中在0～f₆之间，频率分布范围较广，存在较多的低频振动信号。7 MPa时的低频振动信号较5 MPa时更多。

6 振动特性分析 6.1 振动源分析

由于低频振动衰减较慢，对振动源分析析时主要关心低频振动。从图8以及图11(b)可看出，f₁以下的低频振动主要来源于旋转机械1以及航行器器壁振动。

旋转机械1等流体机械的振动主要有转子的不平衡带来的振动以及流致振动。流致振动主要由流体与壁面的结构耦合、回流、脱流、空化、流动失稳、流体对壁面的冲击等造成。从低频振动信号分布时间来看，转子转速较高，其造成振动的频率也较高，由此可知低频部分主要是由流致振动引起的。航行器器壁振动是低频振动的又一来源，主要由流体的冲击以及振动传递造成。

对f₂～f₄的中频振动，各子系统对该频段的振动均有贡献，但在该频段内各系统较强成分的有差别。对f₅以上的高频振动，振动源主要为旋转机械1以及航行器器壁。

6.2 传递特性分分析

实际系统中旋转机械1和减速器连接，从图9可看出，在减速器上未检测到相应的低频信号，主要是由于旋转机械1低频振动较弱。减速器上未检测到f₅以上的高频信号，主要与高频信号的衰减速度有关。

从图11可以看出，连接处的低频振动主要由航行器器壁传递而来。其中7 MPa下低频振动最强。连接处高频振动一方面由旋转机械2传递而来，另一方面航行器器壁振动传递而来。

7 结　语

本文首先基于变分模态分解和希尔伯特变换得到了助推系统中旋转机械1、减速器、旋转机械2、航行器连接处振动信号的时频分布情况。以振动信号时频分布情况为依据，进一步分析了助推系统的振动源构成情况、振动传递特性。

旋转机械1的振动频率最强的成分均分布在f₃附近，最强振动频率之间差距较小在60 Hz内。旋转机械1是低频振动的重要来源，但是旋转机械1低频振动较弱，传播距离有限。减速器振动最强的频率成分集中在f₃附近，工况1和工况3的最强振动频率值相差较小。

旋转机械2振动最强的频率成分集中在f₄附近，不同工况之间相差50～591 Hz。旋转机械2与航行器连接处振动频率较强的成分集中在f₃～f₄附近。不同工况之间相差281～846 Hz。连接处2振动频率最强的成分主要集中在f₃～f₄之间，器壁的振动是低频振动的另一个来源。结合旋转机械1和器壁低频振动的传递路径上的测点振动信号来看，器壁的低频振动强度比旋转机械1的强。

同一子系统不同工况下对应的最强振动频率值有差别，但是总体上差别不大，最强的振动频率值差值在10～846 Hz之间。对同一子系统，助推压力增大时对应的最强振动频率值不一定增大。从振动源构成的角度讲，低频振动控制将重点放在旋转机械1和航行器器壁。从振动传递特性角度讲，需将减振降噪设计应用于如连接处这样振动传递效率较高的地方。