CN112254798A - 一种预报海洋矢量声场的方法、系统及介质 - Google Patents

Abstract

本发明公开一种预报海洋矢量声场的方法、系统及介质，该方法步骤包括：S1.获取待预报海洋区域的现场测量数据以及声源参数信息，建立柱坐标系水声Helmholtz方程，经过变换后得到深度方程；S2.求解深度方程，获得各介质层分界面的声压核函数与垂直振速核函数；S3.根据声压核函数计算声压，根据垂直振速核函数使用基于Hankel反变换积分式垂直方向导数计算垂直振速，以及基于Hankel反变换积分式水平方向导数计算水平振速；S4.根据求解出的声压、振速矢量，计算出待预报海洋区域的水声传播损失与声强矢量。本发明能够基于海洋环境测量数据实现振速矢量的预报，同时能够提高振速矢量预报的精度。

Description

一种预报海洋矢量声场的方法、系统及介质

技术领域

本发明涉及水声场探测技术领域，尤其涉及一种预报海洋矢量声场的方法、系统及介质。

背景技术

声波是目前水下通信、海洋环境与目标探测的主要手段，在军事与经济领域均有重要应用价值。水下声波的接收传感器一般称为水听器，传统水听器为标量水听器，只能测量声场中的标量参数(如声压)。目前较为先进的水听器是矢量水听器，既能够测量标量参数，也可测量声场中的矢量参数(如质点振动速度，简称振速)，对提高水下探测系统的性能具有重要意义。

由于海底地形、海洋环境(主要为声速、密度，具有时变特性)、声源频率与位置等因素都对水声传播具有重要影响，在水听器布放选址、阵列形状优化、接收信号处理与分析等过程中，都需要对海洋声场进行预报，即结合动态变化的海洋环境进行声场数值模拟，并将预报结果融入到水听器信号处理过程中，因此海洋声场预报对提高水听器应用水平、增强水下探测系统性能具有重要作用。

但是传统的水声模型(如波数积分法、简正波法等)仅能预报标量声场，而若要获得矢量声场，需要在声场内布置较密的网格点，在求解出各点声压后，再采用有限差分法计算各点的声压导数，最后将声压导数转换成振速。上述传统的采用有限差分法计算声场矢量(振速)的方法，存在以下不足：

(1)声场网格间距受到有限差分法的限制。为了保证导数差分计算的正确性需要使网格间距足够小(声压值在相邻网格点区间内变化不大)，在每个方向上的网格间距一般都要取参考波长的若干分之一(如取十分之一)，在仅需要计算水听器周围少数空间位置振速的情况下，该限制条件将对声场矢量的计算量与计算方法灵活性产生不利影响。

(2)有限差分法计算声压导数将引入额外的数值误差。实际应用中一般采用有限阶数的有限差分法(如二阶精度)，在网格间距为有限值时(不是趋于零)，差分计算导数的过程将引入有限差分法的耗散与色散误差。

发明内容

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种预报海洋矢量声场的方法、系统及介质，能够基于海洋环境测量数据实现振速矢量的预报，并基于直接进行水平波数积分的振速计算方法可提高振速矢量预报的精度。

为解决上述技术问题，本发明提出的技术方案为：

一种预报海洋矢量声场的方法，步骤包括：

S1.获取待预报海洋区域的现场测量数据以及声源参数信息，建立水平分层海洋环境下的柱坐标系水声Helmholtz方程，并经过积分变换后得到以声压核函数为变量的深度方程；

S2.采用传递函数矩阵法求解所述深度方程，获得各介质层分界面的声压核函数与垂直振速核函数；

S3.根据所述声压核函数使用Hankel反变换积分式计算声压，根据所述垂直振速核函数使用基于Hankel反变换积分式垂直方向导数计算垂直振速，以及基于Hankel反变换积分式水平方向导数并将零阶Bessel函数的导数变换成一阶Bessel函数得到水平振速积分式以计算水平振速，得到振速矢量；

S4.根据步骤S3求解出的声压、振速矢量，计算出待预报海洋的水声传播损失与声强矢量。

进一步的，所述步骤S1的步骤包括：

S11.按照下式建立所述柱坐标系水声Helmholtz方程：

其中，P(r,z)为频率域相对声压，ρ为声传播介质密度，k为波数且k＝2πf/c，f为声源频率，c为介质声速，r为水平方向的坐标，z为垂直或深度方向的坐标，z_s为声源深度、δ为狄拉克函数；

S12.对水声传播介质在深度方向上划分为N层，并将每层内近似为均匀介质；在划分的每一层内，对所述柱坐标系水声Helmholtz方程进行Hankel变换，以将(r,z)空间的声压P(r,z)转换到(k_r,z)空间，即为：

其中，φ(k_r,z)为声压核函数，k_r为水平波数，J₀为Bessel函数；

对所述柱坐标水声Helmholtz方程两边同时作积分

可得到深度方程为：

进一步的，所述步骤S2中求解深度方程时，包括将声压核函数与垂直振速核函数形成联合矢量的步骤，具体步骤包括：

S211.在任意不含声源的介质层内，声压核函数的通解形式为：

其中，φ(k_r,z)为声压核函数，k_z为垂直波数且

A⁺(k_r)表示向下传播的项，A^-(k_r)表示向上传播的项，w(k_r,z)为z方向的垂直振速核函数，且满足：

其中，Γ＝k_z/(ρω)，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

则将所述声压核函数与所述垂直振速核函数形成联合矢量为：

S212.根据步骤S211形成的所述联合矢量得到第m层的上界面z_m-1处的联合矢量为：

以及第m层的下界面z_m处的联合矢量为：

S213.将步骤S212得到的第m层的上、下界面联合矢量表达式联立、消去

后，得到所述联合矢量由下至上传递公式为：

v_m(k_r,z_m-1)＝M_m(k_r)v_m(k_r,z_m)

其中，M_m(k_r)为第m层介质由下至上的传递矩阵，若令第m层厚度为h_m＝z_m-z_m-1，则M_m(k_r)的表达式为：

以及由上至下的传递公式及传递矩阵为：

进一步的，所述步骤S2中获取各介质层分界面的声压核函数与垂直振速核函数具体步骤包括：

S221.根据在计算域上边界面声能量只可能向上传播、向下传播的项为零，即

使用下标“0”表示上边界，得到上边界处的声压核函数φ₀(k_r,z₀)、垂直振速核函数w₀(k_r,z₀)分别为：

以及上边界处的联合矢量为：

其中，w₀(k_r,z₀)为上边界的垂直振速核函数，矢量(1，B₀)^T为上边界声矢量，其中

根据在计算域下边界面声能量只可能向下传播、向上传播的项为零，即

使用下标“N”表示下边界，得到下边界处的联合矢量为：

其中，w_N(k_r,z_N)为下边界的垂直振速核函数，φ_N(k_r,z_N)为下边界的声压核函数，(1，B_N)^T为下边界声矢量，其中

S222.将声矢量从下边界z_N向声源深度z_s传递，即从下边界z_N开始，逐层传递计算声矢量，直至声源深度z_s，其中声源深度紧下方的联合矢量计算式为：

式中，w_N(k_r,z_N)为待定的下边界垂直振速核函数，声源深度紧下方声矢量

S223.将声矢量从上边界z₀向声源深度z_s传递，其中声源深度紧上方传递得到的联合矢量计算式为：

式中，w₀(k_r,z₀)为待定的上边界垂直振速核函数，声源深度紧上方声矢量

S224.根据声源界面条件计算出下、上边界垂直振速核函数，其中在声源界面满足：

即v_s+1(k_r,z_s)-v_s(k_r,z_s)＝Δv(k_r,z_s)，展开得到：

则求解得到下、上边界垂直振速核函数分别为：

S225.根据计算出的下、上边界垂直振速核函数w_N+1(k_r,z_N)与w₀(k_r,z₀)计算各介质层界面位置的声压核函数φ(k_r,z)、垂直振速核函数w(k_r,z)的值。

进一步的，所述步骤S3中求解声压时，具体Hankel反变换积分式为：

将所述Hankel反变换积分式中的水平波数k_r进行离散后得到声压离散计算式为：

其中，Δk_r为水平波数步长且Δk_r＝2π/(r_maxn_w)，r_max为声场最大水平距离，n_w为Bessel函数在一个2π振荡周期内最小采样点数，k_r,n为离散的水平波数且k_r,n＝nΔk_r-iε_k，i为虚数单位，ε_k为复偏移量且ε_k＝3Δk_r/(2πlog₁₀e)，M为离散的水平波数的最大索引号且M＝k_max/Δk_r，k_max为截止波数。对声场各位置(r,z)处采用所述声压的离散计算式进行水平波数积分计算，得到对应各位置(r,z)处的声压。

进一步的，所述步骤S3中计算振速矢量时，求解垂直振速的步骤包括：

S311.对Hankel反变换积分式在深度z方向求导数，得到第一变换式为：

其中，P(r,z)为频率域相对声压，φ(k_r,z)为声压核函数，k_r为水平波数，r为水平方向的坐标，z为竖直或深度方向的坐标，w(k_r,z)为垂直振速核函数，J₀为零阶Bessel函数，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

S312.质点的振速矢量V(r,z)＝(V_r(r,z),V_z(r,z))中的垂直振速分量V_z(r,z)根据所述第一变换式进行积分计算得到：

并对水平波数k_r进行离散，得到垂直振速的离散积分计算式为：

S313.对声场各位置(r,z)处采用所述垂直振速的离散积分计算式进行水平波数积分计算，得到对应各位置(r,z)处的垂直振速。

进一步的，所述步骤S3中计算振速矢量时，求解水平振速的具体步骤包括：

S321.对Hankel反变换积分式在深度r方向求导数，并将零阶Bessel函数的导数转换成一阶Bessel函数，即

得到第二变换式为：

其中，P(r,z)为频率域相对声压，φ(k_r,z)为声压核函数，k_r为水平波数，r为水平方向的坐标，z为竖直或深度方向的坐标，J₀为零阶Bessel函数，

为一阶Bessel函数；

S322.根据所述第二变换式进行积分计算得到振速矢量中的水平振速分量V_r(r,z)为：

其中，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

并对水平波数k_r进行离散，得到水平振速的离散积分计算式为：

S323.对声场各位置(r,z)处采用所述水平振速的离散积分计算式进行水平波数积分计算，得到对应各位置(r,z)处的水平振速。

进一步的，所述步骤S4中，还包括根据步骤S3求解出的声压计算传播损失，形成声场传播损失标量云图，以及根据求解出的声压P(r,z)、振速矢量V(r,z)＝(V_r(r,z),V_z(r,z))按照下式计算声场各位置(r,z)的时均声强矢量；

其中，I(r,z)为时均声强矢量，即在垂直于振速方向的平面上，单位面积上通过的声能量的时间平均值，

分别表示复数声压P(r,z)、振速V(r,z)的共轭复数，r为水平方向的坐标，z为竖直或深度方向的坐标。

一种预报海洋矢量声场的系统，包括处理器以及存储器，所述存储器用于存储海洋环境数据、声源参数与计算机程序，所述处理器用于执行所述计算机程序，所述处理器用于执行所述计算机程序，以执行上述方法。

一种计算机可读存储介质，该计算机可读存储介质上存储有海洋环境数据、声源参数与被编程或配置以执行上述预报海洋矢量声场的方法的计算机程序。

与现有技术相比，本发明的优点在于：本发明预报海洋矢量声场的方法、系统及介质，获取到待预报海洋区域的现场测量数据以及声源参数信息后，在使用波数积分法水声模型计算声压的基础上，通过对Hankel反变换积分式求导，获得垂直振速与水平振速的波数积分式，可以直接对水平波数进行离散并积分计算声场任意点的振速矢量，突破了传统方法因差分计算声压导数求振速引起的对声场网格点间距的限制，同时可避免有限差分法的数值误差，从而提高振速矢量预报精度。

附图说明

图1是本实施例预报海洋矢量声场的方法的详细实现流程示意图。

图2是本发明在具体应用实施例中实现海洋矢量声场预报方法的系统结构示意图。

图3是本发明在具体应用实施例中得到的声场传播损失标量云图。

图4是本发明在具体应用实施例中得到的声强矢量流线图。

具体实施方式

以下结合说明书附图和具体优选的实施例对本发明作进一步描述，但并不因此而限制本发明的保护范围。

本实施例具体以绝对硬海底等速波导海洋环境下的声场预报为例，各项参数具体为：海水密度均匀ρ_w＝1.0g/cm³，水体声速均匀c_w＝1500m/s，海底水平且海深z_N＝100m、深度z方向步长dz＝1m，r方向最大求解距离为r_max＝1000m、步长Δr＝1m。声源频率f＝100Hz，声源深度z_s＝25m，上边界(海面z₀＝0m)取压力释放边界条件(声压为零)、下边界(海底z_N＝100m)取绝对硬边界条件(声压z向导数为零)。

如图1所示，本实施例预报海洋矢量声场的方法的详细步骤包括：S1.获取待预报海洋区域的现场测量数据以及声源参数信息，建立水平分层海洋环境下的柱坐标系水声Helmholtz方程，并经过积分变换后得到以声压核函数为变量的深度方程；

S2.采用传递函数矩阵法求解深度方程，获得各介质层分界面的声压核函数与垂直振速核函数；

S3.根据所述声压核函数使用Hankel反变换积分式计算声压，根据所述垂直振速核函数使用基于Hankel反变换积分式垂直方向导数计算垂直振速，以及基于Hankel反变换积分式水平方向导数并将零阶Bessel函数的导数变换成一阶Bessel函数得到水平振速积分式以计算水平振速，得到振速矢量；

S4.根据步骤S3求解出的声压、振速矢量，计算出待预报海洋的水声传播损失与声强矢量。

使用波数积分法水声模型可以计算出声压，同时在使用波数积分法的基础上，利用Hankel反变换积分式在垂直方向求导，可得到垂直振速积分式，进而可求解出垂直振速，同时考虑Bessel函数具有特性：

其中J₀为零阶Bessel函数，J′₀为J₀的导数，J₁为一阶Bessel函数，即可以将零阶Bessel函数转换为一阶Bessel函数，则在Hankel反变换积分式的基础上，在水平方向求导后利用Bessel函数的上述特性，又可以得到水平振速积分式，进而可以求解出水平振速。本实施例利用上述特性，在获取到待预报海洋区域的现场测量数据以及声源参数信息后，使用波数积分法水声模型计算声压，同时通过对Hankel反变换积分式求导，分别获得垂直振速与水平振速的波数积分式，其中利用垂直振速核函数积分计算出垂直振速，利用Bessel函数的特性以及Hankel反变换积分式水平方向导数求解出水平振速，得到振速矢量，从而可实现水声传播损失与声强矢量的预报。本实施例上述方法，基于直接进行水平波数积分的振速计算方式，可以直接计算获得声场任意点的振速矢量，突破了传统方法因差分计算声压导数求振速引起的对声场网格点间距的限制，同时可避免有限差分法的数值误差从而提高振速矢量预报精度，进而提升矢量水听器的应用技术水平。

本实施例步骤S1中首先获取待预报海洋区域的现场测量数据以及声源参数信息，其中海洋区域的现场测量数据包括海洋深度、声速、密度等数据(本实施例中分别为海深1000米、水中声速均匀为1500m/s，全场密度均为1.0g/cm³)，声源参数信息包括声源频率、位置等数据(本实施例中声源频率100Hz、深度25米)，根据获取的数据建立水平分层海洋环境下的柱坐标系水声Helmholtz方程，再对该柱坐标系水声Helmholtz方程进行Hankel变换，得到以声压波数核函数为变量的深度方程。

本实施例中，步骤S1具体的步骤包括：

S11.按照下式建立柱坐标系水声Helmholtz方程：

其中，P(r,z)为频率域相对声压，ρ为声传播介质密度，k为波数且k＝2πf/c，f为声源频率，c为介质声速，r为水平方向的坐标，z为垂直或深度方向的坐标，z_s为声源深度、δ为狄拉克函数；

S12.对水声传播介质在深度方向上划分为N层，并将每层内近似为均匀介质；在划分的每一层内，对柱坐标系水声Helmholtz方程进行Hankel变换，以将(r,z)空间的声压P(r,z)转换到(k_r,z)空间，即为：

其中，φ(k_r,z)为声压核函数，k_r为水平波数，J₀为Bessel函数；

对柱坐标水声Helmholtz方程两边同时作积分

可得到深度方程为：

本实施例具体对水声传播介质在竖直方向上划分为N＝100层，每层厚度为dz＝1m，且每层内为均匀介质，在每一层内，对柱坐标系水声Helmholtz方程按照上式(2)进行Hankel变换，并对柱坐标水声Helmholtz方程两边同时作积分

得到以声压波数核函数为变量的深度方程(3)。

本实施例步骤S2中求解深度方程时，包括将声压核函数与垂直振速核函数形成联合矢量的步骤，具体步骤包括：

S211.在任意不含声源的介质层内，声压核函数的通解形式为：

其中，φ(k_r,z)为声压核函数，k_z为垂直波数且

指数项

分别具有向上、下传播的物理意义，A⁺(k_r)表示向下传播的项，A^-(k_r)表示向上传播的项，w(k_r,z)为z方向(垂直方向)的垂直振速核函数，且满足：

其中，Γ＝k_z/(ρω)，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

将声压核函数与垂直振速核函数形成联合矢量为：

S212.根据步骤S211形成的联合矢量得到第m层的上界面z_m-1处的联合矢量为：

以及第m层的下界面z_m处的联合矢量为：

S213.将步骤S212得到的第m层的上、下界面联合矢量表达式联立、消去

后，得到联合矢量由下至上传递公式为：

v_m(k_r,z_m-1)＝M_m(k_r)v_m(k_r,z_m) (9)

其中，M_m(k_r)为第m层介质由下至上的传递矩阵，若令第m层厚度为h_m＝z_m-z_m-1，则M_m(k_r)的表达式为：

以及由上至下的传递公式及传递矩阵分别为：

按照上述步骤，可以建立由声压核函数与垂直振速核函数组成的联合矢量，同时得到联合矢量由下至上以及由上至下的传递公式，便于后续对深度方程进行求解。

本实施例步骤S2中求解深度方程，获取各介质层分界面的声压核函数与垂直振速核函数的具体步骤包括：

S221.根据在计算域上边界面声能量只可能向上传播、向下传播的项为零，即

使用下标“0”表示上边界，得到上边界处的声压核函数φ₀(k_r,z₀)、垂直振速核函数w₀(k_r,z₀)分别为：

以及上边界处的联合矢量为：

其中，w₀(k_r,z₀)为上边界的垂直振速核函数，矢量(1，B₀)^T为上边界声矢量，其中

本实施例具体

ρ₀表示海面上方介质密度，在绝对软(压力释放)边界条件下，海面上方可视为真空，即ρ₀＝0g/cm³；

同理，根据在计算域下边界面声能量只可能向下传播、向上传播的项为零，即

使用下标“N”表示下边界，得到下边界处的联合矢量为：

其中，w_N(k_r,z_N)为下边界的垂直振速核函数，φ_N(k_r,z_N)为下边界的声压核函数，(1，B_N)^T为下边界声矢量，

本实施例具体

ρ_N+1分别表示海底下方介质密度，在绝对硬边界条件下，海底下方可视为密度非常大的岩石，即取ρ_N+1＝10⁹⁹g/cm³；

S222.将声矢量从下边界z_N向声源深度z_s传递，即利用矢量由下至上传递公式，从下边界z_N开始，逐层传递计算声矢量，直至声源深度z_s，具体z_N＝100m、z_s＝25m，其中声源深度紧下方的联合矢量计算式为：

式中，w_N(k_r,z_N)为待定的下边界垂直振速核函数，声源深度紧下方声矢量

S223.将声矢量从上边界z₀向声源深度z_s传递，具体z₀＝0m，其中声源深度紧上方传递得到的联合矢量计算式为：

式中，w₀(k_r,z₀)为待定的上边界垂直振速核函数，声源深度紧上方声矢量

S224.根据声源界面条件计算出下、上边界垂直振速核函数，其中在声源界面满足：

即v_s+1(k_r,z_s)-v_s(k_r,z_s)＝Δv(k_r,z_s)，展开得到：

则求解得到下、上边界垂直振速核函数分别为：

S225.根据计算出的下、上边界垂直振速核函数w_N+1(k_r,z_N)与w₀(k_r,z₀)计算各介质层界面位置的声压核函数φ(k_r,z)、垂直振速核函数w(k_r,z)的值。

通过上述步骤，使用传递函数矩阵法先从下边界声矢量逐层向上传递求解各层声矢量，直至求解出声源界面紧下方声矢量；再从上边界声矢量逐层向下传递求解各层声矢量，直至求出声源界面紧上方声矢量；在声源界面处，利用声源条件建立紧下方与紧上方联合矢量方程并求解出下、上边界的垂直振速；最后将下、上边界的垂直振速乘以各层界面保存的声矢量，即可得到声压核函数与垂直振速核函数。

本实施例具体先利用各层界面求解出的声压核函数，根据Hankel反变换积分式求解出声压；然后利用各层界面求解出的垂直振速核函数，根据Hankel反变换积分式的垂直方向导数求解出垂直振速；再通过将零阶Bessel函数的导数转换成一阶Bessel函数，根据Hankel反变换积分式的水平方向导数求解出水平振速，详细如下所示。

本实施例步骤S3中求解声压时，具体Hankel反变换积分式为：

将Hankel反变换积分式中的水平波数k_r进行离散后得到的声压离散计算式为：

其中，Δk_r为水平波数步长且Δk_r＝2π/(r_maxn_w)，r_max为声场最大水平距离(本实施例取r_max＝3000m)，n_w为Bessel函数在一个2π振荡周期内最小采样点数(本实施例取n_w＝10)，k_r,n为离散的水平波数且k_r,n＝nΔk_r-iε_k，i为虚数单位，ε_k为复偏移量以用于防止深度方程求解过程出现奇异，ε_k＝3Δk_r/(2πlog₁₀e)，M为离散的水平波数的最大索引号且M＝k_max/Δk_r，k_max为截止波数，本实施例具体取k_max＝20k₀，其中海水波数k₀＝2πf/1500。对声场各位置(r,z)处采用所述声压的离散计算式进行水平波数积分计算，得到对应各位置(r,z)处的声压。

本实施例步骤S3中计算振速矢量时，求解垂直振速的具体步骤包括：

S311.对Hankel反变换积分式在垂直z方向求导数，得到第一变换式为：

其中，P(r,z)为频率域相对声压，φ(k_r,z)为声压核函数，k_r为水平波数，r为水平方向的坐标，z为竖直或深度方向的坐标，w(k_r,z)为垂直振速核函数，J₀为零阶Bessel函数，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

S312.质点的振速矢量V(r,z)＝(V_r(r,z),V_z(r,z))中的垂直振速分量V_z(r,z)根据第一变换式进行积分计算得到：

并采用与上述声压计算相同的方法对水平波数k_r进行离散，得到垂直振速的离散积分计算式为：

S313.对声场各位置(r,z)处采用垂直振速的离散积分计算式进行水平波数积分计算，得到对应各位置(r,z)处的垂直振速。

本实施例步骤S3中计算振速矢量时，求解水平振速的步骤包括：

S321.对Hankel反变换积分式在水平r方向求导数，并将零阶Bessel函数的导数转换成一阶Bessel函数，即

得到第二变换式为：

其中，P(r,z)为频率域相对声压，φ(k_r,z)为声压核函数，k_r为水平波数，r为水平方向的坐标，z为竖直或深度方向的坐标，J₀为零阶Bessel函数，J₁为一阶Bessel函数；

S322.根据第二变换式进行积分计算得到振速矢量中的水平振速分量V_r(r,z)为：

其中，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

并采用与上述声压计算相同的方法对水平波数k_r进行离散，得到水平振速的离散积分计算式为：

S323.对声场各位置(r,z)处采用水平振速的离散积分计算式进行水平波数积分计算，得到对应各位置(r,z)处的水平振速。

通过上述步骤，通过对Hankel反变换积分式进行垂直、水平两个方向求导获得垂直振速与水平振速的积分式，采用与声压积分式相同的方法对水平波数进行离散，可以求解出垂直振速与水平振速，从而实现振速矢量V(r,z)＝(V_r(r,z),V_z(r,z))的预报。

本实施例步骤S4中，还包括步骤S3求解出的声压计算传播损失，形成声场传播损失标量云图，根据声压计算传播损失(标量)的计算式具体为：

TL(r,z)＝-20log₁₀|P(r,z)| (30)

以及根据求解出的声压P(r,z)、振速矢量V(r,z)＝(V_r(r,z),V_z(r,z))按照式(31)计算声场各位置(r,z)的时均声强矢量；

其中，I(r,z)为时均声强矢量，即在垂直于振速方向的平面上，单位面积上通过的声能量的时间平均值，

分别表示复数声压P(r,z)、振速V(r,z)的共轭复数。

在具体应用实施例中得到的声场传播损失标量云图如图3所示，以及由计算得到的声场声强矢量流线(背景色为传播损失)如图4所示，由图4可见，传播损失小(能量高)的区域声强矢量线相对密集、传播损失大(能量低)的区域声强矢量线相对稀疏。

如图2所示，在具体应用实施例中应用本实施例上述方法时，结合高性能计算机工作站，在高性能计算机工作站的数据存储介质中加载能够实现本实施例上述预报海洋矢量声场的方法功能的程序模块，由高性能计算机工作站接收包括海洋深度、声速、密度等现场测量数据，以及包括声源频率、位置等声源参数信息，经过上述预报海洋矢量声场的方法步骤后，生成海洋矢量声场预报图形图像，进一步还可提供给后续进行矢量水听器声信号处理与分析。

本实施例还提供一种预报海洋矢量声场的系统，包括处理器以及存储器，存储器用于存储海洋环境数据、声源参数与计算机程序，处理器用于执行计算机程序，处理器用于执行计算机程序，以执行上述预报海洋矢量声场的方法。本实施例系统具体可采用如图2所示结构，高性能计算机工作站配置有上述处理器及存储器。

此外，本实施例还提供一种计算机可读存储介质，该计算机可读存储介质上存储有海洋环境数据、声源参数与被编程或配置以执行上述提高预报海洋矢量声场的方法的计算机程序。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

上述只是本发明的较佳实施例，并非对本发明作任何形式上的限制。虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明。因此，凡是未脱离本发明技术方案的内容，依据本发明技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均应落在本发明技术方案保护的范围内。

Claims (10)

1.一种预报海洋矢量声场的方法，其特征在于，步骤包括：

S1.获取待预报海洋区域的现场测量数据以及声源参数信息，建立水平分层海洋环境下的柱坐标系水声Helmholtz方程，并经过积分变换后得到以声压核函数为变量的深度方程；

S2.采用传递函数矩阵法求解所述深度方程，获得各介质层分界面的声压核函数与垂直振速核函数；

S3.根据所述声压核函数使用Hankel反变换积分式计算声压，根据所述垂直振速核函数使用基于Hankel反变换积分式垂直方向导数计算垂直振速，以及基于Hankel反变换积分式水平方向导数并将零阶Bessel函数的导数变换成一阶Bessel函数得到水平振速积分式以计算水平振速，得到振速矢量；

S4.根据步骤S3求解出的声压、振速矢量，计算出待预报海洋的水声传播损失与声强矢量。

2.根据权利要求1所述的预报海洋矢量声场的方法，其特征在于，所述步骤S1的步骤包括：

S11.按照下式建立所述柱坐标系水声Helmholtz方程：

其中，P(r,z)为频率域相对声压，ρ为声传播介质密度，k为波数且k＝2πf/c，f为声源频率，c为介质声速，r为水平方向的坐标，z为垂直或深度方向的坐标，z_s为声源深度、δ为狄拉克函数；

S12.对水声传播介质在深度方向上划分为N层，并将每层内近似为均匀介质；在划分的每一层内，对所述柱坐标系水声Helmholtz方程进行Hankel变换，以将(r,z)空间的声压P(r,z)转换到(k_r,z)空间，即为：

其中，φ(k_r,z)为声压核函数，k_r为水平波数，J₀为Bessel函数；

对所述柱坐标水声Helmholtz方程两边同时作积分

可得到深度方程为：

3.根据权利要求1所述的预报海洋矢量声场的方法，其特征在于，所述步骤S2中求解深度方程时，包括将声压核函数与垂直振速核函数形成联合矢量的步骤，具体步骤包括：

S211.在任意不含声源的介质层内，声压核函数的通解形式为：

其中，φ(k_r,z)为声压核函数，k_z为垂直波数且

A⁺(k_r)表示向下传播的项，A^-(k_r)表示向上传播的项，w(k_r,z)为z方向的垂直振速核函数，且满足：

其中，Γ＝k_z/(ρω)，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

则将所述声压核函数与所述垂直振速核函数形成联合矢量为：

S212.根据步骤S211形成的所述联合矢量得到第m层的上界面z_m-1处的联合矢量为：

以及第m层的下界面z_m处的联合矢量为：

S213.将步骤S212得到的第m层的上、下界面联合矢量表达式联立、消去

后，得到所述联合矢量由下至上传递公式为：

v_m(k_r,z_m-1)＝M_m(k_r)v_m(k_r,z_m)

其中，M_m(k_r)为第m层介质由下至上的传递矩阵，若令第m层厚度为h_m＝z_m-z_m-1，则M_m(k_r)的表达式为：

以及由上至下的传递公式及传递矩阵为：

4.根据权利要求3所述的预报海洋矢量声场的方法，其特征在于，所述步骤S2中求解所述深度方程，获取各介质层分界面的声压核函数与垂直振速核函数的具体步骤包括：

S221.根据在计算域上边界面声能量只可能向上传播、向下传播的项为零，即

使用下标“0”表示上边界，得到上边界处的声压核函数φ₀(k_r,z₀)、垂直振速核函数w₀(k_r,z₀)分别为：

以及上边界处的联合矢量为：

其中，w₀(k_r,z₀)为上边界的垂直振速核函数，矢量(1，B₀)^T为上边界声矢量，其中

根据在计算域下边界面声能量只可能向下传播、向上传播的项为零，即

使用下标“N”表示下边界，得到下边界处的联合矢量为：

其中，w_N(k_r,z_N)为下边界的垂直振速核函数，φ_N(k_r,z_N)为下边界的声压核函数，(1，B_N)^T为下边界声矢量，其中

S222.将声矢量从下边界z_N向声源深度z_s传递，即从下边界z_N开始，逐层传递计算声矢量，直至声源深度z_s，其中声源深度紧下方的联合矢量计算式为：

式中，w_N(k_r,z_N)为待定的下边界垂直振速核函数，声源深度紧下方声矢量

S223.将声矢量从上边界z₀向声源深度z_s传递，其中声源深度紧上方传递得到的联合矢量计算式为：

式中，w₀(k_r,z₀)为待定的上边界垂直振速核函数，声源深度紧上方声矢量

S224.根据声源界面条件计算出下、上边界垂直振速核函数，其中在声源界面满足：

即v_s+1(k_r,z_s)-v_s(k_r,z_s)＝Δv(k_r,z_s)，展开得到：

则求解得到下、上边界垂直振速核函数分别为：

S225.根据计算出的下、上边界垂直振速核函数w_N+1(k_r,z_N)、w₀(k_r,z₀)计算各介质层界面位置的声压核函数φ(k_r,z)、垂直振速核函数w(k_r,z)的值。

5.根据权利要求1所述的预报海洋矢量声场的方法，其特征在于，所述步骤S3中求解声压时，具体Hankel反变换积分式为：

将所述Hankel反变换积分式中的水平波数k_r进行离散后得到声压离散计算式为：

其中，Δk_r为水平波数步长且Δk_r＝2π/(r_maxn_w)，r_max为声场最大水平距离，n_w为Bessel函数在一个2π振荡周期内最小采样点数，k_r,n为离散的水平波数且k_r,n＝nΔk_r-iε_k，i为虚数单位，ε_k为复偏移量且ε_k＝3Δk_r/(2πlog₁₀e)，M为离散的水平波数的最大索引号且M＝k_max/Δk_r，k_max为截止波数；

对声场各位置(r,z)处采用所述声压的离散计算式进行水平波数积分计算，得到对应各位置(r,z)处的声压。

6.根据权利要求1～5中任意一项所述的预报海洋矢量声场的方法，其特征在于，所述步骤S3中计算振速矢量时，求解垂直振速的步骤包括：

S311.对Hankel反变换积分式在垂直z方向求导数，得到第一变换式为：

其中，P(r,z)为频率域相对声压，φ(k_r,z)为声压核函数，k_r为水平波数，r为水平方向的坐标，z为竖直或深度方向的坐标，w(k_r,z)为垂直振速核函数，J₀为零阶Bessel函数，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

S312.质点的振速矢量V(r,z)＝(V_r(r,z)，V_z(r,z))中的垂直振速分量V_z(r,z)根据所述第一变换式进行积分计算得到：

并对水平波数k_r进行离散，得到垂直振速的离散积分计算式为：

S313.对声场各位置(r,z)处采用所述垂直振速的离散积分计算式进行水平波数积分计算，得到对应各位置(r,z)处的垂直振速。

7.根据权利要求1～5中任意一项所述的预报海洋矢量声场的方法，其特征在于，所述步骤S3中计算振速矢量时，求解水平振速的具体步骤包括：

S321.对Hankel反变换积分式在水平r方向求导数，并将零阶Bessel函数的导数转换成一阶Bessel函数，即

得到第二变换式为：

其中，P(r,z)为频率域相对声压，φ(k_r,z)为声压核函数，k_r为水平波数，r为水平方向的坐标，z为竖直或深度方向的坐标，J₀为零阶Bessel函数，

为一阶Bessel函数；

S322.根据所述第二变换式进行积分计算得到振速矢量中的水平振速分量V_r(r,z)为：

其中，ρ为声传播介质密度，ω＝2πf为声源振动角频率，f为声源振动频率；

并对水平波数k_r进行离散，得到水平振速的离散积分计算式为：

S323.对声场各位置(r,z)处采用所述水平振速的离散积分计算式进行水平波数积分计算，得到对应各位置(r,z)处的水平振速。

8.根据权利要求1～5中任意一项所述的预报海洋矢量声场的方法，其特征在于，所述步骤S4中，还包括根据步骤S3求解出的声压计算传播损失，形成声场传播损失标量云图，以及根据求解出的声压P(r,z)、振速矢量V(r,z)＝(V_r(r,z),V_z(r,z))按照下式计算声场各位置(r,z)的时均声强矢量；

其中，I(r,z)为时均声强矢量，即在垂直于振速方向的平面上，单位面积上通过的声能量的时间平均值，

分别表示复数声压P(r,z)、振速V(r,z)的共轭复数，r为水平方向的坐标，z为竖直或深度方向的坐标。

9.一种预报海洋矢量声场的系统，包括处理器以及存储器，所述存储器用于存储海洋环境数据、声源参数与计算机程序，所述处理器用于执行所述计算机程序，其特征在于，所述处理器用于执行所述计算机程序，以执行如权利要求1～8中任意一项所述方法。

10.一种计算机可读存储介质，其特征在于，该计算机可读存储介质上存储有海洋环境数据、声源参数与被编程或配置以执行权利要求1～8中任意一项所述预报海洋矢量声场的方法的计算机程序。

Images