Matlab入门：主动声呐模型与匹配滤波测距仿真

简介面向水声通信与信号处理学习者这份 MatLab 代码资源演示了水中声呐模型的简单建立方法覆盖主动声呐与被动声呐的基本原理和仿真实现可帮助理解水下声波传播、目标检测与回波处理等核心环节。资源共 14 个文件主要包括 9 个 .m 脚本和 5 个 .asv 自动备份文件压缩包仅 8KB代码结构紧凑。脚本中可见 main 主程序、sonar 声呐核心模块、init 初始化、distance 距离计算、Arrive 信号到达、Gettheta 角度估计、motion 运动模型及 jd 等模块便于按功能拆解学习。目前已有 3160 人浏览学习。通过阅读和运行这些代码读者能快速搭建一个简易声呐仿真环境掌握发射信号生成、传播衰减、回波接收与处理的编写思路为后续水声建模研究提供可扩展的参考框架。 “声速设成340这模型是给空气用的不是给水用的。”这句话是我上周帮人调代码时说原话。那位兄弟用Matlab写了一个简单的声呐测距模型代码逻辑看着没啥毛病探测距离却怎么都对不上最后把声速参数一改现象当场消失。这几乎是每个刚开始接触声呐仿真的人都会遇到的一幕代码本身不难难的是你没把水声传播的物理逻辑先捋顺。这篇文章就干一件事把一套最简单的主动声呐模型用Matlab完整搭起来。模型覆盖发射、传播、接收、检测四条链路代码不长不需要任何额外工具箱纯基础Matlab就能跑。适合刚接触水声建模、想快速跑通一个能看的声呐demo或者做课程设计、写论文仿真的人参考。我会把每个参数为什么这么设也讲清楚你跑通之后想改频率、加目标、换信号都能自己动手改。1. 声呐模型不是玄学先拆解水声传播的四个物理环节动手写代码之前先把声呐模型在模拟什么物理过程说清楚。一个最简单的主动声呐模型绕不开四个环节声速、传播衰减、目标回波、环境噪声。这四个东西搞明白了Matlab代码其实就是把它们逐行翻译出来。1.1 声速水下声波的“限速牌”很多人建模时顺手写c340这是空气中的声速放在水里整个模型就废了。声波在水中的速度大约是1500 m/s比空气里快4倍多声速错了延迟、波长、距离全部跟着错。严格来说海水声速不是常数会随温度、盐度、深度变化。工程上常用简化公式c ≈ 1473 1.7T - 0.006T² 1.3(S - 35) 0.017D其中T是温度摄氏度S是盐度PSUD是深度米。对于浅水近程探测这种场景直接用1500 m/s往往够用误差在1%以内对几百米的探测距离来说影响不大。真正的高精度测距系统才需要声速剖面数据这个后面讲踩坑时再展开。1.2 传播衰减为什么探测距离不是想多远就多远声波在水中传播会变弱原因主要有两部分。第一部分是几何扩展。声能从一个点源向四周扩散能量要均匀分布在不断变大的球面上幅度随距离的一次方衰减换算成分贝就是20log10(R)R是传播距离米。这个衰减和频率无关任何声波都躲不掉。第二部分是吸收衰减。声波的能量会被水介质本身转化成热能频率越高吸收越厉害。这就是为什么低频声呐能传得远、高频声呐测距近但分辨率高。工程上常用Thorp公式估算吸收系数alpha 0.1f²/(1f²) 40f²/(4100f²) 2.75e-4 f² 0.003注意这里的f单位是kHz算出来的alpha单位是dB/km。建模时还需要除以1000转成dB/m。这个单位换算坑了我很久后面专门写一节。1.3 目标回波延迟、幅度、多普勒主动声呐的工作原理很简单发射机发出声波声波碰到目标反射回来接收机收到回波。回波和发射信号之间差了一个双程传播时间tau 2R/c这就是测距的基础公式。回波的幅度由目标强度Target Strength, TS决定。不同目标差别很大一群鱼可能是-20 dB一块礁石可能接近-10 dB一颗光滑的金属球又不一样。简单模型里你可以给每个目标设一个相对反射系数不用去抠TS的精确物理定义。目标如果运动回波还会带多普勒频移fd 2vf0/c。v是目标相对声呐的径向速度f0是发射频率。用20kHz的发射信号目标以1 m/s速度靠近多普勒约26.7Hz这个量级在窄带系统里是能测出来的。但静态模型可以先不管它。1.4 环境噪声最简单也最容易被忽略的环节水下的噪声来源很多波浪、降雨、船舶、生物活动甚至海水本身的热噪声。对模型来说不需要一个个去建通常用一个高斯白噪声来代表环境噪声基底信噪比SNR直接反映信号和噪声的相对强弱。Matlab里一行awgn(rx, SNR, measured)就搞定了。别小看这一步没有噪声匹配滤波这一整套检测理论就失去意义了。决定模型真实感的关键往往就在SNR怎么定。2. 一个能跑的Matlab主动声呐模型从发射到测距一条龙下面给出完整代码。这段代码我加了不少注释你直接复制就能跑。它在场景上做了三件事生成一个20kHz的CW脉冲作为发射信号模拟三个不同距离的目标回波叠加最后用匹配滤波检测出目标距离。%% 参数设置 clear; clc; close all; c 1500; % 声速 m/s fs 100e3; % 采样率 Hz T 20e-3; % 脉宽 20ms f0 20e3; % 发射频率 Hz %% 发射信号CW脉冲 t 0:1/fs:T-1/fs; N length(t); s sin(2*pi*f0*t); %% 目标与传播参数 targetRanges [50 120 200]; % 多目标距离 m targetTS [1.0 0.5 0.3]; % 相对目标强度反射系数 fc f0 / 1e3; % Thorp公式用kHz alpha 0.1*fc^2/(1fc^2) 40*fc^2/(4100fc^2) ... 2.75e-4*fc^2 0.003; % 吸收系数 dB/km alpha alpha / 1000; % 转成 dB/m %% 构造接收信号多目标回波叠加 maxDelay 2*max(targetRanges)/c; totalLen ceil(maxDelay*fs) N; rx zeros(1, totalLen); for k 1:length(targetRanges) R targetRanges(k); tau 2*R/c; % 双程传播延迟 delaySamples round(tau*fs); Att 20*log10(R) alpha*R; % 总衰减 dB A targetTS(k) * 10^(-Att/20); rx(delaySamples1 : delaySamplesN) ... rx(delaySamples1 : delaySamplesN) A*s; end %% 加噪声 SNR 10; % 信噪比 dB rxNoisy awgn(rx, SNR, measured); %% 匹配滤波检测 mf xcorr(rxNoisy, s); mf mf(N:end); % 截取零延迟及之后的部分 [~, idx] max(mf); rangeEst (idx-1) / fs * c / 2; fprintf(峰值对应距离: %.2f m\n, rangeEst); %% 画图 tAxis (0:totalLen-1)/fs; rangeAxis tAxis*c/2; figure; subplot(2,1,1); plot(rangeAxis, rxNoisy, b); grid on; xlabel(距离 (m)); ylabel(幅度); title(接收信号含噪声); subplot(2,1,2); plot(rangeAxis, mf, r); grid on; xlabel(距离 (m)); ylabel(幅度); title(匹配滤波输出);2.1 发射信号为什么选择20kHz CW脉冲CW脉冲指的就是正弦波在有限时间内的一段sin(2*pi*f0*t)简单直接。选择20kHz这个频率是基于浅水近程探测的常见工程折中。频率越高吸收衰减越大。20kHz在海水里的吸收系数大约是2~3 dB/km用Thorp公式算对于几百米的探测范围影响很小。但频率也不能太低低频就需要大尺寸换能器指向性差还容易撞上环境噪声较多的低频段。所以20kHz左右是很多近程声呐的常用起点代码里想换频率直接改f0就行。脉宽20ms同样有讲究见后面的参数权衡章节。2.2 传播建模延迟、衰减、多目标叠加构造接收信号这一段的思路是对每个目标先算双程延迟再把发射信号按延迟位置放进去乘上衰减系数。衰减系数A的推导过程是先把球面扩展和吸收衰减换算成分贝再用10^(-Att/20)从分贝转回线性幅度。这个转换别搞反了Att是能量衰减幅度衰减要除以20不是除以10。很多初学朋友在这里栽过。循环里用的是叠加操作因为多个目标的回波在物理上就是线性叠加的。rx(delaySamples1 : delaySamplesN)这一段用累加而不是直接赋值这样多个目标回波重叠时也不会互相覆盖。2.3 匹配滤波检测把回波从噪声里捞出来匹配滤波的本质是把接收信号和发射信号的“镜像”做相关运算。当接收信号里存在与发射信号形状相同的回波时相关输出会出现一个尖峰噪声和发射信号不相关输出就不会有尖峰。从频域角度看匹配滤波器就是对信号频带内的成分加强、带外噪声压制。xcorr(rxNoisy, s)算的是互相关返回长度N totalLen - 1的序列。其中第N个元素对应零延迟时刻所以用mf mf(N:end)截取之后mf(k)就对应于延迟(k-1)/fs秒距离就是(k-1)/fs*c/2。这段对齐逻辑是很多人的知识盲区我第5章会专门展开讲。max找到峰值位置后输出距离就是探测到的主要目标距离。注意这里只找了一个最强峰如果要输出多个目标得用门限加峰值搜索后面进阶部分再提。3. 关键参数与工程权衡这些值不能照着抄把代码复制下来跑通很容易但如果你不理解里面的参数一改场景就废。这一节把代码里几个关键参数的设定逻辑讲透。3.1 采样率、频率、脉宽的三角关系采样率fs 100e3是发射频率f0 20e3的5倍。理论上采样定理要求fs大于2倍最高频率但实际做信号处理时我习惯留2.5到5倍的余量尤其是后面还要做匹配滤波、频率分析采样率太紧会让波形失真、峰位偏移。脉宽T 20e-3决定了两个东西信号能量和频带宽度。脉宽越长发射能量越大回波信噪比越高但CW脉冲的带宽约等于1/T脉宽越长带宽越窄距离分辨力越差。你的系统是想要更远的探测距离还是更精细的距离分辨率直接决定脉宽怎么选。代码里20ms对应50Hz带宽这个窄带特性适合做简单测距demo但实际要用它分辨近距离的多个目标就很吃力。3.2 距离分辨力与盲区CW脉冲的天然限制CW脉冲有两个绕不开的硬限制一个是模糊的一个是根本性的。先算一笔账。两个目标的回波如果间隔小于一个脉宽匹配滤波输出就会合成一个峰分不开。这个最小可分辨时间间隔约等于脉宽T换算成距离deltaR c*T/2 1500*0.02/2 15 m。也就是说这套模型里两个目标距离差小于15米就只剩一个峰。盲区更直接发射机发射脉冲的这20ms里接收机收到的全是自己的发射信号什么目标都看不见。这段时间对应的距离c*T/2 15 m所以15米以内的目标完全被“屏蔽”了。这就是为什么实际声呐很少单纯用CW脉冲做测距要么缩短脉宽要么换调频信号后面4.1会说到。3.3 信噪比SNR怎么设才合理代码里我直接给了SNR 10 dB这在demo里无所谓但真实系统里SNR不是凭空拍脑袋定的背后有个简单的水声链路预算SNR SL - 2TL TS - NL其中SL是声源级TL是单程传播损失TS是目标强度NL是环境噪声级。例如声源级200 dB目标距离100m时TL约44 dB球面扩展40 dB 吸收4 dB目标强度10 dB环境噪声60 dB算下来SNR约62 dB。这个数字远高于10 dB说明这个场景下目标很容易检测。你的模型里如果目标是模拟一个很弱的小目标或者距离拉远SNR就会下降检测可能失效。所以改代码时不要只调SNR的数字要想想它对应实际物理场景中的哪一个环节。4. 模型跑通后往真实系统推进的三条路线基础模型跑通之后你会发现它只是一个“玩具”。把它变成接近真实系统的仿真有三条最常见的升级路线按性价比排序。4.1 换用LFM线性调频信号脉冲压缩的威力如果想立刻感受到声呐建模的“高级感”把CW换成LFM线性调频信号是成本最低的一步。LFM信号的频率在脉宽内从低扫到高代码只需要改发射信号那一行B 5e3; % 带宽 5kHz s chirp(t, f0-B/2, T, f0B/2, linear);匹配滤波的代码完全不用动。但效果天差地别LFM的匹配滤波输出主瓣宽度约等于1/B 0.2 ms对应距离分辨力变成了c*0.0002/2 0.15 m比刚才CW的15米好了整整100倍。同时时间带宽积B*T 100匹配滤波后SNR增益也提升100倍相当于20 dB。这就是脉冲压缩技术的核心脉宽越来越长信号能量越积越多但依靠频率调制把主瓣压到足够窄距离分辨力和探测能力两头都占。20ms脉宽、5kHz带宽的LFM在实际声呐里是个相当常规的参数组合。4.2 加入阵列与波束形成测距之外还能测向单个接收换能器只能告诉你目标多远告诉你方向得靠阵列。均匀线阵是最容易理解和实现的阵列结构每个阵元接收到同一回波的时间略有差异把这个时延差补偿掉再叠加就能让某个方向来的信号“同相相加”其他方向的信号互相抵消。代码上就是先算各阵元相对参考阵元的时延差tau_i (i - 1)*d*sin(theta)/c再对各通道做时延补偿后求和。需要注意的是阵元间距d要小于等于lambda/2否则会出现栅瓣拿到一个假方向。20kHz信号波长为7.5cm所以d取3~4cm比较合适。这个扩展对初学来说是个不小的跳跃但它能帮你理解声呐从“测距仪”变成“成像仪”的关键一步。4.3 目标强度与多普勒效应建模简单模型里目标强度我用了相对反射系数但真实的目标强度随方位角变化、随频率变化。如果只是做数值仿真可以给每个目标定一个TS值然后用TS替换掉代码里的targetTS让衰减计算更物理。多普勒效应则是运动场景里绕不开的。设目标以速度v径向运动回波信号会在发射信号基础上产生频率偏移fd 2vf0/c。模拟多普勒时可以直接在回波里乘一个复指数exp(1j*2*pi*fd*t)。20kHz信号、目标速度2 m/s时fd约为53Hz这个偏移虽然小但积累足够长时间后用FFT就能测出来这也是一切多普勒声呐的基础。5. 我在调试这套模型时踩过的几个坑最后这部分是我反复栽跟头之后总结的。模型跑不出来问题大概率出在下面这四个环节里。5.1 匹配滤波的峰值位置偏移这是我见过最多人犯的错。xcorr和filter都能做匹配滤波但结果对齐方式不一样。如果用filter(fliplr(s), 1, rx)做匹配滤波输出峰值的位置并不在回波到达时刻而是滞后了N-1个采样点——这是FIR滤波器本身的固有群延迟。如果你不修正20ms脉宽对应就会产生15米左右的距离误差。我现在的习惯是想省事就用文中代码的xcorr版本想提高计算效率再用filter版本但用filter时记得把峰值位置减掉N-1再换算距离。两种方法的输出对齐方式一定要在注释里写清楚不然隔两天自己看代码都会懵。5.2 衰减公式的单位换算错误Thorp公式算出来是dB/km但代码里目标的距离单位是米。如果你直接用alpha*RR是米衰减就被低估了1000倍远处的目标回波会大得离谱整个模型的动态范围全是乱的。反过来球面扩展20*log10(R)里的R必须以米为单位。有人会把R写成公里于是衰减又从分贝变成了负几十分贝目标信号微弱到完全看不见。我建议把单位换算写在代码里用变量名区分清楚alphaPerMeter、RinMeters让接手的下一个人一眼就能看懂。5.3 仿真长度与边界效应构造接收信号时totalLen必须覆盖最远目标的完整回波也就是至少delaySamples N个点。如果长度不够回波尾巴被截断匹配滤波的能量会损失峰值位置也可能偏移。还有一个隐蔽问题xcorr结果的头部和尾部存在半窗效应前几个采样点的相关值不完全幅度偏小。如果你的目标特别近回波延迟很短峰值出现在头几毫秒这时要小心别把第一组有效样本当成了噪声。我的做法是仿真时在最远目标后面预留一段保护间隔比如再多加几个毫秒的采样点宁可多算一点也别让边界效应污染检测结果。5.4 声速偏差对测距精度的影响常数声速1500在浅水近程场景下够用但如果你把模型推到更严谨的场景这里是个隐患。实际的浅海声速可能在1450到1520之间变化冬季表层甚至可能低于1480。1500和1480的偏差大约1.3%200米处目标会算出2.6米的距离误差对于数十米量级的分辨力系统来说还能接受但如果你后面要接高精度定位、成像就必须引入声速剖面。我的经验是先在常声速下把算法链路调通再做声速剖面敏感性分析。声速从常数改成随深度变化的剖面本质上只是把每一小段的传播时间累加起来代码复杂度会上一个台阶但物理结果会明显更可信。这套模型我前前后后改过很多版从最简单的单目标CW到多目标加LFM再到后来的阵列波束形成。回头看初始那版代码只是个空壳真正让模型变“真”的是那些物理参数和边界条件的细节。如果你刚开始跑别急着加功能先把这段基础代码调通再去碰上面的扩展方向。遇到测距不准优先检查声速、延迟对齐和单位换算三件事这能省掉你多半天的调试时间。本文还有配套的精品资源点击获取