声学工具箱学习随着水声学的发展，计算水下声场越来越重要，很多人都开启了学习如何使用声学工具箱(Acoustics Too

前言

随着水声学的发展，计算水下声场越来越重要，很多人都开启了学习如何使用声学工具箱(Acoustics Toolbox)计算水下声场，本文就介绍了水下声场的计算和声学工具箱使用基础内容。

一、波动方程原理和使用场景

波动方程的推导 $\begin{cases}\nabla^2p-\frac{1}{c^2}\frac{\partial^2p}{\partial t^2}=0\\\frac{\partial\rho^{\prime}}{\partial t}=-\rho_0\nabla\cdot\boldsymbol{v}\\p=\rho^{\prime}c^2\end{cases}\Rightarrow\nabla^2p-\frac{1}{c^2}\frac{\partial^2p}{\partial t^2}=0$ 波动方程由运动方程、连续性方程和物态方程推导出，其中p为声压，也被称为余量压强，c为声速，t为时间，ρ'为介质密度，ρ_0为介质静态密度，v为质点振速。波动方程是一个非常难解的问题，从其表达式上看，它属于二阶多元偏微分方程，且在非齐次和边界条件复杂的条件下的理论求解更加困难，在大多数情况下，难以找到该方程的解析解。科学家们将物理图像应用一个声源在水中发射声波的情况，使用了数值求解的方法，对水下声场进行求解，目前比较热门的数值求解方法为简正波模型（KRAKEN）、射线模型（BELLHOP）和抛物方程模型（RAM）三种。三种方法各有利弊，适用于不同环境条件下的声场计算中。

(1) 简正波原理： $p(r,z)=\sum_{m=1}^\infty\Phi_{m}(r)\Psi_{m}(z)$ 其中， $p(r,z)$ 为 $(r,z)$ 位置上接收器的声压， $\Phi_m(r)$ 为汉克尔函数， $\Psi_{m}(z)$ 为简正波本征函数。将该式带入存在声源的波动方程中。通过求解本征值 $k_{rm}$ 和本征函数 $\Psi_{m}(z)$ 得到声场的解。得到解为： $p(r,z)=\frac{i}{\rho(z_{s})\sqrt{8\pi r}}e^{-i\pi/4}\sum_{m=1}^{\infty}\Psi_{m}(z_{s})\Psi_{m}(z)\frac{e^{ik_{rm}r}}{\sqrt{k_{rm}}}$ 注：详细推导过程见《Computational Ocean Acoustic》Chapter. 5

(2) 射线原理：基本思想:声波能量沿一定路径通过一定形状的几何面向外传播，确定射线空间走向及其强度大小即可完全确定声传播规律 $p(x,y,z,t)=A(x,y,z)e^{i[\omega t-k(x,y,x)\Phi_1(x,y,z)]}$ 代入波动方程，得到程函方程和强度方程 $\frac{\nabla^2A}A-k_0^2\nabla\Phi\cdot\nabla\Phi+k^2=0$ $\nabla^2\Phi+\frac2A\nabla A\cdot\nabla\Phi=0$ 将程函方程进行高频近似，强度方程略去小量，求解得到声场。注：详细推导过程见《Computational Ocean Acoustic》Chapter. 3

(3) 抛物方程原理：抛物方程法通过对三维波动方程简化求解，得到抛物方程： $2ik_0\frac{\partial\varphi}{\partial r}+\frac{\partial^2\varphi}{\partial z^2}+k_0^2(n^2-1)\phi=0$ 对上式简化求解得到递推关系： $p(r,z)=p(r_0,z)\sqrt{\frac{r_0}{r}}\exp[ik_0\left(\frac{k_0^2+3k_{rm}^2}{3k_0^2+k_{rm}^2}\right)(r-r_0)]$ $p(x,y,z,t)=A(x,y,z)e^{i[\omega t-k(x,y,x)\Phi_1(x,y,z)]}$ 上方程为Claerbout抛物方程的解。在实际的算法中，可以利用简正波法求解 $p(r_0,z)$ ，再利用上式递推求解声场。注：详细推导过程见《Computational Ocean Acoustic》Chapter. 6 不同场景下的模型适用度

在一些书中给出了更宽泛的模型和应用场景，但是我是做应用的，我其实推荐在大多数场景下使用射线模型，速度快，宽容性高，可解释性好，能够计算时域到达结构。

二、使用步骤

1. 下载Acoustics Toolbox环境

可以使用声学工具箱软件（Acoustics Toolbox）搭配Matlab来快速使用这些模型，在 Acoustics Toolbox 网站上可以下载这个at文件，一般下载 Binaries文件。推荐下载我整编的工具箱，包含了ram_p等修复模型和函数，at_qp at官方下载

将zip下载到本地后，解压到任意路径，在matlab中打开这个文件夹，注意当前路径的最后一层一定要是atWin10_2020_11_4，这样使用pwd命令才能找到需要依赖的文件夹。打开at_init_matlab.m，新增蓝色框代码，添加.exe路径，并保存路径，下次就可以直接使用工具箱，而不需要重新设置路径。

代码设置路径

代码如下（示例）：

addpath( fullfile(Home, 'windows-bin-20201102' ) );             % AT binaries
Matdir = fullfile( Home, 'Matlab' );          % AT Matlab routines
addpath( Matdir );
savepath;

或者您也可以在Matlab主页中设置路径，添加并包含文件夹选择atWin10_2020_11_4文件夹。 Matlab设置路径添加并包含子文件夹选择at文件夹

2. KRAKEN使用说明

kraken.m的输入和输出

.env和.flp为KRAKEN的环境文件； .prt、.mod和.shd为KRAKEN生成文件；
.prt为KRAKEN生成的注释文件，用来对环境文件的读取和KRAKEN计算过程进行说明，是文本文件可以直接阅读；
.mod文件储存KRAKEN计算出的特征值和特征函数，是二进制文件，无法直接阅读，可以用read_modes函数阅读；
.shd文件为KRAKEN计算出声压数据，是二进制文件，无法直接阅读，可以用read_shd函数读取声压数据，用plotshd画传播损失图。

pekerisK.env（浅海示例）：

'pekeris'		! 标题
500			! 频率
1			! NMEDIA，介质层数
'NVW'		! TOP OPTION，海面选项。第一个字母代表声速剖面的插值方法，第二个字母代表上半边界条件，第三个字母代表衰减系数单位。
0  0.0  100.0		! 海水竖直网格层个数，界面粗糙度，海水深度
    0.0  1500.00  0.0  1.0  0.0  0.0 	! 深度，纵波速度，横波速度，密度，纵波衰减系数、横波衰减系数
    100.0 1500.00  /		! 深度，纵波速度，'/'代表后面几项与上面保持一致
'A' 0.0              	! BOTTOM OPTION，海底选项，A：海底半空间条件，0.0海底粗糙度
    100.0  1800.00 0.0 1.8 0.0 0.0 / ! 海底深度，纵波速度度，横波速度，密度，纵波衰减，横波衰减
1400 1800		! CLOW，CHIGH，相速度计算的范围
0			! RMAX (km) 
1			! NSD 规定声源个数
25.0/			! SD(1:NSD) (m) 规定声源深度
101			! NRD  规定接收器深度取值个数
0 100/		! RD(1:NRD) (m) 规定接收器深度

pekerisK.flp（浅海示例）：

'pekeris'		! 标题
'RA'			! SOURCE OPTION 声源选项，第一个字母中'R'代表点源，'X'代表线源，第二个字母代表模式理论的选择，'C'代表耦合模式理论，'A'代表绝热模式理论。
9999			! M  (number of modes to include) 模式个数
1			! 声速剖面个数
     0.0 100.0 /		! RPROF(1:NPROF) (km)  声速剖面范围
300			! NR 接收器水平个数
0.1 30 /		! RMIN,   RMAX (km) 接收器水平距离范围
1			! NSD 声源个数
   25 /			! SD(1:NSD) 声源深度范围
101			! NRD 接收器竖直个数
  0 100  /		! RD(1:NRD) 接收器深度范围
101			! NRD 接收器竖直个数
    0  /			! RR(1:NRR)   (m) 接收器在水平方向的偏移（倾斜阵）

运行代码

% 运行kraken
kraken pekerisK
kraken('pekerisK') 
% 画声速剖面
plotssp('pekerisK', 'KRAKEN')
% 画简正波模态函数
plotmode('pekerisK.mod', 0)
% 读取声场
[ PlotTitle, PlotType, freqVec, freq0, atten, Pos, pressure ] = read_shd( 'pekerisK.shd' );
% 画传播损失
figure; plotshd pekerisK.shd

运行结果为 kraken计算声场结果

3. BELLHOP使用说明

bellhop.m的输入和输出

.env为BELLHOP的环境文件，将top option和bottom option 加入~读取海底地形.bty、海面形状. ati文件，可以将top option和bottom option 改为F读取海面反射系数.trc和海底反射系数.brc文件，可以将runtype改为*读取声源指向性.sbp文件；
.prt、.ray和.shd等为BELLHOP生成文件；
.prt为BELLHOP生成的注释文件，用来对环境文件的读取和BELLHOP计算过程进行说明，是文本文件可以直接阅读；
.ray为声线文件，可以使用plotray函数画声线图；
.shd文件为BELLHOP计算出声压数据，是二进制文件，无法直接阅读，可以用read_shd函数读取声压数据，利用plotshd画传播损失图。

pekerisB.env（浅海示例）：

'pekeris'		! 标题
500			! 频率
1			! NMEDIA，介质层数
'NVW'		! TOP OPTION选项。第一个字母代表声速剖面的插值方法，第二个字母代			表上半边界条件，第三个字母代表衰减系数单位。
51  0.0  100.0		! 海水分层个数，界面粗糙度，海水深度
    0.0  1500.00  0.0  1.0  0.0  0.0 	! 深度，纵波速度度，横波速度，密度，纵波衰减，横波衰减
    100.0 1500.00  /		! 深度，纵波速，'/'代表后面几项与上面保持一致
'A' 0.0              		! BOTTOM OPTION选项（A：海底半空间条件） 海底粗糙度
    100.0  1800.00 0.0 1.8 0.0 0.0 / ! 海底深度，纵波速度度，横波速度，密度，纵波衰减，横波衰减
1			! NSD 规定声源个数
25.0/			! SD(1:NSD) (m) 规定声源深度
101			! NRD  规定接收器深度取值个数
0 100.0/			! RD(1:NRD) (m) 规定接收器深度
300			! NR 规定接收器水平范围取值个数
0.1  30.0 /		! R(1:NR ) (km)规定接收器水平范围取值
'CG'	  		! Runtype（运行方式）'R/C/I/S/A'
0			! NBeams 规定绘制的声线数量
-90.0 90.0 /		! ALPHA1,2 (degrees) 声源声线出射掠射角范围
0.0  150.0  31.0		! STEP (m), ZBOX (m), RBOX (km)规定步长和声场计算范围（ZBOX需要大于最大海深，RBOX需要大于接收器最大水平距离）

运行代码

% 运行bellhop
bellhop pekerisB
% 画声速剖面
figure
plotssp('pekerisB', 'BELLHOP')
% 读取声场
[ PlotTitle, PlotType, freqVec, atten, Pos, pressure ] = read_shd( 'pekerisB.shd' );
% 画传播损失
figure; plotshd pekerisB.shd
% % 读取.arr文件，得到声场到达结构
% [ Arr, Pos ] = read_arrivals_asc('pekerisB_Arr')
% % 画声线传播轨迹
% plotray('pekerisB_ray')

运行结果为 bellhop计算声场结果

4. RAM使用说明

ram_p.m的输入和输出

.in为RAM的环境文件；
.line、.grid为RAM生成文件；
.line为RAM计算出的声源深度的传播损失，为文本文件；
.grid为RAM计算出的声压数据，为二进制文件。

pekerisR.env（浅海示例）：

pekeris
500.0 40.0 30.0           ! freq zs zr 频率，声源深度、.line文件
30000.0 5.0 1          ! rmax dr ndr 最大距离，距离计算间隔2-5dz，.grid文件的距离采样步长（ndr=1代表全采样）
100.0 0.5 2 100.0       ! zmax dz ndz zmplt 最大深度，深度计算间隔波长/5~12，.grid文件的深度采样步长，输出到.grid文件的最大深度
1500.0 8 1 0.0           ! c0 np ns rs 起始位置声速，Pade展开数，稳定限制数，稳定限制的最大距离
    0.0 100                ! rb zb 海底距离和深度(m)
30000.0 100.0
-1 -1
  0 1500              ! z cw 海水深度和声速
100 1500
-1 -1
0.0 1800              ! z cb 海底深度和声速
-1 -1
0.0 1.8               ! z rhob 海底深度和密度
-1 -1
   0.0  0.0
 900.0  0.0           ! z attn 海底深度和衰减系数
1000.0 10.0
-1 -1
30000.0                !  rp 第二个距离剖面环境参数
   0 1500              z cw
3000 1500
-1 -1
0.0 1700               z cb
-1 -1
0.0 1.5                  z rhob
-1 -1
   0.0  0.5
 900.0  0.5                z attn
1000.0 10.0
-1 -1

运行代码

% 使用ram模型计算声场
ram_p pekerisR

% 画声场传播损失
figure; plotgrid pekerisR.grid

% 画出在接收深度的传播损失
figure; plotline pekerisR.line

% 读取声压等参数
[ PlotTitle,  freqVec, Pos, pressure ]  = read_grid('pekerisR.grid');

运行结果为在这里插入图片描述

5. scooter使用说明

scooter的输入输出与kraken比较类似，但flp与kraken不同 pekerisS.env

'Pekeris problem'
50.0
1
'NVW'
0    0.0  100.0
     0.0  1500.0 0.0  1.0  0.0  0.0
   100.0  1500.0 /
'A'  0.0
   100.0  1800.00 0.0 1.8 0.0 0.0 /
1400.0  1800.0
50.0                    ! RMAX (km)
1                        ! NSD
25.0 /                  ! SD(1:NSD) (m)
101                        ! NRD
0 100 /                 ! RD(1:NRD) (m)

pekerisS.flp

'RP'                        ! 'R/X (coord), Lin/DB, Pos/Neg/Both'
300
0.1  30.0  /            ! RMIN, RMAX, NR

OPTIONS 语法及各参数含义概述这段文字主要是在介绍名为 “OPTIONS”（选项）相关的内容，包括其语法格式以及不同参数所代表的具体含义，用于设定特定程序或计算过程中的相关配置选项。 OPT (1:1) - 坐标相关语法：使用 “OPT”，通过对其第一个元素（OPT (1:1)）进行设置来确定坐标相关的选项。具体选项及描述： 'R'：表示采用圆柱坐标（R - Z），并且在计算时会运用远场近似的汉克尔函数来处理。 'X'：代表使用笛卡尔坐标（X - Z）来进行相关操作。 'H'：意味着运用汉克尔函数来提供精确的变换，不过这样做的速度会稍慢一点，但结果更加准确。需要注意的是，此选项仅在 Matlab 的 “fieldsco.m” 中可用。 OPT (2:2) - 频谱相关语法：同样是通过 “OPT” 的第二个元素（OPT (2:2)）来配置频谱方面的选项。具体选项及描述： 'P'：选择正频谱（这是推荐使用的方式）。 'N'：表示采用负频谱。 'B'：意味着同时包含正频谱和负频谱。从理论上来说，频谱积分应该沿着整个实 k 轴来进行，但实际上负频谱部分仅在近场中贡献比较显著。如果忽略负频谱部分的计算，运行时间会相应减少。 OPT (3:3) - 插值类型相关语法：通过 “OPT” 的第三个元素（OPT (3:3)）来指定插值类型。具体选项及描述： 'O'：代表多项式插值（适用于宽带运行情况，并且是默认的插值方式）。 'A'：指的是帕德（Pade）插值，这种插值方式对于连续波（CW）运行情况可能会产生更好的结果，但可靠性相对较低。同时提到帕德选项（OPT (3:3)）的鲁棒性较差，如果出现下溢或上溢情况，建议使用多项式插值。 OPT (4:4) - 源波束模式相关语法：利用 “OPT” 的第四个元素（OPT (4:4)）来选定源波束模式。具体选项及描述： '*'：表示读取一个源波束模式文件来获取相应的源波束模式信息。 'O'：是全向性（这是默认的源波束模式选项）。

运行结果 scooter结果

kraken和Scooter的模型对比结果如下：

总结

提示：以上就是今天要讲的内容，本文简单介绍了水下声场的简正波、射线、抛物方程模型原理和声学工具箱中KRAKEN、BELLHOP、RAM的使用。