CN107576388B - 一种浅海信道下三维结构声源辐射声场预报方法 - Google Patents

Abstract

本发明公开了一种浅海信道下三维结构声源辐射声场预报方法，属于声源辐射声场预报领域，其特征是首先建立浅海信道下结构声辐射多物理场耦合模型，根据分析频段、计算距离和浅海环境等因素选择海洋声传播理论建立二维轴对称下信道传输函数G，通过方位耦合解拓展为信道三维声场传输函数，然后采用浅海信道多边界声场耦合计算理论建立多物理场耦合模型，计算获取浅海信道下受多边界影响的结构声源表面速度U，最后利用Fourier波叠加法理论预报信道下结构声源在任意场点声场信息P。该算法有效解决了目前在研究复杂浅海信道下三维弹性结构声源辐射声场特性时所遇到的计算量大、耦合物理场多和信道环境复杂等诸多瓶颈性问题。

Description

一种浅海信道下三维结构声源辐射声场预报方法

技术领域

本发明涉及声源辐射声场预报领域，具体涉及一种浅海信道下三维结构声源辐射声场预 报方法。

背景技术

研究弹性结构在浅海信道下的声辐射特性，对开展潜器的辐射噪声实时预报和噪声有效 控制具有极其重要的理论和军事价值，是以后我国水声技术领域长期关注的热点和难点问题 之一。

然而目前关于水下弹性结构多物理场耦合振动与声辐射问题大多考虑为无界或半空间流 体域下近距离声场问题，对于浅海信道下弹性结构多物理场耦合声辐射问题的研究尚不多见。 因为传统数值法(边界元法、有限元法、统计能量法等)将受大网格计算量和复杂海洋环境 和多物理场耦合等因素的严重限制，无法开展相关研究工作；解析解法只能分析部分简单结 构在简单信道环境下的声场问题，对复杂浅海信道下任意弹性结构的声辐射研究也无能为力。 为了解决在研究水下大型复杂结构声源辐射声场所涉及的瓶颈性问题，Koop-mann引入了声 场等效计算的波叠加法，因其适应性强和计算效率高的优势逐步地成为了一种研究弹性结构 声辐射行之有效的新方法，近年来受到了国内学者极大的关注，并广泛应用自由场或半空间 环境下的近场声全息、噪声源识别和结构声辐射预报等研究领域。但对于信道下复杂三维结 构声源辐射波叠加法研究却很少涉及，因为信道下结构声场问题涉及流固耦合、声壳耦合甚 至声与复杂海底边界的声边界耦合等复杂多物理场耦合环境，而且浅海信道下更关注远距离 声辐射问题，在这些耦合条件下采用理论解析解推导难以求解，采用传统数值法(有限元、 边界元和统计能量法等)因计算量大而无法建立数值模型，导致浅海信道下结构表面振动信 息难以获取。为了避免复杂结构和浅海环境的多物理场理论推导问题，也有部分学者提出直 接忽略结构与流体、结构与环境的耦合作用，把结构视为点声源，但这样直接忽略了结构的 近场声辐射特性。以上共同导致了无法从理论解和数值法的角度有效地开展信道下弹性结构 声辐射研究，但它对浅海中水下结构声辐射预报具有极为重要的意义，急需探索一种新的研 究方法来解决浅海信道中结构声辐射预报问题。

发明内容

本发明为开展浅海信道环境下任意弹性结构声源辐射声场预报提供了一种新方法，有效 地解决传统数值法和解析法在研究浅海信道下复杂弹性结构声源辐射声场预报时所遇到的计 算量大、耦合物理场多和信道环境复杂等诸多瓶颈性问题。该发明通过采用信道传输函数、 波叠加法和多边界耦合声场计算理论联合的方法计算预报了浅海信道下弹性结构声源的辐射 声场，该发明方法具有对信道环境和结构适应性强、声场计算稳定性好、计算结果准确性高 和计算效率高等优势，为今后开展浅海信道下弹性结构声源辐射声场预报提供了一种新的研 究途径。

本发明通过联合信道传输函数、波叠加法和多边界耦合声场计算理论进行信道下结构声源 辐射声场特性研究，如图1所示，其特点主要包括以下几个步骤：

一种浅海信道下三维结构声源辐射声场预报方法，其特征在于,包含以下步骤：

步骤一利用简正波理论求解等声速浅海信道下传输函数的表达式：

其中，r、z分别为轴对称坐标下水平方向和竖直方向，z_s为声源位置，z为信道任意场 点位置；

为汉克函数，n为简正波模态数，ξ_n为水平行波方向上的特征值；Z_n(z)为 竖直方向特征值函数也称为格林函数，其解形式由海面、海底边界决定；

步骤二对于低频声场问题采取简正波理论求解三维浅海信道下信道传输函数G；通过格 林函数求解边界条件，格林函数Z_n(z)解形式与海面、海底边界决定，海面通常为Dirichlet 边界，满足的边界条件为：

Z_n(z)|_z＝0＝0 (2)

在各种海底模型交互面上，各模型满足声学边界方程为：

其中，

m₁为海水ρ₁与海底密度ρ₂比，即

k₂为海底波速定义 为

c₂、α₂为海底声速和吸收系数；V为海底反射系数，对于复杂海底界面反射系数可采用“三参数”界面反射模型获取；由确定的海底边界便可确定特征函数(格林函数)Z_n(z)，根据方程(1)便可得到与信道格林函数相关的传输函数；

步骤三信道下点源声场传输函数空间维度变换，把二维浅海信道模型下与格林函数相关 的传输函数变换为三维条件下声场函数，可采取方位耦合解进行空间维度变换；

步骤四通过波叠加法在结构内部设置虚拟源面S、，并在源面上合理布放若干个单极源， 采用浅海多边界声场耦合理论计算获取信道下结构表面受界面影响的振动速度U，通过信道 传输函数G建立起结构内部虚拟源强Q与结构表面振速U的传递矩阵D，满足以下关系：

{U}_M×1＝{D}_M×N·{Q}_N×1 (5)

其中，G为轴对称条件下基于简正波理论推导与格林函数相关的信道传输函数表达式；r_s为结构表面任意一点的方向矢量，r₀为结构内部任意虚拟源的方向矢量；

为法向求导算子；

步骤五根据虚拟源源强Q、表面振速U和传递矩阵D的关系，则虚拟源源强Q则可表示 为：

[Q]＝[D]^-1[U] (7)

步骤六对病态传递函数[D]进行奇异值分解：

其中U，V分布为M阶和N阶方阵，Σ＝diag(σ₁,σ₂,σ₃...σ_n)为包含n个非奇异值σ_i的对 角阵，且σ₁≥σ₂≥σ₃≥...≥σ_n≥0，T为矩阵转置运算；

步骤七采用Tikhonov正则化方法对病态传递矩阵[D]进行处理，利用奇异值分解有 D^TD＝VΣ²V^T，则源强表达式为：

其中，U^δ为结构表面速度矩阵，α为正则化参数；

步骤八利用L曲线进行最优参数α的选取；

步骤九获取虚拟源源强后，便可根据以下公式求解任意场点的声压：

P＝TQ＝TV(Σ²+αI)^-1ΣU^TU^δ (10)

其中，T为虚拟源Q与信道任意场点的单极矩阵，与浅海信道环境下对应的信道声场传 输函数关系为T＝jρωG(|r-r₀|)，j为虚数，ω＝2πf为角频率，f为频率；ρ为海水密度。

所述步骤三的方位耦合解为考虑水平环境变化的影响，采取绝热模式二维图案进行变换， 三维亥姆霍兹方程为：

其中，

为求导算子，ρ(x,y,z)和c(x,y,z)为浅海信道水层密度和声速，p(x,y,z)空间声场 声压，x_s为声源在x轴上的固定位置，ω为角频率，δ(y)为特殊函数；

设解的形式即三维信道传输函数G(x,y,z)为：

其中，Γ_m(x,y,z)为局部模式，m为模数，E_m(y,z)为二维声场函数；

其中，Γ_m(x,y,z)为局部模式，m为模数，E_m(y,z)为二维声场函数；

Γ_m(x,y,z)为局部模式，通过绝热近似，可得到水平折射方程：

利用局部简正波把三维问题消去了x维变为二维，得到一个关于(y，z)新的亥姆霍兹 方程，水平方向的有效折射率可由

控制，根据简正波理论求解式(13)所表示的二维 声场，然后代入式子(12)可求解三维声场传输函数。

所述步骤四的获取结构表面速度U，其特征在于，首次采用采用浅海多边界声场耦合理 论计算获取结构表面速度U，充分考虑流体、结构和边界之间的相互耦合作用，突破多边界 耦合声场环境下结构声振信息获取难题，计算获取结构表面准确的速度信息U，其计算理论 包含以下几个部分：

流固耦合方程，根据Helmholtz方程结合流固耦合边界条件得出结构与流体的耦合有限 元方程为：

其中，F_st、F_at分别为结构、流体介质声的耦合激励载荷,Kc、Mc分别为耦合刚度矩阵、 耦合质量矩阵，且

PML(PerfectlyMethodedLayer，简称PML)技术，进行信道四周无限大边界模拟，在PML吸收层域控制方程：

海面Dirichlet边界，满足的边界条件为：

p_an(x,y,z)|_z＝0＝0 (16)

海底边界，其声与海底耦合建模方法主要有Neumann边界模型，Rayleigh定律，Sommerfeld模型，Cauchy边界和Geoacoustics模型，大致可分为硬海底、液态海底、弹性 海底和多孔海底等，为了区别海底与海水域，下标a和b分别标示流体层和海底层；在各类 海底交互面(z＝H)上建立流体介质与海底物理模型的耦合方程：对于硬质海底，满足：

对于液态海底，满足声压和法向振速连续：

p_a＝p_b(18)

对于各向同性弹性海底，建立声-弹性波耦合方程，即法向上满足位移连续和应力连续， 切向的应力为零，即：

其中w_b、u_b分别表示海底水平方向(r)的位移和垂直方向(z)的位移，λ、μ为拉米 常数，Δ的定义为

对于多孔海底，所涉及声学物理参数较多，建立声-多孔弹 性波耦合方程，满足边界条件为牵引力连续、流体压力连续和流体流量连续。

通过流固耦合方程、PML控制方程和浅海多边界耦合方程联立求解，便可计算局部浅海 信道流体下弹性结构表面受界面和环境影响的准确振动信息，以作为波叠加计算源强Q求解 的输入条件。此外，源强Q的求解为声学反问题，应对传递矩阵D择优选择正则化工具进行 去病态处理后，便可计算获取结构在场点准确的辐射声场，提高辐射声场计算精度和稳定性。

与已有的技术相比，本发明的优势和创新体现在：

1、该发明对浅海信道环境类型和弹性结构声源适应能力强，只要确定信道环境(主要是 确定浅海信道海底边界类型)，便可求解对应信道的声场传输函数G；传输函数G求解理论可 根据信道环境进行灵活调整，当信道环境参数随距离或深度变化、分析频段不同时，可选择 合适的海洋声传播理论对相应信道下的传输函数进行求解，以达到使整个计算过程具有极高 的准确性和灵活性。对于任意形状、尺寸的弹性结构声声源，均可通过试验法和多边界耦合 声场理论获取在考虑浅海信道边界影响作用下的表面振动数据U，进而开展各种浅海信道环 境下任意弹性结构的声辐射特性准确计算。

2、该发明突破了多边界耦合声场计算难题，不同于以往浅海信道下声场研究多为点源或 者忽略多边界影响以及自由场和半空间环境下结构声振问题，信道下结构声振将受信道上下 界面的严重影响，具体表现为界面反射声不仅作为二次激励而且对总声场具有干涉作用，尤 其是海底边界类型多、参数复杂、影响显著，加大了进行信道下结构表面声振信息获取的难 度和工作量。该发明方法为了考虑界面反射声影响，针对信道下低频声辐射，有限元网格量 相对较少，采用有限元法建立多边界下耦合声场计算理论，计算获取结构声振信息，充分考 虑了流体、结构和边界之间的相互作用，所以计算获取多种耦合条件下结构表面速度是准确 的，浅海信道下场点声场信息整个计算过程均充分考虑了界面的影响，而且该发明不受计算 距离限制，该发明突破了多边界耦合声场远距离声辐射计算的难点问题。

3、该发明针对浅海信道下三维结构声源辐射声场计算问题，解决了目前在研究三维信道 声场时所遇到的计算量大，环境复杂和耦合声场多等突出难题。结合了理论解析解法和数值 有限元法建立了贴近实际的三维浅海信道下有限弹性结构声辐射问题，采取三维空间变换技 术对二位轴对称下的信道声场传输函数G进行三维空间变换，进行浅海信道下三维声场信息 的准确获取，以利用波叠加法进行三维信道下弹性结构声辐射预报研究。

4、该发明在进行浅海信道下大型结构声源辐射声场预报时具有计算用时少、效率高、计 算稳定性好的特点，在波叠加计算了过程中，通过对结构离散点和虚拟源进行最优配置，可 以在实际应用中大大降低测点和虚拟源的数目，减少操作工作量，使得在实际工程中更加简 便、高效；采用最小二乘法(LMS)和最小均方误差法(HELS)进行虚拟源的最优配置，达到 在满足所需精度要求下极大提高计算效率，减少在实际应用中的工作量，简化了该发明在实 际应用中的操作过程，进一步提高预报效率和准确性。该发明方法可以快速的计算传统数值 法在相同计算硬件条件下无法模拟或计算复杂声辐射问题，而且可快速准确计算信道下结构 在距离、深度和空间的声场分布特性以及辐射声功率。采用正则化工具曲线对计算过程涉及 的离散非适定性问题进行处理，在低信噪比条件下计算精度明显优于传统处理方法，使得该 发明计算稳定性增强，计算精度提高。

附图说明

图1为本发明信道下弹性结构辐射声场计算方法流程图；

图2为本发明包含的浅海信道下声场传输函数模型；

图3为本发明包含的浅海信道下弹性结构场点辐射声场计算示意图；

图4为本发明包含的浅海信道下多边界耦合声场计算理论模型；

图5为浅海信道下脉动球壳声辐射模型；

图6为采用本发明方法计算信道下脉动球壳辐射声压级随距离变化曲线与理论解的比较 图，频率为50Hz；

图7为采用本发明方法计算信道下脉动球壳辐射声压级随距离变化曲线与理论解的比较 图，频率为60Hz；

图8为采用本发明方法计算信道下脉动球壳辐射声压级随距离变化曲线与信道下脉动球 壳有限元法计算结果的比较图，频率为40Hz；

图9为采用本发明方法计算信道下脉动球壳辐射声压级随距离变化曲线与信道下脉动球 壳有限元法计算结果的比较图，频率为60Hz；

图10为采用本发明方法计算信道下脉动球壳辐射声压级随距离变化曲线与信道下脉动 球壳有限元法计算结果的比较图，频率为80Hz；

图11为采用本发明方法计算信道下脉动球壳辐射声压级随距离变化曲线与信道下脉动 球壳有限元法计算结果的比较图，频率为100Hz；

图12为采用本发明方法计算信道下振动小球辐射声压级随距离变化曲线与信道下脉动 球壳有限元法计算结果的比较图，频率为60Hz；

图13为采用本发明方法计算信道下振动小球辐射声压级随距离变化曲线与信道下脉动 球壳有限元法计算结果的比较图，频率为80Hz；

图14为浅海信道下弹性圆柱壳声辐射模型；

图15为浅海信道下脉动球壳表面数据不包含噪声下本发明方法计算的声辐射曲线与理 论解对比图，频率为80Hz；

图16为浅海信道下脉动球壳表面数据不包含噪声下本发明方法计算的声辐射曲线与理 论解对比图，频率为80Hz；

图17为信噪比10dB条件下通过常规求解源强预报的声辐射曲线与理论解在辐射距离上 的相对误差曲线，频率为80Hz；

图18为信噪比10dB条件下通过本发明方法包含的正则化处理求解源强预报的声辐射曲 线与理论解在辐射距离上的相对误差曲线，频率为80Hz；

图19为信噪比20dB条件下通过常规求解源强预报的声辐射曲线与理论解在辐射距离上 的相对误差曲线，频率为80Hz；

图20为信噪比20dB条件下通过本发明方法包含的正则化处理求解源强预报的声辐射曲 线与理论解在辐射距离上的相对误差曲线，频率为80Hz；

图21为信噪比40dB条件下通过常规求解源强计算预报的声辐射曲线与理论解在辐射距 离上的相对误差曲线，频率为80Hz；

图22为信噪比40dB条件下通过本发明方法包含的正则化处理求解源强计算预报的声辐 射曲线与理论解在辐射距离上的相对误差曲线，频率为80Hz；

图23为信噪比80dB条件下通过常规求解源强计算的声辐射曲线与理论解在辐射距离上 的相对误差曲线，频率为80Hz；

图24为信噪比80dB条件下通过本发明方法包含的正则化处理求解源强计算的声辐射曲 线与理论解在辐射距离上的相对误差曲线，频率为80Hz；

具体实施方式

下面结合理论模型和算例具体说明本发明的实施方式，参见图2，该发明计算理论步骤 简要为：

步骤a、首先建立浅海信道下多物理场耦合模型，通过海洋声传播理论求解二维轴对称 坐标下信道传输函数G，并采用空间变换技术进行传输函数G的三维空间变换。

步骤b、建立浅海信道下多边界多物理场耦合声场计算模型，通过流固耦合、声边界耦 合和PML技术考虑结构、流体和浅海边界的相互作用，获取浅海信道环境下影响下的结构表 面真实的振动速度U。

步骤c、信道下结构辐射声场波叠加法计算分为声学“反”问题和声学“正”问题两部 分：通过步骤a求解和空间变换得到的三维信道传输函数G建立虚拟源强Q与结构表面速度 U的传递矩阵D，并采用正则化方法进行声学“反”问题源强Q的准确求解；声学“正”问题，通过信道三维传输函数G建立源强Q与流体域场点声压P的单极矩阵T，然后便可计算获取任意场点的声压信息P，然后对P进行声学运算，便可获取浅海信道下结构辐射声功率、空间声场指向性分布和辐射声传播损失等声场信息。

本发明方法的应用领域为浅海信道环境下三维结构声源辐射声场问题，本发明的基本理 论思想是弹性结构在浅海信道下任意场点的辐射声场P可以由置于结构内部一系列虚拟源产 生的声场线性叠加进行等效替代，不同于以往自由场条件下的声场等效，本发明方法通过求 解信道传输函数以及空间变换技术、多边界耦合声场计算理论和浅海下波叠加声场快速计算 等首次建立了三维信道环境下弹性结构辐射声场计算理论。该发明与传统分析方法相比，该 发明方法避免了在研究信道下结构声辐射时所遇到的有限元大量网格划分而无法计算、边界 元奇异积分计算结果失效以及统计能量法浅海模型环境无法建立等瓶颈性难题，而且通过合 理优化等效源配置和正则化处理，进一步提高了该发明方法的计算效率、稳定性和精度。改 发明方法的理论模型主要由信道传输函数、信道波叠加法和多边界耦合声场计算理论三部分 组成，具体理论细节如下：

1.信道传输函数

建立如图2所示的浅海信道模型，把信道下点源在任意场点的声场传输函数G(r,z)表示 成n个简正模式之和，令

在二维轴对称条件下各方向上满足的方程为：

Z_n(z)＝Asin(k_zz)+Bcos(k_zz) (1)

其中，Z_n(z)为竖直方向特征值函数也称为格林函数，A、B为待定常数，由边界条件和 归一化条件确定；R_n(r)为水平方向的特征值函数，

为汉克函数，k_z和ξ_n分别为竖 直方向和水平方向的特征值，满足关系为

格林函数Z_n(z)解形式与海面、海底边界决定，海面通常为Dirichlet边界，满足的边界 条件为：

Z_n(z)|_z＝0＝0 (3)

海底是一个反射和散射边界，典型海底边界模型有Neumann边界模型，Rayleigh定律， Sommerfeld模型，Cauchy边界和地声模型等。在各种海底模型交互面上，各模型满足声学边 界方程为：

其中，

m₁为海水ρ₁与海底密度ρ₂比，即

k₂为海底波速定义 为

c₂、α₂为海底声速和吸收系数；g(ξ_n)由各海底模型声学特性决定，V为海底反射系数，对于复杂海底界面反射系数可采用“三参数”界面反射模型获取，由各种确定的海底边界便可确定特征函数(格林函数)Z_n(z)，根据方程(2)求解出函数R_n(r)，然 后带入方程

便可得到信道传输函数。

2.信道下波叠加法

波叠加法的基本理论是结构体的辐射声场可以由置于辐射体内部等效源产生的声场叠加 代替，如图4所示，信道下波叠加计算的基本思路为通过数值计算域计算获取多边界声场耦 合环境下结构表面速度U，然后应用波叠加法结合信道格林函数G计算结构辐射声场P。振动 结构体置于均质流体密度为ρ₁，流体声速为c₁。振动体的表面记为S'，n为表面的外法向矢 量。P为外部流体中的一点。Ω为结构内部连续分布的等效声源体，其虚拟边界S作为与结 构共形的等效源面，Q为S上一点。Q和P间的距离记为r。经过离散处理，可由N个简单声 源可以构建复杂振源的表面某一点的法向速度u_n为：

式中，q_i是虚拟面上第i小段的源强，r_s为S结构表面一点的位移矢量，r_0i为S'虚拟面各点对应的位移矢量，G(|r_s-r_0i|)为信道下与格林函数相关的信道传输函数，

为求导 算子。

把上式写成矩阵形式：

{U}_M×1＝{D}_M×N·{Q}_N×1 (7)

结构表面法向速度矩阵U是已知的，D为源点和场点的传递函数，一旦配置好等效源， 传递函数

可以通过相应浅海信道环境下的传输函数得到，便可计算虚拟源 的源强矩阵Q：

[Q]＝[D]^-1[U] (8)

其中，T为单极矩阵，与浅海信道环境下对应的声场传输函数成正比即T＝jρωG(|r-r₀|)。 求解结构内部虚拟源强度Q首先需要通过数值法或实验法获取结构在信道环境下的表面振 动信息，其过程属于声学反问题，在波叠加法数值计算中属于离散不适定问题，即解非唯一 性或不连续依赖于结构表面测点数据，小的测点误差带来计算结果的巨大振荡导致数值计算 失稳，但可以通过Tikhonov正则化以及L曲线后验参数选取法来消除不适定问题，以提高计 算精度。

首先对病态传递函数[D]进行奇异值分解：

其中U和V分布为M阶和N阶方阵，Σ＝diag(σ₁,σ₂,σ₃...σ_n)为包含n个非奇异值σ_i的对 角阵，且σ₁≥σ₂≥σ₃≥...≥σ_n≥0。

采用Tikhonov正则化方法进行传递矩阵[D]病态处理，利用奇异值分解有 D^TD＝VΣ²V^T，则式子(8)的正则化方程为：

其中，U^δ为结构表面测点速度矩阵，α为正则化参数，采用Hansen等人提出的L曲线准则来预估正则化参数α。

写出源强矩阵的平滑范数和残余范数的表达式分别为：

绘出平滑范数||Q^α,δ||₂对残余范数||DQ^α,δ-U^δ||₂的单调递减对数曲线，该曲线称为L曲线。 L曲线清楚地反映了最小化残余范数和最小化解的“大小”之间的一种平衡关系，其水平部 分为过正则化起重要作用，垂直部分为欠正则化起重要作用。在L曲线由竖直到水平的转角 处为最佳平衡关系，因为在此处解有最小的残值余量和比较少的震荡项。

通过L曲线选择最优正则化参数，然后利用Tikhonov正则化对传递矩阵处理获取辐射体 结构内部一系列的虚拟源源强后，便可计算任意一点的声压：

P＝TQ＝TV(Σ²+αI)^-1ΣU^TU^δ (13)

3.多边界耦合声场理论

从信道下弹性结构波叠加计算理论可知，声场计算精度取决于源强求解的准确性程度， 而源强求解是由结构表面速度U进行计算获取，所以结构表面速度获取的准确性将决定整个 声场计算过程的精度。不同于自由场和半空间环境下结构声阵问题，信道下结构声振将受信 道上下界面的严重影响，具体表现为界面反射声不仅作为二次激励而且对总声场具有干涉作 用，尤其是海底边界类型多、参数复杂、影响显著，加大了进行信道下结构表面声振信息获 取的难度和工作量。为了考虑界面反射声影响，可采用数值法和试验法获取结构表面的振动 数据，因试验法成本高，且针对信道下低频声辐射，有限元网格量相对较少，所以采用有限 元法建立多边界下结构声辐射局域数值模型是可行的，本节将充分考虑流体、结构和边界之 间的相互作用，数值计算获取多种耦合条件下结构的表面速度。

建立如图5所示的浅海信道下弹性结构声辐射多边界声场耦合计算模型，其中Ω_s为弹性 结构内部流体介质域，Ω_a为浅海海水介质域，四周Ω_p为完美匹配层(PerfectlyMethoded Layer简称PML)域，PML宽度为h_PML，内外边界为Γ_I、Γ_O，Γ为流固耦合边界，Γ_l、Γ_h分别 为声场与浅海海面与海底交互作用边界。根据对应边界下的连续条件，建立结构与流体，结 构与边界声耦合作用方程。

在稳定简谐激励海水流体介质域下的频域声场表达式即Helmholtz方程为：

其中，k＝ω/c＝2πf/c为波数；ω＝2πf为角频率，f为频率(Hz)；

为拉普拉斯算子，有

ρ_a为海水密度。

对Helmholtz方程进行权重积分并结合高斯理论，可写出声学有限元方程为：

(K_a+jωC_a-ω²M_a){p_i}＝{F_i} (15)

其中，{F_i}为声学激励，质量矩阵M_ij、刚度矩阵K_ij和阻尼矩阵C_ij均为n×n阶矩阵。下 标a为声学系统矩阵。类似声学有限元方程推导，对于弹性结构，其有限元振动方程为：

(K_s+jωC_s-ω²M_s){u_i}＝{F_si} (16)

其中，M_s、K_s和C_s分别为结构网格上没有受到约束(位移u_i)部分的刚度矩阵、质量矩 阵和阻尼矩阵；{F_si}为结构上的激励载荷。

在流固耦合面Γ上，结构表面法向速度v_sn与外部流体介质法向速度v_an连续：

则联立式15、式16和式17，可写出结构与流体的耦合方程为

其中，F_st、F_at分别为结构、流体介质声的耦合激励载荷，K_c、M_c分别为耦合刚度矩阵、 耦合质量矩阵，且

不同于自由场下结构声辐射数值计算，信道下结构声场计算涉及辐射声与多边界交互作 用的复杂问题。本发明将采用PML技术和边界耦合条件对浅海信道四周边界和上下界面进行 模拟，充分考虑边界对结构声辐射的影响作用。信道四周无限大边界(xoy平面)采用PML (PerfectlyMethodedLayer，简称PML)技术处理，在PML吸收层域Ω_p的方程与有限元求 解域Ω_a满足的方程不同，PML通过对控制方程增加吸收系数转换为吸收层的控制方程，为了 简化方程描述，令x轴为x₁轴、y轴为x₂轴，利用分离变量可写出频域下的PML方程：

其中，σ_i为与空间变量x_i相关的吸收系数，v_i，p_i为匹配层域的速度和声压

联立式(19)-(20)，便可求得：

计算流体域包括介质域Ω_a和PML域Ω_p，两者不同在于PML域有声吸收σ_i，且 σ_i(x)＝σ_i(-x)，当|x_i|＜X_i，σ_i(x)＝0，当X_i＜|x_i|＜X_i'，σ_i(x)有值，引入空间吸收系数：

PML方程与Helmholtz方程可统一写为：

PML的作用是在边界上增加几层网格来吸收声学量，其厚度h_PML需大于最大波长λ_max的十 五分之一即h_PML＞λ_max/15，通过式(23)对PML网格层进行声吸收，最终在PML内外网格层 达到的目的为：

p_a＝p_ponΓ_I (24)

v_an＝v_pnonΓ_I (25)

p_p＝0onΓ_O (26)

海面平整边界通常为Dirichlet边界，满足的边界条件为：

p_an(x,y,z)|_z＝0＝0 on Γ_H (27)

海底边界大致可分为硬海底、液态海底、弹性海底和多孔海底等，为了区别海底与海水 域，下标a和b分别标示流体层和海底层。在各类海底交互面(z＝H)上建立流体介质与海底 物理模型的耦合方程：对于最典型硬质海底，满足：

对于液态海底，满足声压和法向振速连续：

p_a＝p_b on Γ_L (29)

对于各向同性弹性海底，建立声-弹性波耦合方程，即法向上满足位移连续和应力连续， 切向的应力为零，即：

其中w_b、u_b分别表示海底水平方向(r)的位移和垂直方向(z)的位移，λ、μ为拉米 常数，Δ的定义为

对于多孔海底，所涉及声学物理参数较多，建立声-多孔弹 性波耦合方程，满足边界条件为牵引力连续、流体压力连续和流体流量连续。

通过有限元数值法建立的流固耦合方程、多边界耦合方程和PML控制方程，便可计算局 部浅海信道流体环境下下弹性结构表面准确的振动信息，以作为信道下结构远距离声辐射波 叠加计算源强求解的输入条件，便可计算获取结构在场点准确的辐射声场。

以下将从具体算例进行该发明方法的准确性验证、计算高效性说明和低信噪比条件下声 场计算的稳定性测试，也进一步说明该发明方法的实施方式、理论和实际应用价值。

1.发明方法准确性检验

算例1：浅海信道下脉动球壳声辐射验证

建立如图6所示有限元模型，信道深度为H＝50m，四周为无限大域(有限元软件采用 PML模拟)；上边界为Dirichlet软边界，下边界为Neumann边界，点源在信道下的传输函数 满足的边界条件为：

Z_n(z)|_z＝0＝0 (34)

求解点源在该信道下场点的信道声场函数为：

其中，其中，r、z分别为轴对称坐标下水平方向和竖直方向，z₀为声源位置，z为信道 任意场点位置；

为汉克函数，n为简正波模态数，ξ_n为水平行波方向上的特征值；

为简正波阶次n的最大可取正整数，k_zn为简正波特征值。

根据脉动球与点源的等效关系为

建立轴对称条件下的信道脉动球壳声辐射 模型，如图6所示，球壳半径r＝3m，球壳中心深度为h₁＝25m，在忽略上下界面反射声对脉 动球壳的二次作用情况下，脉动球壳的辐射声场等效为点源产生的声场。

为了达到一定计算精度，在球壳均匀布放201测点，并提取每个测点的振动数据，内部 虚拟源面与结构面共形且源面半径r₀与结构半径r_s满足r₀≤0.76r_s，虚拟源点数与结构离散点数 之比τ控制在0.6-1范围内，利用式(8)求解源强Q，然后根据式(13)即源强Q与场点单 极矩阵T作用求解场点声场信息。

图7和图8分别用该方法计算了f＝50Hz、f＝60Hz频率下脉动球壳的辐射声压级随距离变 化的曲线(每个距离上场点深度为30m，相连两个场点之间距离1m)，并与等效点源的解析解 进行了对比，对比结果表明在进行浅海信道下脉动球壳声辐射计算中，本发明的信道下结构 声辐射算法与理论解析法式(36)计算结果吻合很好。

利用有限元数值分析软件建立了点源在在本算例所模拟信道下声传播模型，分别计算了 点源在40Hz、60Hz、80Hz和100Hz频率下声压级随距离的曲线，点源深度为25m，场点深度 30m。并与本发明方法计算结果，如图9、图10、图11和图12所示。对比结果表明在进行浅 海信道下脉动球壳声辐射计算中，本发明的信道下结构声辐射算法与有限元法计算结果吻合 很好。

算例2：浅海信道下振动球壳声辐射验证

建立如图6所示有限元模型，模型参数与算例1相同，球壳半径为r＝1.2m，球壳中心深 度为h₁＝25m，材料为钢，顶端施加F_z＝1000N垂直向下的简谐力，结构边界为自由边界。通 过数值法计算、提取振动球壳在该信道下表面的振动数据，结构表面均匀离散120个测点， 虚拟源点数与结构离散点数之比τ控制在0.6-1范围内，内部共形虚拟源面上的单极源的数目 为81个，源面的半径r₀与结构半径r_s满足r₀≤0.76r_s取虚拟源面的半径为0.6m。利用该发明方 法计算给出了频率为60Hz、80Hz声压级随距离(各场点深度为30m，相连两个场点之间距离 1m)的变化曲线，并与有限元计算结果进行对比如图13和图14所示，对比结果表明本文开 发的波叠加法程序与有限元结果基本吻合

2.发明方法高效性检验

算例3：浅海信道下大型圆柱壳声辐射特性计算时间测试

建立如图15所示的浅海信道结构声辐射多物理耦合模型，浅海信道环境与算例1相同。 结构为圆柱壳，半径为a＝3m，长度为l＝30m，圆柱壳几何中心深度为h₁＝25m，在中间施加 径向简谐点力F_r＝1000N，结构边界为自由边界。结构表面255点离散点，虚拟源数目208 点。了说明改发明方法的高效性，先采用了有限元法计算了本算例中所模拟的浅海环境下圆 柱壳声辐射特性，但由于弹性结构和浅海环境尺寸大，在相同场点计算密度下，60Hz频率下 的最大计算距离范围只能为0-100m，而且网格计算量为千万量级，最大计算距离范围上限大 于100m时，CPU＝256GB的大型计算机已经无法计算。但本发明方法是完全可以计算任意距离 范围下的浅海信道结构声辐射特性，如表1所示，对本算例的模型进行了不同距离范围的计 算时间测试，频率为60Hz，该发明方法不仅可以计算任意距离范围下的声辐射特性，而且计 算用时少、效率高，当辐射距离范围取到0-100km条件下，计算时间仅为157min，说明信道 下结构辐射距离完全可以扩展为计算任意距离范围。表2，在相同距离范围下(0-5km)，对 不同分析频率下的计算的时间进行了测试，当频率上升到100Hz时，计算时间为69min，说 明该发明方法可高效快速计算不同频率下的声辐射问题。

表1

表2

3.发明方法稳定性和计算精度检验

算例4：低信噪比条件下浅海信道弹性结构辐射声场计算精度验证

根据式(13)可知，在利用该发明方法进行信道下结构声辐射声场计算时，大体可以分为两 个部分：声学反问题和声学正问题。声学反问题是指利用表面振动信息U进行源强Q的逆运 算，然后利用多个虚拟点源的辐射声场叠加进行声场正问题求解即计算场点声压P大小。所 以源强求解的精度至为重要，将影响整个计算过程，而声学反问题涉及离散不适定问题，即 解非唯一性或不连续依赖于结构表面测点数据。而结构表面振动数据U不可避免包含误差， 这将导致在求解源强过程中出现矩阵病态，严重影响波叠加法声场计算的准确度。

以浅海信道下脉动球壳为例，根据脉动球与点源的等效关系为

建立轴对称 条件下的信道脉动球壳声辐射模型，如图6所示，信道参数与算例1相同；球壳半径r＝3m， 球壳中心深度为h₁＝25m，因为球壳结构尺寸远小于信道尺寸，可忽略上下界面反射声对脉动 球壳的二次作用，则脉动球壳的辐射声场等效为点源产生的声场。根据简正波理论求解该信 道下的信道传输G，便可求解与传输函数相关的信道下点源声场解析解(式(36))，等效为 信道脉动球壳声场解析解。在脉动球壳表面振动数据U不包含误差前提下，通过该发明方法 对信道下脉动球壳声辐射特性进行了计算，结构离散点数和虚拟源配置方式与算例相同，计 算了60Hz频率下辐射声压级随距离的变化曲线，并与信道下脉动球壳声辐射解析解进行了对 比，如图16所示，可以看出，在结构表面振动数据不包含误差前提下，该发明方法计算值与 解析解完全吻合。

为了更好描述振动数据误差对信道下结构辐射声场计算的精度影响程度，定义相对误差 η(n)，平均相对误差λ公式分别为：

其中，p_w(n)为波叠加法计算的声压，p_t(n)为理论解计算的声压，n表示在距离方向上的 场点离散个数，N为场点离散总数。

图17、图18、图19和图20为在脉动球表面振动数据U中分别施加包含噪声信噪比10dB、 20dB、40dB、80dB的高斯白噪声下，通过传统广义伪逆求解源强Q计算结果与理论解在整个 辐射距离的相对误差曲线图，在低信噪比SNR＝10dB的条件下，相对误差曲线在某些传播距离 上已经达到了40％左右，其他信噪比条件下计算结果整体误差也比较高，已经无法满足计算 精度要求。

图21、图22、图23和图24为在表面振动数据U分别包含噪声信噪比10dB、20dB、40dB、 80dB高斯白噪声下通过本发明方法包含的Tikhonov+L曲线正则化方法对噪声进行正则化处 理后，计算结果与理论解在整个辐射距离的相对误差曲线图，可以看出，经过本发明方法包 含的Tikhonov+L曲线正则化方法对振动数据误差进行处理后，明显提高了弹性结构在整个辐 射距离上声辐射计算精度。为了清楚地描述Tikhonov+L曲线正则化方法在低信噪比条件下对 结构声辐射计算稳定性和计算精度的巨大作用，表3列举了在不包含高斯白噪声和包含信噪 比为10dB、20dB、40dB、80dB的高斯白噪声下，经正则化对高斯白噪声进行处理前后，计算 结果与理论解的相对误差大小，可以看出，在包含高斯白噪声和含高信噪比噪声的情况下， 通过传统方法计算结果与采用Tikhonov+L曲线正则化计算精度基本一致。在低信噪比条件 下，利用本发明方法的正则化处理后的计算误差明显小于常规方法求解的声辐射误差。说明 该发明方法能够很好地抑制结构表面振动数据误差对计算结构稳定性的扰动，极大地提高低 信噪比条件下结构声辐射计算的稳定性和精度。

表3

Claims (1)

1.一种浅海信道下三维结构声源辐射声场预报方法，其特征在于,包含以下步骤：

步骤一 利用简正波理论求解等声速浅海信道下传输函数的表达式：

其中，r、z分别为轴对称坐标下水平方向和竖直方向，z_s为声源位置，z为信道任意场点位置；

为汉克函数，n为简正波模态数，ξ_n为水平行波方向上的特征值；Z_n(z)为竖直方向特征值函数也称为格林函数，其解形式由海面、海底边界决定；

步骤二 对于低频声场问题采取简正波理论求解三维浅海信道下信道传输函数G；通过格林函数求解边界条件，格林函数Z_n(z)解形式与海面、海底边界决定，海面通常为Dirichlet边界，满足的边界条件为：

Z_n(z)|_z＝0＝0 (2)

在各种海底模型交互面上，各模型满足声学边界方程为：

其中，

m₁为海水ρ₁与海底密度ρ₂比，即

k₂为海底波速定义为

c₂、α₂为海底声速和吸收系数；V为海底反射系数，对于复杂海底界面反射系数可采用“三参数”界面反射模型获取；由确定的海底边界便可确定特征函数(格林函数)Z_n(z)，根据方程(1)便可得到与信道格林函数相关的传输函数；

步骤三 信道下点源声场传输函数空间维度变换，把二维浅海信道模型下与格林函数相关的传输函数变换为三维条件下声场函数，可采取方位耦合解进行空间维度变换；

步骤四 通过波叠加法在结构内部设置虚拟源面S、，并在源面上合理布放若干个单极源，采用浅海多边界声场耦合理论计算获取信道下结构表面受界面影响的振动速度U，通过信道传输函数G建立起结构内部虚拟源强Q与结构表面振速U的传递矩阵D，满足以下关系：

{U}_M×1＝{D}_M×N·{Q}_N×1 (5)

D＝▽_nG(|r_s-r₀|) (6)

其中，G为轴对称条件下基于简正波理论推导与格林函数相关的信道传输函数表达式；r_s为结构表面任意一点的方向矢量，r₀为结构内部任意虚拟源的方向矢量；▽_n为法向求导算子；

步骤五 根据虚拟源源强Q、表面振速U和传递矩阵D的关系，则虚拟源源强Q则可表示为：

[Q]＝[D]^-1[U] (7)

步骤六 对病态传递函数[D]进行奇异值分解：

其中U，V分布为M阶和N阶方阵，Σ＝diag(σ₁,σ₂,σ₃...σ_n)为包含n个非奇异值σ_i的对角阵，且σ₁≥σ₂≥σ₃≥...≥σ_n≥0，T为矩阵转置运算；

步骤七 采用Tikhonov正则化方法对病态传递矩阵[D]进行处理，利用奇异值分解有D^TD＝VΣ²V^T，则源强表达式为：

其中，U^δ为结构表面速度矩阵，α为正则化参数；

步骤八 利用L曲线进行最优参数α的选取；

步骤九 获取虚拟源源强后，便可根据以下公式求解任意场点的声压：

P＝TQ＝TV(Σ²+αI)^-1ΣU^TU^δ (10)

其中，T为虚拟源Q与信道任意场点的单极矩阵，与浅海信道环境下对应的信道声场传输函数关系为T＝jρωG(|r-r₀|)，j为虚数，ω＝2πf为角频率，f为频率；ρ为海水密度；

所述步骤三的方位耦合解为考虑水平环境变化的影响，采取绝热模式二维图案进行变换，三维亥姆霍兹方程为：

其中，▽为求导算子，ρ(x,y,z)和c(x,y,z)为浅海信道水层密度和声速，p(x,y,z)空间声场声压，x_s为声源在x轴上的固定位置，ω为角频率，δ(y)为特殊函数；

设解的形式即三维信道传输函数G(x,y,z)为：

其中，Γ_m(x,y,z)为局部模式，m为模数，Ε_m(y,z)为二维声场函数；

其中，Γ_m(x,y,z)为局部模式，m为模数，Ε_m(y,z)为二维声场函数；

Γ_m(x,y,z)为局部模式，通过绝热近似，可得到水平折射方程：

利用局部简正波把三维问题消去了x维变为二维，得到一个关于(y，z)新的亥姆霍兹方程，水平方向的有效折射率可由

控制，根据简正波理论求解式(13)所表示的二维声场，然后代入式子(12)可求解三维声场传输函数；

所述步骤四的获取结构表面速度U，其特征在于，首次采用采用浅海多边界声场耦合理论计算获取结构表面速度U，充分考虑流体、结构和边界之间的相互耦合作用，突破多边界耦合声场环境下结构声振信息获取难题，计算获取结构表面准确的速度信息U，其计算理论包含以下几个部分：

流固耦合方程，根据Helmholtz方程结合流固耦合边界条件得出结构与流体的耦合有限元方程为：

其中，F_st、F_at分别为结构、流体介质声的耦合激励载荷,Kc、Mc分别为耦合刚度矩阵、耦合质量矩阵，且

PML(Perfectly Methoded Layer，简称PML)技术，进行信道四周无限大边界模拟，在PML吸收层域控制方程：

海面Dirichlet边界，满足的边界条件为：

p_an(x,y,z)|_z＝0＝0 (16)

海底边界，其声与海底耦合建模方法主要有Neumann边界模型，Rayleigh定律，Sommerfeld模型，Cauchy边界和Geoacoustics模型，大致可分为硬海底、液态海底、弹性海底和多孔海底等，为了区别海底与海水域，下标a和b分别标示流体层和海底层；在各类海底交互面(z＝H)上建立流体介质与海底物理模型的耦合方程：对于硬质海底，满足：

对于液态海底，满足声压和法向振速连续：

p_a＝p_b (18)

对于各向同性弹性海底，建立声-弹性波耦合方程，即法向上满足位移连续和应力连续，切向的应力为零，即：

其中w_b、u_b分别表示海底水平方向(r)的位移和垂直方向(z)的位移，λ、μ为拉米常数，Δ的定义为

对于多孔海底，所涉及声学物理参数较多，建立声-多孔弹性波耦合方程，满足边界条件为牵引力连续、流体压力连续和流体流量连续。

Images