简介:针对水中声呐模型构建需求,这份Matlab代码包提供了一套从零开始的简实现方案,适合水声通信、信号处理方向的初学者与科研人员快速上手。模型覆盖主动声呐与被动声呐的基本工作流程,包括声速剖面计算、发射脉冲生成、传播衰减模拟、回波接收及目标检测等环节,代码结构清晰,便于逐段理解原理并改造复用。压缩包共14个文件,以9个.m脚本为主干,分别对应主程序、距离/角度解算、运动仿真、初始化等模块,另有5个.asv自动保存文件可作调试参考,整体仅8KB,小巧轻量。目前已有3160人学习下载,配套代码可帮助读者在Matlab中直观复现声呐探测过程,掌握信号生成、滤波与目标识别的基础方法,为进一步研究水声通信与声呐算法打下实践基础。
1. 项目概述与建模思路
我一直觉得,很多刚接触水下声学的人对“声呐模型”这四个字有莫名的畏惧感,好像必须要有深潜器、水听器阵列、一堆硬件才能动手。其实用Matlab建立一版简化的水中声呐模型,完全是可以在宿舍里完成的事情。我自己就是从一段几十行的代码开始,把一个主动声呐从发射、传播、回波到检测的完整链路跑通的。这篇博文我会直接把整套思路、物理公式和代码拆开讲清楚,适合刚入门水声工程、做水下机器人感知、或者单纯想用Matlab做信号仿真的人参考。
这里说的“声呐模型”,核心解决的是这样一个问题:在水下环境中,一个声源发出一段脉冲信号,信号经过水介质传播,碰到目标后反射回来,接收端收到一个微弱且被噪声污染的回波,我们怎么根据回波的时延推断目标距离,怎么判断目标是否存在。“简单建立”则意味着我们在模型里做必要的理想化简化,不追求信道多径、海面海底边界反射等复杂因素,先跑通主干流程,再去逐步加复杂度。
1.1 水中声呐模型到底在模拟什么
声呐(Sonar)这个词来源于Sound Navigation and Ranging,本质上就是通过声波在水下的传播特性来探测目标。主动声呐的工作流程可以用四步概括:发射换能器把电信号转成声信号向水中辐射;声波在海水里向前传播,同时因为扩展和吸收产生衰减;遇到目标后一部分声波反射回来,这个回波强度跟目标本身的反射能力有关;接收换能器把声信号转回电信号,经过放大、滤波、时延估计,最终给出目标的距离和方位。
这四步放到Matlab里,每一步都能用一段代码对应。发射信号可以用一个脉冲波或者线性调频信号表示;传播衰减可以用球面扩展加海水吸收模型估算;目标回波可以等效成发射信号的延迟、缩放版本;接收端则通过匹配滤波或者相关运算来抑制噪声、检测回波。这个过程并不涉及复杂的偏微分方程求解,也不需要用有限元把整个声场网格化,而是用“射线声学+声呐方程”的思路做系统级仿真,这对理解声呐原理和验证算法来说已经够用了。
1.2 建模时做了哪些简化
既然是“简单建立”,就必须明确划掉哪些物理过程。我的第一版模型做了五个理想化假设:海水是均匀介质,声速恒定取1500m/s,不考虑温盐深变化;声波按球面波扩展,不存在波导效应;目标当作单个点目标处理,忽略目标形状和姿态;传播路径只有直达路径,没有海底海面反射造成的多径效应;目标和声呐平台都是静止的,不考虑多普勒频移。
这些假设看起来“很不真实”,但恰恰是这种简化让问题变得可分析、可调试。如果你一开始就把多径、随机信道、阵列波束全塞进模型里,代码跑出来的结果有问题时,你根本分不清是发射端的问题、传播模型的问题还是检测算法的问题。先做一版理想模型,把每一行代码和每一步物理过程对应起来,后续再逐步替换模块,这是我比较推荐的学习路径。
2. 声呐模型的理论基础:从声呐方程到传播损失
写代码之前,先把背后的公式理解透。Matlab代码本质上就是把声呐方程和信号处理流程翻译成程序语言,公式理解了,代码只是表达形式的问题。
2.1 主动声呐方程逐项拆解
主动声呐方程是系统设计的核心工具,它把声源级、传播损失、目标强度、噪声级和接收指向性指数统一到一个等式中,用来衡量接收端的信噪比。方程写成:
SNR = SL - 2TL + TS - (NL - DI)
逐项解释一下。SL是声源级,单位dB,代表声源辐射声强的对数表示,声源级越高,信号能传得越远。TL是单程传播损失,声波从声源到目标是一趟,从目标反射回来又是一趟,所以方程里是2TL。TS是目标强度,反映目标反射声波的能力,一个半径1米左右的水下目标,典型TS值在-20dB到-10dB之间。NL是环境噪声级,海洋里风浪、生物、航运都会产生噪声。DI是接收指向性指数,物理含义是接收阵相比全向接收能抑制多少环境噪声。
举个例子,假设SL=210dB,TL=60dB,TS=-15dB,NL=70dB,DI=20dB,代入方程得到SNR=210-120-15-(70-20)=25dB。这个值是正的,说明回波能从噪声里被检测出来。如果算出来是负数,检测就比较困难了。在我们的信号级仿真里,不会直接用这个方程算最终结果,但可以用它来估算参数设置得是否合理,比如发射功率够不够、目标距离是否在检测范围内。
2.2 传播损失与Thorp吸收公式
传播损失TL包含两部分:几何扩展损失和介质吸收损失。几何扩展在我们假设的球面波条件下,等于20倍的对数距离,也就是TL_geo = 20log10(R),这里R的单位是米,结果单位是dB。更精确的声呐方程还会区分球面扩展和柱面扩展,但在简化模型里球面扩展就够了。
吸收损失跟声波频率关系很大。低频声波在水里传播损失小,能传很远的距离,这也就是为什么远程声呐通常用几百赫兹到几千赫兹的频率。高频声波分辨力好,但衰减快,只适合短距离高精度测量。Thorp公式是工程上常用的海水吸收经验公式,以kHz为频率单位,给出吸收系数alpha,单位为dB/km。公式长这样:
alpha = 0.11 * f^2 / (1 + f^2) + 44 * f^2 / (4100 + f^2) + 2.75e-4 * f^2 + 0.003
公式第一项和括号里的分式用于描述低频段硼酸和硫酸镁的弛豫吸收,最后那个常数项代表纯水吸收。计算时注意f是频率,单位kHz,算出来的alpha是dB/km,要换算成dB/m需要除以1000。举个例子,f=20kHz时,第一项约0.11400/401≈0.1097,第二项约44400/4500≈3.911,第三项0.11,最后加0.003,alpha≈4.13dB/km。也就是说20kHz的声波每传播一公里强度要衰减4.13dB。
于是单程传播损失写成TL = 20log10(R) + alpha * R / 1000。这样一个式子就把几何扩展与吸收都考虑了。注意这里alpha用dB/km,R是米,所以第二项要除以1000换算成km。
3. Matlab代码实现与逐段讲解
接下来进入正题:把上述物理模型转成Matlab代码。我的目标不是写一个复杂的功能包,而是用最直白的方式呈现主流程,每一行代码都有明确对应。整套代码分为三块:参数初始化与发射信号生成、目标回波模拟与噪声叠加、匹配滤波检测与绘图。
3.1 参数初始化与发射信号生成
我先定义仿真参数。载频选择20kHz,这个频率在水下算是中高频,适合几百米范围内的探测场景,吸收衰减可以接受,波长也足够短,能保证分辨率。采样率设为200kHz,这样每个载波周期有10个采样点,既能保住波形细节,又不会让数据量过大。脉宽取5ms,对应的时间带宽积在单频脉冲情况下比较小,为了检测分辨率更好,后续可以升级成线性调频信号。
% 声呐模型仿真参数设置 c = 1500; % 水下声速,单位m/s fs = 200e3; % 采样率,单位Hz fc = 20e3; % 发射信号中心频率,单位Hz T = 5e-3; % 脉冲宽度,单位s R_true = 100; % 目标真实距离,单位m SNR_dB = 10; % 接收端信噪比,单位dB(回波信号与噪声功率比) % 发射信号:单频矩形脉冲 t = 0 : 1/fs : T - 1/fs; tx_signal = sin(2 * pi * fc * t); tx_power = mean(tx_signal.^2);这里t是脉冲持续时间内的时间轴,从0到T,步进1/fs。tx_signal用正弦函数生成单频脉冲,mean(tx_power)用于后续计算噪声功率时做基准。为什么不直接用cos而用sin?其实没有本质区别,但sin(0)=0能让发射信号从零开始,避免在仿真起始时刻出现电流跳变,虽然这个影响在回波检测里基本可以忽略,但养成分段信号从零开始的习惯总没坏处。
3.2 目标回波模拟与噪声叠加
信号传播到目标再反射回来,距离是2倍R_true,所以回波时延用2R/c计算,100米距离对应约0.1333秒。回波幅度用传播损失来决定:强度衰减因子是10的负TL/10次方,幅度衰减因子则是强度衰减因子的平方根。同时按SNR_dB设置噪声功率,确保回波在噪声中处于合理的可见程度。
% 计算单程传播损失TL R_km = R_true / 1000; f_khz = fc / 1e3; alpha = 0.11*f_khz^2/(1+f_khz^2) + 44*f_khz^2/(4100+f_khz^2) + 2.75e-4*f_khz^2 + 0.003; TL = 20*log10(R_true) + alpha * R_km; % 回波时延和衰减 delay = 2 * R_true / c; % 搜索往返时延 n_delay = round(delay * fs); % 换算成采样点数 amp_loss = 10^(-(2*TL)/20); % 双程衰减对应的幅度系数 % 构建接收信号时间轴,时间长度要覆盖回波 recv_len = n_delay + length(tx_signal) + 2000; % 尾部多留2000点 rx_signal = zeros(1, recv_len); rx_signal(n_delay+1 : n_delay+length(tx_signal)) = amp_loss * tx_signal; % 加入高斯白噪声,使回波信噪比近似为SNR_dB signal_power = mean((amp_loss*tx_signal).^2); noise_power = signal_power / (10^(SNR_dB/10)); noise = sqrt(noise_power) * randn(1, recv_len); rx_signal_noisy = rx_signal + noise;这里的核心是amp_loss = 10^(-(2*TL)/20)。为什么括号里是2TL而不是TL?因为声波经历的是双程传播,传播损失要算两次,然后用20除以是因为我们要的是幅度衰减而不是功率衰减,功率衰减因子是10^(-(2TL)/10),幅度要开平方,所以变成10^(-(2TL)/20)。这个细节很容易算错,我一开始就是只用了单程TL,导致回波幅度虚高,检测距离被严重高估。
3.3 匹配滤波检测与结果绘图
接收数据准备好了,接下来就是检测。匹配滤波本质上是让接收信号与发射信号的共轭翻转序列做卷积,当接收信号中出现与发射信号相似的回波时,卷积输出会出现一个明显的峰值,峰值位置对应回波时延。
% 匹配滤波处理 mf_output = filter(fliplr(tx_signal), 1, rx_signal_noisy); mf_output = mf_output / max(abs(mf_output)); % 搜索峰值并估计距离 [peak_val, peak_idx] = max(abs(mf_output)); est_delay = (peak_idx - 1) / fs; est_range = est_delay * c / 2; fprintf('真实距离: %.2f m\n', R_true); fprintf('估计距离: %.2f m\n', est_range); fprintf('距离误差: %.2f m\n', abs(est_range - R_true)); % 绘图 figure('Position', [100 100 1200 800]); subplot(3,1,1); plot((0:length(tx_signal)-1)/fs*1000, tx_signal); xlabel('时间 (ms)'); ylabel('幅度'); title('发射信号'); xlim([0 T*1000]); subplot(3,1,2); plot((0:length(rx_signal_noisy)-1)/fs*1000, rx_signal_noisy); xlabel('时间 (ms)'); ylabel('幅度'); title('带噪接收信号'); xline(delay*1000, 'r--', '真实回波时刻'); subplot(3,1,3); plot((0:length(mf_output)-1)/fs*1000, mf_output); xlabel('时间 (ms)'); ylabel('归一化输出'); title('匹配滤波输出'); xline(delay*1000, 'r--', '真实回波时刻'); grid on;用filter函数实现匹配滤波时,fliplr(tx_signal)是对发射信号做时间翻转,这一行是整个检测的核心。如果信号是单频脉冲,匹配滤波的效果其实和自相关差不多,但遇到线性调频信号时,匹配滤波能获得脉冲压缩增益,峰值会明显更尖锐,这也是为什么实际声呐里更常用LFM信号的原因。
4. 运行结果与参数影响分析
代码写完直接跑,默认参数是100米距离、10dB信噪比、20kHz载频。我实际跑完这版代码,匹配滤波峰值出现在0.1333秒附近,换算成距离是99.98米左右,误差很小,在几厘米量级。这个误差主要来自时延量化,因为回波时延要取整成采样点,你取整损失的那点时间换算成距离就是量化误差。
4.1 默认参数下的检测结果解读
从三张图能看得很清楚。第一张图里发射信号就是一个5毫秒的20kHz正弦波。第二张图里,如果不告诉你回波在哪,肉眼基本看不到100米距离处那个微弱的回波,因为幅度经过双程传播衰减后已经非常小,而且被高斯白噪声淹没了。第三张匹配滤波输出在0.1333秒处出现一个明显的峰值,其他位置的噪声被有效抑制,这就是匹配滤波对白噪声的抑制能力和对已知信号的积累增益。
这里要强调一个概念:匹配滤波不是“把噪声去掉”了,而是把分散在整个脉冲时间里的信号能量集中到一个峰值点上。单频脉冲本身没有频率调制,所以匹配滤波输出的峰值宽度大约是1/T,也就是200Hz带宽对应的时域分辨率,5毫秒脉宽对应的主瓣宽度在毫秒量级,这决定了两个距离相近的目标能不能被分辨开。
4.2 改距离、噪声、目标强度会有什么变化
为了验证模型行为是否合理,我做了三组参数实验。第一组把距离从50米拉远到500米,保持其他参数不变。结果是:距离越远,匹配滤波峰值越低,因为双程传播损失快速增加。50米时信噪比绰绰有余,250米时峰值已经明显变矮,500米时在10dB信噪比设置下勉强可见。这说明在不加时间增益控制或者脉冲积累的情况下,单频脉冲声呐的探测距离是有限的。
第二组调整SNR_dB,从5dB改到20dB。噪底明显下降导致峰值检测稳定性大幅提升,但这个“SNR_dB”在代码里指的是回波信号功率与噪声功率的比值,它是在已知目标距离和衰减后反推噪声功率得到的,相当于“上帝视角”设置。实际系统中你无法预知信噪比,只能通过积累时间和带宽来控制。
第三组把目标距离改成350米,同时把目标强度由默认的-20dB改成-10dB,你会发现回波幅度明显增强。这说明目标反射特性对检测影响很大,同样的声呐系统探测大鱼群和探测小鱼群的性能可能差很多,这个现象在实际渔业声呐中非常明显。
我把这三组实验的关键观察整理成一个表:
| 实验变化 | 固定参数 | 观察到的趋势 | 物理原因 |
|---|---|---|---|
| 距离50m到500m | SNR=10dB,fc=20kHz | 峰值逐渐降低,500m时接近噪底 | 双程球面扩展+吸收损失增大 |
| SNR 5dB到20dB | R=100m,fc=20kHz | 噪底降低,峰值更突出 | 噪声功率减少,检测更稳定 |
| TS -20dB到-10dB | R=350m,SNR=10dB | 回波幅度明显提升 | 目标反射能力增强,回波强度升高 |
这个模型的行为和物理直觉一致,说明简化版的声呐方程和信号级模型搭得是自洽的。做仿真最重要的不是追求代码复杂,而是先让你的模型输出符合物理规律。
5. 常见问题与排查技巧实录
写这版代码的时候,我自己踩了不少坑,也帮别人查过不少类似问题。这里整理几个典型的报错和参数陷阱,如果你照着代码跑出奇怪的结果,优先检查这几个环节。
5.1 典型报错与参数陷阱速查
| 现象 | 常见原因 | 排查方法 |
|---|---|---|
| 匹配滤波峰值出现在0时刻 | 接收信号长度不够,回波还没进来就截断了 | 检查recv_len是否足够覆盖n_delay加发射信号长度 |
| 估计距离总是偏大或偏小固定误差 | delay取整导致时延量化偏差 | 确认n_delay用的是round不是floor,可改用interp1做亚采样精度估计 |
| 回波被噪声完全淹没,看不到峰值 | noise_power计算错误(符号或倍率) | 检查signal_power是否用衰减后信号计算,SNR公式是否用了10而不是20 |
| alpha算出来是负数或异常大 | 频率单位或系数搞错 | Thorp公式频率必须用kHz,alpha单位是dB/km |
| 绘制图形时报错维度不匹配 | 时间轴长度和信号长度不一致 | 统一用length()而不是硬编码数字 |
5.2 实测中容易忽略的三个细节
第一个细节是噪声功率的单位域问题。很多人在算noise_power时直接拿发射信号功率除以10^(SNR/10),但发射信号还没有经过双程衰减,如果拿原始信号功率去配噪声功率,加入噪声后回波根本看不见,因为你在用“未衰减的发射信号”去定义信噪比。正确的做法是先算衰减后的回波信号功率,再以此配置噪声功率,这样SNR_dB才真正代表接收端回波的信噪比。
第二个细节是时延不一定落在整数采样点上。100米距离在1500m/s声速下时延是0.133333秒,乘以200kHz采样率等于26666.6个采样点,取整后丢失0.6个采样点的信息,对应约2厘米的距离误差。这个误差在小范围高精度测距场景里不能忽略。想进一步降低量化误差,可以对匹配滤波输出做抛物线插值,或者用更高采样率、更高载频。
第三个细节是max(abs(mf_output))搜索范围内的边界效应。接收信号尾部多留的2000个采样点并不是随意写的,如果尾部留得太少,匹配滤波输出还没完全衰减到尾部边界,峰值搜索可能把边界处的截断尖峰误判成目标。我建议尾部至少留一个脉冲宽度的长度,也就是T*fs个点,这样才能让滤波器状态充分衰减。
6. 从简单模型到实用模型的扩展方向
一个能跑的简化模型只是起点,实际工程声呐系统远比我这个demo复杂。我在做完基础版本之后,按以下几个方向逐步加功能,每条路都能看到明显的效果变化。
6.1 在Matlab里做二维/三维波束
把单个接收换能器升级成均匀线阵或者平面阵,就能在Matlab里仿真波束形成。核心思路是把每个阵元的接收信号按不同时延对齐相加,某个方向的信号因为同相叠加被增强,其他方向的信号因为相位错开被抑制,这就是数字波束形成的基本原理。Matlab的Phased Array System Toolbox可以直接用phased.UCA、phased.URA这些对象建阵,但自己手写一个简单的延迟求和波束形成器,反而更能理解这个过程的本质。
6.2 接入工具箱和真实数据
如果要处理真实的水声数据,或者做更复杂的水声信道仿真,认真推荐几个方向。MATLAB的Phased Array System Toolbox里提供了一系列声呐系统对象,包括信号源、信道、接收机全链路,适合系统级验证。UCL和CMRE等机构发布过一些公开的水声实验数据,把真实采集的WAV文件读进Matlab,套上这段匹配滤波流程做检测,你会立刻体会到真实信道的“恶意”——多径、起伏、瞬态干扰都会让检测变难。这时候再回头优化模型,你对每个处理模块的理解就会更深入。
我个人在实际搭建声呐仿真模型时的体会是:先让代码跑通,再让它跑对,最后才让它跑快。第一个版本不要追求模型复杂,能检测出100米处一个简单回波,你已经理解了主动声呐的完整链路。之后每加一个模块,都要在固定距离、固定目标下做对比测试,这样任何新引入的错误都能被及时发现。如果你照着这版代码跑出了类似的结果,下一步不妨试着把发射信号换成线性调频信号,对比一下匹配滤波输出峰值的形状变化,这个实验做完,你对声呐信号处理的理解会上一个台阶。
本文还有配套的精品资源,点击获取