简介：面向水声通信与海洋声学领域的研究人员及工程师，基于Kraken模型开发，用于计算声传播损耗并支持水平阵与垂直阵各阵元接收波形分析，适用于声纳性能评估、水声信道研究及阵形优化。压缩包共19个文件，约713KB，主体为3个Matlab脚本与测试数据，另有12张结果图、说明文档、公式文档和任务图示，目录划分明确。已有1108人学习下载。阅读README和示例代码，可快速掌握海洋环境参数设置、声压级衰减曲线绘制以及阵列指向性分析方法；result目录下的传播损耗、深度距离图像与LFM波形图，能直观对比不同阵型的接收波形差异。这套资源覆盖从环境参数建模、传播计算到阵列输出分析的完整流程，适合相关课程实践或二次开发，便于深入理解不同海洋条件下的传播特性。
1. Kraken 计算传播损耗：一个简正波工具如何变成波形发生器
在水声工程里，Kraken 这类简正波模型常被当成“算损耗的黑匣子”，除了画 TL 曲线之外似乎没有别的用途。实际上，配合 Matlab 做频域合成之后，它完全可以输出每个接收阵元上带相位信息的复声压，再通过逆傅里叶变换还原成时域接收波形。这对做水听器阵列信号仿真的人来说非常实用：不需要铺设整条试验链路，在代码里就能评估不同声速剖面、海底沉积层参数和阵元深度对接收波形的拉扯。这篇笔记面向两种人：一是刚接触 Kraken、想知道环境文件怎么配、传播损耗怎么算的入门者；二是已经在跑 TL 但想把声场数据进一步加工成阵元波形的熟手。你会看到环境文件配置、Matlab 调用的完整思路，以及我踩过的几个滤波器级参数误区。
2. 搭好 Kraken 的“翻译层”：环境文件格式与 Matlab 调用方式
2.1 环境文件里最容易被忽略的三个量
Kraken 的输入核心是环境文件，后缀通常为.env。这个文件描述的是声速剖面、介质分层、海底沉积层和半空间参数。刚上手的人容易把注意力放在声速剖面的形状上，却忽略了三个同样致命的量：频率、介质层数和界面类型。
频率直接影响简正波数量，Kraken 在运行时自动把计算得到的模态数写进日志文件，如果你指定的频率低于该环境下的截止频率，模态数可能会掉到 1 甚至 0，最后算出来的传播损耗大幅偏低或干脆接近 0。介质层数NMedia指的是水体和沉积层的分层数，不是阵元数，搞混之后环境文件的读取位置会对不上，程序直接报错。界面类型在水体表面通常是理想压力释放界面，海底则要看是硬海底还是弹性海底，写错后低频段的传播损耗会偏差好几 dB。
2.2 用 Matlab 生成环境文件而不是手写
我一般不会手工去编辑.env文件，因为 Fortran 格式对列宽有要求，小数点位置错一位，Kraken 就会读出荒谬的声速梯度值。更稳妥的做法是在 Matlab 里写一个小函数，把参数传进去，按固定格式打印出环境文件。下面这段代码来自我常用的生成函数，只保留最关键的段落：
function write_env(freq, depth, ssp, envbase) % write_env 生成 Kraken 环境文件 % freq: 计算频率 (Hz) % depth: 水深 (m) % ssp: N×2 矩阵，[深度(m), 声速(m/s)] % envbase: 输出文件名前缀 fid = fopen([envbase, '.env'], 'w'); fprintf(fid, '''%s''\n', [envbase, '.env']); fprintf(fid, '%8.2f ! freq\n', freq); fprintf(fid, '1 ! NMedia\n'); % 单层水介质 fprintf(fid, '0 ! Nsd\n'); % 源深度由后文指定 fprintf(fid, '0 ! Nrd\n'); % 接收深度由 field 程序指定 fprintf(fid, ' %8.1f\n', depth); % 水深 fprintf(fid, ' %d\n', size(ssp, 1)); % 声速剖面采样点数 % 声速剖面 for k = 1:size(ssp, 1) fprintf(fid, ' %8.1f %8.2f\n', ssp(k, 1), ssp(k, 2)); end fprintf(fid, ' %d\n', 0); % 无分层沉积层 % 海底沉积层参数：声速、衰减、密度 fprintf(fid, ' 1600.0 0.5 1.5\n'); % 半空间参数 fprintf(fid, ' 0.0\n'); fclose(fid); end这个函数把声速剖面按[深度, 声速]传进来，然后逐行写入环境文件。NMedia固定为 1 表示把水层当作单一介质，沉积层参数直接跟在水层后面。这里最关键的是%8.1f和%8.2f这样的宽度控制，Kraken 读取环境文件时按固定列宽扫描，如果输出宽度不足，会把后续数据挤到错误的位置上。在写这个函数之前，我用手工编辑的方式翻过几次车，典型现象是明明声速剖面正确，算出来的模态却全部集中在表面，后来才发现是小数点位没对齐。
2.3 在 Matlab 里直接调用 Kraken 可执行文件
环境文件生成后，接下来就是调用 Kraken 和 field 两个程序。在 Windows 系统下，我习惯用 Matlab 的system函数直接执行，这样能把整个链路放进同一个脚本里。代码结构通常是这样：
envbase = 'mymodel'; write_env(100, 200, ssp, envbase); % 运行 Kraken 计算本征值和本征函数 [status, logout] = system(['kraken ', envbase, ' > kraken.log']); if status ~= 0 error('kraken 运行失败，查看 kraken.log'); end % 运行 field 计算声场 [status, logout] = system(['field ', envbase, ' > field.log']); if status ~= 0 error('field 运行失败，查看 field.log'); endsystem的返回值能帮你快速判断程序是否崩溃，而不是等到读取输出文件时才发现里面全是空数据。Kraken 计算的是频域本征函数，field 程序则基于这些本征函数算出各个接收点的复声压。注意field不一定需要重新跑kraken，只要环境文件没变，模态文件可以直接复用，这在扫频计算中能省下大量时间。
3. 传播损耗计算：从单条曲线到距离-深度切片
3.1 读懂.shd文件输出与归一化约定
Kraken 的 field 程序输出是.shd文件，里面保存的是复声压，单位不是 dB，而是带幅度和相位的线性值。很多人第一次读取时直接用幅度画图，发现数值小到离谱，这是因为输出经过了球面扩展归一化。传播损耗 TL 的通用定义是：
% 将复声压转换为传播损耗 p = read_shd('mymodel.shd', rvec, zrec); TL = -20 * log10(abs(p));也就是说，Kraken 输出的声压已经默认你拿它和自由场球面扩展做比较。如果某个接收点的声压幅度是 0.001，对应的传播损耗是 60 dB，这个量纲在后续做波形合成时很容易出错。之前我给某实验室做模拟项目 X 的数据处理时，就因为直接用abs(p)做波形幅度，导致接收信噪比整体偏离，排查了很久才发现归一化这一步漏掉了。
3.2 用 Matlab 自动批量扫频并绘制距离-深度切片
实际工程里不会只算一个频点，往往需要在带宽内扫几十个频率。如果每个频率都手动改环境文件、重新调用 Kraken 和 field，效率非常低。我会写一个批量循环脚本，把频率列表作为外层循环，内部逐次更新环境文件并记录现场数据。下面是一个典型的calc_tl_slice函数，它返回距离-深度切片上的传播损耗矩阵。
function [TL, rvec, zrec] = calc_tl_slice(envbase, rvec, zrec, freq) % 计算单频点下的距离-深度传播损耗切片 % rvec: 1×M 水平距离向量 (m) % zrec: 1×N 接收深度向量 (m) % freq: 计算频率 (Hz) write_env(freq, 200, ssp, envbase); system(['kraken ', envbase, ' > kraken.log']); system(['field ', envbase, ' > field.log']); % 读取声场数据 p = read_shd([envbase, '.shd'], rvec, zrec); % 避免 log10(0)，将极小值替换为量程下限 p_abs = abs(p); p_abs(p_abs < 1e-12) = 1e-12; TL = -20 * log10(p_abs); end调用时，rvec和zrec需要分别设定距离网格和深度网格，如果用等间距向量，绘图前需要做一次meshgrid。我一般把距离范围设为 1 到 10 公里，步长 50 米，接收深度从 5 米到 150 米，间隔 2 米。这里有个经验值：如果接收深度间隔太粗，后续画切片时会看到明显的条纹伪影；但间隔太细会让.shd文件体积暴涨，field 计算时间也跟着翻倍。
3.3 切片图里的模态干涉条纹怎么看
距离-深度切片画出后，你会看到典型的模态干涉条纹。这些条纹不是噪声，而是不同简正波之间相位叠加的结果。条纹的疏密主要由声速剖面的梯度决定，负梯度环境下条纹往往更密，因为高阶模态更容易被束缚在声道内。看到条纹后，第一步不是调色标，而是检查几条特征曲线是否平滑：如果在某个深度上出现突然断裂的边界，往往是接收深度超出了模态函数的有效覆盖范围，而不是物理现象。这个判断标准我后面还会用到，它在阵元波形计算中更关键。
4. 把声场转成各阵元接收波形：频域合成与逆傅里叶变换
4.1 为什么 Kraken 默认给不了时域波形
Kraken 在单个频率上输出的复声压，本质上是环境对点源的格林函数采样。你要的是宽带信号通过水声信道后的时域波形，那就得把一个频带内所有频率点的复响应都算出来，再乘以源信号的频谱，最后逆傅里叶变换。这个原理不难，但操作上有个容易被忽视的约束：频率分辨率必须足够密，否则时域波形会出现时间混叠。Kraken 本身不含自动扫频模块，所以这个铺频点的工作只能在 Matlab 里完成。
4.2 逐频点计算复声压并合成阵元波形
下面这段代码是我的常用方案，针对一条八元垂直线阵，阵元深度从 20 米到 34 米，间距 2 米，频带从 50 Hz 到 200 Hz，频率步长 0.5 Hz。每次循环只更新环境文件里的频率，然后调用 Kraken 和 field，最后读取该频点下所有阵元的复声压，存进复数矩阵。
% 参数定义 fs = 1000; % 采样率 (Hz) T_dur = 2; % 信号时长 (s) t = (0:fix(fs*T_dur)-1)' / fs; % 时间向量 % 源信号：简单脉冲 fc = 100; % 中心频率 bw = 30; % 带宽 src = exp(-((t - 0.5) / 0.2).^2) .* sin(2*pi*fc*t); % 频率轴 fvec = 50 : 0.5 : 200; % 频带与步长 z_array = 20 : 2 : 34; % 阵元深度 [rvec, zrecv] = meshgrid(1000, z_array); % 距离 1 km，深度取阵元位置 % 预分配复声压矩阵 H = zeros(length(fvec), length(z_array)); for jf = 1:length(fvec) write_env(fvec(jf), 200, ssp, 'broadband'); system(['kraken broadband > /dev/null']); system(['field broadband > /dev/null']); % 读取当前频点的复声压 p_f = read_shd('broadband.shd', 1000, z_array); H(jf, :) = p_f(:).'; end这段代码里，H的每一行代表一个频率点，每一列代表一个阵元深度。读出来之后先别急着画图，检查read_shd返回的数据是否包含正确的相位信息。Kraken 输出的复声压默认带相位，但某些封装程序在读取时只保留了幅度，导致后续逆傅里叶变换丢失所有时延信息，波形会变成完全对称的零相位脉冲。这个错误非常隐蔽，因为幅度谱看起来完全正常，但时域波形完全不是预期的样子。
4.3 逆傅里叶变换：时间窗和频域补零的取舍
得到频域响应H之后，源信号的频谱要左乘传输函数。为了保证时域波形长度足够长，不会因循环卷积产生首尾混叠，我习惯对频域数据做补零，把频域长度扩展到与源信号 FFT 长度一致。具体操作如下：
% 源信号频谱 NFFT = length(src); S = fft(src, NFFT); f_fft = (0:NFFT-1)' / NFFT * fs; % 插值传输函数到 FFT 频点 H_interp = interp1(fvec, H, f_fft(1:length(fvec)), 'linear', 0); % 注意矩阵维度：H_interp 是 NFFT × N 向量，需要按频率展开 H_full = zeros(NFFT, length(z_array)); for k = 1:length(fvec) H_full(k, :) = H(k, :); end % 频率响应与信号频谱相乘 Y = H_full .* S; % 逆傅里叶变换并取实部 r_wave = real(ifft(Y, NFFT, 1));这里有几个参数要特别解释。interp1的插值把稀疏频率点的传输函数插到 FFT 的密集频点上，如果频带外的插值全置 0，实际上相当于做了一个频带限制滤波器，带外噪声被直接滤掉。H_full在k=1:length(fvec)范围内填充，其余频率点保持 0，这样做的好处是源信号在频带外的能量不会被错误放大。补零后ifft得到的时域波形长度是 2 秒，但有效信号只有脉冲附近的部分，末尾会出现一些低幅度的尾响应，这是多径到达的体现，不是数值误差。
4.4 阵元间波形差异与信道的“指纹”
八元阵输出波形后，你会看到不同阵元之间的到达时刻存在微小偏移，这本质上就是不同深度的信道冲激响应差异。如果你只关心传播损耗，可能意识不到这些相位差的存在，但一旦进入波束形成阶段，相位差就是决定阵列增益的关键。所以建议每次合成完阵元波形，先把相邻阵元的互相关峰值时间差打印出来，检查是否在同一直线阵几何关系可解释的范围内。如果时间差跳变超过一个周期，大概率是read_shd读取的深度索引错位，而不是信道本身的问题。
5. 常见问题与排查：Kraken 与 Matlab 结合时的典型坑
5.1 模数少到离谱：频率和声速剖面不匹配
现象：日志文件里显示的模态数只有 1 或 2，传播损耗曲线在远距离处异常平滑，衰减速率明显快于球面扩展。
原因：Kraken 自动选择的模态数与频率直接相关，当你给出的频率低于该环境下基阶模态的截止频率时，能够传播的模态数就极少。另一个常见原因是声速剖面里存在明显的负梯度，但负梯度层深度过浅，导致能量被束缚在表面附近，深部模态无法激发。
解决：先提高频率测试，看模态数是否上升。如果频率提高后模态数显著增加，说明原始频率确实接近截止频率。如果是负梯度的问题，需要检查声速剖面数据是否来源于实测 CTD 插值，插值方法不当可能把细微梯度放大。我处理模拟项目 X 时遇到过这种情况，最后把声速剖面插值深度增加了一层，问题就消失了。
5.2 读取到的阵列深度尺寸与预期不一致
现象：read_shd返回的声压矩阵维度不是预期值，比如你指定了 8 个阵元，返回的却是 7 行或 9 行。
原因：.shd文件里的接收深度向量由 field 程序的输入参数决定，而不是环境文件里的Nrd。环境文件里的Nrd=0只是为了占位，实际接收深度必须通过 field 程序的参数文件或命令行指定。如果你在环境文件里把Nrd写成了 8，而 field 程序实际只用了 5 个接收深度，就会出现维度不匹配。
解决：在调用field之前，先确认现场参数文件里指定的接收深度向量与你的z_array完全一致。一个稳妥的做法是从.shd文件中直接读取实际的接收深度列表，再和你 Matlab 里的z_array做一次去重检查。
5.3 逆傅里叶变换后出现剧烈的前置脉冲
现象：合成出的阵元波形在信号到达前就已经有大幅度的振荡，而且振荡幅度与主信号相当。
原因：这是频域补零不当导致的。如果你把H_full中超出频带的部分全部置零，相当于在频谱上做了一个矩形窗，矩形窗的时域旁瓣会产生剧烈的振荡。更隐蔽的原因是源信号中心频率fc与 Kraken 计算的频带中心不一致，导致频带边缘出现陡峭的幅度跳变。
解决：在频带边缘做 5% 的过渡带衰减，而不是直接硬截断。我常用的办法是用一个汉宁窗对H_full的幅频响应做平滑过渡，只在过渡带内施加衰减，频带内保持原始值。这样能有效压低前置脉冲，波形看起来也干净得多。
5.4 海底衰减参数单位不一致导致损耗突变
现象：某个频段内传播损耗突然出现不连续台阶，而且台阶位置随频率移动。
原因：Kraken 环境文件中的海底衰减系数单位有dB/λ和dB/m两种表示方式，如果输入时没有留意单位，Kraken 内部会按dB/λ解析，结果低频段衰减被严重放大，高频段反而正常。这会让扫频计算出的损耗曲线出现频率依赖性异常。
解决：统一按dB/λ或统一按dB/m做换算，换算关系是α_dB_per_m = α_dB_per_wavelength × f / c_bottom，其中c_bottom是海底声速。我一般会在生成环境文件的函数里把单位换算写死，这样即使换一个人来跑他的数据也不会搞混。
6. 最后补一个实用习惯：核对阵元深度与模态覆盖范围
在阵列波形合成这条链路里，我最看重的是验证。这里分享一个固定动作：在正式批量计算之前，先用单频点算一条垂直模态分布图，检查你刚才指定的阵元深度是否落在模态函数的有效覆盖范围内。具体做法是在环境文件不变的前提下，把接收深度设为从上到下连续变化的向量，跑一次 field 程序，然后把各深度的复声压幅度画成深度方向的曲线。如果某个深度的声压幅度接近零，说明这个位置处于模态函数的零点附近，你在这里布阵元，算出来的波形信噪比会非常差。
这个验证动作不是可有可无。以前我直接拿一条 20 到 34 米的阵元深度去批量扫频，结果波束形成出来的目标方向图完全没有主瓣，排查了很久，才发现这段深度恰好落在某一阶模态的零点带附近，所有阵元接收到的都是微弱信号。把阵元深度避开零点带之后，方向图立刻恢复正常。从那以后，我每次做阵列波形仿真之前，都会强制走一遍“单频点垂直声压扫描”，先看深度分布再定阵元位置，这个习惯帮我省下了一整周的返工时间。希望这个习惯也能帮到你。
本文还有配套的精品资源，点击获取
