0 引　言

在水下航行器减振降噪工程中，如何准确地在机械噪声系统中识别出起主要作用的噪声源，并从频率和空间域上查明航行器各主要噪声源的空间分布和频率特征，定量识别各主要噪声源对航行器辐射噪声的贡献大小，从而有针对性地采取有效减振降噪措施，成为安静型水下航行器研制的一项关键技术。工程上圆柱壳是一种常见壳体结构，水下航行器的主体结构可近似为圆柱壳体结构，例如潜艇舱段和鱼雷。壳体振动辐射噪声是水下航行器的主要噪声源，为提高水下舰艇的作战和隐蔽性能，圆柱壳体的减振降噪成为一项长期研究的重要工程^[1]。Xia等^[2]对充满空气刚性连接的圆柱壳的空气中噪声传播和辐射进行了建模，辐射源表征为单极子，讨论了不同厚度和长度的圆柱形壳体对远场辐射声压的影响，以及不同激励位置对辐射声场的影响。Lee等^[3]通过使用几何关系将刚性无限障板中的方形活塞的归一化辐射阻抗表示为四重表面积分，并以数值方法对辐射阻抗做了对比。Greenspon等^[4]针对无限长刚性圆柱壳计算得到了圆柱壳上矩形活塞的辐射声压和辐射阻抗表达式，并且针对典型情况得到了数值结果，对圆柱壳上活塞的相互辐射阻抗进行比较。胡世猛等^[5]对多源激励下的有限弹性圆柱壳辐射声场进行了分析，在激励力源总大小不变的条件下，得到了将激励源分散为多个激励源可以降低辐射噪声的结论。刘涛^[6]对水下有限长弹性圆柱壳共振进行了理论分析，从Donnell薄壳理论出发，推导出了水下简支圆柱壳对称激励和非对称激励下的辐射声场公式，分析了对称模态和非对称模态的特性，讨论了考虑互辐射阻抗和不考虑互辐射阻抗对圆柱壳辐射声场的影响。石焕文等^[7]基于圆柱壳受迫振动方程推导了声辐射阻抗表达式，利用数值法计算了不同频率下圆柱壳的辐射声场指向性，研究了结构损耗因子对结构辐射声功率和辐射效率的影响。

本文提出一种基于振声传递矩阵的广义逆波束形成方法^{[8 − 9]}，根据规则障板表面活塞辐射声场的解析表达式建立目标结构表面声源与接收基阵之间的声传递矩阵，代替传统点源球面波扩展传递矩阵，克服了传统点源球面扩展声传播模型对于结构噪声源识别精度不高的局限性，解决了复杂结构声源声传播模型的精细化表征问题；利用基于振声传递矩阵的广义逆波束形成方法实现了对圆柱壳结构噪声源的识别^{[10 − 12]}；理论分析与仿真结果表明，该方法具有更高的空间分辨率，有效提高圆柱壳结构声源的识别定位性能。

1 圆柱壳建模及声辐射特性 1.1 圆柱壳测量模型

圆柱壳测量模型如图1所示，水下有限长弹性圆柱壳的两端镶嵌于半无限长刚性圆柱障板中，边界条件简支（任意时刻法向、切向位移恒为0，且法向力矩、切向斜率恒为0）。在笛卡尔坐标系中，圆柱壳轴线与Z轴重合，壳体底端中心与坐标系原点重合。阵列面采用平面阵，阵元数为M，壳面均匀划分为N个表面共形活塞，第n个结构表面活塞的几何中心坐标为$ \left( {{x_n},{y_n},{z_n}} \right) $，第m个阵元的空间坐标位置为$ \left( {{x_m},{y_m},{z_m}} \right) $。共形活塞面到平面阵的垂直距离为$ {r_{nm}} $。

1.2 辐射声功率、均方振速及辐射效率

结构声辐射特性的研究对于结构辐射噪声控制、结构声学性能优化有着重要的意义。在声学领域，通常用辐射阻抗、辐射声功率、均方振速以及声辐射效率等特性参数来描述结构的辐射特性。辐射声功率是表征声源特性的重要物理量，反映了声源向外辐射声能的能力大小，与声波传播距离无关，其定义为声源在单位时间内向外辐射的声能，一般不能直接测量得到。除辐射声功率外，声辐射效率也是衡量结构声辐射性能的重要参数。反映了声源向外辐射声能的效率，与一般效率定义类似，声辐射效率定义为声源在单位时间向外辐射的声能与结构振动总能量的比。根据辐射声功率、均方振速及辐射效率的定义，通过推导可得有限长弹性圆柱壳的辐射声功率公式：

$ P\left(\omega\right)=\frac{S}{4}{\mathrm{Re}}\left\{\sum\limits_{\alpha=0}^1\sum\limits_{n=0}^{\infty}\sum\limits_{m=1}^{\infty}\frac{1}{\varepsilon_n}\dot{W}_{nm}^{\alpha}Z_{nmm}\left(\dot{W}_{nm}^{\alpha}\right)^*\right\}。$

(1)

式中：S为圆柱壳外表面积；$ {\varepsilon _n} $为Neumann函数；Re为取实部；*为共轭。

表面均方振速$ < v_n^2 > $：

$ < v_n^2 > =\frac{1}{4}\sum\limits_{\alpha=0}^1\sum\limits_{n=0}^{\infty}\sum\limits_{m=1}^{\infty}\frac{1}{\varepsilon_n}\dot{W}_{nm}^{\alpha}\left(\dot{W}_{nm}^{\alpha}\right)^*。$

(2)

辐射效率$ \sigma \left( \omega \right) $：

$ \sigma(\omega)=\frac{Re\left(\displaystyle\sum\limits_{\alpha=0}^1\displaystyle\sum\limits_{n=0}^{\infty}\displaystyle\sum\limits_{m=1}^{\infty}\frac{1}{\varepsilon_n}\dot{W}_{nm}^{\alpha}Z_{nmm}\left(\dot{W}_{nm}^{\alpha}\right)^*\right)}{\rho c\displaystyle\sum\limits_{\alpha=0}^1\displaystyle\sum\limits_{n=0}^{\infty}\displaystyle\sum\limits_{m=1}^{\infty}\frac{1}{\varepsilon_n}\dot{W}_{nm}^{\alpha}\left(\dot{W}_{nm}^{\alpha}\right)^*}。$

(3)

式中：$\rho $、$c$分别为介质密度和声速。

2 基于传递矩阵的波束形成方法 2.1 广义逆波束形成原理

广义逆波束形成技术是一种基于阵列接收信号处理的噪声源识别定位技术。它通过阵列接收数据处理以及声场传递矩阵逆向求解算法，可以得到噪声源声压分布。广义逆波束形成算法首先要建立声源到接收阵元的传递矩阵，利用传递矩阵直接正则化求逆与阵列接收声压向量，计算得到声源表面声压分布。

对于圆柱壳结构声源，首先将结构辐射面划分为很多均匀规则共形活塞，将每个活塞看作一个独立辐射声源。圆柱壳结构声源测量模型如图1所示，第$ n $个活塞到第$ m $个阵元的空间距离为 ${r_{mn}} = \sqrt {{{\left( {{x_m} - {x_n}} \right)}^2} + {{\left( {{y_m} - {y_n}} \right)}^2} + {{\left( {{z_m} - {z_n}} \right)}^2}} $，根据广义逆波束形成基本理论，建立声源表面活塞到阵元的传递函数，通过方程可准确描述声场空间信息，可得：

$ \left\{\begin{array}{l}\left[\boldsymbol{A}\right]\cdot x=\boldsymbol{p}，\\ \left[\boldsymbol{B}\right]\cdot v=\boldsymbol{p}。\end{array}\right. $

(4)

式中：$ \left[ {\boldsymbol{A}} \right] $为辐射面到阵列面的点源格林函数传递矩阵；$ \left[ {\boldsymbol{B}} \right] $为辐射面到阵列面的单元辐射法振声传递矩阵；$ \left[ x \right] $为声源表面声压；$ \left[ v \right] $为声源表面法向振速矢量；$ {\boldsymbol{p}} $为M元阵列接收声压。M、N分别为阵元个数和结构辐射面离散活塞个数。为达到更好的识别效果，往往要求辐射面离散活塞尺寸尽可能小，因此活塞个数M远大于阵元个数N。

由于活塞个数大于阵元个数，矩阵$ \left[ {\boldsymbol{A}} \right] $、$ \left[ {\boldsymbol{B}} \right] $的行数小于列数，式(4)为欠定方程组，不能对矩阵$ \left[ {\boldsymbol{A}} \right] $、$ \left[ {\boldsymbol{B}} \right] $直接求逆。在这种情况下，声学响应传递方程是具有不适定性的，无法通过直接对传递矩阵求逆得出满意的声源面成像结果。为此，需要采用结合正则化方法的广义逆方法在消除这种不适定性的同时近似求解上述方程。经典的正则化方法是Tikhonov正则化方法，该方法是基于上式逆问题解的残余范数和解范数之间的联合加权以达到最小的思想，即：

$ q=\text{arg min}\left\{ \left\| p-Gq \right\| _2^2+\varepsilon^2\Omega\left(q\right)^2\right\}。$

(5)

式中：$ {\boldsymbol{q}} $为声源表面物理量；$ \Omega \left( q \right) = {\left\| {{L^{{\text{ - }}1}}q} \right\|_2} $，$\Omega (q)$为离散平滑范数；$ {L^{ - 1}} $为正则化方阵，在传统的广义逆波束形成方法中，是令$ {L^{ - 1}} = {I_N} $，通过最小二乘法结合Tikhonov正则化方法可得解最终声源表面声压矢量和声源法向振速矢量：

$ \left\{\begin{gathered} x={\boldsymbol{A}}^{\mathbf{H}}\left(\boldsymbol{A}\boldsymbol{A}^{\mathbf{H}}\mathbf{+}\varepsilon\boldsymbol{I}\right)^{\mathbf{\mathrm{-1}}}p，\\ v={\boldsymbol{B}}^{\mathbf{H}}\left(\boldsymbol{B}\boldsymbol{B}^{\mathbf{H}}\mathbf{+}\varepsilon\boldsymbol{I}\right)^{\mathbf{\mathrm{-1}}}p 。\end{gathered} \right. \\$

(6)

式中：$ {\text{H}} $表示矩阵的共轭转置；$ {\mathbf{I}} $为对角单位矩阵；$ \varepsilon $表示Tikhonov正则化参数；对于基于广义逆波束形成的声源识别问题，其经验值为$ {\boldsymbol{A}}{{\boldsymbol{A}}^{\mathbf{H}}} $和$ {\boldsymbol{B}}{{\boldsymbol{B}}^{\mathbf{H}}} $的最大特征值1%~10%。通过式(6)就可以重构声源表面声压$ {\boldsymbol{p}} $和表面振速$ v $。

2.2 建立传递矩阵

单元辐射叠加法^[10]建立了结构表面振动到辐射声场之间的传递模型。因此，只需得到声场传递模型中的声场传递函数，再结合结构表面法向振速分布，就可得到辐射声场中任意一点的声压值。声场传递函数决定于与辐射面共形的刚性活塞面辐射声场。根据点源球面波扩展模型，构建自由场格林函数传递矩阵$ \left[ {\boldsymbol{A}} \right] $，第$ \left( {m,n} \right) $个元素的表达式为：

$ {{\boldsymbol{A}}_{mn}} = \frac{{{e^{ - jk{r_{mn}}}}}}{{4{\text{π}} {r_{mn}}}}。$

(7)

根据圆柱障板表面共形活塞辐射声场解析式，构建振声传递矩阵$ \left[ {\boldsymbol{B}} \right] $，第$ \left( {m,n} \right) $个元素的表达式为：

$ \begin{split} {{bm{B}}_{mn}} =& \frac{{2\rho cv{\alpha _0}{L_0}}}{{{{\text{π}} ^2}{r_{mn}}\sin \theta }}{{\text{j}}_0}\left( {k{L_0}\cos \theta } \right){e^{{\text{j}}k\left( {{r_{mn}} - {L_1}\cos \theta } \right)}}\cdot \\ & \sum\limits_{m = - \infty }^\infty {\frac{{{{\left( { - {\text{j}}} \right)}^m}{{\text{j}}_0}\left( {m{\alpha _0}} \right)}}{{{{\text{H}}_m}^\prime \left( {ka\sin \theta } \right)}}} {e^{{\text{j}}m\left( {\phi - {\alpha _1}} \right)}} 。\end{split} $

(8)

式中：a为共形活塞所在圆柱壳周向半径，2$ {L_0} $为活塞轴向长度，2$ {\alpha _0} $为活塞周向夹角；$ {{\boldsymbol{H}}_m}^\prime $表示第一类m阶汉克尔函数的导数；在球坐标系下，原点为共形活塞所在圆柱壳的轴线中心点，极轴与Z轴平行；$ {L_1} $为活塞中心距离直角坐标系原点的轴向距离，$ {\alpha _1} $为活塞中心的水平方位角。

为了比较2种方法识别误差，定义重构误差：

$ \mathrm{error}\text{ = }\frac{ \left\| \mathit{{P}_{{re}}}-\mathit{{P}_{{th}}} \right\| _2}{\left|\mathit{|{P}_{{th}}} \right\| _2}\times100\text{% }。$

(9)

式中：$ \mathit{{P}_{{re}}} $为重构声压；$ \mathit{{P}_{{th}}} $为对应的理论声压；$ {\Vert \cdot\Vert }_{2} $表示矩阵2-范数。

2.3 结构噪声源识别流程

利用GFTM-GIB和VATM-GIB方法，对圆柱壳结构噪声源进行识别。图2为结构噪声源识别示意图，由于GFTM-GIB方法传递矩阵为结构表面声压到声场声压的传递矩阵，VATM-GIB方法为结构表面法向振速和声场声压的传递矩阵，为了方便比较，选取距离结构表面很近的重构面，识别重构面声压分布。首先，利用解析法计算测量阵的辐射声压数据。其次，将测量数据代入广义逆波束形成器，其中GFTM-GIB方法可以直接得到重构面声压分布，VATM-GIB方法先计算得到结构表面法向振速，再利用结构表面法向振速计算重构面声压分布。最后，将这2种方法得到的重构面声压分布与理论值对比，并分析2种方法的识别性能。

3 仿真分析

仿真中选取由钢材料制成的简支圆柱壳作为目标声源，圆柱壳体长度L=2 m，内半径a=0.6 m，厚度为0.006 m，共形活塞的尺寸为$ 0.05 \times 0.015\;{{\text{m}}^2} $（一般要求活塞的尺寸要小于计算波长的1/12，尺寸越小计算结果越精确），满足Neumann边界条件，圆柱壳表面周向和轴向划分为40个单元。钢板的弹性模量为$ 2.1 \times {10^{11}}\;{\text{N}}/{{\text{m}}^2} $，泊松比为0.3，取圆柱壳底面几何中心为坐标原点，壳体表面(a, 0, L/2)处施加简谐激励力，激励钢板向外辐射噪声，激励力幅值为1 N，测量面的大小为$ 2 \times 2{\text{ }}{{\text{m}}^2} $，阵元数目为$ 11 \times 11 $，阵元间距为0.2 m，测量距离为1 m，重构面距离壳面$ 0.01{\text{ m}} $，水中的声速为$ 1\;500{\text{ m/s}} $，水的密度为$ 1\;000{\text{ kg/}}{{\text{m}}^{\text{3}}} $，参考声压为$ 1 \times {10^{ - 6}}{\text{ Pa}} $。

振动模态阶数m、n分别取30阶，频率f取10～800 Hz，利用单元辐射叠加法，仿真圆柱壳在空气和水2种介质中，其结构辐射声功率、均方振速以及辐射效率的频响曲线，仿真结果如图3所示。其中尖峰代表频率达到了模态的本征频率，发生了共振现象，200 Hz在共振频率附近，400 Hz为非共振频率。可以看出，圆柱壳在水中的各阶共振频率相比空气中较低。在非共振频率处，圆柱壳在水中的辐射声功率基本高于空气，而水中的均方振速低于空气。

图4和图5分别为200、400 Hz时圆柱壳表面法向振速分布。通过对比可以看出，VATM-GIB方法计算的表面法向振速分布与理论分布基本一致，说明该方法可以实现对圆柱壳表面法向振速的重构。

图6和图7分别为200、400 Hz时重构面声压分布。通过对比可以看出，VATM-GIB方法可以实现对圆柱壳表面声压的识别重构，而GFTM-GIB方法不能完全识别重构，VATM-GIB方法较GFTM-GIB方法提高了结构噪声源识别定位的空间分辨率。

表1为频率分别为200、400 Hz时，2种方法声压重构误差，可以看出VATM-GIB方法比GFTM-GIB方法重构误差小。

4 结　语

本文基于广义逆波束形成方法，利用平面阵研究了基于传递矩阵的圆柱壳结构噪声源识别，克服了传统点源球面扩展声传播模型对于结构噪声源识别精度不高的局限性，解决了复杂结构声源声传播模型的精细化表征问题，理论仿真表明该方法对圆柱壳结构声源识别具有较好的空间分辨率，实现了对圆柱壳结构表面法向振速和辐射声压的重构，提高了圆柱壳结构声源的识别定位性能。本文可为圆柱壳结构噪声源减振降噪研究提供一定参考，具有较好的应用前景。