基于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 Toolbox（AT）官网下载源码包，或者使用网上流传的、已经附带了 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 文件。

Windows 上建议直接用预编译 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 Hz：

'Pekeris: 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), dtype='<f4') rd = np.frombuffer(f.read(4*nrd), dtype='<f4') rr = np.frombuffer(f.read(4*nrr), dtype='<f4') # 复数声压，单精度交错存储 n = nrr * nrd * nsd raw = np.frombuffer(f.read(8*n), dtype='<f4').reshape(-1, 2) pressure = (raw[:, 0] + 1j*raw[:, 1]).reshape(nrr, nrd, nsd, order='F') 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)

读取思路是：头部文本段以空行为界，接着三个整数表示三个维度的点数，再依次读出源深数组、接收深度数组、距离数组，最后读复压数据并重排成 (距离, 深度, 源深) 的结构。这里 order='F' 表示按 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 kHz，dr 必须缩到 0.03 m 以下，否则 TL 曲线会看到高频抖动，那不是物理，是采样混叠。

参数说明：接收深度范围不要只设到水层顶和底。海底半空间内部 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, cw=1500.0, cb=1700.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, cw=1500.0, cb=1700.0): with open(prefix + '.env', 'w') as f: f.write(env_text) proc = subprocess.run(['./krakenc', '-f', prefix], capture_output=True, text=True) 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(f'estimated: {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 类工具包里最值得自己重构的一层。

本文还有配套的精品资源，点击获取