一种可稳定计算负跃层浅海环境下声场的波数积分方法

于晓林，许伟杰，杨春梅，骆文于

(1.中国科学院声学研究所东海研究站，上海 201800；2.自然资源部第一海洋研究所，山东青岛 266100；3.中国科学院声学研究所声场声信息国家重点实验室，北京 100190)

0 引言

在夏、秋季节，浅海中常存在较强的负跃层。在负跃层海域，海水可以分成三层，最上层以及最下层的声速随深度变化缓慢，两层之间的中间层的声速变化剧烈。浅海负跃层对水下声监测、渔业生产、舰艇航行等有极大的影响[1]，因此对浅海负跃层声传播规律的研究极为重要[2]。

波数积分方法首先是由Pekeris引入到水声学中的[3]，简正波[4-5]方法是通过复围线积分法计算该积分，得到的是不同阶简正波对应的留数之和，而波数积分方法是通过直接对波数进行数值积分求解的，国内外许多学者对波数积分方法进行了相应的研究[6-8]。在某些低频浅海环境下，波数积分算法的计算结果比简正波方法更准确[9]。

文献[10]提出了一个负跃层浅海环境中的波数积分解法，该解法在理论上是正确的，但是在实际的仿真中发现，当水平波数的取值过大时，会造成艾里(Airy)函数的数值溢出，系数矩阵存在数值不稳定的现象。为了解决文献[10]中给出的解法中存在的数值不稳定的问题，本文通过合理地归一化，提出了一种负跃层浅海环境下无条件稳定的波数积分方法。

1 理论推导

1.1 波数积分理论

1.2 负跃层浅海环境下不稳定的波数积分算法

文献[10]给出了计算负跃层浅海环境下的深度格林函数的线性方程组，该线性方程组在理论推导上是没问题的，但是由于未对Airy函数做任何的变换，因此在实际的仿真中，系数矩阵存在数值不稳定现象。下面简单介绍下文献[10]中给出的线性方程组并分析其不稳定的原因，然后通过对Airy函数合理的变换，从而求得稳定的深度格林函数。

负跃层浅海环境中的声速如式(5)所示。

其中：cu为第一层中的声速，cd为第三层中的声速，cb为海底声速，ρw为海水密度，ρb为海底密度，H为海水深度，H1、H2分别为负跃层的起、止深度。

以声源位于第一层为例，各层中的线性方程可以表示如下[11]：

该矩阵方程的系数矩阵里存在趋于无穷大的项Bi(ζ)以及Bi'(ζ)，因此可能存在数值溢出的现象。在仿真中发现，当水平波数取值大于第三层海水波数k3时，仿真计算会提示矩阵接近奇异值、结果可能不准确的警告。可见，由于Airy函数的渐近特性，该深度格林函数的解在数值上是不稳定的，为了解决该问题，需要在系数项里避免趋于无穷大的项Bi(ζ)以及Bi'(ζ)，因此本文提出了一种无条件稳定的波数积分算法。

1.3 负跃层浅海环境下稳定的波数积分算法

由于Airy函数的渐近特性，上述得到的矩阵方程存在数值不稳定的问题，下面给出一种数值稳定的计算方法。负跃层浅海环境中，声源位于第一层、第三层、第二层时的示意图如图1～3所示。

图1 声源处于第一层Fig.1 The sound source in the first layer

1.3.1 声源处于第一层时的矩阵方程

为了避免系数矩阵的系数项是奇异矩阵，需要对Airy函数进行变换，将Airy函数表示成式(22)中的形式[12]：

图2 声源处于第三层Fig.2 The sound source in the third layer

图3 声源处于第二层Fig.3 The sound source in the second layer

该矩阵方程的系数项不存在趋于无穷大的项Bi(ζ)以及Bi'(ζ)，是无条件稳定的，可以求得未知变量，然后通过式(6)～(9)求得声源在第一层时，各个层的深度格林函数，最后利用式(3)可以求得声场的位移势能，从而求得各个层的声场。

1.3.2 声源处于第三层时的矩阵方程

当声源处于第三层时，同样可以得到声源位于第三层时的无条件稳定的矩阵方程：

此时：A、x与式(30)中取值相同，b的取值如下：

通过该矩阵方程，可以求得7个未知变量，进而求得总声场。

1.3.3 声源处于第二层时的矩阵方程

当声源处于第二层时，需要在声源深度处构建一个虚拟界面，同时在该虚拟边界z=zs处添加两个边界条件：

从而得到声源位于第二层的矩阵方程：

此时：

利用该矩阵方程，可以求出9个未知变量，进而求得总声场。

1.4 声场计算

通过前面的推导，可以得到声源位于三个层中的无条件稳定的深度格林函数，再利用式(3)可以得到声场的位移势。在实际的仿真中，对格林函数尖峰的欠采样，容易引入较大的误差。为了解决这一问题，采用复围线积分的方法[11]：

式中：M是采样点数，kmax为数值积分的上限，kmin为数值积分的下限。

2 数值算例

考虑图4的海洋环境，海洋环境为负跃层浅海环境，海水深度为100 m，海水密度为1.0 g·cm-3，0～25 m的声速为1 500 m·s-1，25～45 m的声速由1 500 m·s-1变为1 460 m·s-1，45～100 m的声速为1 460 m·s-1，海底声速为1 800 m·s-1，海底密度为1.8 g·cm-3，海底吸收系数为0.15 dB·λ-1。

图4 负跃层浅海声场计算实例的海洋声速剖面示意图Fig.4 A typical sound speed profile for sound field calculation in the shallow-water with a negative thermoclines

分别将声源置于第一层等声速层以及第二层负跃层中，比较数值不稳定的算法与本文所提出的数值稳定的解法的结果差异：对于数值不稳定算法来说，水平波数最大取值只能取到第三层海水中波数k3附近，其中k3=ω/cd，当水平波数的取值超过k3值时，系数矩阵存在数值不稳定的现象，线性方程求解器在求解时会发生预警；而对于优化后的算法来说，水平波数可以取到一个极大值。在实际的仿真计算中，水平波数的取值应该根据算例的具体情况来确定：理论上，为保证格林函数的收敛性，水平波数的取值应该越大越好，但是为保证计算效率，水平波数取值又应该越小越好。对于本算例来说，1.2k3是一个可以保证格林函数收敛到0值附近的极小取值，因此本文中优化后算法的水平波数一律取到1.2k3。图5是声源位于20 m(声源在第一层)时，各层的深度格林函数。在本算例中，数值不稳定算法的深度格林函数的水平波数最大值只能取到k3(即0.215)，优化后的算法的水平波数最大取值为1.2k3(即0.258)时，就可以保证该算例下深度格林函数的收敛性。图6给出相应的传播损失结果。

图6 声源深度为20 m，不同接收深度处的传播损失计算结果（稳定与不稳定波数积分算法对比）Fig.6 The transmission losses at different depths for the sound source at the depth of 20 m(comparison between stable and unstable wave-number integral algorithms)

图7是声源位于40 m(声源在第二层)时，不同深度的深度格林函数，可以看到，当接收深度为40 m时，数值不稳定算法的深度格林函数不能收敛到0附近，而优化后的算法可以保证深度格林函数的足够收敛性。图8给出了相应的传播损失结果。

图7 声源深度为40 m，不同接收深度处的深度格林函数计算结果（稳定与不稳定波数积分算法对比）Fig.7 The magnitudes of the Green’s function at different depths for the sound source at the depth of 40 m(comparison between stable and unstable wave-number integral algorithms)

图8 声源深度为40 m，不同接收深度处的传播损失的计算结果(稳定与不稳定波数积分算法对比)Fig.8 The transmission losses at different depths for the sound source at the depth of 40 m(comparison between stable and unstable wave-number integral algorithms)

为了进一步验证本文所提方法的准确性，比较本文所提方法的计算结果与简正波模型KRAKENC的计算结果[14]，图9给出声源深度在20 m(声源位于第一层)时各个接收深度的传播损失结果，图10给出声源深度在40 m(声源位于第二层)时的各个接收深度的传播损失结果。结果显示，两种方法的计算结果吻合得较好。由于简正波方法在某些低频浅海环境中的近场计算上存在误差[9]，因此本文所提方法可以作为负跃层浅海环境中声场计算的有效补充。

图9 声源深度为20 m，不同接收深度处的传播损失的计算结果(本文算法与KRAKENC算法对比)Fig.9 The transmission losses at different depths for the sound source at the depth of 20 m(comparison between the proposed algorithm and KRAKENC algorithm)

图10 声源深度为40 m，不同接收深度处的传播损失的计算结果(本文算法与KRAKENC算法对比)Fig.10 The transmission losses at different depths for the sound source at the depth of 40 m(comparison between the proposed algorithm and KRAKENC algorithm)

3 结论

本文通过对Airy函数的处理，得到了负跃层浅海环境下无条件稳定的波数积分方法。该方法保持了系数矩阵方程中的数值稳定性，避免了由于数值不稳定造成的深度格林函数的收敛性不足的问题，能够提高负跃层浅海环境下声场计算的稳定性。

将本文所提方法与简正波模型KRAKENC的计算结果进行比较，结果显示，两种方法有很好的一致性。本文所提的方法在数值上是稳定的，并且可以直接通过稳定的系数矩阵方程求得相应的解析解。因此，本文所提的方法可以作为负跃层浅海环境中声场计算的有效补充。