一种海面起伏环境下的水下低频声场快速计算方法及系统

1.本发明涉及水声物理领域，尤其涉及一种海面起伏环境下的水下低频声场快速计算方法及系统。


背景技术：

2.海面作为水声波导边界的重要组成部分，其不均匀的起伏将对水下声场产生散射作用，调制水下信号特征，限制水下声场的传输能力。考虑现今水下声探测的低频化发展趋势，为实现起伏海面环境下的低频声场预报，评估海面起伏对声呐探测距离的影响程度，建立海面起伏环境下的低频声场快速计算方法具有重要意义。
3.在已有起伏海面下的声场模拟中，常采用海面散射修正下的反射损失来间接实现声场的求解，或直接根据抛物方程等模型实现非水平环境中的声场的求解。然而，间接计算方法受限于海面散射假设条件，往往无法应用于高海况等海面起伏情况，且起伏海面的细节影响无法明确表征；直接计算方法虽然不受海况限制，可直观体现海面起伏导致的声场波动等现象，但低频频段下计算效率较低、内存需求大。为此，需综合考虑起伏海面与低频声传播特征，建立海面起伏环境下的低频声场快速计算方法。


技术实现要素：

4.现有技术存在海面起伏环境下水下低频声场直接计算效率慢，计算内存需求大等缺点，本发明的目的在于克服上述技术缺陷，提出了一种海面起伏环境下的低频声场快速计算方法及系统。
5.为了实现上述目的，一种海面起伏环境下的水下低频声场快速计算方法，所述方法包括：
6.根据算子理论和前向假设，在起伏海面的波导环境下建立声压场的抛物型微分方程；通过高阶有理近似将抛物型微分方程转化为水平递推方程；
7.根据起伏海面环境与声场计算稳定性，对海洋波导垂直方向进行非均匀的离散间隔区域划分；
8.根据划分的波导离散间隔区域，结合水平递推方程计算各非均匀离散点处的算子作用表示；根据划分的波导离散间隔区域和各非均匀离散点处的算子作用表示，构建非均匀离散间隔下的水平递推矩阵；
9.根据非均匀离散间隔下的水平递推矩阵，计算得到海面起伏环境下的水下低频声场。
10.作为上述方法的一种改进，所述根据算子理论和前向假设，在起伏海面的波导环境下建立声压场的抛物型微分方程；具体为：
11.建立二维柱坐标系，以点声源在海平面的投影点为原点，海深方向为z轴的正方向，水体密度与声速为ρw和cw，海底密度、声速和吸收为ρb、cb和αb，海面由均值为且标准差为σ
srf
的高斯分布随机向量z
srf
(r)来表征，为避免出现负数坐标而取的垂直平移量，
时间因子exp(-iωt)的单频点声源位于水中且深度为zs，根据算子理论和波动方程，前向声压场在二维柱坐标系的点(r,z)下所满足的抛物型微分方程为：
[0012][0013]
其中，k0为参考波数；
[0014]
声场计算量u(r,z)与声压p(r,z)之间的关系为：
[0015][0016]
其中，exp(ik0r)为参考波数k0下提取的相位因子，ρ(r,z)为介质密度，c0为参考声速，c(r,z)为介质声速；
[0017]
令ρ＝ρ(r,z)、c＝c(r,z)，算子x满足下式：
[0018][0019]
其中，k(r,z)为介质波数，取值k1(r,z)或k2(r,z)，当k(r,z)＝k1(r,z)时表示水体的介质波数，当k(r,z)＝k2(r,z)时表示海底的介质波数。
[0020]
作为上述方法的一种改进，所述通过高阶有理近似将抛物型微分方程转化为水平递推方程；具体为：
[0021]
从原点开始，沿水平方向进行数值离散，取水平递推步长恒为δr，第l个水平离散点与原点的距离r
l
＝lδr,l＝0,1,2,...,nr，r
max
为设定的最远水平距离，表示向原点方向取整数；
[0022]
第l 1个水平离散点与原点的距离为r
l 1
，根据高阶pad
é
有理近似方法和抛物型微分方程，转化为水平递推形式，具体为：
[0023][0024]
其中，u(r
l 1
,z)为r
l 1
处的声场计算量，u(r
l
,z)为r
l
处的声场计算量，和为第j
p
，阶pad
é
有理近似系数，j
p
＝1,...,n
p
，n
p
为pad
é
有理近似阶数；
[0025]
根据分裂-步进理论，将高阶pad
é
有理近似方法下的n
p
阶乘法运算转化为n
p
步的循环水平递推形式，根据u(r
l
,z)计算u(r
l 1
,z)时，循环水平递推形式为：
[0026][0027]
其中，循环水平递推的输入量由u(r
l
,z)给出，令u0(r
l
,z)＝u(r
l
,z)；根据第一步水平递推形式由u0(r
l
,z)得到第一步输出量u1(r
l
,z)，再将u1(r
l
,z)作为下一步水平递推形式的输入，得到第二步输出量u2(r
l
,z)，以此类推直至得到第n
p
步输出量第n
p
步输出量等于第l 1个水平离散点的声场计算量u(r
l 1
,z)，即
[0028]
作为上述方法的一种改进，所述根据起伏海面环境与声场计算稳定性，对海洋波导垂直方向进行非均匀的离散间隔区域划分；具体为：
[0029]
根据起伏海面环境与声场计算稳定性，对海洋波导垂直方向上进行非均匀的离散间隔区域划分，得到海面密集离散区、过渡离散区和稀疏离散区；其中，所述海面密集离散区的离散间隔为δzs；所述过渡离散区的离散间隔由δzs至aδzs线性递增，所述稀疏离散区的离散间隔为aδzs，a为离散间隔倍数。
[0030]
作为上述方法的一种改进，所述根据划分的波导离散间隔区域，结合水平递推方程计算各非均匀离散点处的算子作用表示；具体为：
[0031]
在第l个水平离散点开始的第j
p
阶pad
é
级数处，根据划分的波导垂直离散间隔区域和水平递推方程，对于z1和两垂直离散点处不考虑算子作用，令这两点处的声场计算量当水平递推方程是左侧计算量，下标当水平递推方程是右侧计算量，下标对于z1和以外的其余各垂直离散点zj,j＝2,3,...,n
z-1处的算子作用，根据推导表示为：
[0032][0033][0034][0035]
其中，δzj＝z
j 1-zj，δzj表示第j个垂直离散点处的离散间隔；ρj＝ρ(r
l
,zj)，表示(r
l
,zj)处的介质波数平方差值，cj＝c(r
l
,zj)，表示(r
l
,zj)处的介质声速，表示第j个垂直离散点处的介质波数平方差值。
[0036]
作为上述方法的一种改进，所述根据划分的波导离散间隔区域和各非均匀离散点处的算子作用表示，构建非均匀离散间隔下的水平递推矩阵；具体为：
[0037]
第l个水平离散点，根据划分的非均匀波导垂直离散间隔区域和各非均匀离散点处的算子作用表示，将每阶pad
é
步骤下的水平递推方程转化为三对角矩阵形式，三对角矩阵形式方程为：
[0038][0039]
其中，j
p
＝1,...,n
p
为pad
é
有理近似的各阶数，系数矩阵顶部元素为：有理近似的各阶数，系数矩阵顶部元素为：系数矩阵底部元素为：
[0040]
z＝z2至海面z
srf
(r
l
)各离散点zj,j＝2,3,...,n
zsrf
处的系数矩阵顶部元素为：系数矩阵底部元素为：系数矩阵底部元素为：
[0041]
海面下至截断边界上方各离散点zj,j＝n
zsrf
1,...,n
z-1处的系数矩阵顶部元素为：
[0042]
[0043][0044][0045]
系数矩阵底部元素为：
[0046][0047][0048][0049]
作为上述方法的一种改进，所述根据非均匀离散间隔下的水平递推矩阵，计算得到海面起伏环境下的水下低频声场；具体为：
[0050]
根据抛物方程自初始场理论或简正波理论，由第1个水平离散点与原点的水平距离为r1，各非均匀垂直离散点zj,j＝1,2,...,nz的声压p(r1,zj)，得到初始声场计算量u(r1,zj)：
[0051][0052]
作为三对角矩阵表征下循环水平递推方程的输入量u0(r1,zj)＝u(r1,zj)，进行n
p
步循环递推求解，得到第2个水平离散点与原点的水平距离为r2的声场计算量并以u(r2,zj)作为下一水平递推输入值进行循环求解，依次求解得到各水平离散点处的声场计算量u(r
l
,zj),l＝3,4,...nr，直到计算至设定的最远水平距离；
[0053]
根据u(r
l
,zj),l＝1,2,...nr,j＝1,2,...,nz，(r
l
,zj)处的声压p(r
l
,zj)为：
[0054][0055]
一种海面起伏环境下的水下低频声场快速计算系统，所述系统包括：声压场方程建立模块、海洋波导区域划分模块、水平递推矩阵建立模块和水下低频声场求解模块；其中，
[0056]
所述声压场方程建立模块，用于根据算子理论和前向假设，在起伏海面的波导环境下建立声压场的抛物型微分方程；通过高阶有理近似将抛物型微分方程转化为水平递推方程；
[0057]
所述海洋波导区域划分模块，用于根据起伏海面环境与声场计算稳定性，对海洋波导垂直方向进行非均匀的离散间隔区域划分；
[0058]
所述水平递推矩阵建立模块，用于根据划分的波导离散间隔区域，结合水平递推方程计算各非均匀离散点处的算子作用表示；根据划分的波导离散间隔区域和各非均匀离散点处的算子作用表示，构建非均匀离散间隔下的水平递推矩阵；
[0059]
所述水下低频声场求解模块，用于根据非均匀离散间隔下的水平递推矩阵，计算得到海面起伏环境下的水下低频声场。
[0060]
与现有技术相比，本发明的优势在于：
[0061]
本发明采用非均匀空间网格下的离散技术，建立了适应于起伏海面的低频声场计算方法，既保证了起伏海面附近的空间精细划分与声场计算精度，又可考虑低频声波波长特征，大幅度降低空间离散数目，提高计算效率，降低计算时间与计算内存。
附图说明
[0062]
图1是本发明实施例1的具有起伏海面的海洋波导环境模型示意图；
[0063]
图2是本发明实施例1的水体声速剖面图；
[0064]
图3是本发明实施例1中声场损失分布计算结果；
[0065]
图4是本发明实施例1中接收深度50m处的传播损失曲线对比图。
具体实施方式
[0066]
本发明公开了一种海面起伏环境下的水下低频声场快速计算方法，该计算方法主要在依据非均匀空间网格下的galerkin离散形式建立适应于海面起伏环境下的抛物方程声场计算模型，模型主要包括：以抛物方程方法和高阶有理近似理论建立海面起伏波导环境下的前向声场水平递推方程；依据算子理论和海面起伏环境下声场快速计算需求，从空间数值离散形式为基本出发点解决海面起伏对声场快速求解的限制，即以较密集的离散步长处理起伏海面及其附近区域，并采用较稀疏的离散步长处理波导其他区域，同时在步长密集与稀疏区域之间设置过渡层来防止计算步长突变导致的计算不稳定的产生，实现对角形式的水平递推矩阵构建，到达声场快速求解的目的；起伏环境下的仿真结果表明，该计算方法在保证海面微弱起伏影响得到充分考虑的同时，又可满足远距离低频、甚低频声场快速计算的需要，克服了起伏海面环境下声场计算时间过长、计算内存过大等问题。
[0067]
下面结合附图和实施例对本发明的技术方案进行详细的说明。
[0068]
实施例1
[0069]
本发明的实施例1提出了一种海面起伏环境下的水下低频声场快速计算方法。该起伏海面的波导环境模型如图1所示：建立二维柱坐标系，以点声源在海平面的投影点为原点，海深方向为z轴的正方向，水体密度与声速为ρw和cw，海底密度、声速和吸收为ρb、cb和αb，海面由均值为且标准差为σ
srf
的高斯分布随机向量z
srf
(r)来表征，为避免出现负数坐标而取的垂直平移量，海底由实际情况给出，时间因子exp(-iωt)的单频点声源位于水中且深度为zs，则计算量u(r,z)所满足的抛物型微分方程为：
[0070][0071]
声场计算量u(r,z)与声压p(r,z)之间的关系为：
[0072][0073]
其中，exp(ik0r)为参考波数k0下提取的相位因子，ρ(r,z)为介质密度，c0为参考声速，c(r,z)为介质声速；
[0074]
令ρ＝ρ(r,z)、c＝c(r,z)，算子x满足下式：
[0075][0076]
其中，k(r,z)为介质波数，取值k1(r,z)或k2(r,z)，当k(r,z)＝k1(r,z)时表示水体的介质波数，当k(r,z)＝k2(r,z)时表示海底的介质波数。
[0077]
求解过程中，取水平离散点r
l
＝lδr(l＝1,2,...,nr)，水平步长恒为δr，一般水平步长δr≤5λ，λ为声波长，根据分裂-步进理论，高阶pad
é
有理近似方法下的n
p
阶乘法运算转化为n
p
步的循环水平递推形式，在根据水平距离r
l
处声场计算量u(r
l
,z)计算r
l 1
处声场计算量u(r
l 1
,z)时，循环水平递推形式为：
[0078][0079]
其中，循环水平递推的输入量由r
l
处声场计算量u(r
l
,z)给出，即u0(r
l
,z)＝u(r
l
,z)；根据第一步水平递推形式由u0(r
l
,z)得到第一步输出量u1(r
l
,z)，再将u1(r
l
,z)作为下一步水平递推形式的输入，得到第二步输出量u2(r
l
,z)，以此类推直至得到第n
p
步输出量第n
p
步输出量等于r
l 1
处声场计算量u(r
l 1
,z)，即,z)，即和为第j
p
阶pad
é
有理近似系数。
[0080]
对于水平递推方程中算子x的作用，常规可在均匀步长下采用有限差分或galerkin(一维有限元)等方法数值离散直接求解。
[0081]
根据起伏海面环境与声场计算稳定性，对海洋波导垂直方向上进行非均匀的离散间隔区域划分，得到海面密集离散区、过渡离散区和稀疏离散区；其中，所述海面密集离散区的离散间隔为δzs；所述过渡离散区的离散间隔由δzs至aδzs线性递增，所述稀疏离散区的离散间隔为aδzs，a为离散间隔倍数。
[0082]
在水平离散点r
l
＝lδr(l＝1,2,...,nr)和第j
p
(j
p
＝1,...,n
p
)阶pad
é
级数处，根据划分的波导垂直离散间隔区域和水平递推方程，除由于z1和两垂直离散点处不考虑
算子作用而令下标或j
p-1，由水平递推方程左侧和右侧计算量决定外，其余各垂直离散点zj,j＝2,3,...,n
z-1处的算子作用根据推导，表示为：
[0083][0084][0085][0086]
其中，δzj＝z
j 1-zj，δzj表示第j个离散点处的离散间隔；ρj＝ρ(r
l
,zj)，表示第j个离散点处的介质波数平方差值；cj＝c(r
l
,zj)，表示第j个离散点处的介质声速。表示第j个离散点处的介质波数平方差值。
[0087]
本实施例对于算子x的作用，采用非均匀空间网格下的离散技术进行处理，并根据海面起伏情况和声场计算稳定性需要，将波导垂直方向上划分为三类：
[0088]
(一)范围[0,max(z
srf
) 5δzs]内取作海面密集离散区，垂直离散间隔为δzs；
[0089]
(二)范围(max(z
srf
) 5δzs,max(z
srf
) 205δzs]内取作过渡离散区，垂直离散间隔由δzs线性递增至aδzs；
[0090]
(三)范围(max(z
srf
) 205δzs,z
max
]内取作稀疏离散区，垂直离散间隔为aδzs；其中，z
max
为计算深度最大值，a为离散倍数，为保证起伏海面作用考虑充分和声场计算精度，一般垂直间隔δzs≤z
srfm
/10且aδzs≤λ/10，λ为声波长。
[0091]
对于各水平离散点r
l
＝lδr(l＝1,2,...,nr)，根据划分的非均匀波导垂直离散间隔区域和各非均匀离散点处的算子作用表示，将每阶pad
é
步骤下的水平递推方程化为三对角矩阵形式，所化成的三对角矩阵形式方程为：
[0092][0093]
其中，j
p
＝1,...,n
p
为pad
é
有理近似的各阶数，系数矩阵顶部元素为：有理近似的各阶数，系数矩阵顶部元素为：系数矩阵底部元素为：
[0094]
z＝z2至海面z
srf
(r
l
)各离散点(j＝2,3,...,n
zsrf
)处的系数矩阵元素为：
[0095][0096]
海面下至截断边界上方各离散点(j＝n
zsrf
1,...,n
z-1)处的系数矩阵元素为：
[0097][0098]
[0099][0100][0101][0102][0103]
对于该三对角矩阵形式的水平递推方程，在抛物型自初始场理论或其他方法给定的初始场下，由于非均匀网格技术的引入使得矩阵维数远小于常规方法，可依据追赶法等基本数值计算方法进行快速求解。
[0104]
根据非均匀空间离散方法，建立易于求解且计算稳定的三对角形式水平递推矩阵方程；对于海面起伏情况，由于非均匀网格技术的采用，导致水平递推矩阵的维数要小于均匀情况很多，使得计算精度在得到保障下的同时，计算时间、计算内存皆较大幅度下降；对于该三对角形式水平递推方程，在抛物型自初始场理论或其他方法给定的初始场下，可依据追赶法等基本数值计算方法进行求解，也可采用多种并行的矩阵方程计算方法进行求解。
[0105]
考虑具有起伏海面和表层声道的深海波导情况，海面由均值为16m且标准差为4m的高斯分布随机向量z
srf
(r)来表示，即垂直平移量为16m，平移前的声源和接收深度皆为50m，平移前的海深恒为5000m，水体声速剖面如图2所示、密度1.0g/cm3，海底声速为cb＝1520 0.13(z-5016)、密度为1.5g/cm3、吸收为0.25db/λ，声场计算频率为200hz，水平步长为δr＝10m，垂直离散间隔为δzs＝0.1m且离散倍数a＝10，水平距离最大值r
max
＝35km，计算深度最大值z
max
＝7000m，声场传播损失分布如图3所示，均匀空间网格与非均匀空间网格下接收深度处的传播损失对比结果如图4所示。
[0106]
由对比结果可以看出，本发明所建立的声场计算方法具有较高的准确性；同时，在相同的计算机条件和矩阵求解方法下，采用编程语言编写，非均匀空间网格情况下的计算时间近乎为均匀网格情况的1/10，内存近乎为均匀网格情况的1/5左右，明显提高了起伏海面环境下低频声场的计算效率。
[0107]
起伏环境下的仿真结果表明，该计算方法在保证海面微弱起伏影响得到充分考虑的同时，又可满足远距离低频、甚低频声场快速计算的需要，克服了起伏海面环境下声场计算时间过长、计算内存过大等问题。
[0108]
实施例2
[0109]
本发明的实施例2提出了一种海面起伏环境下的水下低频声场快速计算系统，所述系统包括：声压场方程建立模块、海洋波导区域划分模块、水平递推矩阵建立模块和水下低频声场求解模块；本系统的具体实现同实施例1的方法。
[0110]
所述声压场方程建立模块，用于根据算子理论和前向假设，在起伏海面的波导环境下建立声压场的抛物型微分方程；通过高阶有理近似将抛物型微分方程转化为水平递推方程；
[0111]
所述海洋波导区域划分模块，用于根据起伏海面环境与声场计算稳定性，对海洋波导垂直方向进行非均匀的离散间隔区域划分；
[0112]
所述水平递推矩阵建立模块，用于根据划分的波导离散间隔区域，结合水平递推方程计算各非均匀离散点处的算子作用表示；根据划分的波导离散间隔区域和各非均匀离散点处的算子作用表示，构建非均匀离散间隔下的水平递推矩阵；
[0113]
所述水下低频声场求解模块，用于根据非均匀离散间隔下的水平递推矩阵，计算得到海面起伏环境下的水下低频声场。
[0114]
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。