KrakenC简正波模型详解：从环境文件到声传播损失计算

简介面向水下声学建模与声传播计算领域的研究者与工程师一份围绕KrakenC与Kraken工具的声场仿真压缩包解决复杂环境中声传播损失的计算问题。压缩包共1个文件类型为MATLAB脚本大小仅772B核心脚本通过调用KrakenC库在MATLAB中完成几何模型导入、物理参数设置、网格构建、声源定义、波动方程求解及传播损失计算并支持结果可视化。已有598人学习下载。该工具特别适用于大规模声学模拟任务可帮助用户快速搭建声场分析流程为水下环境评估、声纳性能分析及噪声控制研究提供实用参考。资源虽小但功能明确适合需要借助KrakenC进行声传播建模的工程技术人员。1. Kraken与KrakenC解决什么声场问题KrakenC 是 Kraken 的 C 语言移植版krakenc_tl.zip 这个包名里的 tl 是 transmission loss 的缩写也就是声传播损耗。Kraken 基于简正波理论用于计算水平分层介质中的声场分布给定声速剖面、海底参数、声源频率和源深它求解特征值问题得到各阶简正波再叠加得到整个距离-深度平面上的传播损耗。相比射线模型它在低频浅海下不会出现焦散线失真相比抛物方程和全波数值解它算一个频点通常只要几秒。因此做声呐性能预报、水声通信评估、地声参数反演的工程师多数会把 KrakenC 当作第一套跑通声传播链路的选择。2. 准备KrakenC环境文件声速剖面、海底参数与频率2.1 环境文件的段落结构与单位约定KrakenC 的输入是一个文本环境文件后缀习惯用 .env。它不读数据库也不识别配置文件所有环境信息按固定顺序写在同一个文件里。下面是一个浅海场景的完整示例Shelf_200m 200.0 1 0.0 200.0 4 0.0 1500.0 0.0 1.0 50.0 1510.0 0.0 1.0 100.0 1525.0 0.0 1.0 200.0 1550.0 0.0 1.0 200.0 1600.0 0.5 1.8按顺序读一下这些行第一行是标题单引号括起来只作为输出文件命名的前缀第二行是声源频率 200.0Hz第三行是介质层数 1表示水柱只建一层第四行是水层深度范围从 0 到 200m第五行是声速剖面采样点个数 4接下来的四行每行是一个采样点的深度、声速、吸收系数和密度第十行是海底半空间的起始深度 200.0m必须与水层底界对上最后一行是半空间的声速、吸收系数和密度。注意这里“深度、声速、吸收、密度”的列顺序是我在这个版本里常用的部分移植版会写成“深度、声速、密度、吸收”拿到新包先看自带样例确认列顺序再套自己的数据。单位问题值得多说一句。声速用 m/s密度用 g/cm³吸收系数用 dB/波长而不是 dB/m。这个 dB/波长在换算上最容易踩坑同样一个 0.5dB/波长的海底吸收在 200Hz 下对应的 dB/m 数值约等于 0.0033在 2kHz 下约等于 0.033频率变了物理含义完全不同。不要把 dB/波长直接当成 dB/m 用否则长距离 TL 会被算得离谱。我把常用字段整理成一张表调参时对着改即可字段示例值单位对结果的影响频率200.0Hz决定截止频率和模态数水层深度范围0.0 200.0m超出SSP范围会截断积分SSP采样点数4个采样过少会丢失声道结构采样点深度0.0m必须单调递增声速1500.0m/s决定模态的相速度和群速度吸收系数0.0dB/波长水层一般取0密度1.0g/cm³影响边界反射系数半空间深度200.0m必须和水层底界一致半空间声速1600.0m/s控制泄漏模态数量半空间吸收0.5dB/波长控制远场TL斜率这套字段的物理背景是简正波理论里的水平分层假设。介质在水平方向无限延伸且性质只随深度变化声压场按深度方向的本征函数展开KrakenC 的核心工作就是数值求解这个本征问题。环境文件里的每一行本质上是在刻画某一层的声速剖面和边界条件。2.2 用Python从CTD数据生成KrakenC环境文件实测数据通常来自CTD或XBT先要用声速公式把温盐深换算成声速得到一组离散的深度-声速数组。直接手写进 .env 经常会因为深度排序、采样点间距等问题让 krakenc 报错更可靠的方式是用脚本生成环境文件。import numpy as np def write_env(fname, freq, z, c, rho1.0, alpha0.0, bottom_c1600.0, bottom_rho1.8, bottom_alpha0.5): # 按深度排序避免CTD数据出现倒序 order np.argsort(z) z np.asarray(z)[order] c np.asarray(c)[order] depth_max z[-1] lines [CTD derived SSP, f{freq:.1f}] lines.append(1) lines.append(f0.0 {depth_max:.2f}) lines.append(str(len(z))) for zi, ci in zip(z, c): lines.append(f{zi:8.2f} {ci:8.2f} {alpha:5.2f} {rho:5.2f}) lines.append(f{depth_max:.2f}) lines.append(f{bottom_c:.2f} {bottom_alpha:.2f} {bottom_rho:.2f}) with open(fname, w) as f: f.write(\n.join(lines) \n) print(f[ok] {fname} generated, {len(z)} samples)调用时直接传入CTD质量控制后的数组z np.array([0, 20, 50, 100, 150, 200]) c np.array([1500.2, 1506.5, 1518.0, 1532.4, 1546.1, 1552.8]) write_env(ctd_200m.env, freq200.0, zz, cc)这个脚本里值得关注的参数有三个。freq直接决定模态的截止频率200Hz 的浅海环境一般能激发出 8 到 15 阶传播模态alpha传 0.0 表示水层没有内部吸收这在深海是合理假设但在高浑浊度浅海里要按实测调整bottom_alpha是海底吸收系数它控制 TL 曲线远场的斜率常见取 0.1 到 0.5 dB/波长。脚本把depth_max默认取成 SSP 最后一个点的深度这是有意的。如果环境文件里水层底界大于 SSP 覆盖范围krakenc 会对声速做恒定外推在海底地形变化大的区域容易产生虚假的泄漏模态。让水层底界和 SSP 末点保持一致可以减少一个控制变量。提示如果手里是 Fortran 版 Kraken 的用户环境文件可以直接沿用。KrakenC 在文件格式上基本兼容唯一要留意的是浮点解析更严格SSP 里出现 NaN 或 Inf 会直接中断而不是像 Fortran 版那样忽略。2.3 浅海与深海场景的SSP采样和海底参数经验值浅海声传播的显著特征是声波频繁与海底相互作用海床参数对 TL 的影响往往大于水体声速剖面。水深小于 200m 时我会把 SSP 采样点加密到每 5m 一个点因为在温跃层附近声速梯度变化剧烈采样太稀会把模态的垂直结构算偏。海底参数按沉积物类型取砂质海底取声速 1700m/s、密度 2.0g/cm³、吸收 0.5dB/波长泥质海底取声速 1550m/s、密度 1.5g/cm³、吸收 0.2dB/波长。不确定时取两者的中间值或偏大吸收得到的声呐作用距离估计偏保守工程上更安全。深海场景里声道轴附近的 SSP 采样点要加密。声道轴通常在 800 到 1200m 之间这里的声速最小模态能量被约束在声道内可以传播上千公里。如果采样点在声道轴附近太少krakenc 会把声道轴位置算偏TL 曲线在远场出现不该有的波动。一个通用的判断原则是最高计算频率对应的波长 λSSP 采样间隔不要大于 λ/4。400Hz 对应约 3.75m 的波长SSP 采样间隔最大取 0.9m 才稳妥。高频计算时手工加密不可取用脚本按此规则插值最省事。3. 用krakenc与fieldc跑通声传播计算TL导出的完整过程3.1 krakenc_tl.zip里的程序组成与验证顺序拿到的 KrakenC 发布包一般会包含下面这些可执行文件。名字带不带 .exe 取决于平台Linux 下通常没有后缀文件名命令行作用krakenc./krakenc envfile解算简正波特征值和模态函数输出 .modfieldc./fieldc envfile用模态叠加计算声场输出 .shdshdfile./shdfile xxx.shd把二进制 .shd 导出为 ASCII 文本plotssp./plotssp envfile绘制声速剖面快速检查 SSPplotmodes./plotmodes envfile绘制模态函数检查截断这些程序在 Fortran 版 Kraken 中的名字是 kraken.exe 和 field.exeKrakenC 版改为 krakenc 和 fieldc用法一致。拿到包后不要直接跑自己的数据先用发布包自带的例子验证一遍。跑通一个 sample.env确认 .mod 和 .shd 能正常生成再替换成自己的环境文件。这套动作能把“程序本身有问题”和“我的环境文件有问题”分开省掉大量排查时间。如果包里没有可执行文件需要用 CMake 或 Makefile 编译编译时开启优化选项简正波求解在浮点运算密集时会快很多。3.2 用krakenc解算简正波命令和输出解读假设第 2 章生成的环境文件是 ctd_200m.env第一步是解算简正波# 解算简正波输出 ctd_200m.mod ./krakenc ctd_200m.env终端输出会依次显示频率、介质层数、SSP 采样点数以及一列特征值。每个特征值对应一阶模态的水平波数实部是水平传播常数虚部表示沿距离方向的衰减。如果输出里看不到任何模态说明该频率低于截止频率或者海底参数导致模态全部截止先检查环境文件再谈调参。运行成功后工作目录多出 ctd_200m.mod。这个文件保存了模态函数在深度网格上的采样值是 fieldc 的输入之一。注意 krakenc 只依赖 .env不依赖声源位置所以换源深、换接收深度时不需要重新解模态直接重跑 fieldc 即可这是简正波模型的核心优势。3.3 用fieldc计算传输损耗源深与接收网格设置fieldc 默认从 .env 读取源深用内置的默认接收网格计算声场。如果要用自定义的源深和接收范围在 .env 末尾追加一段参数SrcDepth, 30.0 RMax, 50.0 NProf, 201 NDist, 251含义依次是源深 30m最大接收距离 50km垂向接收点 201 个水平距离点 251 个。追加后重新跑 krakenc 和 fieldc输出的 .shd 就对应这套网格。# 重新解算模态并计算声场 ./krakenc ctd_200m.env ./fieldc ctd_200m.envfieldc 的计算量随模态数和网格点数线性增长201×251 的网格在普通笔记本上几秒内完成优势非常明显。提示部分版本会校验接收网格必须落在环境文件定义的水层范围内如果报 “RMAX exceeds water depth”要么减小 RMax要么把水层深度范围调大。3.4 解析.shd并绘制TL等值线.shd 是二进制文件文本编辑器打开是一堆乱码。先用发布包里的 shdfile 工具导出# 把二进制声场文件导出为文本 ./shdfile ctd_200m.shd ctd_200m_tl.txt导出的文本里每一行对应一个距离点行首是距离值后面是该距离上各个深度位置的 TL。用 Python 画等值线图最直观import numpy as np import matplotlib.pyplot as plt # 读入导出的ASCII文件第一列是距离其余列是TL raw np.loadtxt(ctd_200m_tl.txt, dtypefloat) r_km raw[:, 0] tl raw[:, 1:] n_depth tl.shape[1] plt.figure(figsize(8, 6)) plt.pcolormesh(r_km, np.arange(n_depth), tl.T, shadingauto) plt.colorbar(labelTL (dB)) plt.xlabel(Range (km)) plt.ylabel(Depth index) plt.gca().invert_yaxis()等值线图画出来后第一眼看近场 0 到 1km 的 TL 是否随距离平滑下降。如果出现锯齿状跳跃通常是模态截断不足、SSP 采样点太少或者半空间参数不匹配。把这三个方向逐个排除大多数异常都能定位到具体原因。4. 读懂KrakenC输出模态数、海底吸收与调参方向4.1 .mod和.shd文件里存了什么.mod 文件的物理内容是模态函数。用 plotmodes 把它画出来可以看到第 1 阶模态能量集中在哪个深度、第 10 阶模态在海底附近有没有残余能量。.shd 文件存的是声压或 TL 的二维网格排列顺序一般是距离优先、深度其次。KrakenC 输出 TL 时默认按 -20·log10(|p|) 换算但个别编译版本保留复数压力场读数据时要先确认包里的 README 对 .shd 格式的说明。4.2 模态数截断怎么判断算少了模态截断是简正波模型误差的主要来源。物理上200Hz 的浅海环境往往有十几阶传播模态阶数越高衰减越大。KrakenC 求解特征值时会自动丢弃衰减过大的模态但用户仍可以限制最大模态数在 .env 末尾追加一行NMODES, 20判断截断是否合适最简单的办法是用 plotmodes 看最高阶模态。如果最高阶模态在水层底部仍有明显幅度说明更高阶模态可能还有能量模态数需要增加。更定量的方式是从 .mod 文件里读特征值检查虚部import numpy as np # 读取.mod文件每个模态记录水平波数的实部和虚部 data np.fromfile(ctd_200m.mod, dtypenp.float64) N len(data) // 2 for m in range(N): kr data[2 * m] ki data[2 * m 1] if -ki 0.1: print(fmode {m:2d}: k_r{kr:8.4f} /m, alpha{-ki:.5f} dB/m)这里的虚部单位是 dB/m表示模态在传播方向上的衰减系数。经验是传播距离 10km 以内衰减小于 0.1dB/m 的模态都保留超过 50km 的远场只保留衰减小于 0.02dB/m 的模态。保留过多远距离衰减模态会把数值误差放大反而不利于远场分析。浅海环境里模态数太少TL 曲线在近场会出现明显的振荡。如果一个 200Hz 的浅海案例计算出的传播模态少于 5 阶先怀疑海底声速是不是设得过高或者频率低于实际截止频率。提示模态截断不是越多越好。模态数过多时高阶模态的数值噪声会被叠加进总场使 TL 曲线出现高频抖动。如果抖动集中在某个深度区间优先检查该深度附近是否有强梯度界面。4.3 海底参数对TL的敏感性一个可复现的对照海底参数的取值对 TL 的影响往往比 SSP 大。复制环境文件把海底吸收系数从 0.5 改成 0.05其他不动# 复制环境文件并精确替换最后一行海底吸收参数 cp ctd_200m.env ctd_200m_lowalpha.env sed -i s/^1600\.0 0\.5 1\.8$/1600.0 0.05 1.8/ ctd_200m_lowalpha.env ./krakenc ctd_200m_lowalpha.env ./fieldc ctd_200m_lowalpha.env这里的 sed 正则精确匹配最后一行海底参数避免误改水层里的数字。跑完后用 3.4 节的 Python 代码比对两个 .shd 的远场 TL会看到低频段差异小中高频段 50km 处 TL 可能差 5 到 10dB。这说明当实测海底参数不确定时把吸收系数设偏大得到的声呐作用距离是保守估计工程上更稳。这种做法在做地声反演之前特别有用用对照实验量化参数不确定性的影响范围。5. 用KrakenC做扫频声场分析批量计算与结果对比5.1 批量修改频率并逐个计算TL的脚本声呐系统的频带通常覆盖多个倍频程比如 25Hz 到 800Hz。手动改环境文件跑二三十个频点不现实用一个 bash 循环即可# 依次计算6个频点的模态和声场 for f in 25 50 100 200 400 800; do # 把环境文件第二行的频率替换为当前值 sed s/^200\.0$/${f}/ ctd_200m.env ctd_f${f}.env ./krakenc ctd_f${f}.env ./fieldc ctd_f${f}.env mv ctd_f${f}.shd shd_f${f}.shd mv ctd_f${f}.mod mod_f${f}.mod done这个循环里sed只替换第二行的频率值正则^200\.0$保证不会波及声速剖面里的数字。每个频率的 .mod 和 .shd 分别改名保存避免下一轮覆盖。需要注意不同频率的模态数差异很大800Hz 可能有几十阶模态25Hz 可能只有 1 到 2 阶这是正常现象。如果某个频率点算出的模态数为 0说明低于截止频率该频点直接跳过。5.2 不同频率下的TL曲线对比与频散分析批量计算完成后取固定深度位置画 TL-距离曲线可以快速看清单频传播特性。下面这段代码读取多个 .shd取第 100 个深度索引对应约 100m 深度import numpy as np import matplotlib.pyplot as plt def load_tl(fname, nprof201, ndist251): 读取二进制.shd返回二维TL矩阵形状为(ndist, nprof) d np.fromfile(fname, dtypenp.float32) return d.reshape(ndist, nprof) plt.figure(figsize(8, 5)) r_km np.linspace(0, 50, 251) # 取第100个深度索引约对应100m深度 for f in [50, 200, 800]: tl load_tl(fshd_f{f}.shd) plt.plot(r_km, tl[:, 100], labelf{f} Hz) plt.xlabel(Range (km)) plt.ylabel(TL (dB)) plt.legend() plt.gca().invert_yaxis() plt.grid()读出来的 TL 曲线一般符合“低频衰减慢、高频衰减快”的规律。如果某一频率的 TL 在远场异常回升大概率是该频率模态截断不足回到 4.2 节的检查步骤。多频扫频除了对比损耗还能提取频散特性不同频率的模态传播速度不同处理宽带通信信号时特别重要。把多个频率的群速度画在一张图上可以看到声道内的频散曲线。频散曲线对水声通信的作用很直接时延扩展等于群速度随频率变化的斜率乘以传播距离用这组数据可以估计不同频段的码间干扰程度给均衡器设计提供输入。本文还有配套的精品资源点击获取