★ FEATURED ARTICLE 简介Kraken是一款基于Matlab的水声传播计算工具面向水声通信和海洋声学领域的科研人员与工程师用于计算声波在水中的传播损耗并分析水平阵、垂直阵各阵元的接收波形及阵列方向图帮助评估和优化声纳探测性能。压缩包共19个文件大小约713KB内含3个Matlab源程序、12张仿真结果图、1个海洋环境参数数据文件以及README说明、公式文档和任务示意图片覆盖从环境参数设置、传播损耗计算到阵元波形绘制的完整流程。目前已有1108人学习浏览。源码带有详细注释便于理解算法逻辑和自定义修改结果图直观呈现声场分布、传播损耗随距离和深度的变化以及不同阵列配置下的接收波形配套文档则提供了使用步骤与注意事项适合水声工程初学者快速上手也可供研究人员调整参数开展环境适应性实验。1. Kraken到底在算什么一个被当成“黑匣子”的简正波工具箱拿到一个新海域的声学环境第一件让我心里踏实的事就是用 Kraken 把传播损耗和各阵元接收波形先算出来。很多人一上来就调波束形成、匹配场算法结果阵列增益、检测概率全是虚的——因为信道没摸清。Kraken 是水声领域常用的简正波声场计算程序它把海水看作水平分层波导用模态叠加近似求解声传播输出的是复声压。这个复声压既能画成传播损耗TL图也能通过多频点合成还原成每个阵元上的时域接收波形。适合做水声通信、目标被动定位、阵列信号处理的从业者尤其刚接手浅海环境的人用它打底比直接上射线模型更稳。它的核心价值不是“算个数字”而是把信道的关键物理量——模态数、相速度、群速度、深度依赖的声压分布——一次给清楚。你后面写波形仿真也好做匹配场拷贝向量也好全都从这一套输出里来。下面我按自己常用的落地路径从环境文件写到阵元波形合成把参数和坑一次说透。2. 先搭环境文件这部分错了后面算什么都白搭2.1 环境文件的骨架频率、SSP、海底那几行必须配对我在某公司做模拟项目X时第一次跑Kraken就是被环境文件格式卡住的。它不像常规配置文件那样有括号、有注释而是按行列排列的纯文本少一行、多一个空格都可能让程序直接退出。我一般会先写一个最小可用的模板把骨架立住再往里填参数。KRAKEN ! 运行类型KRAKEN表示普通简正波 50.0 ! 计算频率(Hz) 1 ! 介质层数 0.0 100.0 ! 层顶深度(m)和层底深度(m) 1500.0 0.0 0.0 ! SSP点声速(m/s)、纵波吸收(dB/λ)、密度(g/cm3) 1520.0 0.0 0.0 ! 第二个SSP点声速随深度变化 100.0 1.6 0.8 ! 海底半空间声速、密度、衰减(dB/λ) 0.0 ! 海面粗糙度标准差(m)平滑海面填0这里每行位置都有讲究。第一行“KRAKEN”指定用简正波模型而不是抛物方程或射线模型第二行的频率决定了后续模态展开的截断阶数第三行介质层数告诉程序往下读几层SSP。关键坑在于SSP点的深度范围必须和层顶、层底深度完全一致。上面例子里层底是100米SSP最后一个点的深度也需要落在100米否则程序会认为剖面没有闭合直接报错。海底半空间那行容易被忽略。它给的是半空间顶部的声速、密度和纵波衰减系数。密度单位是g/cm3衰减单位是dB/λ不是dB/m。很多人直接把实验室测的衰减值填进去换算关系没做导致TL曲线在高频段往下掉得异常快你还会以为是多途干涉造成的其实只是海底参数没配对。2.2 用分层声速剖面逼近真实温盐数据实际海域的声速剖面来自CTD或历史水文数据通常只有几十个深度点直接塞进Kraken也能跑但简正波方法对SSP的线性分段很敏感点太少会在不连续处产生虚假模态。我的做法是先做深度方向插值把剖面加密到每1~2米一个点然后平滑掉明显毛刺。import numpy as np from scipy.interpolate import interp1d depth_orig np.array([0, 10, 20, 50, 80, 100], dtypefloat) c_orig np.array([1500, 1498, 1495, 1490, 1485, 1488], dtypefloat) f_ssp interp1d(depth_orig, c_orig, kindlinear, bounds_errorFalse, fill_valueextrapolate) depth_new np.arange(0, 101, 1.0) c_new f_ssp(depth_new) # 核心平滑对声速做3点滑动平均去掉微尺度抖动 c_smooth np.convolve(c_new, np.ones(3)/3, modesame)插值之后再平滑是因为简正波方程的算子里包含了声速梯度的二阶项剖面局部震荡会激发一些物理上不应存在的伪模态。平滑时的窗口别开太大3点是常态深度100米左右水层用5点也问题不大超过5点会把真实的跃层结构磨平算出来的声传播特征就失真了。最后写入环境文件时我习惯把层顶深度写成0层底深度写成最大测量深度SSP点按深度升序排列。这里有个容易翻车的细节如果声速剖面数据是从水面到海底的但环境文件里层顶从10米开始那水面这一层就被程序当成了半空间近场传播损耗会乱掉你对着实测数据怎么也对不上。2.3 运行 Kraken 生成模态文件命令与常规检查环境文件准备好之后运行其实就一条命令。我一般在项目目录里新建一个子目录做仿真把环境文件命名为envfile方便后面field和Kraken默认去读。Kraken带不带后缀都能识别但最好别带中文路径和空格容易出现读取异常。kraken envfile跑完会在同目录生成一个envfile.mod文件里面存放了模态特征值、相速度、群速度和模态函数的离散采样值。紧接着用文本模式打开环境文件调高输出信息级观察模态截断是否合理。grep Mode envfile.prt我通常会看.prt里列出的模态数量。浅海50Hz110米水深一般会有20~40阶模态如果频率升高到500Hz模态数量会涨到上百阶。模态太少的直接后果就是后面field时TL曲线不光滑看起来像锯齿其实是高阶模态被截断了能量没有算进去。还有一个常规检查是看相速度谱如果相速度集中在某个窄带内说明环境比较理想如果相速度出现跳变优先检查SSP是否有异常梯度而不是盲目增加模态数。到这里环境文件这一步才算真正站住了。3. 把传播损耗从复声压里提出来field 程序与接收网格3.1 flp 文件怎么定接收网格深度步进、距离步进、相干与非相干kraken只负责算模态真正把它铺成空间声场的是field程序。field 需要读取一个接收场文件一般叫flp里面写接收深度、接收距离和输出选项。这个文件的格式同样很死板写错一个数字就会导致输出全为NaN。20 ! 接收深度点数 5.0 100.0 ! 接收深度范围(m) 50 ! 接收距离点数 10.0 5000.0 ! 接收距离范围(m) A ! 输出选项A相干B非相干C只用幅度接收深度范围最好落在水层内不要跑到海底半空间里去半空间内声压通常不是简正波展开的覆盖范围。距离网格从10米开始是因为点源近场在1米内会出现奇点TL在那儿失去意义打印成图就是一根向负无穷掉的线。用10米作为起点既保留了近场结构又不会让奇点污染色标范围。选项里A表示相干叠加即所有模态按复声压求和保留了相位信息B表示非相干叠加相位被丢掉只按能量求和。匹配场、波形合成必须用A能量探测、平均传播损失分析用B就够了。这一步选错后面波形永远不对因为你在源头就把相位抹掉了。3.2 从 pressure 文件读复声压并画 TL 图field跑完之后除了.prt里的TL列表还会生成一个.pressure文件。这个文件在常规配置下是ASCII文本按“接收深度优先、距离次之”的顺序存放复声压的实部和虚部。我习惯直接用一个Python脚本读进来重新排成二维矩阵。import numpy as np import matplotlib.pyplot as plt data np.loadtxt(envfile.pressure) nrd 20 nrr 50 # 列格式深度(m)、距离(m)、实部、虚部(Pa1m参考) pressure data[:, 2] 1j * data[:, 3] P pressure.reshape(nrd, nrr) # 传播损耗20log10(|p|)参考1m处点源声压 TL -20.0 * np.log10(np.abs(P) 1e-12) depths data[:, 0][:nrd] ranges data[:, 1][::nrr] Z, R np.meshgrid(depths, ranges, indexingij) plt.figure(figsize(8, 5)) plt.pcolormesh(R, Z, TL, cmapviridis, shadingauto) plt.colorbar(labelTL (dB)) plt.ylabel(Depth (m)) plt.xlabel(Range (m)) plt.title(Transmission Loss computed by Kraken)这里有个新手经常踩的坑np.loadtxt读进来的行顺序是深度循环在里还是外不同版本和不同field编译选项下不一样。我的经验是先打印前几行的深度列看它是不是从5米跳到下一个深度再决定reshape的维度顺序而不是盲信某个固定格式。上面代码按“深度在外层、距离在内层”排和field常规输出一致但如果你的.prt表头写的是深度循环在内就把reshape参数反过来TL图会立刻合理很多。3.3 单频TL曲线的三处坑近场奇点、去归一化、dB参考很多人拿到pressure文件就开始画TL画完却和文献对不上。我遇到过三个反复出现的坑。第一个是近场奇点。距离小于1米时Kraken的简正波解包含点源项的奇异性20log10(|p|)会急剧向下这不是真实物理而是模型本身在近场的退化。处理办法是绘图和后续分析都从10米之后取数近场部分直接丢掉。第二个是去归一化。Kraken输出的声压默认是1m参考距离的点源辐射声压归一化后的结果也就是说当距离刚好1米时pressure接近1.0幅度为0dB。若你想得到绝对声压级就必须额外加上声源级SL比如140dB re 1μPa然后在频域合成时叠加到对应频点。漏掉这一步后面算接收波形会整体偏弱几十dB信噪比完全失真。第三个是dB参考混淆。传播损耗通常定义为20log10(|p|)功率量纲才是10log10(|p|^2)。两者数值等价但如果你在代码里将幅度乘以2再取对数就变成功率损耗除以2结果图会整体下移3dB的倍数看起来像多途衰落实际只是对数公式写错了。我一般只看20log10(|p|)以pressure文件的幅度为准避免二次换算带来额外的乘除错误。4. 从单频声压到各阵元接收波形宽带合成的完整链路4.1 为什么不能直接拿 TL 当波形相位信息全没了TL图只给了幅度而水下信道的多途结构恰恰藏在相位里。同一个接收点直达声和海面反射会形成干涉条纹相位差异决定了对齐增强还是抵消。直接用TL拟合一个衰减包络去模拟接收信号你等于把干涉、时延、模式频散全丢了做出来的阵元波形会非常“干净”——但这反而说明它是错的。真实接收波形应该是幅度起伏剧烈、带有明显尾部的瞬态响应。所以想要得到各阵元接收波形正确链路是对频带内每个频点用Kraken算出该频点在该阵元位置的复声压含幅度和相位然后用这个复声压充当信道传递函数H(f)与发射信号频谱S(f)相乘最后做逆傅里叶变换回到时域。听上去不复杂但真正落地时涉及采样率、频点数和时间窗的配套设计下面逐个说。4.2 建立阵元坐标垂直阵、水平阵、倾斜阵怎么排阵元波形合成前先把阵元的几何坐标用脚本生成好不要手写。常见的接收阵有两种垂直阵和水平阵。垂直阵关注深度维度的模态干涉结构阵元间距通常取关心频带内最高频率对应波长的一半水平阵则沿距离方向排列更关心水平方向的波束指向性。import numpy as np f_max 300.0 # 频带上限Hz c0 1500.0 # 参考声速m/s wavelength c0 / f_max spacing wavelength / 2.0 # 半波长布阵 n_elements 8 depth_array 100.0 - np.arange(n_elements) * spacing # 垂直阵从海底往水面排 range_array np.zeros(n_elements) # 阵元坐标表后续写入field的接收深度/距离范围 for i in range(n_elements): print(f{i1:02d} depth{depth_array[i]:.2f}m range{range_array[i]:.2f}m)如果做水平阵保持深度不变距离方向按等间距递增如果做倾斜阵就是深度和距离同时变化。关键在于场计算程序通常按规则的深度-距离网格输出阵元位置不一定恰好落在网格线上。我一般把Kraken的接收网格设得足够密再用双线性插值把阵元处的复声压取出来而不是直接取最近网格点。最近邻取值的误差在高频时会被相位抖动放大合成波形会出现人为的频移。4.3 多频点 IFFT 合成时域波形采样率、频点数和时间窗怎么定真正合成波形时我是对每个频点单独跑一次field或者一次批量读取然后按下面流程综合处理。import numpy as np # 频带设计 f1, f2 100.0, 300.0 # 频带上下限Hz delta_f 0.5 # 频点间隔Hz freqs np.arange(f1, f2 delta_f, delta_f) # 发射信号谱以线性调频为例 T_dur 0.1 # 信号时长s t np.arange(0, T_dur, 1.0 / (2 * f2)) k (f2 - f1) / T_dur tx_signal np.sin(2 * np.pi * (f1 * t 0.5 * k * t**2)) S_f np.fft.fft(tx_signal, n4096) S_f_full np.zeros(len(freqs), dtypecomplex) for ix, f in enumerate(freqs): idx int(round(f * 4096 / (2 * f2))) S_f_full[ix] S_f[idx] # 用前面算好的各频点复声压做信道传递 H_f np.zeros(len(freqs), dtypecomplex) for ix, f in enumerate(freqs): # 实际使用时从pressure文件取对应频点的插值结果 H_f[ix] get_pressure_at(f, depth_target, range_target) R_f S_f_full * H_f # 接收频谱 rx_waveform np.fft.irfft(R_f, n4096) # 按实际采样率整理 fs 2 * f2 time_axis np.arange(len(rx_waveform)) / fs这段代码里最容易调错的是delta_f。频点间隔不是越小越好它决定信道仿真的最大无混叠时延T_max 1 / delta_f。如果信道多途时延扩展到100ms而delta_f0.5Hz则T_max2s足够如果你图省事把delta_f放到5HzT_max0.2s多途尾部会折叠进主脉冲波形就会混叠失真。我一般先跑一次单频域脉冲响应量出最大时延扩展再倒推delta_f取它的1/2作为安全裕量。另外要注意H_f必须是复数值且Kraken输出的幅度谱实际上已经包含1m归一化。若想得到绝对声压级这里面还要叠加一个SL的幅值因子A 10^(SL/20)把H_f乘上A即可。发射信号谱和数据采集系统的频响如果不想建模至少保留带通滤波别把带外噪声混进IFFT。5. Kraken 使用避坑清单我从翻车现场捡回来的 5 条经验5.1 现象field 输出大面积 NaN.prt里出现无穷大这个翻车我在模拟项目X刚启动时遇到过。现象是运行正常但pressure文件里每隔几行就会冒NaNTL图出现一条条白色裂缝。后来排查发现flp里的接收深度写到110米而环境文件里介质层底只有100米110米落在了半空间里。简正波模态在半空间内的采样点不是标准定义域求值就溢出。解决把flp的接收深度范围收敛在水层内部我习惯留5%的余量。比如100米水深接收深度下限设到95米避免触碰半空间边界。如果你真的需要海底附近的声压就把环境文件的层底和海底半空间参数重做而不是强行让field去算界面上奇异的点。5.2 现象提高频率后模态数量骤降TL曲线震荡加剧有次从100Hz切到500Hz发现.mod文件里模态数减半。看.prt提示海底衰减项随频率剧增高阶模态在海底反射时损失太大无法形成稳定的驻波结构程序自动截断。震荡是因为参与叠加的模态太少声场不光滑。解决先检查环境文件里海底衰减参数是否用了dB/λ为单位。同样的衰减值在频率升高后换算成dB/m会成倍增长导致海底“变软”。我一般把半空间衰减的物理含义换算成Q值或等效沙质海底参数而不是机械调低衰减系数。另一个补救是提高介质分层精度若海底附近声速梯度很陡模态会聚集在水层内而不进入半空间模态截断会自然稳定。5.3 现象合成波形拖尾严重尾部还有周期重复这是delta_f没选好的典型症状。我把频点间隔设成5Hz做完IFFT后主脉冲后面跟着一串周期成分。起初以为是信道多途后来把时延轴拉长一看重复周期恰好是1/50.2s这就是时域混叠0.2s外的能量被折叠进来。解决先算信道脉冲响应的最大时延扩展再把delta_f取它的1/2以下。浅海低声速海底环境下2000米范围内的多途时延扩展经常到100ms量级delta_f建议不大于1Hz。同时把IFFT点数加大到4倍以上给拖尾留出足够空间别让尾部截断导致突变。5.4 现象两个阵元的波形幅度对但互相关峰总是斜的有个阵元波形合成任务做匹配滤波时主峰位置总是偏几个采样点。查来查去问题出在pressure文件读取时把复声压的绝对值存了下来相位完全丢失。做IFFT时只用幅度谱等效于对所有多途分量强制同相位叠加时延信息被改造成幅度信息相关峰自然不正。解决直接从pressure读实部、虚部合成复数数组任何过程都不要对复数取abs()再存储。这看起来像低级错误但数据分析流程长了之后很容易在某一步顺手取模。我在管线里加了断言检查如果H_f是实数类型直接抛异常从机制上杜绝。5.5 现象换电脑后 Kraken 运行报错 “invalid mode file”一次从服务器迁移到Windows本机环境文件不变运行Kraken时却报读取环境文件失败。检查后发现是文件编码问题在Windows记事本另存时默认带了BOM头Kraken按Fortran风格逐列读取时BOM字符占了第一行第一个位置整个文件解析错位。解决用支持强制无BOM的编辑器重新另存或者统一用VS Code设置编码为UTF-8 without BOM。另外检查文件行尾符Unix环境用LFWindows下用CRLF可能导致读取异常建议全部转成LF再上传运行环境。这事听起来很玄学但跨平台跑水声工具箱时属于高频事故。6. 一个进阶技巧用模态相位直接做阵元波形预白化如果你已经走通了单频点合成波形这条路下一步可以把Kraken输出里的模态特征值直接拿来用而不只是等最终pressure文件。具体技巧是从.mod文件读出每阶模态的本征值k_l和深度模态函数φ_l(z)对每个阵元位置预先构造阵元响应向量。这样做的好处是匹配场或者宽带波束形成时拷贝向量可以直接写成v [exp(j k_l * r_n) * φ_l(z_n)]的线性组合形式省去每个频点重复跑field的开销。我在处理8元垂直阵、32个频点、50个搜索网格时用这个办法把拷贝向量计算时间从几分钟压到几秒。前提是环境水平不变阵元在同一个距离覆盖范围内如果距离跨越很大或者地形起伏明显简正波的水平不变假设就不成立了别硬套。做这个预白化时我习惯顺手做一次自检拿一个单频点的field输出和直接用模态叠加公式手算的结果对比偏差控制在0.5dB以内才继续用。这个自检成本很低但能挡掉大多数模态截断和参数换算问题。我现在接手任何新环境第一件事永远是先跑一版Kraken看模态数分布再决定后面用简正波还是射线模型。这个习惯帮我绕开了不少后端的“假信号”问题希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
