ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

基于KrakenC的水声传输损失仿真与参数调试

基于KrakenC的水声传输损失仿真与参数调试 简介面向水下声学建模与仿真学习者这份压缩包聚焦 KrakenC 与 Kraken 在声场计算和声传播中的应用解决利用开源水声工具快速评估传播损失的需求。包内仅含一个 krakenc_tl.m 脚本体积 772B属于轻量级 MATLAB 脚本通过调用 KrakenC 接口即可完成从几何模型导入、介质参数设置、网格构建到声源定义、波动方程求解及传播损失计算的可视化全流程。对于希望掌握 KrakenC 调用方法、理解声传播损失计算步骤的初学者或工程师具有直接参考价值。已有 598 人学习下载脚本虽短但串联了水声建模的关键环节可作为二次开发的模板便于替换几何模型和物理参数后快速迁移到不同水下环境场景也适合与 Kraken 原始版本对照体会二者在计算效率与接口设计上的差异。1. 从 krakenc_tl 这个文件名说起水声信道仿真绕不开 Kraken 这个简正波模型。krakenc_tl.zip 拆开看基本就是 KrakenC 编译好的可执行文件加一条计算传输损失TL的调用链路_tl 即 Transmission Loss。评估声纳作用距离、做水声通信链路预算、分析浅海声传播损失都要把声速剖面、水深、海底声学参数写成环境env文件让 KrakenC 先求模态再叠加声场、输出传播损失。Kraken 是 Fortran 原版KrakenC 是 C 重写版两者读同一套 env 文件后者在批处理和跨平台分发上更省事。这篇按实际跑 Kraken 的顺序写选型、env、field、参数调试最后给一个模式截断自检技巧。2. Kraken 与 KrakenC声场模型、编译器差异和 Acoustics Toolbox 选型2.1 Kraken 的简正波模型声场为什么用模态叠加而不是射线浅海波导里声波在自由海面和海底之间反复反射形成稳定的干涉图形。这种场景用射线模型算到达结构会碎成一大堆本征射线近距离还能数得清距离超过几十个水深之后基本没法看。简正波normal mode的思路完全不同把波动方程在深度方向上离散成本征值问题得到一组只由环境剖面决定的模态 $\psi_m(z)$ 和对应水平波数 $k_{r,m}$任意点声压写成这些模态的叠加$$ p(r,z) \frac{i}{\rho(z_s)}\sum_m \psi_m(z_s)\psi_m(z) H_0^{(1)}(k_{r,m}r) $$远场再把 Hankel 函数换成渐近形式就得到每阶模态按柱面扩展传播、带各自衰减系数的清晰物理图像。Kraken 与其他简正波模型的差异在于求解本征值问题的方式它把深度区间用非均匀网格离散用有限差分法求解局部本征值问题。这个做法对连续声速剖面、多层介质、弹性海底都能给出可靠的数值本征函数所以二十多年来一直是水声传播计算的事实基准。说明一下为什么是模态叠加而不是直接数值积分直接求解二维 Helmholtz 方程比如用抛物方程法要做逐距离步进每步都要处理一个大矩阵简正波法把深度和距离解耦先一次性算完深度方向的本征值后面所有距离上的声场只是对模态做解析组合。对需要扫描几百个距离点、几十个接收深度的 TL 任务简正波法快得多。代价是环境必须近似为水平分层介质海底地形变化剧烈时要在多个距离段分别建 env 文件拼接。2.2 KrakenC 与 Kraken 的差异C 版本在批处理里的优势Kraken 的 Fortran 原版可执行名通常是 kraken.exe和 C 重写版krakenc.exe在声学内核上保持一致读同一份 env 文件生成同名 .mod 模态文件送到 field 或 fieldc 里做后处理。也就是说你完全可以把 krakenc 算出的 .mod 拿给 Fortran 版的 field 用两者没有格式冲突。选 C 版主要是工程原因不需要在目标机器上装 gfortran 运行时交叉编译和静态链接都干净在 Python 里用 subprocess 调子进程做批量参数扫描时C 程序的启动开销和退出码行为也更可控。很多流传的 krakenc_tl.zip 这类压缩包实际就是把 krakenc.exe 和 fieldc.exe 加上几个示例 env 打成一个最小工具集命名里的 _tl 强调默认输出直接落到传输损失。需要注意KrakenC 不是把 Fortran 源码逐行翻译成 C而是按相同算法重写。个别版本的 C 实现在复数数学库依赖、模式筛选阈值Mode Cutoff 判断和输出精度上有细微差别。同一组环境参数kraken 和 krakenc 跑出来的模式数应该一致但相位速度的小数位可能差几位。批量扫描时我习惯固定用 krakenc不混着用避免把版本差异误当成物理结果。2.3 获取和编译Linux 与 Windows 环境对比常见做法是从 Acoustics ToolboxAT官网下载源码包或者使用网上流传的、已经附带了 Makefile 的编译产物。拿到源码后在 Linux 上先看根目录的 Makefile 和 Makefile_Linux 等文件cd acoustics_toolbox make -f Makefile_Linux如果没有现成的 Makefile单独编 C 版核心程序也很直接。KrakenC 通常由主程序加若干工具源文件组成标准做法是gcc -O2 -o krakenc krakenc.c readat.c modsediment.c -lm gcc -O2 -o fieldc field.c field_io.c -lm两个命令都依赖 POSIX 数学库-lm 必须带上编译成功与否最直接的验证是运行时不带参数./krakenc正常会打印一行用法说明提示需要 -f 指定前缀如果直接段错误或报缺失共享库先检查是否漏了某个 .c 文件。程序语言输入输出用途kraken / krakencFortran / C*.env*.mod, *.prt求解简正波模态与本征值field / fieldcFortran / C*.flp *.mod*.shd模态叠加计算声压场与 TLWindows 上建议直接用预编译 exe并开一个控制台把 exe 所在目录加进 PATH。如果自己用 MinGW 编注意 AT 源码是老式 C89 风格在较新 GCC 的高告警级别下会产生大量 warning但一般不影响产出。我一贯的做法是编完立刻用仓库自带的某个 .env 样例跑一遍把生成 .mod 的时间和日志留档后面换机器时才有个基准。提示编译报 undefined reference 时优先检查源文件列表是否完整以及 -lm 是否放在命令末尾。3. 用 KrakenC 算声场env 文件结构与最小可运行配置3.1 env 文件的字段拆解频率、声速剖面与海底参数Kraken 的 env 文件是一切的起点。它的通用结构是第一行标题然后频率、层数、上边界类型接着是层内材料和声速剖面SSP最后是半空间参数。不同 AT 版本之间这个文件的解析顺序会有细微移动最稳的办法是拿压缩包里自带的 .env 样例改而不是从空白创建。为了讲清字段含义这里给一个常见的文本格式 Pekeris 波导环境水深 100 米频率 50 HzPekeris: 100m isovelocity water 50.0 ! 频率, Hz 1 ! 水层半空间算 1 层 0 ! 上边界 0压力释放海面 1500.0 0.0 0.0 ! 层介质: 声速 1500 m/s, 密度 1.0, 衰减 0 SSP 1500.0 0.0 ! 声速剖面: 水深 0m 处声速 1500.0 100.0 ! 声速剖面: 水深 100m 处声速 BOTTOM 1700.0 2.0 0.5 ! 半空间: 声速 1700, 密度 2.0, 衰减 0.5 dB/λ频率决定模态个数和传播特性50 Hz 在这个水深下大约只能容纳 3 阶束缚模态换成 500 Hz 就是 30 阶量级不同频率的行为差距很大。上边界 0 对应声压释放海面绝大多数浅海场景都这么设。声速剖面在水层内部按深度给出若干采样点KrakenC 会线性插值到它的内部网格如果剖面是负梯度声速随深度减小模态会向下偏折模式数明显变化这是后面排查声场异常时最先怀疑的对象。海底半空间参数里的衰减以 dB/λ每波长分贝为单位而不是 dB/m这个单位换算错误是新手最容易踩的坑同样写 0.5在 50 Hz 和 5 kHz 下的物理衰减相差 100 倍。提示AT 不同发行版对 env 文件的具体解析可能有微调批量处理时先把自带样例跑通再复制成模板比对着文档盲改省时间。3.2 运行 krakenc 生成模态文件命令与日志解读env 文件准备好之后运行 KrakenC 的主程序。AT 系列程序统一用前缀约定-f 后面的名字不带扩展名程序自己去找对应的 .env 文件并生成同前缀的 .mod 和 .prt。以当前目录的 pekeris.env 为例./krakenc -f pekeris正常运行会在终端打印介质参数、网格点数和每阶模态的摘要行。跑完检查产物ls -la pekeris.mod pekeris.prt head -30 pekeris.prt产物中 .mod 是二进制模态文件后续 field 程序直接读它.prt 是文本摘要包含本征值、相位速度、衰减系数等信息。我一般不看终端回显而是直接 head 一下 .prt因为这里记录了每阶模态的衰减项和水平波数能直观看出某几阶模态是否被海底吸收压得很低。跑不起来时报错基本三类找不到文件-f 前缀写得和文件名不一致读入 env 的字段类型不匹配比如把密度那一列填了字符串还有一类是 0 个模态被找到这种通常是频率太低或者声速剖面使所有模式都截止。第一类问题看文件名第二类检查 env 每一行的列数第三类则是 2.1 节里说的物理问题需要通过调频率或检查剖面解决。3.3 声场后处理field 程序与 krakenc_tl 的分工krakenc 只负责求模态不直接给传播损失。真正把模态叠加成声场并输出 TL 的是 field或 C 版 fieldc。如果压缩包里有一个叫 krakenc_tl 的脚本或可执行它的作用通常就是把这两步串起来调 krakenc 算模态再调 fieldc 按接收网格输出传输损失。field 程序读一个 .flp 文件作为参数文件。常见字段顺序如下Pekeris TL field 50.0 ! 频率, Hz 1 ! 声源个数 20.0 ! 声源深度, m 51 ! 接收深度个数 0.0 100.0 ! 接收深度范围: 0~100m, 51 个点 0.0 10000.0 101 ! 距离范围: 0~10km, 101 个点我把这个文件命名成 pekeris.flp然后执行# 某些版本可执行名是 field.exe对不上时用 ls 确认 ./fieldc -f pekeris输出产物是 pekeris.shd二进制声场文件头部是文本参数段后面是复数声压数组。拿到 .shd 后可以用 Python 读取并画 TL 切片。下面这段代码是读取 AT 二进制声场文件并提取 TL 的简便写法import numpy as np import struct def read_shd(filename): with open(filename, rb) as f: # 头部文本段以空行为界 while True: line f.readline() if line.strip() b: break # 三个维度点数 nsd, nrd, nrr struct.unpack(3i, f.read(12)) # 依次读源深、接收深度、距离数组 sd np.frombuffer(f.read(4*nsd), dtypef4) rd np.frombuffer(f.read(4*nrd), dtypef4) rr np.frombuffer(f.read(4*nrr), dtypef4) # 复数声压单精度交错存储 n nrr * nrd * nsd raw np.frombuffer(f.read(8*n), dtypef4).reshape(-1, 2) pressure (raw[:, 0] 1j*raw[:, 1]).reshape(nrr, nrd, nsd, orderF) return rr, rd, sd, pressure rr, rd, sd, p read_shd(pekeris.shd) tl -20*np.log10(np.abs(p) 1e-12) np.save(pekeris_tl.npy, tl)读取思路是头部文本段以空行为界接着三个整数表示三个维度的点数再依次读出源深数组、接收深度数组、距离数组最后读复压数据并重排成 (距离, 深度, 源深) 的结构。这里 orderF 表示按 Fortran 列优先重排AT 的输出是 Fortran 风格数组如果换一个 AT 版本发现切片形状对不上把 order 参数改掉即可。4. 传输损失计算与参数调试krakenc_tl 链路里的 3 个关键参数4.1 接收网格与距离采样flp 文件的排布方式flp 文件决定 TL 输出的形状。源深个数、接收深度个数、距离点数这三个值直接影响 .shd 的数据量和计算时间。以 3.3 节的示例为例输出是 101 个距离乘 51 个深度再乘 1 个源深的复数矩阵单精度存储不到 0.5 MB跑起来没有压力。但如果把距离点数加到 5000、深度加到 200存储就是 8 MB 量级加上 field 在每一距离点对几十阶模态求和耗时会上来且肉眼已经不可能从曲线里看出更多信息。距离采样有一个经验下限相邻距离点的相位差不要超过 π否则两点的干涉条纹是混叠的。经验做法是让距离步长小于水层中最小声速对应波长的 1/10即dr c_min / (10 * f)100 米水深、50 Hz、声速 1500 m/s 时dr 取 3 m 就足够密如果频率升到 5 kHzdr 必须缩到 0.03 m 以下否则 TL 曲线会看到高频抖动那不是物理是采样混叠。flp 字段含义常见取值NSD声源深度个数1 到几十SD源深数组或单个值布放在声轴附近NRD接收深度个数10 到几百RD接收深度范围覆盖水层NRR/RR距离点数与范围按 dr 反推参数说明接收深度范围不要只设到水层顶和底。海底半空间内部 Kraken 也返回声压如果关心海底穿透把 RD 上界放到海底以下几十米浅海信道仿真有时特意看海底折射路径这个数据是有用的。4.2 模式数截断模态数量不够时的表现与快速估算模态数mode count是整个计算里最需要人工干预的参数。模式太少远场会丢掉高阶模态的干涉TL 曲线在中距离出现不自然的平滑模式太多泄漏模态的数值噪声会把远场搅乱。KrakenC 的 env 文件里可以显式限制最大模式数不设时程序按内部阈值判断但这个判断结果经常和物理直觉不一致。对于等声速水层加液态半空间这种最简波导束缚模态数存在解析上限。束缚模态要求水平波数大于海底波数可以得到近似公式$$ N_{\max} \approx \left\lfloor 2 f D \sqrt{\frac{1}{c_w^2} - \frac{1}{c_b^2}} 0.5 \right\rfloor $$其中 $D$ 是水深$c_w$ 是水中声速$c_b$ 是海底声速。这个公式不用背可以直接放进脚本里做检查。下面是一个最小 Python 实现import numpy as np def max_bound_modes(freq_hz, depth_m, cw1500.0, cb1700.0): 估算 Pekeris 波导的束缚模态数量上限 if cb cw: raise ValueError(海底声速必须大于水中声速才存在束缚模态) return int(2 * freq_hz * depth_m * np.sqrt(1.0/cw**2 - 1.0/cb**2) 0.5) print(max_bound_modes(50.0, 100.0)) # 输出约 3参数说明公式里的 0.5 是修正项来自模态本征值的相角边界频率、深度的单位必须分别是 Hz 和 m声速用 m/s。如果代码算出来是 3而你的 krakenc 日志里出现了 6 阶模态多出来的那几阶基本是泄漏模态在 TL 远场会被海底衰减消耗掉但近场会有异常振荡。这时候就该在 env 里显式把模式数上限设为 3再对比一次结果。4.3 三个必调参数频率、海底吸收与源深第一个必调参数是频率。env 文件里的频率不是算个样看看的摆设它决定模态数和每阶模态的水平波数。同一条声速剖面50 Hz 和 500 Hz 算出的 TL 干涉结构完全不同。做宽带仿真时正确做法不是把某几个频率堆进一个 env而是用脚本按频点批量生成 env 文件逐个跑 krakenc 和 fieldc最后把各频点的 TL 加权合并。第二个必调参数是海底吸收系数。它的单位是 dB/λ与常见的 dB/m 换算关系是alpha_dB_per_m alpha_dB_per_lambda * f / cb例如半空间声速 1700 m/s吸收 0.5 dB/λ在 50 Hz 下约等于 0.0147 dB/m在 5 kHz 下约等于 1.47 dB/m。远场 TL 曲线的斜率对海底吸收极其敏感调参时先保持声速和密度不变只扫 α看哪条曲线的远端衰减接近实测。如果扫了 α 仍然偏快再怀疑海底密度而不是把原因归结到水深或声速剖面上。第三个必调参数是源深。源深通过本征函数在源位置的幅度调制每阶模态的激发强度。源放在本征函数反节点处该模态被强烈激发放在节点附近该模态几乎不出现。实际效果是源深差 1 米近场 TL 可能差 10 dB 以上。field.flp 里的源深数组一次可以填入多个值KrakenC 的模态不用重算只有叠加阶段重复执行这是一个很划算的批量扫描方式。调试时遇到锯齿状抖动先减模式数远端衰减异常先改 α近场剧烈振荡先看距离采样。和我 2.2 节说的一样任何调整都要固定可执行文件版本kraken 和 krakenc 的结果不要混在一起对比。5. 模式截断自检与 krakenc_tl 曲线的自动化验证5.1 用一条命令自动生成 TL 并检查模态数手动改 env、跑 krakenc、看 .prt 的流程只适合调试。批量环境扫描时我一般把三件事写进同一个脚本用 4.2 节的公式估算 $N_{\max}$调用 krakenc 生成模态然后对比实际模态数与估算值的差距。一个实用的检查函数是这样import subprocess, sys def run_tl(prefix, env_text, freq, depth, cw1500.0, cb1700.0): with open(prefix .env, w) as f: f.write(env_text) proc subprocess.run([./krakenc, -f, prefix], capture_outputTrue, textTrue) if proc.returncode ! 0: print(proc.stderr) sys.exit(1) n_est max_bound_modes(freq, depth, cw, cb) # 从 .prt 里数模态行数不同版本行格式略有不同 cnt 0 for line in open(prefix .prt): if line.strip().startswith(Mode): cnt 1 print(festimated: {n_est}, computed: {cnt}) if cnt n_est 2: print(warning: possible leaky modes, check cutoff settings) # env_text 由调用方按模板生成这里示意 env Pekeris 50Hz 50.0 1 0 ... run_tl(scan01, env, 50.0, 100.0)逻辑说明进程调用失败时直接看 stderr模式数从 .prt 的文本中按行计数如果换的 AT 版本里 .prt 的行首不是 Mode改成统计以数字开头的行即可。这里设置cnt n_est 2的容差是因为泄漏模态偶尔会多出 1 到 2 阶超过这个数就值得人工介入。5.2 用 TL 曲线斜率做快速冒烟测试自动化流程里比看数值更省事的是看曲线形状。浅海 TL 在近场大约按球面扩展衰减每十倍距离 20 dB远场转成柱面扩展每十倍距离 10 dB并叠加模态衰减。固定一个深度切片用下面的代码快速判断import numpy as np def check_slope(r, tl_slice, r_range): 返回 TL(dB) 对 log10(r) 的斜率用于粗判扩展规律 mask (r r_range[0]) (r r_range[1]) coeff np.polyfit(np.log10(r[mask]), tl_slice[mask], 1) return round(float(coeff[0]), 2)如果返回的斜率在 -15 到 -5 之间说明场落在柱面扩展区域基本可信如果斜率大于 5 或小于 -40说明数据有问题优先排查模式截断和海底吸收这两个参数。这个检查不依赖具体深度也不依赖精确的声源级归一化可以作为 CI 里每个 env 变更的冒烟测试。以上两段技巧配合 4.2 节的公式足以在参数扫描时把物理异常和代码异常区分开。把这个检查固化成脚本后每次修改 env 文件只需要跑一次输出的 estimated/computed 两个数就能定位问题出在模式截断还是海底参数上这也是网上那些 krakenc_tl.zip 类工具包里最值得自己重构的一层。本文还有配套的精品资源点击获取
返回列表