简介Kraken水声传播计算资源面向海洋声学研究与水声通信开发者基于Matlab实现传播损耗计算和水平、垂直阵元接收波形模拟解决声场建模与阵列性能评估需求。整套资源共19个文件以Matlab脚本.m、仿真结果图.jpg/.png、公式文档.docx、声速剖面数据.mat和说明文档.md为主压缩包整体仅713KB携带分享都很方便。目前已有1108人浏览学习常被用作水声信号处理与通信系统设计的参考实现。代码覆盖垂直阵与水平阵两种布阵方式可计算传播损耗随距离、深度的变化并绘制声场与接收波形图附带的公式文档和任务说明有助于理解Kraken边界模型及算法逻辑。通过修改声源频率、环境参数或阵列配置可探索不同水文条件下的传播特性也适合水声物理、阵列信号处理等方向的课程设计与科研入门。1. 用Kraken计算传播损耗和阵元波形先确认它适不适合你的场景大约半年前我在某个浅海项目里对接过一次实测数据声速剖面基本水平换能器阵元排成一条直线但我拿射线模型算出来的传播损耗和测量差了接近10 dB。后来换成了Kraken这套基于简正波理论的Matlab工具包修改环境文件和接收深度设置后传播损耗曲线才跟实测对得上。Kraken的核心思路是把声场分解为一组简正波再按深度和距离叠加它在浅海、低频、水平分层介质中比射线模型稳定得多。它能给你每个频点的复声压、传播损耗也能通过传递函数继续合成自定义频带下的阵元接收波形。这份资源适合做水声信道仿真、阵列信号处理前端、目标深度估计的同学但如果你要处理高频复杂海底地形它就不太合适建议先认清这个边界。2. 环境与输入文件准备跑通一次正演要过的三关很多第一次拿到Kraken工具包的人会卡在第一步因为Kraken本体是Fortran写的一堆可执行文件Matlab只是外层驱动和数据后处理。这里按我的整理习惯把准备阶段拆成环境路径、环境文件、调用方式三个部分全部走一遍以后后面算传播损耗和阵元波形就有共同的工作基础了。2.1 工具箱路径接进Matlab下载下来的Kraken资源包解压以后首先确认里面包含哪些东西。常见的目录结构里会有一个Kraken子目录放着可执行文件还有一个Matlab子目录装着field.m、plotshd.m、readshd.m这些函数。第一步是把这些路径加入Matlab搜索目录并验证Kraken可执行文件能不能被调用。我一般会写一个setup脚本避免每次重启Matlab都重新设置一遍% setup_kraken.m base fileparts(mfilename(fullpath)); addpath(fullfile(base, Kraken)); addpath(fullfile(base, Matlab)); % 确认可执行文件存在 if ispc exe fullfile(base, Kraken, kraken.exe); else exe fullfile(base, Kraken, kraken); end if exist(exe, file) disp(Kraken executable found); else error(Kraken executable missing, check installation); end % 保存路径下次打开自动生效 savepath;这段代码的逻辑很简单把两个子目录放进搜索路径然后检查对应平台的可执行文件是否存在。在Windows下Kraken工具包一般给出kraken.exe在Linux或macOS下是不带后缀的二进制文件。这里有个很多人忽略的坑如果Matlab安装在Windows下但工具包是从某高校Linux服务器拷贝过来的里面只有Linux二进制即使路径设置对了也跑不起来。所以确认exist(exe, file)为true是第一步。savepath这行建议保留但要注意Matlab可能因为权限问题无法写入pathdef.m这时会输出警告不影响本次运行只是下次要重新执行。路径设置完成后建议在Matlab里直接运行!kraken试试看终端是否输出版本或帮助信息。如果提示找不到命令检查当前工作目录是否在Kraken子目录下或者把子目录临时加入系统PATH。很多同学喜欢把Kraken的exe复制到Matlab当前目录这能跑通但不利于维护环境文件后面我会解释为什么环境文件应该统一管理在一个工作目录里。2.2 环境文件生成后续所有计算的源头Kraken运行起来不像Matlab那样直接在命令行传参数它要读取一个外部的环境文件也就是.env文件。这个文件保存了声速剖面、海底参数、频率范围、声源深度、接收深度等一切物理设置。环境文件的第一行通常是标题行随便写点备注文字第二行开始是频率和运行标志。比如我常用这样一个浅海环境文件Shallow Water Test 50 0 0 0 SVP 0.0 1500.0 50.0 1500.0 SSP 0.0 0.0 50.0 1.0 BOT 1500.0 1.7 0.5这里逐段解释第一行是注释标题必须用单引号括起来Kraken会忽略它但要求它存在。第二行的第一个数50表示计算频率为50 Hz后面三个整数分别是声源深度标志、接收深度标志、上方半空间标志一般写0表示使用默认行为某些版本还会把Nx和NxM等参数放进来具体要看资源包里的样本文件格式。第三行开始是介质描述SVP后跟着声速剖面每行是一个深度值和对应的声速这里表示0米水深处声速1500 m/s50米水深处也是1500 m/s也就是等声速层。然后是海底层我用BOT关键字后面接海底纵波速度1500 m/s、密度1.7 g/cm³、衰减系数0.5 dB/波长。这个环境文件格式和版本的对应关系是新手最容易翻车的地方。工具包里如果自带样例环境文件比如类似balck.env或者sbcx.env这样的演示文件强烈建议先打开样例确认每一行的含义再修改成自己的声速剖面。因为我见过某一次在资源包中解压出来的环境文件模板使用了SSP关键字而某文档版本里用的是SVP差一个字母就会导致程序识别失败。遇到这种情况不要panic去工具包的样例里复制最接近的结构不要自己从零写。2.3 在Matlab里调用而不是去命令行折腾环境文件准备好以后常见做法是在Matlab里用system调用Kraken可执行文件。我一般把环境文件和工作目录放在一起运行前用cd切换过去% run_kraken_case.m workdir fullfile(pwd, case01); if ~exist(workdir, dir) mkdir(workdir); end % 将环境文件和Matlab脚本放在同一目录 copyfile(shallow.env, workdir); cd(workdir); % 调用Kraken生成声场文件 if ispc status system(..\Kraken\kraken.exe shallow.env); else status system(../Kraken/kraken shallow.env); end if status ~ 0 error(Kraken run failed, check .env file); end % Kraken输出的声场文件默认是 shallow.mod用于后续field计算 disp(Kraken completed);这段代码解决了两个问题一是避免手工在终端敲命令二是把环境文件集中管理方便后续批量修改参数。status返回值很关键很多运行失败在Matlab里不会直接抛异常只会在终端输出错误信息。所以一定要检查status如果非零优先打开环境文件检查关键字拼写和数字格式。Kraken对格式非常敏感比如某个数字后面多一个空格、某个图层行数不对都会导致解析中断。参数方面system调用里的命令格式是kraken 环境文件名环境文件名不带引号但路径中有空格时需要格外小心Windows下会用短路径或者用双引号包裹整个命令。我一般会把工作目录保持为无空格路径这样能省掉很多转义问题。运行结束以后工作目录下会出现.mod文件这是后续field.m读取的核心文件。跑通这一遍以后才算真正建立起正演计算的工作流。下一步开始接触传播损耗的具体计算这部分才是资源里最有价值的地方。3. 传播损耗计算从声压到TL曲线的关键参数环境文件和模型程序都跑通后计算传播损耗是一件相对轻松但细节很多的事。Kraken输出的是一堆频域复声压我们要把它转换成工程上通用的传播损耗曲线同时还要理解深度、距离、频率这些参数到底怎么影响结果。3.1 传播损耗的定义与单位约定传播损耗Transmission LossTL是水声学里最常用的量定义为某一点声强与离声源1米处参考声强的比值取20倍对数[ TL -20\log_{10}\left(\frac{|p(r,z)|}{|p_0|}\right) ]在Kraken的默认约定中声源强度归一化因此field.m返回的声压绝对值可以直接代进这个公式。但这里有个单位陷阱某些版本的field.m返回的是复声压Pa某些版本的帮助文件里说是相对值实际测试后我发现按上面的公式计算出的TL在近距离约等于20log10(1/r)这是球面扩散的正常结果。如果算出来在近距离是0 dB或负值多半是参考声压归一化方式不同这时不要硬改公式先检查Kraken版本对应的文档。距离和深度参数的概念也要提前理清距离指水平距离单位是米深度指接收点水深单位也是米声源深度是单个值。Kraken在计算时是在距离-深度网格上给出复声压所以要先定义我们要在哪些距离点和深度点上取值。3.2 用field.m读取声场并绘制TL剖面% compute_TL.m % 先运行Kraken得到shallow.mod后再执行本脚本 % 读取环境文件中的频率和介质信息 freq 50; % 与.env文件中一致 sd [25]; % 声源深度单位m rd linspace(0, 50, 51); % 接收深度单位m rs linspace(0, 2000, 201); % 水平距离单位m % 调用field计算复声压 [pressure, freqOut] field(shallow, freq, sd, rd, rs, GS); % 计算传播损耗 TL -20 * log10(abs(pressure) eps); % 绘制深度-距离二维TL图 imagesc(rs, rd, TL); xlabel(Range (m)); ylabel(Depth (m)); title(Transmission Loss at 50 Hz); colorbar; shading interp;field是Kraken工具包里最核心的函数。它的第一个参数是环境文件重名输入但注意不需要扩展名比如环境文件是shallow.env这里就写shallow。第二个参数是频率必须和环境文件一致如果不一致很多版本不会报错而是自动重新运行Kraken这会导致时间浪费。第三个参数sd是声源深度可以是一个标量第四个参数rd是接收深度可以是数组第五个参数rs是距离采样也是数组。最后一个选项GS表示输出Greens function格式这是声压的常规定义。计算TL时加eps是为了防止log10(0)导致负无穷。这段代码输出的二维矩阵第一维对应rd第二维对应rs绘制imagesc时旋转一下坐标轴让距离在横轴。你会发现TL在近距离呈现明显的干涉条纹这是两种简正波相干叠加形成的如果条纹看起来太密或者太宽直接怀疑接收深度网格和距离网格的分辨率。3.3 关键参数对结果的影响计算传播损耗时三个参数几乎决定了结果的所有形态。第一个是接收深度rd的采样间隔。在浅海环境中某一阶简正波的垂直波长大约在几米这个量级所以rd间隔至少要小于最短垂直波长的一半否则深度方向会出现混叠。比如50 Hz声波在1500 m/s水中的波长为30米垂直模态的尺度通常比30米小我习惯用1米或者更细的网格采样代价是矩阵变大、计算时间变长。第二个是距离rs的最大值和步长。Kraken的传播损耗在近距离有强烈的球面扩散和模态干涉在远距离趋于柱面扩散。所选取的最大距离要根据实际声场确定一般来说在浅海环境中取水深的10到20倍距离能看到明显的模态衰减特征。距离步长影响横轴的平滑程度如果步长小于十分之一波长会算得慢但曲线光滑步长太粗则干涉条纹会被直接跳过看起来像被低通滤波了。第三个是频率freq的选择。Kraken的简正波模型在低频表现最好一般低于几百赫兹。如果你要算20 kHz以上的高频声场简正波数量会膨胀到几百阶计算量剧增且不稳定。这里有个实用判断标准看输出TL曲线是否有连续性如果相邻距离点的TL跳动超过20 dB大概率是模态数不够此时应在环境文件里增加最大模态数或在调用field时让程序自动调整而不是继续用默认设置。传播损耗这一块跑通后才能进入阵元接收波形的计算。阵元波形不是简单读一下TL就行需要把频域传递函数拿去做时域合成这一步是很多人从算场到阵列信号处理的转折点。4. 阵元接收波形合成频域传递函数到时域卷积算完传播损耗后我们拿到的是一张深度-距离上的TL图。如果要模拟一个阵列各阵元实际接收到的时域波形你需要做一次从频域到时域的搬运Kraken给出的每个频点的复声压本质上是信道的频域传递函数H(f)。把发射信号的频谱乘以H(f)再逆傅里叶变换回时域就得到该阵元的接收波形。4.1 为什么不能直接用TL推波形TL是标量强度信息它丢掉了相位。而水声信道中相位包含多径传播的时延信息两个阵元之间波形错位几毫秒可能就对应着几十米的路径差丢掉相位就无法保留。因此阵元波形合成必须用复声压而不是TL。Kraken的field函数正好返回复数压力场这就保留了相位信息。实现时要注意参考距离归一化不同版本的Kraken可能把声源强归一化在不同的距离上这会造成波形幅值缩放但对波形形状和相对时延影响很小。4.2 逐阵元提取频域传递函数假设阵列共有N个阵元每个阵元位于深度数组element_depths中的不同位置水平距离分别是element_ranges。我的做法是把每个阵元的接收深度和距离组成独立的参数对分别调用field提取对应点的复声压构成这个阵元的频域传递函数H_n(f)。% compute_element_H.m % 假设工作目录下shallow.mod已经生成 freqs 30:0.5:80; % 关心的频带注意预留过采样空间 N 8; element_depths linspace(10, 40, N); element_ranges ones(1, N) * 500; % 水平距离500m H zeros(N, length(freqs)); for n 1:N % 每个阵元单独提取复声压 [p, ~] field(shallow, freqs, 25, element_depths(n), element_ranges(n), GS); H(n, :) p(:).; end这里频率数组freqs我取了30到80 Hz步长0.5 Hz。步长不是越大越好0.5 Hz对应最大可观测时间窗口为2秒即1/0.52秒如果多径时延超过2秒则会产生时间混叠。所以需要先估算多径最长传播路径差用路径差除以声速得到最大时延再反推频率步长。如果最大时延是1秒频率步长必须小于1 Hz我习惯留2倍余量取0.5 Hz。调用field时freqs直接传向量这比逐频点循环快得多。field内部的实现会对每个频率分别计算简正波频率点越多计算越慢同时内存占用线性增长。8个阵元和101个频点这个规模大概几百毫秒可以接受。真正慢的是环境文件里楼层多、模态多的情况。4.3 合成时域接收波形的代码拿到H_n(f)后下一步是设计发射信号频谱S(f)再相乘后做IFFT。发射信号我常用线性调频信号或者衰减正弦脉冲下面是LFM信号合成接收波形的完整示例% synth_waveform.m fs 1000; % 时域采样率Hz T 1.2; % 信号时长s t 0 : 1/fs : T - 1/fs; f0 35; f1 75; % LFM起止频率Hz s chirp(t, f0, T, f1, linear); % 把信号变换到频域频点与H对齐 S fft(s, length(freqs)); S_pad interp1(linspace(0, fs, length(S)), S, freqs, linear, 0); % 每个阵元的接收频谱 Y H .* S_pad; % IFFT回到时域 y ifft(Y, length(freqs), 2); % 绘制第1和第4阵元的波形 subplot(2,1,1); plot(t, real(y(1,:))); title(Element 1 received waveform); subplot(2,1,2); plot(t, real(y(4,:))); title(Element 4 received waveform);这段代码里的核心是把信号频谱映射到H的频率轴上。fft得到的是频率范围从0到fs的频谱而freqs只截取了30到80 Hz一段所以用interp1做线性插值。插值外的分量补0等效于带通滤波。ifft在第二个维度上执行得到每个阵元的时域序列长度等于length(freqs)101这不像一个完整的时域信号因为频率步长限制导致时间轴很短。要得到完整的波形需要频率点覆盖整个信号带宽且频率间隔足够小。更稳妥的做法是把频率数组设置成覆盖信号带宽且同时满足fs和时长要求。比如fs1000T1.2秒则频率轴长度应为fs*T1200点对应的频率分辨率为1/T≈0.833 Hz。所以freqs应该设置为从0到fs的1200个频点而不是我之前例子里的30:0.5:80。field接受0频率可能出问题实际处理时把freqs的视频范围设为0到fs但只计算非零频点0频点直接置0即可。工程上我会先把freqs定义成完整频率轴然后对低频段调用field计算高频段用对数插值或者直接补零这样能兼顾计算速度和精度。4.4 频率轴设计是波形正确与否的分水岭很多人在这一步翻车是因为没有理解频域采样和时域窗口的关系。Kraken本身是频域模型field返回的频率点是离散的频率点间隔决定了波形最大可分辨时间长度。假如只算了5个频点IFFT出来的时域波形会被极大的周期化误差污染看起来像多个脉冲叠加其实只是混叠。我一般会用下面的规则确定频率数组频率范围要覆盖发射信号的主要能量带宽频率步长要满足最大时延分辨频率点数即IFFT长度。以fs2000、T2秒为准频率轴是0到2000步长0.5点数2000。field对2000个频点逐个计算每个频点生成一组模态这个计算量在低频浅海场景是可接受的大约几十秒。如果阵元数量多、频率点又密可以考虑先算出一个深度-距离网格上的场再用插值获取阵元位置的值避免重复计算。这一阶段的成功标志是某个阵元波形能明显看到直达波和一次海面反射波到达两个脉冲时延差和几何计算一致。看到这个结果说明Kraken的相位信息和你的时域合成流程都对上了。5. 避坑清单最容易让结果翻车的5个细节Kraken用起来说简单也简单但如果你没有踩过它的几个典型坑某些结果会看起来很合理却经不起推敲。以下就是我实际运行中遇到最多的问题每一条都按现象、原因、解决三步记录。5.1 距离方向出现高密度干涉条纹远看像一团乱麻现象绘制的TL深度-距离图中干涉条纹在近距离密集到无法分辨远距离条纹变粗但中间有大量不连续跳变。原因接收深度或距离网格的分辨率不够导致某些高阶模态在网格之间未被正确取样。解决把rd网格间隔从2米改到1米rs的步长同时改小一般能消除毛刺。如果毛刺仍然存在考虑在field调用中增加模态计算数量或在环境文件中把最大模态数写大一点。5.2 海底参数用的密度值导致TL整体偏移现象算出来的全场TL比相同条件下的解析解大或者小6 dB左右曲线形状基本正常。原因海底密度单位错位。Kraken默认海底密度单位是g/cm³也就是相对水的密度比比如1.7表示1.7倍水密度。如果你误填成kg/m³的1700程序会按1.7进行归一化但实际计算中密度比变为1700阻抗失配严重能量折射方式改变。解决检查环境文件中的BOT段落确保密度在1.0到2.5之间同时确认衰减系数0.5是dB/波长而不是dB/m。5.3 高频条件下运行时间剧增且结果不收敛现象freq从200 Hz提高到500 Hz计算时间从几秒变成几分钟而且TL曲线在不同距离上跳动剧烈。原因简正波数量随频率线性增长高频时总模态数超过程序默认上限截断误差变大。解决这不是Kraken脚本的bug而是模型本身的适用范围限制。改用更高频的射线模型或减小环境文件中海水的层数以降低模态数还可以尝试在field中设置更小的深度网格。如果必须在高频下用简正波检查环境文件里的最大模态数设置通常可以通过增加该值提高精度但计算时间和内存消耗直线上升。5.4 阵元波形出现周期性假目标现象合成的阵元接收波形中出现规律性重复的脉冲间隔恒定像是一串回波但实际环境不可能产生这么规律的反射。原因频域采样率过低或频点数不足导致IFFT后时域发生周期延拓假脉冲来自相邻周期的泄漏。解决把频率步长减小使最大可观测时间超过环境多径最大时延的2倍增加频点覆盖范围使IFFT点数足够包含全部能量。另一个常见因素是发射信号上采样不足把发射信号fs设置成信号带宽的2倍以上避免在频域插值时产生伪峰。5.5 环境文件改了声速剖面后程序直接报错退出现象修改某个样例环境文件的声速剖面数据保存后运行Kraken终端输出End of file reached while reading layer 2程序终止。原因声速剖面的数据行数与文件头声明的层数不一致或者关键字SSP和SVP用错。解决打开环境文件数一遍声速数据行数确保与声明中的层数匹配修改声速时直接覆盖不要删除中间行避免破坏格式。这里有个小技巧每次修改环境文件后先复制一份保留原样再在副本上改动出了问题可以立即对比差异不用靠记忆回退。以上这5条覆盖了环境、网格、单位、频率轴、格式五类最常见的雷区。如果你算出来的结果看起来像那么回事但就是和预期对不上先对照这个清单检查一遍往往比调参数更有效。6. 进阶校验技巧用解析解给Kraken结果兜底Kraken算出来的东西再漂亮也要有个独立验证手段。我的习惯是构造一个Pekeris波导的解析解与Kraken的数值结果做对比。Pekeris模型是均匀水层加均匀半空间海底它的传播损失可以用简正波解析公式直接计算不需要跑数值软件。正因为在同一个物理条件下两边应该一致这个方法能揭露出环境文件里细微的设置错误。构造解析解的方式并不复杂。以水深50米、水层声速1500 m/s、海底声速1600 m/s、密度1.8 g/cm³为例在50 Hz频率下前几阶简正波的水平波数可以写为特征方程的解。我只在距离较远、单模态占优的地方对比TL斜率。若Kraken给出的远距离TL曲线斜率接近每百米衰减3到5 dB同时绝对值和解析解相差在1 dB以内说明模型参数可信。这个对比过程只需几十行Matlab代码但它让我避开了很多次因单位错误导致的整体偏移。批处理阵元波形时还有两个实用的工程技巧。第一优先把field计算得到的复声压网格保存成Mat文件后续所有阵元波形都从Mat文件里插值提取而不是重新跑Kraken。这样当你需要修改发射信号参数时不需要重新计算声场省掉大量重复时间。第二当阵元数超过16、频点超过1000时内存消耗会显著上升此时尽量在单次field调用中返回整个深度-距离网格然后用interp2提取阵元位置的值避免在循环内重复调用field。我实测过一个24阵元、800频点的场景循环调用field需要8秒改成网格插值后压缩到0.2秒。从那以后我每次拿到一个新的Kraken环境文件都强制自己先做一遍Pekeris解析解比对再进入阵元波形合成。这个习惯让很多潜在的单位和网格问题在早期就暴露而不是等波形成形后才回头找原因。Kraken工具包的功能边界和这些校验手段配合使用才能让计算传播损耗和阵元接收波形这件事做得可复现、可信赖希望帮到你。本文还有配套的精品资源点击获取