0 引　言

水下航行器辐射噪声信号研究涉及领域广泛，它是水声对抗的基础，也是被动式声探测和识别设备的主要信号源。在水声领域的研究过程中需要大量信号样本来进行特性的对比、统计、测试。然而实测水下航行器噪声需要场地，经费，人力物力以及时间成本，且参数不可控，很难获得实验所需要的大量信号数据。因此，水下航行器噪声信号的模拟仿真具有重要的现实意义和作战意义。

目前，对于水下航行器辐射噪声的生成主要方法是分析船舶噪声的功率谱特性，对照实测信号的谱特征，对噪声功率谱进行构造^[1]。从水下航行器辐射噪声功率谱中各噪声分量的产生机理出发，将水下航行器噪声信号拆分成线谱，连续谱以及调制谱的仿真^[2]。对于连续谱信号，将高斯白噪声通过设计好的FIR滤波器获得符合功率谱形状的连续谱^[3-4]。调制谱主要通过demon谱分析获得。何正耀等^[5]对信号进行功率谱分析和demon谱分析，研究了舰船辐射噪声功率谱的声学特征。

由于水声信道的作用，信号在从声源到达接收点的过程中经过不同路径传播，从而产生衰减，相位变化以及幅频变化。为了准确获取信号所携带的信息，需要对水声信道的传播特性展开研究。涂峰^[6]分析了Bellhop的数学原理，利用Matlab对水声信道的本征路径模型进行了仿真。李庆来^[7]基于射线理论对分层非均匀浅海声场进行数值预报。孙博等^[8]和丁浩等^[9]利用虚源法建立浅海水声信道模型，对于信道的多普勒频移采用的不同的数学模型进行仿真。

为了开发极地地区的丰富资源，近年来，多个国家开展了对极地区域的研究。由于极地地区特殊的水文环境，加上冰面的覆盖，形成了极地特殊的水下声场环境。因此对冰面展开研究是必要的。李贞刚^[10]对光滑界面和粗糙界面的冰面不同特性进行了阐述，杨威^[11]和明平寿^[12]都基于射线理论结合B-T冰水界面模型进行了冰下声传播仿真计算。

本文主要研究水下航行器辐射噪声的建模以及信号经过海洋信道传播后的声学特性。在水下航行器辐射噪声模拟方面，本文采用任意形状频率响应滤波器设计方法，设计FIR滤波器，同时用经验公式模拟调制谱分量以及线谱分量，使仿真获得的噪声信号具有较高逼真度。然后分别利用Bellhop方法和虚源法模拟水声信道传播后的水下航行器辐射噪声，2种方法计算结果基本一致。运用B-T冰水模型计算出冰面反射系数，导入到Bellhop声学工具箱中进行冰下水声信道仿真，开展了辐射噪声信号在冰下水声信道中的传播特性研究。

1 水下航行器辐射噪声特性

水下航行器的辐射噪声由螺旋桨噪声、机械噪声以及水动力噪声组成。螺旋桨噪声产生于航行器外，当航速较低时，螺旋桨未发生空化现象，此时机械噪声最高，其次是螺旋桨噪声；当航速较高时，会发生空泡现象，此时螺旋桨噪声就会成为水下航行器噪声的最主要成分。机械噪声在水下航行器噪声的低频段表现为单频线谱，其频率与机械设备的基频和谐频有关。水动力噪声是指水下航行器与周边流体发生相对运动时，流体作用于外壳产生的噪声，随着航速的升高水动力噪声淹没到螺旋桨噪声中。

2 水下航行器辐射噪声及水下信道仿真方法

基于上文对水下航行器噪声各噪声成分的分析，将噪声模拟分为对连续谱噪声，线谱噪声以及调制谱噪声的模拟。水下航行器辐射噪声，其时域波形可描述为：

$ S(t) = [1 + a(t)] \times {g_x}(t) + {g_L}(t)。$

(1)

式中：$ a(t) $为调制谱时域波形；$ {g_x}(t) $为连续谱分量的时域波形；$ {g_L}(t) $为线谱时域波形。

2.1 连续谱信号模拟方法

螺旋桨连续谱噪声的功率谱级在低频段随着频率增加而增加，在高频段以−6 dB/Oct衰减，会在10～1000 Hz之间形成一个峰值。本文对水下航行器连续谱噪声的仿真基本思路为构造一个与期望的辐射噪声的频率响应和相位响应一致的FIR滤波器，将符合高斯分布的白噪声通过该滤波器，获得与期望噪声谱形状相同的连续谱曲线。本文构建FIR2滤波器（构成简单，原理清晰）。对于水下航行器的连续谱，只需确定谱峰频率$ {f_0} $，功率谱上升斜率$ {a_1} $(dB/Oct)，功率谱下降斜率$ {a_2} $(dB/Oct)，功率谱峰值声压谱级$ S{L_0} $^[13]，就能构建出完整的连续谱曲线模型如下式：

$ SL = \left\{ \begin{gathered} S{L_0} - {a_1}\lg \left(\frac{{{f_0}}}{f}\right),\quad {f_L} \leqslant f \leqslant {f_0} ，\\ S{L_0},\quad f = {f_0} ，\\ S{L_0} + {a_2}\lg \left(\frac{f}{{{f_0}}}\right),\quad {f_0} \leqslant f \leqslant {f_H} 。\\ \end{gathered} \right. $

(2)

式中：$ {f_0} $为谱峰频率；$ S{L_0} $为峰值声压谱级；$ {a_1} $为功率谱上升斜率；$ {a_2} $为功率谱下降斜率。

对船舶连续谱噪声进行估计后，调用Matlab的FIR2函数设计任意形状频率响应的FIR滤波器，提取其中的谱峰频率$ {f_0} $，峰值声压谱级$ S{L_0} $，其余频率区间内的谱级按倍频程取样，作为滤波器幅频响应曲线的特征，其中理想频率以$ f/{f_{\max }} $来表示，在峰值频率的两侧按倍频程取样。理想幅值为$ {10^{SL/20}}/{m_{\max }} $。将白噪声通过按上述方法形成的滤波器，就能得到水下航行器噪声连续谱。

2.2 线谱信号模拟方法

线谱分量主要集中在低频部分，它是能量集中的体现，具有很高的辨识度。可以采用周期性正弦信号来进行仿真，其数学表达式如下：

$ {g_L}(t) = \sum\limits_{k = 1}^N {{A_k}\sin \left(2{\text π} {f_k}t + {\varphi _k}\right)}。$

(3)

式中：$ {A_k} $为线谱分量的幅值；$ {f_k} $为信号频率；$ {\varphi _k} $为相位。

线谱频率为$ {f_k} = m \times n \times s $，其中$ m $为谐波次数，$ n $为螺旋桨叶片数，$ s $为螺旋桨转速。根据统计学实测，线谱一遍比连续谱高出10～25 dB。

2.3 调制谱信号模拟方法

水下航行器辐射噪声拥有鲜明的节奏感，这是因为螺旋桨在水中旋转一周，每个叶片产生的噪声对连续信号具有调制作用。使得叶频整数倍上的频谱强度远大于轴频整数倍的频谱强度。这种调制信号可以通过DEMON谱分析提取，目前大多学者认为调制包络是具有随机幅度，相同形状，相同周期的脉冲性随机过程，数学模型为：

$\begin{split} a(t) =&\sum\limits_{k = 1}^N {{A_1}{e^{ - pk}}\sin (2{\text π} k{f_1}t + {\varphi _1})} +\\ &\sum\limits_{i = 1}^M {{A_2}{e^{ - qi}}\sin (2{\text π} i{f_2}t + {\varphi _2})} 。\end{split}$

(4)

式中：$ {f_1} $为轴频；$ {f_2} $为叶频；$ {\varphi _1} $及$ {\varphi _2} $为随机起伏相位；q和p为衰减系数；$ i{f_2} $为叶频谐波频率；$ k{f_1} $为非叶频及叶频谐波处的轴频谐波频率。

2.4 基于Bellhop模型的信道模拟方法

为了研究水声信道对于接收信号产生的干扰以及畸变，本文利用模拟获得的水下航行器辐射噪声与Bellhop模型相结合；首先配置Bellhop的海洋环境参数，设置海水表面为真空，插值方式选择线性插值；水底环境选择为声学半空间，平坦地形；计算信号的输出幅度和到达时间，以及传播路径。声速剖面采用典型深海Munk声速剖面。

同时，为了证明仿真的可靠性，利用文献[9]的虚源法和Xiao模型进行水声信道数值模拟。Xiao模型能反映Rice信道的各种物理性质，其接收信号数学模型为：

$ r(t) = \sum\limits_{m = 1}^M {{\zeta _m}(t)x(t - {\tau _m})}，$

(5)

$ {\zeta _m}(t) = {u_{cm}}(t) + j{u_{sm}}(t) ，$

(6)

${u_{cm}}(t) =\frac{{{C_m} \left[ {\sqrt K \cos ( 2{\text π} {f_d}t\cos {\theta_m} + {\varphi _m} ) + \displaystyle\frac{1}{{\sqrt N }} \displaystyle\sum\limits_{n = 1}^N {\cos ( 2{\text π} {f_d}t \cos {\alpha _n} + {\varphi _n} )} } \right]}}{{\sqrt {1 + K} }} ，$

(7)

$ {u_{sm}}(t) = \frac{{{C_m}\left[ {\sqrt K \sin (2{\text π} {f_d}t\cos {\theta _m} + {\varphi _m}) + \displaystyle\frac{1}{{\sqrt N }}\displaystyle\sum\limits_{n = 1}^N {\sin (2{\text π} {f_d}t\cos {\alpha _n} + {\varphi _n})} } \right]}}{{\sqrt {1 + K} }} 。$

(8)

式中：$ {C_m} $为第m条本征声线的衰减系数；$ {\theta _m} $为第m条本征声线的入射角；$ {\varphi _m} $为第m条本征声线的相位。

虚源法把每一根声线等效为一个虚源所射出的直达声线，接收点处的信号便是这些声线的总和。其中主要分量称为本征声线。各本征声线的衰减系数为：

$ \left\{ {\begin{array}{*{20}{l}} {{c_{n1}} = \displaystyle\frac{{V_1^nV_{2n1}^n}}{{{R_{n1}}}}} ，\\ {{c_{n2}} = \displaystyle\frac{{V_1^nV_{2n2}^{n + 1}}}{{{R_{n2}}}}}，\\ {{c_{n3}} = \displaystyle\frac{{V_1^{n + 1}V_{2n3}^n}}{{{R_{n3}}}}} ，\\ {{c_{n4}} = \displaystyle\frac{{V_1^{n + 1}V_{2n4}^{n + 1}}}{{{R_{n4}}}}} 。\end{array}} \right. $

(9)

式中：$ {V_1} $为海面反射系数；$ {V_{2ni}} $为海底反射系数。

$ \left\{ {\begin{array}{*{20}{l}} {{V_1} = - \sqrt {\displaystyle\frac{{1 + 0.1{{({f_0}/{f_1})}^2}}}{{1 + {{({f_0}/{f_1})}^2}}}} } ，\\ {{V_{2ni}} = \displaystyle\frac{{a\cos {\theta _{ni}} - \sqrt {{b^2} - {{\sin }^2}{\theta _{ni}}} }}{{a\cos {\theta _{ni}} + \sqrt {{b^2} - {{\sin }^2}{\theta _{ni}}} }}} ，\\ {{f_1} = \displaystyle\frac{{378}}{{{w^2}}}} ，\\ {{\theta _{ni}} = \arcsin \displaystyle\frac{{{R_0}}}{{{R_{ni}}}}} 。\end{array}} \right. $

(10)

式中：w为风速；a为海底密度与海水密度的比值；b为海水和海底折射率比。

各本征声线的路径长度为：

$ \left\{ {\begin{array}{*{20}{l}} {{R_{n1}} = \sqrt {{r^2} + {{\left(2nH + z - {z_0}\right)}^2}} } ，\\ {{R_{n2}} = \sqrt {{r^2} + {{(2nH + z + {z_0})}^2}} }，\\ {{R_{n3}} = \sqrt {{r^2} + {{[(2n + 2)H - z - {z_0}]}^2}} }，\\ {{R_{n4}} = \sqrt {{r^2} + {{[(2n + 2)H - z + {z_0}]}^2}} } 。\end{array}} \right. $

(11)

式中：z为声源深度；$ {z_0} $为接收器深度；H为海深；n为反射次数。

冰水界面是研究极地环境声学绕不开的一部分，海冰的大小，形状，分布通常是不规律的。冰水界面分为光滑界面与粗糙界面，粗糙界面呈犬齿状，冰的下表面分布不均匀，也有冰脊的情况。本文通过编写冰面文件导入Bellhop模型，即可对海冰环境下的声传播进行研究。

为了获得不规则海冰的反射系数，Burke和Twersky提出了B-T冰层声反射模型^[14]。模型中将冰层下表面看作自由或刚性表面下，随机分布了许多半椭圆凸起。$ \eta $和$ \xi $分别为半椭圆凸起的长半轴和短半轴，冰脊尺寸为$ \rho = \eta /\xi $，当$ \rho = 1 $时冰脊为半圆形，当$ \rho = 0 $时冰脊垂直，当$ \rho \to \infty $时为平面。反射系数R如下：

$ R = {\left| {\frac{{1 + Z}}{{1 - Z}}} \right|^2} = 1 + 4{Re} \left( {\frac{Z}{{{{\left| {1 - Z} \right|}^2}}}} \right) ，$

(12)

$ Z = \frac{n}{{k\cos {\varphi _0}}}{A_f}。$

(13)

式中：n为单位长度内含有半椭圆凸起的平均个数；$ {\varphi _0} $为入射角；k为波数。

半椭圆的长轴表示为$ a = \max (\xi ,\eta ) $，当$ ka \ll 1 $时，即$ 2{\text π} a/\lambda \ll 1 $时，冰脊相对声波波长为小尺寸冰脊，Z由式(14)计算获得；当$ ka \gg 1 $时，即$ 2{\text π} a/\lambda \gg 1 $时，冰脊相对声波波长为大尺寸冰脊，Z由式(15)计算获得。

$\begin{split} Z =& - \frac{n}{{\cos {\varphi _0}}}\Biggr[\frac{{{{\text π} ^2}{k^3}}}{{16}}(2{\xi ^2}{\eta ^2} + {\xi ^2}{(\xi + \eta )^2}{\sin ^2}{\varphi _0}) +\\ & {{i}}\frac{{{\text π} k}}{2}(\xi \eta - \xi (\xi + \eta ){\sin ^2}{\varphi _0})\Biggr] ，\end{split}$

(14)

$\begin{split} Z =& - \frac{{nl}}{{\cos {\varphi _0}}}\left[\cos (\varphi 0 - \gamma ) + \frac{{\tan {\varphi _0}\sin (2kl\sin \gamma \cos {\varphi _0})}}{{2kl}}\right] +\\ &\frac{{n\eta }}{2}{\left(\frac{{{{i}}{\text π} }}{{k\xi \cos {\varphi _0}}}\right)^{1/2}}{e^{ - {{i2k}}\xi {\text{cos}}{\varphi _0}}} 。\\[-10pt] \end{split}$

(15)

式中，$ \gamma $和l的函数式如下：

$ \tan \gamma = {\rho ^{ - 2}}\tan {\varphi _0} ，$

(16)

$ {l^2} = {\xi ^2}(1 + {\rho ^4}{\cot ^2}{\varphi _0}){(1 + {\rho ^2}{\cot ^2}{\varphi _0})^{ - 1}}。$

(17)

3 仿真结果及分析 3.1 海洋信道下水下航行器辐射噪声传播仿真

连续谱模拟：由Ross经验公式计算获得式(2)中的各参数，以国外某型号水下航行器公开的参数进行仿真。螺旋桨直径$ D = 1.6\ {\text{m}} $；进速系数$ {J_p} = 0.47 $；伴流分数$ \omega = 0.21 $；螺旋桨转速为245 r/min；空泡评价$ {\sigma _i} = 3.6 $；螺旋桨为双四叶螺旋桨。计算获得的各参数如表1所示。将特征参数代入FIR滤波器的设计，白噪声信号时域谱如图1所示，将白噪声通过设计出的滤波器，得到连续谱如图2所示。

线谱模拟：本文设定线谱频率为16、32、48、64 Hz，对于随机线谱按照一般情况在1000 Hz以内设计2条。获得的线谱时域信号如图3所示。

调制谱模拟：根据轴频及叶频，调制谱选择为4.1、8.2、12.3、28.7 Hz，频域信号如图4所示

最终利用本节数学模型将3种信号分量整合，获得完整的水下航行器辐射噪声，如图5所示。

水声信道中声传播模拟：对深度为4000 m的海域进行Bellhop声场仿真计算，输入信号为第二章模拟的水下航行器辐射噪声信号，海底为声弹性半空间。海水中的声速梯度采用典型海洋声速梯度接收器深度为60 m，声源深度设置为200 m。计算其接收信号的功率谱，如图6所示。

利用基于虚源法的数值计算方法，计算相同海洋环境下的接收信号，与Bellhop的结果进行对比验证，结果如图7所示，从图中可以看出，接收信号衰减后的幅度基本一致，信号谱峰的趋势也一致，在高频段幅值有些许差距，这是因为虚源法模型中无法精准的考虑声速梯度，且其中包含了随机相位，入射角的选取也与Bellhop后台程序有所差异。验证了本文模拟的可靠性。

3.2 B-T冰水模型信道仿真

设冰脊参数为：$ n = 4/ {\rm{km}} $；$ \xi = 50\ {\rm m} $；$ \eta = 100\ {\rm m} $，结合B-T模型计算，获得反射系数的近似。反射系数如图8所示。

可知，冰水界面的反射系数在掠射角1°～10°的区间内从1逐渐趋近于0，并在10°～12°之间获得最低值，最终趋近于0.4。

3.3 冰下声信道仿真

将3.1节获得的水下航行器辐射噪声信号作为输入信号。分别设置声源深度为60、200和500 m，海深4000 m；水平距离100 km；接收器深度为60 m。声速剖面采用北极典型极地声速剖面。结合Bellhop模型计算声线轨迹，传播损失和经过信道传播的接收点信号。

将3.2节计算获得的反射系数，角度和相位写入到Bellhop中的.trc文件，该文件主要包含了上边界的反射系数。将冰层下表面的半椭圆起伏写入到.ati文件，改变上边界的形状。

由图9可知，在声源处于同一深度的情况下，有冰界面的声线比无冰界面的声线稀疏。这是由于冰水上边界的存在使声线反射次数增加，传播损失也随之增加，到达远方的声线数量减少。由于极地海域声速剖面整体呈上升趋势，所以声线都向上弯曲。

从图10可以看出，无冰和有冰界面的轮廓几乎一致。无冰界面的声线在深海区域能量较低，而有冰界面中，声线由于不规则冰面的反射，声线能量能够到达深海区域；从水平方向来看，前50 km有冰界面的声传播损失增长较慢，但在后50 km传播损失增长幅度大于无冰界面。观察可知，声源位于500 m时的浅海声传播损失小于声源位于60 m和声源位于200 m时的声传播损失，这是因为随着声源深度的增加，受到冰水上边界的影响减弱。但传播损失的增长幅度更大，在深海处到达100 dB。

从图11可知，无冰情况下，接收点信号的功率谱集中在20～30 dB，声源深度的变化对接收信号幅度产生了5 dB左右的变化。这与传播轨迹的变化最后到达接收点的声线不同有关。同时，发现有冰界面下，接收点信号功率谱降低20～30 dB，功率谱初始幅度都集中在0 dB，峰值为10 dB。这是因为冰面的存在增加了声线的传播路径和反射次数，到达接收点的声线数量比无冰情况少。

4 结　语

本文首先对水下航行器辐射噪声的连续谱、线谱、调制谱进行建模，之后将模拟获得的辐射噪声信号与Bellhop模型相结合，研究海洋信道对声传播的调制特征。结合B-T冰水模型模拟冰下信道，计算分析了不同声源深度的接收点声信号特性。结果表明：

1）经过海洋信道传播后的接收点信号功率谱与声源处的辐射噪声功率谱趋势一致，幅值根据水平传播距离的变化产生不同程度的降低；

2）由于冰水界面的反射系数较低、声传播损失较大，有冰界面的声线比无冰界面的声线稀疏，导致部分声线在传播途中消失，声线无法远距离传播；

3）由于冰面的不规则界面对声线的多角度反射，声线被散射到海底，使得声能量扩散在海洋空间中；有冰界面下的声传播损失的增长幅度大于无冰界面下声传播损失的增长幅度，在远距离时传播损失能够到达100 dB；

4）经过水平距离为100 km的无冰界面信道传播后，辐射噪声功率谱幅值下降约100 dB，当冰水界面存在时，增加了声线的传播路径和反射次数，功率谱幅值将额外降低20～30 dB。