☰
Matlab主动声呐建模:LFM信号与匹配滤波仿真实践
2026/10/22 8:32:26 网站建设 项目流程

简介:面向水声通信与信号处理方向的初学者、课程设计学生及刚接触声呐仿真的研究者,提供了一套基于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 * R

n是扩展因子,球面扩展取2;R是传播距离,单位米;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/(1+f_khz^2) + ... 40*f_khz^2/(4100+f_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 Hz,Nyquist要求采样率至少4500 Hz,10k采样率留出了足够余量,还能顺带容纳一些后续扩展的高频分量。如果采样率设得太低,混叠会把回波信号弄得面目全非,匹配滤波结果自然也是乱的。

仿真总时长设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都找不到在哪。这个模型后续还能往波束成形、目标多普勒估计、主动声呐脉冲串处理方向扩展。你如果也在做水声方向的仿真,可以把它当个起点,按自己的频率和场景参数改一改,跑出来的结果会带你去发现很多有意思的新问题。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询