声学工具箱学习
提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档
文章目录
- 前言
- 一、波动方程原理和使用场景
- 二、使用步骤
- 1. 下载Acoustics Toolbox环境
- 2. KRAKEN使用说明
- 3. BELLHOP使用说明
- 4. RAM使用说明
- 5. scooter使用说明
- 总结
前言
提示：这里可以添加本文要记录的大概内容：
随着水声学的发展，计算水下声场越来越重要，很多人都开启了学习如何使用声学工具箱(Acoustics Toolbox)计算水下声场，本文就介绍了水下声场的计算和声学工具箱使用基础内容。
提示：以下是本篇文章正文内容，下面案例可供参考
一、波动方程原理和使用场景
波动方程的推导
{
∇
2
p
−
1
c
2
∂
2
p
∂
t
2
=
0
∂
ρ
′
∂
t
=
−
ρ
0
∇
⋅
v
p
=
ρ
′
c
2
⇒
∇
2
p
−
1
c
2
∂
2
p
∂
t
2
=
0
\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
⎩⎪⎨⎪⎧∇2p−c21∂t2∂2p=0∂t∂ρ′=−ρ0∇⋅vp=ρ′c2⇒∇2p−c21∂t2∂2p=0
波动方程由运动方程、连续性方程和物态方程推导出，其中p为声压，也被称为余量压强，c为声速，t为时间，ρ’为介质密度，ρ_0为介质静态密度，v为质点振速。
波动方程是一个非常难解的问题，从其表达式上看，它属于二阶多元偏微分方程，且在非齐次和边界条件复杂的条件下的理论求解更加困难，在大多数情况下，难以找到该方程的解析解。
科学家们将物理图像应用一个声源在水中发射声波的情况，使用了数值求解的方法，对水下声场进行求解，目前比较热门的数值求解方法为简正波模型（KRAKEN）、射线模型（BELLHOP）和抛物方程模型（RAM）三种。三种方法各有利弊，适用于不同环境条件下的声场计算中。
(1) 简正波原理：
p
(
r
,
z
)
=
∑
m
=
1
∞
Φ
m
(
r
)
Ψ
m
(
z
)
p(r,z)=\sum_{m=1}^\infty\Phi_{m}(r)\Psi_{m}(z)
p(r,z)=∑m=1∞Φm(r)Ψm(z)
其中，
p
(
r
,
z
)
p(r,z)
p(r,z)为
(
r
,
z
)
(r,z)
(r,z)位置上接收器的声压，
Φ
m
(
r
)
\Phi_m(r)
Φm(r)为汉克尔函数，
Ψ
m
(
z
)
\Psi_{m}(z)
Ψm(z)为简正波本征函数。
将该式带入存在声源的波动方程中。通过求解本征值
k
r
m
k_{rm}
krm和本征函数
Ψ
m
(
z
)
\Psi_{m}(z)
Ψm(z)得到声场的解。
得到解为：
p
(
r
,
z
)
=
i
ρ
(
z
s
)
8
π
r
e
−
i
π
/
4
∑
m
=
1
∞
Ψ
m
(
z
s
)
Ψ
m
(
z
)
e
i
k
r
m
r
k
r
m
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}}}
p(r,z)=ρ(zs)8πrie−iπ/4∑m=1∞Ψm(zs)Ψm(z)krmeikrmr
注：详细推导过程见《Computational Ocean Acoustic》Chapter. 5
(2) 射线原理：
基本思想:声波能量沿一定路径通过一定形状的几何面向外传播，确定射线空间走向及其强度大小即可完全确定声传播规律
p
(
x
,
y
,
z
,
t
)
=
A
(
x
,
y
,
z
)
e
i
[
ω
t
−
k
(
x
,
y
,
x
)
Φ
1
(
x
,
y
,
z
)
]
p(x,y,z,t)=A(x,y,z)e^{i[\omega t-k(x,y,x)\Phi_1(x,y,z)]}
p(x,y,z,t)=A(x,y,z)ei[ωt−k(x,y,x)Φ1(x,y,z)]
代入波动方程，得到程函方程和强度方程
∇
2
A
A
−
k
0
2
∇
Φ
⋅
∇
Φ
+
k
2
=
0
\frac{\nabla^2A}A-k_0^2\nabla\Phi\cdot\nabla\Phi+k^2=0
A∇2A−k02∇Φ⋅∇Φ+k2=0
∇
2
Φ
+
2
A
∇
A
⋅
∇
Φ
=
0
\nabla^2\Phi+\frac2A\nabla A\cdot\nabla\Phi=0
∇2Φ+A2∇A⋅∇Φ=0
将程函方程进行高频近似，强度方程略去小量，求解得到声场。
注：详细推导过程见《Computational Ocean Acoustic》Chapter. 3
(3) 抛物方程原理：
抛物方程法通过对三维波动方程简化求解，得到抛物方程：
2
i
k
0
∂
φ
∂
r
+
∂
2
φ
∂
z
2
+
k
0
2
(
n
2
−
1
)
ϕ
=
0
2ik_0\frac{\partial\varphi}{\partial r}+\frac{\partial^2\varphi}{\partial z^2}+k_0^2(n^2-1)\phi=0
2ik0∂r∂φ+∂z2∂2φ+k02(n2−1)ϕ=0
对上式简化求解得到递推关系：
p
(
r
,
z
)
=
p
(
r
0
,
z
)
r
0
r
exp
[
i
k
0
(
k
0
2
+
3
k
r
m
2
3
k
0
2
+
k
r
m
2
)
(
r
−
r
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(r,z)=p(r0,z)rr0exp[ik0(3k02+krm2k02+3krm2)(r−r0)]
p
(
x
,
y
,
z
,
t
)
=
A
(
x
,
y
,
z
)
e
i
[
ω
t
−
k
(
x
,
y
,
x
)
Φ
1
(
x
,
y
,
z
)
]
p(x,y,z,t)=A(x,y,z)e^{i[\omega t-k(x,y,x)\Phi_1(x,y,z)]}
p(x,y,z,t)=A(x,y,z)ei[ωt−k(x,y,x)Φ1(x,y,z)]
上方程为Claerbout抛物方程的解。在实际的算法中，可以利用简正波法求解
p
(
r
0
,
z
)
p(r_0,z)
p(r0,z)，再利用上式递推求解声场。
注：详细推导过程见《Computational Ocean Acoustic》Chapter. 6
在一些书中给出了更宽泛的模型和应用场景，但是我是做应用的，我其实推荐在大多数场景下使用射线模型，速度快，宽容性高，可解释性好，能够计算时域到达结构。
二、使用步骤
1. 下载Acoustics Toolbox环境
可以使用声学工具箱软件（Acoustics Toolbox）搭配Matlab来快速使用这些模型，在 Acoustics Toolbox 网站上可以下载这个at文件，一般下载 Binaries文件。
推荐下载我整编的工具箱，包含了ram_p等修复模型和函数，at_qp
将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文件夹。
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
运行结果为
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')
运行结果为
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’：是全向性（这是默认的源波束模式选项）。
运行结果
kraken和Scooter的模型对比结果如下：
总结
提示：以上就是今天要讲的内容，本文简单介绍了水下声场的简正波、射线、抛物方程模型原理和声学工具箱中KRAKEN、BELLHOP、RAM的使用。
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
77
0- 0
已为社区贡献1条内容
所有评论(0)