CN112254797A - 一种提高海洋声场预报精度的方法、系统及介质 - Google Patents
一种提高海洋声场预报精度的方法、系统及介质 Download PDFInfo
- Publication number
- CN112254797A CN112254797A CN202011083585.XA CN202011083585A CN112254797A CN 112254797 A CN112254797 A CN 112254797A CN 202011083585 A CN202011083585 A CN 202011083585A CN 112254797 A CN112254797 A CN 112254797A
- Authority
- CN
- China
- Prior art keywords
- sound
- depth
- vector
- acoustic
- layer
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H3/00—Measuring characteristics of vibrations by using a detector in a fluid
- G01H3/04—Frequency
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H3/00—Measuring characteristics of vibrations by using a detector in a fluid
- G01H3/10—Amplitude; Power
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
本发明公开一种提高海洋声场预报精度的方法、系统及介质,该方法步骤包括:步骤S1.获取待测海洋水声场的现场测量数据以及声源的参数信息,建立水平分层海洋环境下的柱坐标系水声Helmholtz方程,经过变换后得到深度方程;步骤S2.建立声场上、下边界的声矢量;步骤S3.将声矢量分别从上、下边界向中间深度传递;步骤S4.在中间深度建立声矢量方程,并求解上、下边界的垂直振速,计算出各层声压波数核函数;步骤S5.对声压波数核函数进行水平波数积分,得到接收深度的声压值;步骤S6.计算接收深度的传播损失曲线,实现待测海洋声场的预报。本发明能够基于海洋环境测量数据实现声场预报,同时能够提高海洋声场预报的精度。
Description
技术领域
本发明涉及水声场探测技术领域,尤其涉及一种提高海洋声场预报精度的方法、系统及介质。
背景技术
声波能够在海水介质中远距离传播,是目前水声通信、海水与海底环境探测、水下目标定位的主要信息载体,在军事与经济领域均有重要应用价值。由于水下声传播特性与海水、沉积层的声速与密度分布密切相关,一般需要建立适当的水声模型以描述声压变量与环境介质声学参数、声源参数之间的关系,并借助计算机进行求解,最后通过分析传播损失随水平距离的变化规律获取水声传播特性。水声传播的物理过程受波动方程控制,当声场处于稳态时可将时域内的波动方程采用Fourier变换转为频域内空间三维Helmholtz方程,每次只针对特定频率的声场进行计算,使声场求解难度有所降低。由于Helmholtz方程属于椭圆型方程,采用有限差分、有限元等直接数值离散方法将形成大型线性方程组,并需要采用迭代法求解,这对于追求时效性的水声应用来讲其计算量、存储量仍然过大,因此实际中经常通过各类假设对Helmholtz方程进行简化,衍生出波数积分法、简正波法、抛物方程法、射线法等水声模型。
由于海水中温度、盐度分布主要与深度相关,一般情况下在水平方向变化较慢,在海底地形平坦的区域可良好地近似为水平分层介质(声速、密度不随水平距离变化),满足波数积分法的计算条件。波数积分法未对Helmholtz方程进行简化,通常被称为“精确解”,已广泛应用于海洋声场仿真、海洋声学参数与水下目标定位反演,特别是其提供的参考解可用于校核简正波、抛物方程、射线法等其他水声模型的计算精度,具有重要的理论与应用价值。
传统波数积分模型采用传递函数矩阵法求解深度方程时会遇到计算不稳定问题,原因是当水平波数大于介质波数以后,传递矩阵中包含的指数项使部分介质层分界面处波数核函数的幅值越来越大,直至超过程序变量(提前设定精度)的存储位数限制,从而导致计算崩溃并中断。波数积分过程中声矢量的传递策略也会对计算稳定性有重要影响,传统方法中传递策略为从上、下两个边界分别向声源深度传递,在声源深度利用声源条件建立联合矢量方程,通常实现的具体步骤即为:从下边界声矢量逐层向上传递求解各层声矢量,直至求解出声源界面紧下方声矢量;从上边界声矢量逐层向下传递求解各层声矢量,直至求出声源界面紧上方声矢量;在声源界面处,利用声源条件建立紧下方与紧上方联合矢量方程并求解。但是上述传统方法在深海条件下使用时会存在很大的不足,特别是当声源距离海面较近时,下边界声矢量向上传递至声源深度的距离很远,计算稳定性差(水平波数超过介质波数后很快就出现计算中断),导致声源附近深度的声压波数积分计算误差较大,进而降低声场的预报精度。
发明内容
本发明要解决的技术问题就在于:针对现有技术存在的技术问题,本发明提供一种提高海洋声场预报精度的方法、系统及介质,能够基于海洋环境测量数据实现声场预报,同时能够提高海洋声场预报的精度。
为解决上述技术问题,本发明提出的技术方案为:
一种提高海洋声场预报精度的方法,步骤包括:
步骤S1.获取待预报海洋水声场的现场测量数据以及声源的参数信息,建立水平分层海洋环境下的柱坐标系水声Helmholtz(亥姆霍兹)方程,经过积分变换后得到以声压核函数为变量的深度方程;
步骤S2.根据步骤S1得到的所述深度方程建立声场上、下边界的声矢量;
步骤S3.将步骤S2建立的所述声矢量分别从上、下边界向中间深度传递,以逐层传递计算所述声矢量;
步骤S4.根据步骤S3传递计算的所述声矢量在中间深度建立联合矢量方程,并求解上、下边界的垂直振速,计算出各层声压波数核函数;
步骤S5.对步骤S4计算出的声压波数核函数进行水平波数积分,得到接收深度的声压值;
步骤S6.由步骤S5计算得到的对应不同接收深度的声压值得到接收深度的传播损失曲线,根据所述接收深度的传播损失曲线实现海洋声场的预报。
进一步的,所述步骤S1的具体步骤包括:
S11.按照下式建立所述柱坐标系水声Helmholtz方程:
其中,P(r,z)为频率域相对声压,ρ为声传播介质密度,k为介质波数,r水平方向的坐标,z为竖直或深度方向的坐标,zs为声源深度、δ为狄拉克函数;
S12.对水声传播介质在深度方向上划分为N层,并将每层内近似为均匀介质;在划分的每一层内,对所述柱坐标系水声Helmholtz方程进行Hankel变换,以将(r,z)空间的声压P(r,z)转换到(kr,z)空间,即为:
其中,φ(kr,z)为声压核函数、kr为水平波数,J0为Bessel函数;
进一步的,所述步骤S2具体采用四倍精度程序变量建立所述声场上、下边界的声矢量,具体步骤包括:
S21.将振速核函数与步骤S1得到的所述声压核函数形成联合矢量:
以及上边界处的联合矢量为:
S23.将步骤S22得到的上边界的振速w0、下边界的振速wN、上边界声矢量(1,B0)T、下边界声矢量(1,BN)T均设为四倍精度变量,分别建立声场上、下边界的联合矢量。
进一步的,所述步骤S3中,包括将声矢量从下边界zN向中间深度zh传递的步骤S31,所述步骤S31的具体步骤包括:
S311.将所述声压核函数与所述振速核函数形成联合矢量并转换为矩阵形式得到:
S312.根据步骤S311转换得到的矩阵形式,得到第m层的上界面zm-1处的联合矢量为:
以及第m层的下界面zm处的联合矢量为:
vm(kr,zm-1)=Mm(kr)vm(kr,zm)
其中,Mm(kr)为第m层介质由下至上的传递矩阵,若令第m层厚度为hm=zm-zm-1,则Mm(kr)的表达式为:
S314.使用所述矢量传递公式,从下边界zN开始,逐层传递计算所述声矢量,直至中间深度zh,中间深度由下方传递得到的联合矢量为:
S315.如果声源深度zs处于计算域下边界zN与中间深度zh之间,即zh<zs<zN,则在跨过声源深度之后增加声源条件,此时中间深度由下方传递得到的联合矢量计算式为:
其中,ρs表示声源深度的介质密度,角频率ω=2πf。
进一步的,所述步骤S3中,还包括将声矢量从上边界z0向中间深度zh传递的步骤S32,所述步骤S32的具体步骤包括:
S321.获取第m层介质由上至下的联合矢量传递公式为:
以及中间深度由上方传递得到的联合矢量为:
S322.如果声源深度zs处于计算域上边界z0与中间深度zh之间,即z0<zs<zh,则在跨过声源深度之后需要增加声源条件,此时中间深度由上方传递得到的联合矢量计算式为:
进一步的,所述步骤S4的步骤包括:
S41.根据步骤S31得到的由下方传递得到的联合矢量vh+1(kr,zh)与步骤S32得到的由上方传递得到的联合矢量vh(kr,zh)在中间深度相等,构建二元一次方程组:
vh+1(kr,zh)=vh(kr,zh)
求解所述二元一次方程组,求解得到所述上边界振速w0(kr,z0)与所述下边界振速wN(kr,zN);
其中,h<m≤N;
如果声源深度zs处于计算域下边界zN与中间深度zh之间,即zh<zs<zN,则越过声源深度之后,将所述φm(kr,zm-1)计算式转换为:
其中,1≤m≤h;
如果声源深度zs处于计算域上边界z0与中间深度zh之间,即z0<zs<zh,则在跨过声源深度之后,将所述φm(kr,zm)的计算式转换为:
进一步的,所述步骤S5中,具体利用Hankel反变换的声压积分式对声压波数核函数进行离散形式的水平波数积分,得到所述接收深度的声压值,步骤包括:
S51.Hankel反变换公式为:
其中,P(r,z)为频率域相对声压,φ(kr,z)为声压核函数,kr为水平波数,r为水平方向的坐标,z为竖直或深度方向的坐标,zs为声源深度;
将所述Hankel反变换公式中的水平波数kr进行离散,得到声压离散式为:
其中,Δkr为水平波数步长且Δkr=2π/(rmaxnw),rmax为声场最大水平距离,nw为Bessel函数在一个2π振荡周期内最小采样点数,kr,n=nΔkr-iεk为离散的水平波数,i为虚数单位,εk为复偏移量且εk=3Δkr/(2πlog10e),M=kmax/Δkr为离散的水平波数的最大索引号,kmax为最大截止波数;
S52.在声场任意位置点(r,z)处采用所述声压离散式进行水平波数积分,得到各位置点的声压值。
进一步的,采用预估-试算法确定所述最大截止波数,具体步骤包括:首先设定一个初始预估值,然后按照式kr,n=nΔkr-iεk逐渐增大水平波数,其中n=1,2,3......,M,若当n=j时出现计算崩溃并中断,则可确定最大截止波数kmax=(j-1)Δkr-iεk。
一种提高海洋声场预报精度的系统,包括处理器以及存储器,所述存储器用于存储海洋环境数据、声源参数与计算机程序,所述处理器用于执行所述计算机程序,所述处理器用于执行所述计算机程序,以执行上述提高海洋声场预报精度的方法。
一种计算机可读存储介质,该计算机可读存储介质上存储有海洋环境数据、声源参数与被编程或配置以执行上述提高海洋声场预报精度的方法的计算机程序。
与现有技术相比,本发明的优点在于:
1、本发明通过建立得到声场上、下边界的声矢量后,在波数积分过程中采用分别从上、下边界向中间深度传递的方式逐层传递计算声矢量,同时在中间深度建立联合矢量方程,相比于传统波数积分方法中从上、下两个边界分别向声源深度传递声矢量并在声源深度利用声源条件建立联合矢量方程,更适用于深海环境声场预报,尤其是当声源距离海面较近时(下边界声矢量向上传递至声源深度的距离很远),能够大幅增强海洋声场预报过程中深度方程的求解稳定性,进而大幅提升水声场波数积分模型的计算精度,从而可以有效提高海洋声场预报的精度,减少声场预报的误差。
2、本发明通过采用分别从上、下边界向中间深度传递的方式逐层传递计算声矢量,同时在中间深度建立联合矢量方程,不仅能够增强深度方程的求解稳定性,还能够显著提高最大截止波数,即可使水平波数取值范围更大,从而使声压积分值更精确,提高声场预报精度。
3、本发明进一步通过采用四倍精度变量定义声矢量,可以减小计算机舍入误差,进一步提高海洋声场的预报精度。
附图说明
图1是本实施例提高海洋声场预报精度的方法的实现流程示意图。
图2是本发明在具体应用实施例中实现提高海洋声场预报精度的方法的系统结构示意图。
图3是在具体应用实施例中分别采用传统方法与本发明方法的声矢量传递计算方式对比示意图。
图4是本发明在具体应用实施例中得到的接收深度的传播损失曲线。
具体实施方式
以下结合说明书附图和具体优选的实施例对本发明作进一步描述,但并不因此而限制本发明的保护范围。
本实施例具体以Munk波导海洋环境下的声场预报为例,海水密度均匀ρw=1.0g/cm3,水体声速与深度z相关:cw(z)=1500[1+0.00737(ζ-1+exp(-ζ))],其中参数ζ=2(z-1300)/1300,水体声速单位m/s;沉积层密度均匀且与水体相同ρsed=1.0g/cm3,沉积层声速均匀csed=1600m/s;海底水平且海深zN=5000m、深度z方向步长dz=1m,r方向最大求解距离为rmax=100km、步长Δr=10m。声源频率f=50Hz,声源深度zs和接收器深度zr均为1000m,上边界(海面z0=0m)取压力释放边界条件Pz=0=0,下边界(海底zN=5000m)取Sommerfeld辐射边界条件
如图1所示,本实施例提高海洋声场预报精度的方法的步骤包括:
步骤S1.获取待预报海洋水声场的现场测量数据(本实施例为Munk波导算例的水深、声速与密度等数据)以及声源的参数信息,建立水平分层海洋环境下的柱坐标系水声Helmholtz方程,经过积分变换后得到以声压核函数为变量的深度方程;
步骤S2.根据步骤S1得到的深度方程建立声场上、下边界的声矢量;
步骤S3.将步骤S2建立的声矢量分别从上、下边界向中间深度传递,以逐层传递计算声矢量;
步骤S4.根据步骤S3传递计算的声矢量在中间深度建立联合矢量方程,并求解上、下边界的垂直振速,计算出各层声压波数核函数;
步骤S5.对步骤S4计算出的声压波数核函数进行水平波数积分,得到接收深度的声压值;
步骤S6.由步骤S5计算得到的对应不同接收深度的声压值得到接收深度的传播损失曲线,根据所述接收深度的传播损失曲线实现海洋声场的预报。
本实施例上述方法,通过建立得到声场上、下边界的声矢量后,在波数积分过程中采用分别从上、下边界向中间深度传递的方式逐层传递计算声矢量,同时在中间深度建立联合矢量方程,相比于传统波数积分方法中需从上、下两个边界分别向声源深度传递声矢量以及在声源深度利用声源条件建立联合矢量方程,更适用于深海环境声场预报,尤其是当声源距离海面较近时(下边界声矢量向上传递至声源深度的距离很远),能够大幅增强海洋声场预报过程中深度方程的求解稳定性,进而大幅提升水声场波数积分模型的计算精度,从而可以有效提高海洋声场预报的精度,减少声场预报的误差。
本实施例步骤S1中首先获取待预报海洋水声场的现场测量数据以及声源的参数信息,其中海洋水声场的现场测量数据包括海洋深度、声速、密度等数据(本实施例中分别为海深5000米、水中声速满足Munk声速剖面公式、沉积层声速1600m/s,全场密度均为1.0g/cm3),声源的参数信息包括声源频率、位置等数据(本实施例中声源频率50Hz、深度1000米),根据获取的数据建立水平分层海洋环境下的柱坐标系水声Helmholtz方程。
本实施例步骤S1的具体步骤包括:
S11.按照下式建立柱坐标系水声Helmholtz方程:
其中,P(r,z)为频率域相对声压,ρ为声传播介质密度,k=2πf/c为介质波数,f为声源频率,c为介质声速,r为水平方向的坐标,z为竖直或深度方向的坐标,zs为声源深度、δ为狄拉克函数。
S12.对水声传播介质在深度方向上划分为N层,并将每层内近似为均匀介质;在划分的每一层内,对柱坐标系水声Helmholtz方程进行Hankel变换,以将(r,z)空间的声压P(r,z)转换到(kr,z)空间,即为:
其中,φ(kr,z)为声压核函数、kr为水平波数,J0为Bessel函数;
本实施例中具体频率f=50Hz,c为介质声速,水中声速cw(z)=1500[1+0.00737(ζ-1+exp(-ζ))]、其中ζ=2(z-1300)/1300,沉积层中csed=1600m/s;r与z分别为水平与深度方向坐标;zs为声源深度且具体取zs=1000m,δ为狄拉克函数;对声场传播介质在竖直方向上划分为N=5000层,每层厚度为dz=1m,且每层内近似为均匀介质。
本实施例步骤S2具体采用四倍精度程序变量建立声场上、下边界的声矢量,具体步骤包括:
S21.将振速核函数与步骤S1得到的声压核函数形成联合矢量:
其中,φ(kr,z)为声压核函数(简称核函数);由于在任意不含声源的介质层内,声压核函数具有通解形式:
以及上边界处的联合矢量为:
S23.将步骤S22得到的上边界的振速w0、下边界的振速wN、上边界声矢量(1,B0)T、下边界声矢量(1,BN)T均设为四倍精度变量,建立声场上、下边界的联合矢量。
传统深度方程求解方法使用双倍精度程序变量,计算速度较快,如果提高变量精度会引起计算量增大以及计算时间增加的问题。本实施例通过将w0、wN、(1,B0)T、(1,BN)T设为四倍精度变量,可以减小计算机舍入误差、提高海洋声场的预报精度,并通过采用高性能计算机克服计算量增加带来的不利影响。
在本实施例中,在计算机程序中将w0、wN、(1,B0)T、(1,BN)T均设为四倍精度变量,具体:
上式中ρ0表示海面上方介质密度,在绝对软边界条件下,海面上方视为真空,因此ρ0=0。
本实施例步骤S3中,包括将声矢量从下边界zN向中间深度zh传递的步骤S31,h表示中间深度,即半水深,具体zN=5000m、zh=2500m,步骤S31的具体步骤包括:
S311.将步骤S21得到的声压核函数与振速核函数形成联合矢量并转换为矩阵形式得到:
S312.根据步骤S311转换得到的矩阵形式,得到第m层的上界面zm-1处的联合矢量为:
以及第m层的下界面zm处的联合矢量为:
vm(kr,zm-1)=Mm(kr)vm(kr,zm) (13)
其中,Mm(kr)为第m层介质由下至上的传递矩阵,若令第m层厚度为hm=zm-zm-1(本实施例中各层厚度相等,具体hm=dz=1m),则Mm(kr)的表达式为:
S314.使用矢量传递公式,从下边界zN开始,逐层传递计算声矢量,直至中间深度zh,中间深度由下方传递得到的联合矢量为:
S315.如果声源深度zs处于计算域下边界zN与中间深度zh之间,即zh<zs<zN,则在跨过声源深度之后增加声源条件,此时中间深度由下方传递得到的联合矢量计算式为:
其中,ρs表示声源深度的介质密度,角频率ω=2πf。
本实施例中声源不处于下边界深度zN与中间深度zh之间,而是处于上边界深度z0与中间深度zh之间。
本实施例步骤S3中,还包括将声矢量从上边界z0向中间深度zh传递的步骤S32,步骤S32的具体步骤包括:
S321.按照与步骤S31相同的方式,获取第m层介质由上至下的联合矢量传递公式为:
以及中间深度由上方传递得到的联合矢量为:
S322.如果声源深度zs处于计算域上边界z0与中间深度zh之间,即z0<zs<zh,则在跨过声源深度之后需要增加声源条件,此时中间深度由上方传递得到的联合矢量计算式为:
其中具体取h=2500、s=1000。
通过上述步骤,即可将声矢量从下边界zN向中间深度zh传递,以及将声矢量从上界面z0向中间深度zh传递,实现分别从上、下边界向中间深度传递计算声矢量,相比于传统方法从上、下两个边界分别向声源深度传递,更适用于深海环境声场预报,尤其是当声源距离海面较近时(下边界声矢量向上传递至声源深度的距离很远),能够大幅增强深度方程求解的稳定性,并提高波数积分模型的精度,从而提高海洋声场预报的精度。
本实施例中,步骤S4的具体步骤包括:
S41.根据步骤S31得到的由下方传递得到的联合矢量vh+1(kr,zh)与步骤S32得到的由上方传递得到的联合矢量vh(kr,zh)在中间深度相等,构建二元一次方程组:
vh+1(kr,zh)=vh(kr,zh) (20)
上述联合矢量方程中包含声压(第一分量)与振速(第二分量)两个方程以及w0(kr,z0)与wN(kr,zN)两个待定未知数,则按照求解二元一次方程组的方法,即可求解得到上边界振速w0(kr,z0)与下边界振速wN(kr,zN);本实施例中声源位于海面与中间深度之间,上式的具体形式为:
S42.在中间深度下方,使用步骤S41求解得到的下边界振速wN(kr,zN)与各层声矢量从下至上逐层计算各层界面的声压核函数φ,如其中第m(h<m≤N)层在上边界zm-1处的φm(kr,zm-1)计算式为:
在其他实施例中,如果声源深度zs处于计算域下边界zN与中间深度zh之间,即zh<zs<zN,则越过声源深度之后,将φm(kr,zm-1)计算式转换为:
S43.在中间深度上方,使用步骤S41求解得到的上边界振速w0(kr,z0)与步骤S32保存的各层声矢量从上至下逐层计算各层界面的声压核函数φ,如其中第m(1≤m≤h)层在下边界zm处的φm(kr,zm)计算式为:
在本实施例中,声源深度zs处于计算域上边界z0与中间深度zh之间,即z0<zs<zh,因此在跨过声源深度之后,将φm(kr,zm)的计算式转换为:
本实施例步骤S5中,具体利用Hankel反变换的声压积分式对声压波数核函数进行离散形式的水平波数积分,得到接收深度的声压值,具体步骤包括:
S51.Hankel反变换公式为:
其中,P(r,z)为频率域相对声压,φ(kr,z)为声压核函数,kr为水平波数,r水平方向的坐标,z为竖直或深度方向的坐标,zs为声源深度;
将Hankel反变换公式中的水平波数kr进行离散,得到声压离散式为:
其中,Δkr为水平波数步长且Δkr=2π/(rmaxnw),rmax为声场最大水平距离,nw为Bessel函数在一个2π振荡周期内最小采样点数(一般取10),kr,n=nΔkr-iεk为离散的水平波数,i为虚数单位,εk为复偏移量以防止深度方程求解过程出现奇异,εk=3Δkr/(2πlog10e),M=kmax/Δkr为离散的水平波数的最大索引号,kmax为截止波数。
S52.在声场任意位置点(r,z)处采用声压离散式进行水平波数积分,得到各位置点的声压值。
上述kmax配置时,从提高波数积分精度的角度kmax应尽量大,但受计算机变量精度的限制以及传递函数矩阵算法的影响,当kmax增大到一定值后,会引起声压核函数幅值过大而超出提前设定的变量精度存储位数的限制,使计算崩溃并中断。本实施例通过采用上述分别从上、下边界向中间深度传递的方式逐层传递计算声矢量,同时在中间深度建立并求解联合矢量方程,能够显著提高深度方程求解稳定性的同时,还能够显著提高最大截止波数,即可使水平波数取值范围更大,使声压积分值更精确,从而大幅提高声场预报精度。
为验证本发明的有效性,分别采用双精度变量、在声源深度建立联合矢量方程的传统方法与本实施例上述采用四倍精度变量、在中间深度建立联合矢量方程的方法进行海洋声场预报,声矢量传递计算方式对比如图3所示,采用双精度变量、在声源深度建立联合矢量方程的传统方法(对应图3中方法①)时,可取得的最大截止波数为kmax=1.295kref,其中参考波数kref=2πf/1500=0.2094;而按照本实施例上述方法采用四倍精度变量、在中间深度建立联合矢量方程(对应图3中方法②),可取得的最大截止波数为kmax=21.699kref。即相比于传统方法,本发明方法可取得更大的截止波数,也即本发明方法可使水平波数取值范围更大,进而使海洋声场波数积分方法的计算结果更精确。
本实施例通过采用增强深度方程求解稳定性的声矢量采用四倍精度变量与传递计算方式,稳定性较传统方法大幅提高,截止波数kmax也获得相应增加,但也带来了难以准确判断最大截止波数的困难。本实施例进一步采用预估-试算法确定最大截止波数,即首先设定一个充分大的初始预估值,在本实施例中,具体可取预估值为其中kref=2πf/1500为参考波数;然后逐渐增大水平波数kr,n=nΔkr-iεk,其中Δkr为水平波数步长且Δkr=2π/(rmaxnw),rmax为声场最大水平距离,nw为Bessel函数在一个2π振荡周期内最小采样点数(一般取10),i为虚数单位,n=1,2,3......,M,εk为复偏移量且εk=3Δkr/(2πlog10e),若当n=j时出现计算崩溃并中断为止,即可确定得到最大截止波数kmax=(j-1)Δkr-iεk。
本实施例中,步骤S5中计算指定接收器位置处声场的声压值及传播损失值时,首先根据接收器的深度zr确定所在深度上的声压值P,然后利用下式(27)计算出接收深度上的传播损失值:
TL(r,zr)=-20log10|P(r,zr)| (27)
本实施例中,步骤S5进一步根据声压计算出传播损失并绘制传播损失随水平距离变化曲线,由该传播损失曲线即可实现海洋声场的预报。在具体应用实施例中将zr=1000m的所有水平网格点上的传播损失按上式计算并绘制得到的传播损失曲线如图4所示。
如图2所示,在具体应用实施例中应用本实施例上述方法时,结合高性能计算机工作站,在高性能计算机工作站的数据存储介质中加载能够实现本实施例上述提高海洋声场预报精度的方法功能的程序模块,由高性能计算机工作站接收包括海洋深度、声速、密度等现场测量数据,以及包括声源频率、位置等声源的参数信息,经过上述提高海洋声场预报精度的方法步骤后,生成海洋声场预报图形图像,进一步还可提供给后续进行水听器声信号处理与分析。
本实施例进一步提供一种提高海洋声场预报精度的系统,包括处理器以及存储器,存储器用于存储海洋环境数据、声源参数与计算机程序,处理器用于执行所述计算机程序,以执行上述提高海洋声场预报精度的方法。本实施例系统具体可采用如图2所示结构,高性能计算机工作站配置有上述处理器及存储器。
本实施例还提供一种计算机可读存储介质,该计算机可读存储介质上存储有海洋环境数据、声源参数与被编程或配置以执行上述提高海洋声场预报精度的方法的计算机程序。
本领域内的技术人员应明白,本申请的实施例可提供为方法、系统、或计算机程序产品。因此,本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
上述只是本发明的较佳实施例,并非对本发明作任何形式上的限制。虽然本发明已以较佳实施例揭露如上,然而并非用以限定本发明。因此,凡是未脱离本发明技术方案的内容,依据本发明技术实质对以上实施例所做的任何简单修改、等同变化及修饰,均应落在本发明技术方案保护的范围内。
Claims (10)
1.一种提高海洋声场预报精度的方法,其特征在于,步骤包括:
步骤S1.获取待预报海洋水声场的现场测量数据以及声源的参数信息,建立水平分层海洋环境下的柱坐标系水声Helmholtz方程,经过积分变换后得到以声压核函数为变量的深度方程;
步骤S2.根据步骤S1得到的所述深度方程建立声场上、下边界的声矢量;
步骤S3.将步骤S2建立的所述声矢量分别从上、下边界向中间深度传递,以逐层传递计算所述声矢量;
步骤S4.根据步骤S3传递计算的所述声矢量在中间深度建立联合矢量方程,并求解上、下边界的垂直振速,计算出各层声压波数核函数;
步骤S5.对步骤S4计算出的声压波数核函数进行水平波数积分,得到接收深度的声压值;
步骤S6.由步骤S5计算得到的对应不同接收深度的声压值得到接收深度的传播损失曲线,根据所述接收深度的传播损失曲线实现海洋声场的预报。
2.根据权利要求1所述的提高海洋声场预报精度的方法,其特征在于,所述步骤S1的具体步骤包括:
S11.按照下式建立所述柱坐标系水声Helmholtz方程:
其中,P(r,z)为频率域相对声压,ρ为声传播介质密度,k为介质波数,r为水平方向的坐标,z为竖直或深度方向的坐标,zs为声源深度、δ为狄拉克函数;
S12.对水声传播介质在深度方向上划分为N层,并将每层内近似为均匀介质;在划分的每一层内,对所述柱坐标系水声Helmholtz方程进行Hankel变换,以将(r,z)空间的声压P(r,z)转换到(kr,z)空间,即为:
其中,φ(kr,z)为声压核函数、kr为水平波数,J0为Bessel函数;
3.根据权利要求1所述的提高海洋声场预报精度的方法,其特征在于,所述步骤S2具体采用四倍精度程序变量建立所述声场上、下边界的声矢量,具体步骤包括:
S21.将振速核函数与步骤S1得到的所述声压核函数形成联合矢量:
以及上边界处的联合矢量为:
S23.将步骤S22得到的上边界的振速w0、下边界的振速wN、上边界声矢量(1,B0)T、下边界声矢量(1,BN)T均设为四倍精度变量,分别建立声场上、下边界的联合矢量。
4.根据权利要求3所述的提高海洋声场预报精度的方法,其特征在于,所述步骤S3中,包括将声矢量从下边界zN向中间深度zh传递的步骤S31,所述步骤S31的具体步骤包括:
S311.将所述声压核函数与所述振速核函数形成联合矢量并转换为矩阵形式得到:
S312.根据步骤S311转换得到的矩阵形式,得到第m层的上界面zm-1处的联合矢量为:
以及第m层的下界面zm处的联合矢量为:
vm(kr,zm-1)=Mm(kr)vm(kr,zm)
其中,Mm(kr)为第m层介质由下至上的传递矩阵,若令第m层厚度为hm=zm-zm-1,则Mm(kr)的表达式为:
S314.使用所述矢量传递公式,从下边界zN开始,逐层传递计算所述声矢量,直至中间深度zh,中间深度由下方传递得到的联合矢量为:
S315.如果声源深度zs处于计算域下边界zN与中间深度zh之间,即zh<zs<zN,则在跨过声源深度之后增加声源条件,此时中间深度由下方传递得到的联合矢量计算式为:
其中,ρs表示声源深度的介质密度,角频率ω=2πf。
6.根据权利要求5所述的提高海洋声场预报精度的方法,其特征在于,所述步骤S4的步骤包括:
S41.根据步骤S31得到的由下方传递得到的联合矢量vh+1(kr,zh)与步骤S32得到的由上方传递得到的联合矢量vh(kr,zh)在中间深度相等,构建二元一次方程组:
vh+1(kr,zh)=vh(kr,zh)
求解所述二元一次方程组,求解得到所述上边界振速w0(kr,z0)与所述下边界振速wN(kr,zN);
其中,h<m≤N;
如果声源深度zs处于计算域下边界zN与中间深度zh之间,即zh<zs<zN,则越过声源深度之后,将所述φm(kr,zm-1)计算式转换为:
其中,1≤m≤h;
如果声源深度zs处于计算域上边界z0与中间深度zh之间,即z0<zs<zh,则在跨过声源深度之后,将所述φm(kr,zm)的计算式转换为:
7.根据权利要求1~6中任意一项所述的提高海洋声场预报精度的方法,其特征在于,所述步骤S5中,具体利用Hankel反变换的声压积分式对声压波数核函数进行离散形式的水平波数积分,得到所述接收深度的声压值,步骤包括:
S51.Hankel反变换公式为:
其中,P(r,z)为频率域相对声压,φ(kr,z)为声压核函数,kr为水平波数,r为水平方向的坐标,z为竖直或深度方向的坐标,zs为声源深度;
将所述Hankel反变换公式中的水平波数kr进行离散,得到声压离散式为:
其中,Δkr为水平波数步长且Δkr=2π/(rmaxnw),rmax为声场最大水平距离,nw为Bessel函数在一个2π振荡周期内最小采样点数,kr,n=nΔkr-iεk为离散的水平波数,i为虚数单位,εk为复偏移量且εk=3Δkr/(2πlog10e),M=kmax/Δkr为离散的水平波数的最大索引号,kmax为最大截止波数;
S52.在声场任意位置点(r,z)处采用所述声压离散式进行水平波数积分,得到各位置点的声压值。
8.根据权利要求7所述的提高海洋声场预报精度的方法,其特征在于,采用预估-试算法确定所述最大截止波数,具体步骤包括:首先设定一个初始预估值,然后按照式kr,n=nΔkr-iεk逐渐增大水平波数,其中n=1,2,3......,M,若当n=j时出现计算崩溃并中断,则可确定最大截止波数kmax=(j-1)Δkr-iεk。
9.一种提高海洋声场预报精度的系统,包括处理器以及存储器,所述存储器用于存储海洋环境数据、声源参数与计算机程序,所述处理器用于执行所述计算机程序,其特征在于,所述处理器用于执行所述计算机程序,以执行如权利要求1~8中任意一项所述方法。
10.一种计算机可读存储介质,其特征在于,该计算机可读存储介质上存储有海洋环境数据、声源参数与被编程或配置以执行权利要求1~8中任意一项所述提高海洋声场预报精度的方法的计算机程序。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011083585.XA CN112254797B (zh) | 2020-10-12 | 2020-10-12 | 一种提高海洋声场预报精度的方法、系统及介质 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011083585.XA CN112254797B (zh) | 2020-10-12 | 2020-10-12 | 一种提高海洋声场预报精度的方法、系统及介质 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112254797A true CN112254797A (zh) | 2021-01-22 |
CN112254797B CN112254797B (zh) | 2022-07-12 |
Family
ID=74242394
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011083585.XA Active CN112254797B (zh) | 2020-10-12 | 2020-10-12 | 一种提高海洋声场预报精度的方法、系统及介质 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112254797B (zh) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113108897A (zh) * | 2021-04-23 | 2021-07-13 | 自然资源部第三海洋研究所 | 一种基于非均匀风关声源的海洋环境噪声场预报方法 |
CN113190984A (zh) * | 2021-04-21 | 2021-07-30 | 中国海洋大学 | 水下声场模型bellhop并行实现方法 |
CN113630686A (zh) * | 2021-08-13 | 2021-11-09 | 南京工程学院 | 一种基于模式识别的高强度Helmholtz声源设计方法 |
CN114264721A (zh) * | 2021-11-26 | 2022-04-01 | 海鹰企业集团有限责任公司 | 新型声速处理系统 |
CN114510848A (zh) * | 2022-04-20 | 2022-05-17 | 自然资源部第一海洋研究所 | 一种海上风电场水下噪声计算方法、软件和测量装置 |
CN116341408A (zh) * | 2023-03-13 | 2023-06-27 | 中国人民解放军国防科技大学 | 一种距离相关水声学应用的海洋模式数据坐标转换方法 |
CN117249894A (zh) * | 2023-11-16 | 2023-12-19 | 自然资源部第一海洋研究所 | 一种水下远场声传播在海底透射厚度的诊断方法 |
Citations (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060080418A1 (en) * | 2004-10-13 | 2006-04-13 | Wayne State University | Farfield analysis of noise sources |
WO2010003837A1 (en) * | 2008-07-08 | 2010-01-14 | Brüel & Kjær Sound & Vibration Measurement A/S | Reconstructing an acoustic field |
US20150326966A1 (en) * | 2013-07-01 | 2015-11-12 | The University Of North Carolina At Chapel Hill | Methods, systems, and computer readable media for source and listener directivity for interactive wave-based sound propagation |
CN109341845A (zh) * | 2018-08-15 | 2019-02-15 | 中国人民解放军92942部队 | 一种海洋环境稳态声场空间实时仿真的方法及装置 |
CN109443516A (zh) * | 2018-12-25 | 2019-03-08 | 西北工业大学 | 一种基于噪声场垂直振速信号的海底声速被动获取方法 |
CN109489796A (zh) * | 2018-09-01 | 2019-03-19 | 哈尔滨工程大学 | 一种基于单元辐射法的水下复杂结构辐射噪声源定位识别与声辐射预报方法 |
CN109657257A (zh) * | 2017-10-12 | 2019-04-19 | 中国船舶重工集团公司第七六研究所 | 一种海洋声信道中复杂结构辐射声场计算方法 |
CN110135052A (zh) * | 2019-05-12 | 2019-08-16 | 哈尔滨工程大学 | 浅海信道下弹性结构辐射声场的计算方法 |
CN110750934A (zh) * | 2019-11-01 | 2020-02-04 | 哈尔滨工程大学 | 深海弹性结构与环境耦合声辐射预报方法 |
CN110864802A (zh) * | 2019-11-28 | 2020-03-06 | 中国舰船研究设计中心 | 基于虚拟声源波叠加的舰壳声纳平台区自噪声预报方法 |
CN111639429A (zh) * | 2020-05-29 | 2020-09-08 | 中国人民解放军国防科技大学 | 基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质 |
-
2020
- 2020-10-12 CN CN202011083585.XA patent/CN112254797B/zh active Active
Patent Citations (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060080418A1 (en) * | 2004-10-13 | 2006-04-13 | Wayne State University | Farfield analysis of noise sources |
WO2010003837A1 (en) * | 2008-07-08 | 2010-01-14 | Brüel & Kjær Sound & Vibration Measurement A/S | Reconstructing an acoustic field |
US20150326966A1 (en) * | 2013-07-01 | 2015-11-12 | The University Of North Carolina At Chapel Hill | Methods, systems, and computer readable media for source and listener directivity for interactive wave-based sound propagation |
CN109657257A (zh) * | 2017-10-12 | 2019-04-19 | 中国船舶重工集团公司第七六研究所 | 一种海洋声信道中复杂结构辐射声场计算方法 |
CN109341845A (zh) * | 2018-08-15 | 2019-02-15 | 中国人民解放军92942部队 | 一种海洋环境稳态声场空间实时仿真的方法及装置 |
CN109489796A (zh) * | 2018-09-01 | 2019-03-19 | 哈尔滨工程大学 | 一种基于单元辐射法的水下复杂结构辐射噪声源定位识别与声辐射预报方法 |
CN109443516A (zh) * | 2018-12-25 | 2019-03-08 | 西北工业大学 | 一种基于噪声场垂直振速信号的海底声速被动获取方法 |
CN110135052A (zh) * | 2019-05-12 | 2019-08-16 | 哈尔滨工程大学 | 浅海信道下弹性结构辐射声场的计算方法 |
CN110750934A (zh) * | 2019-11-01 | 2020-02-04 | 哈尔滨工程大学 | 深海弹性结构与环境耦合声辐射预报方法 |
CN110864802A (zh) * | 2019-11-28 | 2020-03-06 | 中国舰船研究设计中心 | 基于虚拟声源波叠加的舰壳声纳平台区自噪声预报方法 |
CN111639429A (zh) * | 2020-05-29 | 2020-09-08 | 中国人民解放军国防科技大学 | 基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质 |
Non-Patent Citations (3)
Title |
---|
EMERSON DE SOUSA COSTA 等: "Underwater Acoustics Modeling in Finite Depth Shallow", 《INTECNOPEN》 * |
刘巍 等: "声场波数积分最大截止波数自动选取算法", 《国防科技大学学报》 * |
刘巍 等: "水声波数积分模型稳定性与积分算法改进", 《2020中国西部声学学术交流会论文集》 * |
Cited By (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113190984A (zh) * | 2021-04-21 | 2021-07-30 | 中国海洋大学 | 水下声场模型bellhop并行实现方法 |
CN113190984B (zh) * | 2021-04-21 | 2022-10-14 | 中国海洋大学 | 水下声场模型bellhop并行实现方法 |
CN113108897A (zh) * | 2021-04-23 | 2021-07-13 | 自然资源部第三海洋研究所 | 一种基于非均匀风关声源的海洋环境噪声场预报方法 |
CN113630686A (zh) * | 2021-08-13 | 2021-11-09 | 南京工程学院 | 一种基于模式识别的高强度Helmholtz声源设计方法 |
CN114264721A (zh) * | 2021-11-26 | 2022-04-01 | 海鹰企业集团有限责任公司 | 新型声速处理系统 |
CN114264721B (zh) * | 2021-11-26 | 2024-04-19 | 海鹰企业集团有限责任公司 | 新型声速处理系统 |
CN114510848A (zh) * | 2022-04-20 | 2022-05-17 | 自然资源部第一海洋研究所 | 一种海上风电场水下噪声计算方法、软件和测量装置 |
CN116341408A (zh) * | 2023-03-13 | 2023-06-27 | 中国人民解放军国防科技大学 | 一种距离相关水声学应用的海洋模式数据坐标转换方法 |
CN116341408B (zh) * | 2023-03-13 | 2024-05-28 | 中国人民解放军国防科技大学 | 一种距离相关水声学应用的海洋模式数据坐标转换方法 |
CN117249894A (zh) * | 2023-11-16 | 2023-12-19 | 自然资源部第一海洋研究所 | 一种水下远场声传播在海底透射厚度的诊断方法 |
CN117249894B (zh) * | 2023-11-16 | 2024-04-05 | 自然资源部第一海洋研究所 | 一种水下远场声传播在海底透射厚度的诊断方法 |
Also Published As
Publication number | Publication date |
---|---|
CN112254797B (zh) | 2022-07-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112254797B (zh) | 一种提高海洋声场预报精度的方法、系统及介质 | |
WO2020228547A1 (zh) | 基于倒置式多波束回声仪的声速剖面反演方法 | |
CN106595834B (zh) | 一种获得深海大深度声场水平纵向相关性的方法 | |
CN112254798B (zh) | 一种预报海洋矢量声场的方法、系统及介质 | |
CN107179535A (zh) | 一种基于畸变拖曳阵的保真增强波束形成的方法 | |
JP6718098B2 (ja) | 位置推定装置及び方法 | |
CN106093849B (zh) | 一种基于测距和神经网络算法的水下定位方法 | |
CN110399680B (zh) | 一种浅海弹性结构辐射声场计算方法 | |
CN110132281B (zh) | 一种基于询问应答模式的水下高速目标高精度自主声学导航方法 | |
CN107016159A (zh) | 本征值确定方法及装置 | |
KR102082263B1 (ko) | 수중위치 추정 시스템 및 방법 | |
CN109579850A (zh) | 基于对水速度辅助惯导的深水智能导航方法 | |
Tu et al. | Applying the Chebyshev–Tau spectral method to solve the parabolic equation model of wide-angle rational approximation in ocean acoustics | |
JP2012202941A (ja) | 水中物体までの水平距離を算出するための水平距離算出システム及び水平距離算出方法 | |
CN110046374B (zh) | 一种基于高斯型声束的声场计算方法 | |
CN116106875A (zh) | 岸基阵坐标联合校准方法、系统、电子设备及存储介质 | |
CN114286279A (zh) | 基于全局矩阵耦合简正波的三维声场模型的声场获取方法 | |
Jian et al. | Effect of mesoscale eddies on underwater sound propagation | |
CN113218421B (zh) | 北斗拒止条件下捷联惯导系统鲁棒自适应动态对准方法 | |
CN108267743A (zh) | 基于拟合插值的快速迭代水下定位方法 | |
CN116608864A (zh) | 一种通信时延影响下基于因子图的auv协同定位方法 | |
CN113126029B (zh) | 适用于深海可靠声路径环境的多传感器脉冲声源定位方法 | |
Guo et al. | A Dynamic‐Weighted Attenuation Memory Extended Kalman Filter Algorithm and Its Application in the Underwater Positioning | |
CN103839104B (zh) | 一种海浪有效波高反演模型建模方法 | |
CN114139105A (zh) | 一种基于多项式拟合的快速声线追踪算法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |