CN111639429B - 基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质 - Google Patents

Abstract

本发明公开了一种基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质，本发明方法包括建立柱坐标系下声学传播Helmholtz方程的简化控制方程；将简化控制方程完成z坐标向x坐标的自变量变换；对简化控制方程实施切比雪夫正变换，将物理空间映射到谱空间上，形成谱空间中离散化后的线性方程组；求解谱空间中离散化后的线性方程组，得到谱空间上的解；对谱空间上的解实施切比雪夫反变换，将谱空间上的解映射回物理空间；完成x坐标向z坐标的自变量反变换获得u(r，z)并求出声压场p；计算传播损失。本发明适合在使用较少离散网格点的条件下获得更高的计算精度，能够充分发挥硬件平台的计算性能，显著加快数值模拟的速度。

Description

基于切比雪夫多项式谱的水下声场数值模拟方法、系统及 介质

技术领域

本发明涉及水下声学传播技术，具体涉及一种基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质。

背景技术

声波是水下远程传播信息的主要载体。在地、空、天领域中得到广泛应用的其他通信手段，往往无法有效用在水下，例如：电磁波在水中会迅速吸收衰减，无法实现远距离传输。充分认识与正确掌握声波在水下的传播规律，是十分重要且必要的。以海洋为例，覆盖地球表面积约71％的海洋中含有丰富的能源、矿产和生物资源，对海洋资源的保护和开发离不开对水下各类目标的探测识别、定位通信等技术，这就要求对作为探测主要载体的声波在水下的传播、以及特定声场的分布等规律有着清晰的认识与掌握。此外，掌握水下声场规律，对于内河、湖泊、近海等水域上的农林渔业生产等也具有重要支撑作用。

由于水体、水面、水底等环境复杂多变，声波在水下并不是直线传播的，通常有着复杂的传播路径，声能量在空间上也呈现复杂的分布状态。从水下声场的原理上看，特定频率的声波在水下传播时，满足Helmholtz偏微分方程，若再辅以一定的边界条件，则可获得整个空间上的声场解。然而，在绝大多数现实条件下，要获得数学上的解析解是不可能的。

传统上大多通过数值计算的方法获得近似解，目前采用数值计算方法求解Helmholtz方程，目前有直接求解和简化求解两大类方法。

(1)直接求解法直接针对原始Helmholtz方程，采用各种数值离散化手段，完成数值求解。直接求解法很少在实际应用中使用，原因有三：一是Helmholtz方程是一类椭圆型方程，实践中的声学问题很难给出其定解所需要的所有边界条件，二是后期形成的大型代数方程组很难收敛求解，三是对于传播距离较远的海洋环境，求解需要较细的离散网格，造成了巨大的存储和计算开销，从而限制了其实用化。相对而言，简化求解法是现实中普遍使用的方法，具有简化易理解、计算量适中等优点，但由于其引入的额外假设条件互不相同，因此适用范围、近似误差也大相径庭，使用时需仔细选择。

(2)简化求解法则是针对不同的水声环境条件，对原始的Helmholtz方程进行二次简化，形成更容易数值求解的新方程。简化求解法是目前实践中广泛使用的一类方法。根据二次简化后方程的不同，简化求解法又可进一步划分为简正波方法、射线方法、抛物化方程法、快速场方法及各种复合方法等。简正波方法建立在与距离无关的假设基础之上，通过求解与深度有关的方程，获得一组振动模式，并按照将声源深度加权、将每个模式的贡献累加起来构成总声场。射线方法建立在高频近似的基础之上，只适用于分析高频声场问题，对低频及一些焦散问题处理困难，且对于非水平分层介质求解过程较烦琐。快速场方法按简正波近似方法来分离波动方程参数，并用快速傅里叶变换算法计算最终结果。抛物化方程法通过忽略反向声波、引入远场近似等处理，将原始Helmholtz方程简化成一个单向传播的抛物方程。不论直接求解法还是各种简化求解法，都需要借助数值离散化手段，将原始的微分或积分方程转化为计算机上方便求解的代数方程组，数值离散化手段包括有限差分法、有限体积法和有限元法等。其中，目前在水下声场计算中最为普及的是有限差分法。

但是，在数值离散化手段方面，有限差分法是目前主流的离散化方法，具有易于理解、实施简单的优点，其主要不足包括：一是精度不高，常用的二阶差分格式难以满足某些声场变化剧烈区域的精度需求，采用更高阶差分格式可以弥补些差距，但却带来计算过程更复杂、需补充额外边界条件等缺点；二是采用等距网格的方式难以同时顾及物理量变化平滑和变化剧烈的不同区域，为了提高计算空间分辨率，往往需要缩小网格间距，从而急剧增加存储和计算开销，声场计算速度无法满足应用需求。

谱(切比雪夫谱)方法是在上世纪六、七十年代发展起来的一种数值求解偏微分方程的方法，至今已成为继有限差分法和有限元法之后又一种重要的数值方法。谱方法把问题的解在一组光滑基函数下进行展开，形成解的谱展开式，再根据这个展开式和原方程来求解展开式系数的方程组。切比雪夫多项式谱方法的最大优势在于它的精度可直接由级数展开式的项数来决定，因而常被视为具有“无穷阶”收敛性的方法。若取切比雪夫多项式作为光滑基函数，则称其为切比雪夫多项式谱方法，或简称为切比雪夫谱方法。切比雪夫多项式谱方法在计算流体力学、最优控制理论、非线性动力系统、地震波研究和数值天气预报等问题中得到了广泛应用。因此将切比雪夫多项式谱(简称切比雪夫谱)方法引入到水下声场模拟，有望解决现有水下声场模拟的问题。

发明内容

本发明要解决的技术问题：针对现有技术的上述问题，提供一种基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质，本发明实现了将切比雪夫多项式谱应用到求解水声场模拟的计算过程中的数值离散环节中，可在使用较少离散网格点的条件下获得更高的计算精度，能够充分发挥硬件平台的计算性能，显著加快数值模拟的速度。

为了解决上述技术问题，本发明采用的技术方案为：

一种基于切比雪夫多项式谱的水下声场数值模拟方法，该方法包括：

1)针对具体水声环境参数，建立柱坐标系下声学传播Helmholtz方程的简化控制方程；

2)将简化控制方程完成z坐标向x坐标的自变量变换；

3)对简化控制方程实施切比雪夫正变换，将物理空间映射到谱空间上，形成谱空间中离散化后的线性方程组；

4)求解谱空间中离散化后的线性方程组，得到谱空间上的解；

5)对谱空间上的解实施切比雪夫反变换，将谱空间上的解映射回物理空间；

6)完成x坐标向z坐标的自变量反变换，获得水平方向距离为r、深度为z位置处的缩比分量函数u的函数值u(r，z)，并利用该函数值u(r，z)求出声压场p；

7)计算指定接收器位置处声场的声压值及传播损失值。

可选地，步骤1)中的简化控制方程为下式所示的方程：

上式中，i为虚数单位，k₀为参考波数，Δr为r方向离散网格点的间距，L是Padé级数的阶数，α_j为Padé级数的分子多项式的第j项系数，β_j为Padé级数的分母多项式的第j项系数，X为算子，u^(m)＝u(r_m，z)为缩比分量函数u在水平方向距离为r_m、深度为z位置处的函数值，u^(m+1)＝u(r_m+1，z)为缩比分量函数u在水平方向距离为r_m+1、深度为z位置处的函数值，r_m＝m·Δr为r方向第m个离散网格点的坐标，r_m+1＝(m+1)·Δt为r方向第m+1个离散网格点的坐标，m＝0，1，2，…；其中算子X的函数表达式如下式所示：

上式中，n是介质折射率，算子X的函数表达式由第一个算子项、第二个算子项两者之和构成。

可选地，步骤3)的详细步骤包括：

3.1)确定数值计算对切比雪夫反变换的函数表达式中的无穷项级数近似的截断项数N；

3.2)给定r＝0处的初场边界条件：使用经典的高斯声源法、Greene声源法、简正波声源、镜像声源、自启动型声源中的一种方法给定r＝0处的u⁽⁰⁾值，其中r为水平方向距离，u⁽⁰⁾为缩比分量函数u在水平方向距离r＝0处的函数取值；根据r＝0处的u⁽⁰⁾值，利用切比雪夫正变换的函数表达式，计算得到谱空间上的N+1项谱系数

其中N为无穷项级数近似的截断项数，

分别为缩比分量函数u在水平方向距离r为0处的N+1项谱展开系数；设定水平方向距离上的离散网格点计数变量m的值为0；

3.3)获得自变量变换后的算子X中第一个算子项在谱空间中表示的矩阵D，获得自变量变换后的算子X中第二个算子项在谱空间中表示的矩阵E，将矩阵D和矩阵E相加，得到算子矩阵X；

3.4)计算Padé级数的分子多项式的第j项系数α_j和Padé级数的分母多项式的第j项系数β_j，基于算子矩阵X装配谱空间上线性方程组的左端系数矩阵A及右端系数矩阵B：

上式中，i为虚数单位，k₀为参考波数，Δr为r方向离散网格点的间距，L是Padé级数的阶数，α_j为Padé级数的分子多项式的第j项系数，β_j为Padé级数的分母多项式的第j项系数，I为单位矩阵，X为算子矩阵；

3.5)加入合适的水面、水底边界条件，修正左端系数矩阵A及右端系数矩阵B，从而得到最终在谱空间中离散化后的线性方程组：

上式中，A为左端系数矩阵，

为缩比分量函数u在水平方向距离为r_m+1处的谱展开系数，B为右端系数矩阵，

为缩比分量函数u在水平方向距离为r_m处的谱展开系数。

可选地，所述获得自变量变换后的算子X中第一个算子项在谱空间中表示的矩阵D的详细步骤包括：

3.3.1A)取(N+1)×(N+1)阶的方阵D，将其所有元素初始化为零值；

3.3.2A)设方阵D的行下标及列下标均取从0到N的非负整数；对于其第k行的第j列元素D[k][j]而言，若k与j的奇偶性相反且j＞k，则重新置D[k][j]的值为2j，其中k＝0，1，…，N；j0，1，…N；

3.3.3A)将矩阵D中第1行的每个元素乘以0.5，所得结果原地保存于矩阵D；

3.3.4A)按照矩阵相乘的规则，计算乘积矩阵D×D，并将结果原地保存在矩阵D中；

3.3.5A)将矩阵D中每个元素乘以4/H²，所得结果仍原地保存于矩阵D，从而得到第一个算子项在谱空间中表示的矩阵D。

可选地，所述获得自变量变换后的算子X中第二个算子项在谱空间中表示的矩阵E的详细步骤包括：

3.3.1B)记物理空间上的函数

为v(x)，利用切比雪夫正变换的函数表达式求出相应的谱空间系数

k＝0，1，…，N；

3.3.2B)取(N+1)×(N+1)阶的方阵E，将其所有元素初始化为零值；

3.3.3B)设方阵E的行下标及列下标均取从0到N的非负整数；对于其第k行的第j列元素E[k][j]而言：若0≤k-j≤N，则将E[k][j]在原值基础上增加

若0≤j-k≤N，则将E[k][j]在原值基础上增加

若0≤j+k≤N，则将E[k][j]在原值基础上增加

若上述条件同时满足多个，则每次都执行相应的增加操作，其中k＝0，1，…，N；j＝0，1，…N；

3.3.4B)将矩阵E中所有元素乘以0.5，所得结果原地保存于矩阵E，从而得到第二个算子项在谱空间中表示的矩阵E；

可选地，步骤3.5)中修正左端系数矩阵A及右端系数矩阵B的步骤包括：

3.5.1)利用切比雪夫反变换的函数表达式，将边界条件

转换成谱空间的方程式

此处a_k＝(-1)^k，其中

为缩比分量函数u在水平方向距离为r_m、深度为0处的取值，N为无穷项级数近似的截断项数，k为缩比分量函数u的展开系数的序号，

为缩比分量函数u在水平方向距离为r_m处的第k个展开系数；

3.5.2)修正左端系数矩阵A第N-1行的第k个元素取为a_k，右端系数矩阵B第N-1行的第k个元素取值为0；

3.5.3)利用切比雪夫反变换的函数表达式，将边界条件

转换成谱空间的方程式

此处a_k＝1；

3.5.4)修正左端系数矩阵A第N-1行的第k个元素取为a_k，右端系数矩阵B第N-1行的第k个元素取值为0。

可选地，步骤4)的详细步骤包括：

4.1)使用标准的线性方程组求解算法解得谱空间中离散化后的线性方程组中第m+1个距离上的谱系数解

且在求解方程组的计算区域利用众核CPU和GPU进行协同加速计算；

4.2)判断水平方向距离达到预定范围r_max是否成立，若成立则跳转执行步骤5)；否则，令m＝m+1，跳转执行步骤4.1)继续求解谱空间中离散化后的线性方程组。

可选地，步骤5)的详细步骤包括：针对每个m＝1，2，…利用切比雪夫反变换将谱空间的定义域x∈[-1，1]离散成指定分辨率的网格，第j个格点的坐标为x_j，利用z_j＝((1-x_j)H)/2将谱空间上的定义域变换为物理空间上的定义域，得到各深度对应的声压，其中z_j为将谱空间上的坐标x_j解映射回物理空间后的坐标，H为目标海域的深度。

此外，本发明还提供一种基于切比雪夫多项式谱的水下声场数值模拟系统，包括计算机设备，该计算机设备被编程或配置以执行所述基于切比雪夫多项式谱的水下声场数值模拟方法的步骤，或者该计算机设备的存储器上存储有被编程或配置以执行所述基于切比雪夫多项式谱的水下声场数值模拟方法的计算机程序。

此外，本发明还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行所述基于切比雪夫多项式谱的水下声场数值模拟方法的计算机程序。

和现有技术相比，本发明具有下述优点：

本发明实现了将切比雪夫多项式谱应用到求解水声场模拟的计算过程中的数值离散环节中，可在使用较少离散网格点的条件下获得更高的计算精度，能够充分发挥硬件平台的计算性能，显著加快数值模拟的速度。

附图说明

图1为本发明实施例方法的基本流程示意图。

图2为本发明实施例得到的全空间场中传播损失分布图

图3为本发明实施例得到的在深度为1000m处的传播损失曲线图。

具体实施方式

下文将以Munk声速剖面的理想流体波导为例，对本发明进行进一步的详细说明。本例中海水密度均匀，恒为ρ≡1.0g/cm³，声速取为Munk声速剖面，海底水平且海深H为5000米，边界条件取海面海底压力释放条件，即p(z＝０)＝p(z＝H)＝0。声源和接收器均位于z_s＝z_r＝1000米深处，声源角频率为ω＝2×π×50Hz。取谱截断阶数为N＝200，r方向步长为Δr＝10m，最大求解距离为r_max＝100km，Padé级数的项数取为L＝4，以自启动场作为初始条件。

如图1所示，本实施例基于切比雪夫多项式谱的水下声场数值模拟方法包括：

1)针对具体水声环境参数，建立柱坐标系下声学传播Helmholtz方程的简化控制方程；

2)将简化控制方程完成z坐标向x坐标的自变量变换；

3)对简化控制方程实施切比雪夫正变换，将物理空间映射到谱空间上，形成谱空间中离散化后的线性方程组；

4)求解谱空间中离散化后的线性方程组，得到谱空间上的解；

5)对谱空间上的解实施切比雪夫反变换，将谱空间上的解映射回物理空间；

6)完成x坐标向z坐标的自变量反变换获得水平方向距离为r、深度为z位置处的缩比分量函数u的函数值u(r，z)，并利用该函数值u(r，z)求出声压场p；

7)计算指定接收器位置处声场的声压值及传播损失值。

本实施例中，步骤1)中的简化控制方程为下式所示的半离散形式的方程：

上式中，i为虚数单位，k₀为参考波数，Δr为r方向离散网格点的间距(步长)，L是Padé级数的阶数，α_j为Padé级数的分子多项式的第j项系数，β_j为Padé级数的分母多项式的第j项系数，X为算子，u^(m)＝u(r_m，z)为缩比分量函数u在水平方向距离为r_m、深度为z位置处的函数值，u^(m+1)＝u(r_m+1，z)为缩比分量函数u在水平方向距离为r_m+1、深度为z位置处的函数值，r_m＝m·Δr为r方向第m个离散网格点的坐标，r_m+1＝(m+1)·Δr为r方向第m+1个离散网格点的坐标，m＝0，1，2，…；其中算子X的函数表达式如下式所示：

上式中，n是介质折射率，算子X的函数表达式由第一个算子项、第二个算子项两者之和构成。

其中，式(1)所示的方程的推导过程如下：

在柱坐标系下，刻画声波传播规律的原始齐次Helmholtz方程可表示为：

上式中，p＝p(r，z)为水平方向距离为r、深度为z位置处的声压值，k₀＝ω/c₀是参考波数，ω为声源频率，n＝c₀/c是介质折射率，c＝c(z)是声速，c₀为预先确定的参考声速。要计算的最远水平方向距离记作r_max。本实施例中c(z)按照Munk声速剖面模型取值，其一般形式为：

上式中，参数ε＝0.00737，中间变量

参考声速取为c₀＝1500m/s。

引入分解变换：

上式中，u(r，z)为水平方向距离为r、深度为z位置处的缩比分量函数u的函数值，r为水平方向距离。引入分解变换后，再基于远场假设及Padé级数近似原理，则上述Helmholtz方程可简化为抛物方程：

上式中，i是虚数单位，k₀为参考波数，r为水平方向距离，Δr为r方向离散网格点的间距(步长)，L是Padé级数的阶数，α_j为Padé级数的分子多项式的第j项系数，β_j为Padé级数的分母多项式的第j项系数，X为算子，见式(2)，u为缩比分量函数。

步骤2)完成z坐标向x坐标的变换的目的在于方便后续使用切比雪夫谱方法。本实施例中，步骤2)的详细步骤包括：针对简化控制方程引入线性变换x＝1-2z/H，其中H为目标海域的深度，将原始问题的定义域从z∈[0，H]变换到x∈[-1，1]，形成r-x坐标系下新的微分方程，该微分方程与式(1)的形式相同但含义上有差别：u^(m)＝u(r_m，x)，且算子X的函数表达式如下式所示：

上式中，H为目标海深，k₀是参考波数，n是介质折射率，算子X的函数表达式由第一个算子项、第二个算子项两者之和构成。

切比雪夫变换的定义如下：切比雪夫变换将一个定义在x∈[-1，1]上的光滑函数f(x)，与离散谱空间上的一个无穷维向量

相对应。当使用有限的N维进行近似时，二者的关系用数学式表示为：

上式中，f(x)为待展开的函数，

为函数f(x)经切比雪夫多项式展开后得到的第k个系数，T_k(x)为第k阶切比雪夫多项式，N为无穷项级数近似的截断项数，d_k为Gauss-Lobatto节点常系数，f(x_j)为函数f(x)在x＝x_j处的值，T_k(x_j)为第k阶切比雪夫多项式在x＝x_j处的值，ω_j为权重系数。其中，式(8)为切比雪夫正变换，式(7)为切比雪夫反变换。

T_k(x)和T_k(xj)的函数表达式为：

T_k(x)＝cos(kcos^-1(x)) (9)

表示第k阶切比雪夫多项式，k＝0，1，2，…。

权重系数ω_j的取值为：当1≤j≤N-1时，

当j＝0或j＝N时，

切比雪夫正变换中的离散节点取Gauss-Lobatto节点为

j＝0，1，2，…。d_k的取值为：当1≤k≤N-1时，

当j＝0或j＝N时，d_k＝π。

本实施例中，步骤3)的详细步骤包括：

3.1)确定数值计算对切比雪夫反变换的函数表达式[式(7)]中的无穷项级数近似的截断项数N；通常可根据声速、介质密度等物理参数在深度方向上的光滑程度进行选择，默认取100；实践中也可采用试算的办法逐步调整修正此参数的取值；

3.2)给定r＝0处的初场边界条件：使用经典的高斯声源法、Greene声源法、简正波声源、镜像声源、自启动型声源中的一种方法给定r＝0处的u⁽⁰⁾值，其中r为水平方向距离，u⁽⁰⁾为缩比分量函数u在水平方向距离r＝0处的函数取值；根据r＝0处的u⁽⁰⁾值，利用切比雪夫正变换的函数表达式，计算得到谱空间上的N+1项谱系数

其中N为无穷项级数近似的截断项数，

分别为缩比分量函数u在水平方向距离r＝0处的N+1项谱展开系数；设定水平方向距离上的离散网格点计数变量m的值为0；

3.3)获得自变量变换后的算子X中第一个算子项在谱空间中表示的矩阵D，获得自变量变换后的算子X中第二个算子项在谱空间中表示的矩阵E，将矩阵D和矩阵E相加，得到算子矩阵X；

3.4)计算Padé级数的分子多项式的第j项系数α_j和Padé级数的分母多项式的第j项系数β_j，基于算子矩阵X装配谱空间上线性方程组[式(1)]的左端系数矩阵A及右端系数矩阵B：

上式中，i为虚数单位，k₀为参考波数，Δr为r方向离散网格点的间距，L是Padé级数的阶数，α_j为Padé级数的分子多项式的第j项系数，β_j为Padé级数的分母多项式的第j项系数，I为单位矩阵，X为算子矩阵；

3.5)加入合适的水面、水底边界条件，修正左端系数矩阵A及右端系数矩阵B，从而得到最终在谱空间中离散化后的线性方程组：

上式中，A为左端系数矩阵，

为缩比分量函数u在水平方向距离为r_m+1处的谱展开系数，B为右端系数矩阵，

为缩比分量函数u在水平方向距离为r_m处的谱展开系数。

本实施例中，获得自变量变换后的算子X[式(6)]中第一个算子项(二阶导数算子)在谱空间中表示的矩阵D的详细步骤包括：

3.3.1A)取(N+1)×(N+1)阶的方阵D，将其所有元素初始化为零值；

3.3.2A)设方阵D的行下标及列下标均取从0到N的非负整数；对于其第k行的第j列元素D[k][j]而言，若k与j的奇偶性相反且j＞k，则重新置D[k][j]的值为2j，其中k＝0，1，…，N；j0，1，…N；

3.3.3A)将矩阵D中第1行的每个元素乘以0.5，所得结果原地保存于矩阵D；

3.3.4A)按照矩阵相乘的规则，计算乘积矩阵D×D，并将结果原地保存在矩阵D中；

3.3.5A)将矩阵D中每个元素乘以4/H²，所得结果仍原地保存于矩阵D，从而得到第一个算子项在谱空间中表示的矩阵D。

本实施例中，获得自变量变换后的算子X[式(6)]中第二个算子项在谱空间中表示的矩阵E的详细步骤包括：

3.3.1B)记物理空间上的函数

为v(x)，利用切比雪夫正变换[式(8)]的函数表达式求出相应的谱空间系数

k＝0，1，…，N；

3.3.2B)取(N+1)×(N+1)阶的方阵E，将其所有元素初始化为零值；

3.3.3B)设方阵E的行下标及列下标均取从0到N的非负整数；对于其第k行的第j列元素E[k][j]而言：若0≤k-j≤N，则将E[k][j]在原值基础上增加

若0≤j-k≤N，则将E[k][j]在原值基础上增加

若0≤j+k≤N，则将E[k][j]在原值基础上增加

若上述条件同时满足多个，则每次都执行相应的增加操作，其中k＝0，１，…，N；j＝0，１，…N；

3.3.4B)将矩阵E中所有元素乘以0.5，所得结果原地保存于矩阵E，从而得到第二个算子项在谱空间中表示的矩阵E；

本实施例中，步骤3.5)中修正左端系数矩阵A及右端系数矩阵B的步骤包括：

3.5.1)利用切比雪夫反变换的函数表达式[式(7)]，将边界条件

转换成谱空间的方程式

此处a_k＝(-1)^k，其中

为缩比分量函数u在水平方向距离为r_m、深度为0处的函数取值，N为无穷项级数近似的截断项数，k为缩比分量函数u的展开系数的序号，

为缩比分量函数u在水平方向距离为r_m处的第k个展开系数；

3.5.2)修正左端系数矩阵A第N-1行的第k个元素取为a_k，右端系数矩阵B第N-1行的第k个元素取值为0；

3.5.3)利用切比雪夫反变换的函数表达式，将边界条件

转换成谱空间的方程式

此处a_k＝1；

3.5.4)修正左端系数矩阵A第N-1行的第k个元素取为a_k，右端系数矩阵B第N-1行的第k个元素取值为0。

本实施例中，利用切比雪夫正变换[式(8)]时，在区间[-1，1]上取N＝200个Gauss-Lobatto节点，200个节点的位置为

j＝0，1，2，…200。

本实施例中，步骤4)的详细步骤包括：

4.1)使用标准的线性方程组求解算法解得谱空间中离散化后的线性方程组中第m+1个距离上的谱系数解

且在求解方程组的计算区域利用众核CPU和GPU进行协同加速计算；本实施例中，利用众核CPU和GPU进行协同加速计算具体是指在求解方程组的计算区域，使用OpenMP编译指导指令，对通用众核计算设备(即众核CPU)上的计算进行并行加速；使用OpenACC编译指导指令，将关键计算部分分载到异构GPU加速设备上，实现协同加速计算，充分挖掘通用众核处理器及GPU加速器异构平台的计算潜力，显著加快了计算速度，减少计算时间；

4.2)判断水平方向距离达到预定范围r_max是否成立，若成立则跳转执行步骤5)；否则，令m＝m+1，跳转执行步骤4.1)继续求解谱空间中离散化后的线性方程组。本实施例中，重复求解线性方程组

次，直到水平方向距离达到最大求解距离r_max＝100km为止。

本实施例中，步骤5)的详细步骤包括：针对每个m＝1，2，…利用切比雪夫反变换将谱空间的定义域x∈[-1，1]离散成指定分辨率的网格(默认为间距为0.001的2001个格点，实际操作中可根据需要改变默认设置)，第j个格点的坐标为x_j，利用z_j＝((1-x_j)H)/2将谱空间上的定义域变换为物理空间上的定义域(即实际深度)，得到各深度对应的声压，其中z_j为将谱空间上的坐标x_j解映射回物理空间后的坐标，H为目标海域的深度。本实施例中，使用切比雪夫反变换，将谱空间上的解

映射回物理空间u(r_m，z)。针对每个m＝1，2，…10000，利用切比雪夫反变换将谱空间的定义域x∈[-1，1]离散成间距为0.001的2001个格点。

本实施例中，步骤6)用于完成x坐标向z坐标的自变量反变换获得水平方向距离为r、深度为z位置处的缩比分量函数u的函数值u(r，z)，获得水平方向距离为r、深度为z位置处的缩比分量函数u的函数值u(r，z)后，利用式(4)即可求出声压场p。

本实施例中，步骤7)计算指定接收器位置处声场的声压值及传播损失值时，首先根据接收器的深度z_r确定所在深度上的声压值p，然后利用下式计算出接收深度上的传播损失值TL。

上式中，p₀为距离声源1米处的声压值。本实施例中步骤7)中还包括绘制图形的步骤，得到的全空间场中传播损失分布图如图2所示，在深度为1000m处的传播损失曲线图如图3所示。

综上所述，本实施例方法在水声传播计算领域引入切比雪夫谱方法，引入切比雪夫谱方法需要克服以下技术障碍或技术问题：

1、传统有限差分方法容易从数学上证明其相容性、稳定性和收敛性，但切比雪夫谱方法除对少数简单情形外，大多数情形下不易证明其相容性、稳定性和收敛性。本发明中，主要关注的是介质水平分层的海洋声学问题，它相比于其它与水平方向距离相关的物理问题来说，处理更为简单，更容易满足相容性、稳定性和收敛性的条件。

2、切比雪夫谱方法要求原始声速剖面数据尽量光滑，但传统上声速数据过于稀疏，若插值不当难以满足光滑性，从而无法保证步进迭代的过程中的稳定性。本发明通过增加声速剖面观测值密度、简化声速剖面插值方法加以解决这一问题。

3、切比雪夫谱方法比传统有限差分法更难处理边界条件。本发明中选择使用经典谱方法中被称为Tau的处理方式，将边界条件自动离散为最终方程组中的两个方程，同时配合使用Gauss-Lobatto积分节点进行谱系数的数值计算，从而解决了边界条件问题。

4、有限差分法得到的是稀疏系数方程组，而切比雪夫谱方法则得到是稠密系数方程组，带来了更大的计算量和更多的计算资源消耗，初看起来，这似乎是切比雪夫谱方法的缺点；然而，为了获得相同的计算精度，切比雪夫谱方法方程组阶数比有限差分方法要小得多。综合下来，切比雪夫谱方法反而更有优势：总的求解时间更快、解的精度更高。本发明中默认设置N＝100，实际计算中用户可以根据所研究的范围大小，海洋环境的复杂程度，解变量的精度要求，自主设置截断阶数(方程组的数量)，在需求和计算量之间做出取舍，取得平衡。

相对于传统有限差分方法，本实施例方法具有收敛速度快、计算精度高、计算时间开销小等优点，可在使用较少离散网格点的条件下获得更高的计算精度，同时能够充分发挥硬件平台的计算性能优势，显著加快数值模拟的速度。

此外，本实施例还提供一种基于切比雪夫多项式谱的水下声场数值模拟系统，包括计算机设备，该计算机设备被编程或配置以执行前述基于切比雪夫多项式谱的水下声场数值模拟方法的步骤，或者该计算机设备的存储器上存储有被编程或配置以执行前述基于切比雪夫多项式谱的水下声场数值模拟方法的计算机程序。

此外，本实施例还提供一种计算机可读存储介质，该计算机可读存储介质上存储有被编程或配置以执行前述基于切比雪夫多项式谱的水下声场数值模拟方法的计算机程序。

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

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (10)

1.一种基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于该方法包括：

1)针对具体水声环境参数，建立柱坐标系下声学传播Helmholtz方程的简化控制方程；

2)将简化控制方程完成z坐标向x坐标的自变量变换；

3)对简化控制方程实施切比雪夫正变换，将物理空间映射到谱空间上，形成谱空间中离散化后的线性方程组；

4)求解谱空间中离散化后的线性方程组，得到谱空间上的解；

5)对谱空间上的解实施切比雪夫反变换，将谱空间上的解映射回物理空间；

6)完成x坐标向z坐标的自变量反变换获得水平方向距离为r、深度为z位置处的缩比分量函数u的函数值u(r，z)，并利用该函数值u(r，z)求出声压场p；

7)计算指定接收器位置处声场的声压值及传播损失值。

2.根据权利要求1所述的基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于，步骤1)中的简化控制方程为下式所示的方程：

上式中，i为虚数单位，k₀为参考波数，Δr为r方向离散网格点的间距，L是Padé级数的阶数，α_j为Padé级数的分子多项式的第j项系数，β_j为Padé级数的分母多项式的第j项系数，X为算子，u^(m)＝u(r_m，z)为缩比分量函数u在水平方向距离为r_m、深度为z位置处的函数值，u^(m+1)＝u(r_m+1，z)为缩比分量函数u在水平方向距离为r_m+1、深度为z位置处的函数值，r_m＝m·Δr为r方向第m个离散网格点的坐标，r_m+1＝(m+1)·Δr为r方向第m+1个离散网格点的坐标，m＝0，1，2，…；其中算子X的函数表达式如下式所示：

上式中，n是介质折射率，算子X的函数表达式由第一个算子项、第二个算子项两者之和构成。

3.根据权利要求2所述的基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于，步骤3)的详细步骤包括：

3.1)确定数值计算对切比雪夫反变换的函数表达式中的无穷项级数近似的截断项数N；

3.2)给定r＝0处的初场边界条件：使用经典的高斯声源法、Greene声源法、简正波声源、镜像声源、自启动型声源中的一种方法给定r＝0处的u⁽⁰⁾值，其中r为水平方向距离，u⁽⁰⁾为水平方向距离r＝0时缩比分量函数u的函数取值；根据r＝0处的u⁽⁰⁾值，利用切比雪夫正变换的函数表达式，计算得到谱空间上的N+1项谱系数

其中N为无穷项级数近似的截断项数，

分别为N+1项的谱展开系数；设定水平方向距离上的离散网格点计数变量m的值为0；

3.3)获得自变量变换后的算子X中第一个算子项在谱空间中表示的矩阵D，获得自变量变换后的算子X中第二个算子项在谱空间中表示的矩阵E，将矩阵D和矩阵E相加，得到算子矩阵X；

3.4)计算Padé级数的分子多项式的第j项系数α_j和Padé级数的分母多项式的第j项系数β_j，基于算子矩阵X装配谱空间上线性方程组的左端系数矩阵A及右端系数矩阵B：

上式中，i为虚数单位，k₀为参考波数，Δr为r方向离散网格点的间距，L是Padé级数的阶数，α_j为Padé级数的分子多项式的第j项系数，β_j为Padé级数的分母多项式的第j项系数，I为单位矩阵，X为算子矩阵；

3.5)加入合适的水面、水底边界条件，修正左端系数矩阵A及右端系数矩阵B，从而得到最终在谱空间中离散化后的线性方程组：

上式中，A为左端系数矩阵，

为缩比分量函数u在水平方向距离为r_m+1处的谱展开系数，B为右端系数矩阵，

为缩比分量函数u在水平方向距离为r_m处的谱展开系数。

4.根据权利要求3所述的基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于，所述获得自变量变换后的算子X中第一个算子项在谱空间中表示的矩阵D的详细步骤包括：

3.3.1A)取(N+1)×(N+1)阶的方阵D，将其所有元素初始化为零值；

3.3.2A)设方阵D的行下标及列下标均取从0到N的非负整数；对于其第k行的第i列元素D[k][j]而言，若k与j的奇偶性相反且j＞k，则重新置D[k][j]的值为2j，其中k＝0，1，…，N；j＝0，1，…N；

3.3.3A)将矩阵D中第1行的每个元素乘以0.5，所得结果原地保存于矩阵D；

3.3.4A)按照矩阵相乘的规则，计算乘积矩阵D×D，并将结果原地保存在矩阵D中；

3.3.5A)将矩阵D中每个元素乘以4/H²，所得结果仍原地保存于矩阵D，从而得到第一个算子项在谱空间中表示的矩阵D。

5.根据权利要求3所述的基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于，所述获得自变量变换后的算子X中第二个算子项在谱空间中表示的矩阵E的详细步骤包括：

3.3.1B)记物理空间上的函数

为v(x)，利用切比雪夫正变换的函数表达式求出相应的谱空间系数

3.3.2B)取(N+1)×(N+1)阶的方阵E，将其所有元素初始化为零值；

3.3.3B)设方阵E的行下标及列下标均取从0到N的非负整数；对于其第k行的第j列元素E[k][j]而言：若0≤k-j≤N，则将E[k][j]在原值基础上增加

若0≤j-k≤N，则将E[k][j]在原值基础上增加

若0≤j+k≤N，则将E[k][j]在原值基础上增加

若上述条件同时满足多个，则每次都执行相应的增加操作，其中k＝0，1，…，N；j＝0，1，…N；

3.3.4B)将矩阵E中所有元素乘以0.5，所得结果原地保存于矩阵E，从而得到第二个算子项在谱空间中表示的矩阵E。

6.根据权利要求3所述的基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于，步骤3.5)中修正左端系数矩阵A及右端系数矩阵B的步骤包括：

3.5.1)利用切比雪夫反变换的函数表达式，将边界条件

转换成谱空间的方程式

此处a_k＝(-1)^k，其中

为缩比分量函数u在水平方向距离为r_m、深度为0处的取值，N为无穷项级数近似的截断项数，k为缩比分量函数u的展开系数的序号，

为缩比分量函数u在水平方向距离为r_m处的第k个展开系数；

3.5.2)修正左端系数矩阵A第N-1行的第k个元素取为a_k，右端系数矩阵B第N-1行的第k个元素取值为0；

3.5.3)利用切比雪夫反变换的函数表达式，将边界条件

转换成谱空间的方程式

此处a_k＝1；

3.5.4)修正左端系数矩阵A第N-1行的第k个元素取为a_k，右端系数矩阵B第N-1行的第k个元素取值为0。

7.根据权利要求3所述的基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于，步骤4)的详细步骤包括：

4.1)使用标准的线性方程组求解算法解得谱空间中离散化后的线性方程组中第m+1个距离上的谱系数解

且在求解方程组的计算区域利用众核CPU和GPU进行协同加速计算；

4.2)判断水平方向距离达到预定范围r_max是否成立，若成立则跳转执行步骤5)；否则，令m＝m+1，跳转执行步骤4.1)继续求解谱空间中离散化后的线性方程组。

8.根据权利要求3所述的基于切比雪夫多项式谱的水下声场数值模拟方法，其特征在于，步骤5)的详细步骤包括：针对每个m＝1，2，…利用切比雪夫反变换将谱空间的定义域x∈[-1，1]离散成指定分辨率的网格，第j个格点的坐标为x_j，利用z_j＝((1-x_j)H)/2将谱空间上的定义域变换为物理空间上的定义域，得到各深度对应的声压，其中z_j为将谱空间上的坐标x_j解映射回物理空间后的坐标，H为目标海域的深度。

9.一种基于切比雪夫多项式谱的水下声场数值模拟系统，包括计算机设备，其特征在于，该计算机设备被编程或配置以执行权利要求1～8中任意一项所述基于切比雪夫多项式谱的水下声场数值模拟方法的步骤，或者该计算机设备的存储器上存储有被编程或配置以执行权利要求1～8中任意一项所述基于切比雪夫多项式谱的水下声场数值模拟方法的计算机程序。

10.一种计算机可读存储介质，其特征在于，该计算机可读存储介质上存储有被编程或配置以执行权利要求1～8中任意一项所述基于切比雪夫多项式谱的水下声场数值模拟方法的计算机程序。

Images