
简介面向水声通信与信号处理方向的初学者、课程设计学生及刚接触声呐仿真的研究者提供了一套基于Matlab的简化声呐建模实现围绕主动声呐与被动声呐的基本原理给出可运行、易扩展的代码结构与实验路径。模型主线涵盖声速剖面计算、发射脉冲生成、水中传播衰减模拟、回波接收、信号滤波与目标检测判断等环节涉及Boussinesq方程/水深声速数据、正弦脉冲序列、距离-衰减因子以及匹配滤波等常见处理手段能够帮助读者将抽象的水声理论落到具体代码层面。压缩包共14个文件主体为9个.m脚本另有5个.asv自动备份文件整体仅8KB非常轻量适合在Matlab中逐行阅读和反复调试。已有3160人学习代码按主程序、参数初始化、角度范围限制、距离与到达时间计算等模块拆分结构清晰便于定位关键逻辑无论是用于理解声呐系统工作原理、完成信号处理实验还是作为进一步构建复杂水声模型的基础这套代码都能提供直观的起点。 最近在整理水下目标探测的仿真项目需要一套能快速验证信号处理算法的声呐模型就用Matlab从零搭了一个最简版本。所谓“最简”不是三五行代码糊弄过去而是把主动声呐的完整链路——发射、传播、目标反射、接收、匹配滤波——每个环节都保留但每个环节只做必要精度。这篇文章会给出这套水中声呐模型的完整Matlab实现思路、参数设计方法以及调试过程中踩过的坑。不管你是刚开始接触水声信号处理的学生还是需要在项目里快速验证目标检测算法的工程师这套代码都可以直接拿来当底座。1. 内容整体设计与思路拆解1.1 为什么从单基地主动声呐模型起步声呐系统大类上分主动和被动。被动声呐“只听不说”靠接收目标辐射噪声来探测但目标噪声特性随机性强、待估参数多对入门者很不友好。主动声呐自己发射声波等回声用回波时延推算目标距离物理链路直观每个环节都能在示波器上看到波形非常适合作为仿真建模的起点。主动声呐里又分单基地、双基地和多基地。单基地就是发射机和接收机放在同一个位置声波发出后沿直线到目标再沿同一条路径反射回来几何关系最简单。我们先把这个跑通后面如果把发射端和接收端拆开到不同坐标就自然扩展成双基地模型。绝大多数水声工程教材里的声呐方程推导也都是从单基地开始的所以这个选择不只是为了省事而是贴合整个学科体系的演进路径。在仿真精度上我刻意做了取舍。初始版本不考虑海面海底反射、声速分层、多途效应也不做海底地形建模只保留最基本的直达波路径和球面扩展衰减。这样的好处是一旦结果不对能很快定位是信号处理问题还是环境建模问题不会出现一堆因素搅在一起无从下手的情况。1.2 声呐方程到代码变量的映射关系主动声呐的工作状态可以用简化主动声呐方程描述SNR SL - 2TL TS - (NL - DI)其中SL是声源级TL是单程传播损失TS是目标强度NL是环境噪声级DI是接收指向性指数。对于入门模型DI可以先忽略NL直接用高斯白噪声来模拟。但实际写代码时我并不建议原封不动地在分贝域做全套计算然后反推回波信号幅度。更实用的做法是把每个物理量拆成代码里的数值变量在生成回波时直接折算成幅度因子。我写代码前先列了下面这张映射表照着写能省掉很多返工声呐方程分量物理含义代码变量说明SL声源级SL_dB 200相对参考值可调TL单程传播损失TL 20*log10(R) alpha*R几何扩展加吸收TS目标强度TS_dB -15不同目标差异很大NL环境噪声级noise_var用高斯白噪声近似后续可换有色噪声c水中声速c 1500淡水海水略有差异2R/c往返时延delay 2*R/c单基地最核心的观测量这里把声呐方程翻译成回波幅度时要注意一个常用公式接收信号幅度与传播距离的关系是A_echo A_ref * 10^(-TL/20) * 10^(TS/20)。20出现的原因是把功率域的压力幅度折算到电压/声压幅度域。很多初学者容易在这里丢掉系数导致回波幅度差几十倍。这套模型的目标不是精确预测某型声呐的探测距离而是快速得到一个“链路通、算法能跑、趋势正确”的仿真结果。精度问题等验证完基本逻辑再逐步优化前期过度设计只会让调参变成一场灾难。2. 核心细节解析与实操要点2.1 坐标几何与声线往返时延的计算建模的第一步是定义场景几何。我把发射机和接收机放在同一位置比如src_pos [0, 0, -20]目标放在[200, 80, -15]单位都是米。这里需要注意水下坐标通常是三维的不用把问题强行压到二维平面直接算欧氏距离最省事R norm(target_pos - src_pos); delay 2 * R / c;这段代码看着简单却是整个基带模型的地基。声波往返时延是所有后续信号时间对齐的依据一旦这里算错匹配滤波的峰位就会整体偏移你还会以为是滤波器写错了。有一个特别容易踩的坑空气中声速约340 m/s水中约1500 m/s差出将近4.5倍。仿真时如果习惯性沿用雷达仿真的参数习惯把声速设成光速量级回波时延会小到不可思议最后匹配滤波的峰值位置和理论值对不上时排查半天才发现是物理常量的问题。建议在代码顶部直接用注释固定好水中声速的物理含义避免头脑发热时改错。如果之后要做深海模型声速就不能再用常数了同一个深度的声速受温度和压力影响很明显可以用声速剖面函数比如Munk剖面替代。但初始版本建议先固定常数把基础链路验证完再叠加分层效应。2.2 传播损失球面扩展加吸收系数的实用近似传播损失TL由扩展损失和吸收损失两部分组成工程上常用这个公式TL n * 10 * log10(R) alpha * Rn是扩展因子球面扩展取2R是传播距离单位米alpha是吸收系数单位dB/m。声音在水下从声源向四周均匀扩散时声压按距离反比衰减对应球面扩展的20log10(R)。这个衰减是几何效应和频率无关频率再低也躲不掉。吸收损失则与频率强相关高频声波被水介质吸收得快所以低频声呐能传得更远。粗略仿真时alpha可以直接取常数比如0.001 dB/m。这个值在几kHz频段附近大致可用但如果仿真频率范围很大最好用更精确的Thorp公式alpha 0.1 * f^2 / (1 f^2) 40 * f^2 / (4100 f^2) 2.75e-4 * f^2注意这个公式输出的单位是dB/km要转成dB/m还需要除以1000。工程实现上我建议把吸收计算单独写成一个函数方便以后替换成更精确的模型不用动主流程function alpha waterAbsorption(f) f_khz f / 1000; % Thorp公式用kHz alpha_db_per_km 0.1*f_khz^2/(1f_khz^2) ... 40*f_khz^2/(4100f_khz^2) 2.75e-4*f_khz^2; alpha alpha_db_per_km / 1000; % dB/m end调参时可以把TL打印出来看一眼数值范围。比如目标在200米外alpha取0.001那么TL大约是20log10(200) 0.001200 ≈ 46.2 dB。这个量的数量级心里要有数后面如果回波幅度太小应该先怀疑是不是TL算得太狠而不是急着调信噪比。2.3 发射波形选型LFM线性调频的工程理由初始模型里发射波形有两种常见选择CW单频脉冲和LFM线性调频脉冲。CW的频率恒定信号结构最简单但距离分辨率受脉宽限制想分辨两个距离很近的目标就得缩短脉宽脉宽一短又牺牲了发射能量探测距离就缩水。LFM的特点是发射期间频率线性扫过一段带宽本质上是在时间上给信号“编了号”。接收端通过匹配滤波可以把能量在时间上重新压缩成窄峰从而在长脉宽能量高和大带宽分辨率高之间同时占住。这就是为什么声呐和雷达里普遍优先考虑LFM波形而不是一个简单脉冲。生成LFM信号的Matlab代码非常短fs 10e3; % 采样率 T 0.05; % 脉宽 B 500; % 带宽 f0 2000; % 起始频率 k B / T; % 调频斜率 t 0 : 1/fs : T; tx cos(2*pi*(f0*t 0.5*k*t.^2));这段代码里起始频率f0是2000 Hz扫频范围从2000到2500 Hz瞬时频率随时间线性上升。后面做匹配滤波时这个发射信号会被当作滤波器的参考模板所以它长什么样直接影响整个系统的性能。有没有必要一开始就把f0定到几十kHz没必要。频率越高吸收损失越大200米距离可能还能接受但到几公里距离时高频信号早被水吃干净了。而且采样率和仿真时长也会随频率上升而猛增对小模型来说纯属浪费。3. 实操过程与核心环节实现3.1 参数声明与场景搭建把基础参数集中放在脚本开头方便后续统一调整。下面是我跑通的第一版完整参数声明% 基本环境参数 c 1500; % 水中声速 m/s % 场景几何 src_pos [0, 0, -20]; % 声源/接收机坐标 target_pos [200, 80, -15]; % 目标坐标 R norm(target_pos - src_pos); % 发射信号参数 fs 10e3; % 采样率 T 0.05; % 脉宽 B 500; % 带宽 f0 2000; % 起始频率 k B / T; % 声呐方程相关 SL_dB 200; % 声源级 TS_dB -15; % 目标强度 alpha 0.001; % 吸收系数 dB/m SNR_dB 20; % 接收端信噪比 % 时间轴 t_total 0 : 1/fs : 1.0; % 仿真时长至少覆盖到回波返回 n_total length(t_total);采样率为什么取10 kHz这不是拍脑袋。LFM的最高瞬时频率是f0 B/2 2250 HzNyquist要求采样率至少4500 Hz10k采样率留出了足够余量还能顺带容纳一些后续扩展的高频分量。如果采样率设得太低混叠会把回波信号弄得面目全非匹配滤波结果自然也是乱的。仿真总时长设1秒对应的最大可测距离是c * T_total / 2 750 m远大于目标距离够用。这个时间长度和采样率决定了数组维度大概是10000个点Matlab处理起来非常快几乎可以实时调参。场景搭建这一步值得多花一点时间确认几何坐标。三维坐标里距离不是简单把x分量相减就完事一定用norm函数求欧氏距离。坐标写好后再把R和delay打印出来目测是否符合物理直觉。比如这里的往返时延约2*215/1500 ≈ 0.287秒合理。3.2 回波生成与匹配滤波回波生成是核心环节。思路是先生成发射信号然后按往返时延把它放到接收序列的对应位置再按传播损失和目标强度折算幅度最后叠加高斯白噪声。t_tx 0 : 1/fs : T; tx cos(2*pi*(f0*t_tx 0.5*k*t_tx.^2)); % 传播损失与回波幅度 TL 20*log10(R) alpha*R; A_echo 10^((-TL TS_dB) / 20); % 初始化接收序列 rx zeros(1, n_total); % 计算发射信号在接收序列中的起始位置 delay_samples round(2 * R / c * fs) 1; end_idx delay_samples length(tx) - 1; if end_idx n_total rx(delay_samples : end_idx) A_echo * tx; end % 叠加高斯白噪声 noise_power 10^(-SNR_dB/10) * mean(rx.^2); noise sqrt(noise_power) * randn(1, n_total); rx_noisy rx noise;匹配滤波器本质上是对发射信号的“时间反转 共轭”卷积也就是说用发射信号本身作为模板在接收序列中逐个位置做相关运算。实现上可以用filter函数h fliplr(tx); % 时间反转因为是实信号这里不需要共轭 y filter(h, 1, rx_noisy);这里有一个概念容易混淆匹配滤波器的输出是相关峰峰位置对应回波到达时刻而不是目标距离本身。从峰值索引换算距离时要记得除以2[~, idx] max(abs(y)); estimated_delay idx / fs; estimated_range estimated_delay * c / 2;用max找峰值只能找到最强的那一个如果后续要检测多个目标或者存在明显的旁瓣建议用findpeaks并设置最小峰高阈值避免被噪声毛刺带偏。关于幅度换算我在最初的版本里直接用了A_echo 10^(-TL/20)结果回波小得几乎看不见后来才想起目标强度TS还没加进去。这个负的TS会进一步削弱回波组合起来就是双重衰减。调通以后理解了目标反射不是完美反射大部分能量被目标吸收或散射到别的方向只有一小部分回到接收机目标强度就是描述“这一小部分有多少”的物理量。3.3 结果可视化与目标距离判读模型跑通后我最关心的三个图分别是发射信号、接收回波、匹配滤波输出。用subplot放在同一张图里一眼就能看出信号处理链路哪里出了问题figure; subplot(3,1,1); plot(t_tx, tx); title(发射信号LFM); xlabel(时间/s); subplot(3,1,2); t_rx (0:n_total-1) / fs; plot(t_rx, rx_noisy); title(接收回波含噪声); xlabel(时间/s); subplot(3,1,3); t_y (0:n_total-1) / fs; plot(t_y, y); title(匹配滤波输出); xlabel(时间/s); ylabel(幅度); % 标注目标距离 [peak_val, peak_idx] max(abs(y)); hold on; plot(peak_idx/fs, peak_val, ro); text_str sprintf(估计距离: %.1f m, peak_idx/fs*c/2); text(peak_idx/fs, peak_val, text_str);匹配滤波输出图上最明显的特征就是在目标时延处出现一个尖锐的窄峰。因为LFM的脉压比是TB 0.05500 25理论上峰值宽度约为脉宽的四分之一左右虽然不算特别窄但已经远高于单个CW脉冲的分辨能力。实际看波形时接收回波可能完全淹没在噪声里人眼根本看不出有信号但匹配滤波输出处的峰值比噪声高出一大截。这就是匹配滤波“把信噪比集中到某个时间点”的能力也是整个模型最有演示价值的地方。第一次跑通并看到那个峰的时候基本就能确定整条链路没有大的逻辑问题。4. 常见问题与排查技巧实录4.1 异常结果速查表仿真模型跑起来之后异常结果会层出不穷。我把自己踩过的问题整理成了一张速查表方便大家对照定位现象可能原因检查方向匹配滤波输出全是平的找不到峰值回波幅度太小被噪声完全淹没检查TL是否过大、TS_dB是否负得离谱、SNR_dB是否设置过高峰值位置与理论距离不一致声速设置错误、时延采样点数取整误差、时间轴起点不对核对c值、检查delay_samples计算、确认时间轴从0开始回波波形严重畸变采样率不满足Nyquist、信号频率超过fs/2提高fs或降低f0/B峰值旁边出现很多次级峰LFM带宽太小导致旁瓣较高、噪声功率过大增加B或加窗抑制旁瓣回波幅度小到看不见传播损失计算时忘记加目标强度或衰减系数设置过大打印TL值检查数量级图像里出现时间上的“双峰”匹配滤波器长度导致边界效应或目标回波与噪声叠加后出现假峰用findpeaks配合最小峰间距处理这张表里最常踩的其实就是第一行和第三行。第一行是链路“信号太弱”的典型症状第三行是参数“超限”的典型症状。其他问题多半是这两类的变体。针对峰值位置问题还要提醒一句delay_samples计算用了round取整这个舍入误差在目标很近时会造成不可忽略的距离偏差。比如声速1500、采样率10k一个采样点对应的距离是0.15米单基地里再折半是0.075米。对于大部分入门仿真这精度够了但如果要做高精度测距应该用高精度时间轴或者对匹配滤波输出做插值峰估计而不是简单取整。4.2 三个容易忽略的工程细节第一个细节是实信号与复信号混用的问题。上面代码里发射信号用的cos实信号匹配滤波用fliplr就够。但如果后面把发射信号换成复数形式的解析信号比如tx exp(1j*2*pi*(f0*t 0.5*k*t.^2))那匹配滤波就需要同时做时间反转和共轭也就是h conj(fliplr(tx))。漏掉共轭会导致匹配滤波后信号里出现频率加倍的分量虽然看起来不影响主峰但会让整个输出底噪抬高影响弱目标检测。经验教训是明确自己代码里用的是实信号还是复信号然后全程保持一致不要混着来。第二个细节是filter函数的边界效应。Matlab的filter输出长度与输入相同但刚开始的一段和结尾一段是滤波器的暂态响应不是真实的相关峰。我之前在短数据上直接找max结果峰值总是出现在序列开头一度以为回波建模错了。排查后发现是边界效应把暂态大值当成了峰值。解决办法是在找峰值前把前L-1个点和后L-1个点截掉其中L是匹配滤波器长度。第三个细节是噪声功率的标定方式。代码里用noise_power 10^(-SNR_dB/10) * mean(rx.^2)来标定噪声功率这里mean(rx.^2)是回波信号的平均功率这样定义SNR后不同SNR下的回波可检测性才具有可比性。但要注意rx序列只有一个时刻附近有回波大多数时间是零所以mean(rx.^2)本身偏小。严格一点应该用发射信号的功率来标定而不是接收序列的平均功率。我后来改成了signal_power mean(tx.^2)再去折算噪声这样SNR的物理含义更清晰。这个细节在很多教材里不会写但对结果的可靠性影响很大。另外还有一个可视化上的经验如果接收回波里噪声太明显、匹配滤波峰不容易看清可以先把噪声去掉跑一遍确认回波到达时间和峰位都正确再加噪声跑第二遍。两次结果对照能快速区分“模型错误”和“噪声干扰”这是排查链路问题最快的路径。这套模型整体跑通以后我又往里面加了多目标回波、随机起伏噪声和收发分置的扩展。每次加复杂度都要回头修改参数表尤其是传播损失那一段。我的经验是先把主链路用最简单版本写透、跑顺再一步步替换物理模型如果一开始就上射线追踪和声场积分大概率连bug都找不到在哪。这个模型后续还能往波束成形、目标多普勒估计、主动声呐脉冲串处理方向扩展。你如果也在做水声方向的仿真可以把它当个起点按自己的频率和场景参数改一改跑出来的结果会带你去发现很多有意思的新问题。本文还有配套的精品资源点击获取