做水声被动目标识别DEMON谱分析基本上是绕不开的第一步。哪怕现在深度学习的当量越来越大传统特征仍然没有退出战场很多实际系统里DEMON谱仍然是目标分类识别的重要输入之一。原因很简单它用最简单的包络分析手段从宽带的舰船辐射噪声里“挖”出螺旋桨转动的周期性特征进而估计出轴频、叶频甚至叶片数信息量非常大计算量却小到可以在嵌入式平台实时跑。这篇文章围绕Matlab实现DEMON谱分析的完整链路展开从物理机理、仿真信号建模、核心算法实现到关键参数调整和实测数据中常见的“坑”我会连同可直接运行的源码一起给到。适合水声信号处理的初学者打基础也适合刚接手被动声纳目标识别任务的工程师做工程参考源码可以直接拿去改。1. DEMON谱的物理基础为什么“包络”能透露螺旋桨的秘密1.1 舰船辐射噪声中的“隐性节拍”先说机理。舰船在水下航行时辐射噪声由很多成分构成但对于被动声纳关心的中高频段大约几百赫兹到几十千赫兹范围最主要的宽带成分来自螺旋桨空化。螺旋桨叶片高速旋转时叶片梢部压力降低产生空泡空泡破灭的时候辐射出宽带噪声。这个噪声不是平稳白噪声它受到螺旋桨旋转周期调制每个叶片转到某一位置时空化强度不一样于是宽带噪声的幅度会随轴频呈周期性起伏。问题在于这个“起伏”非常微弱而且通常是随机宽带噪声叠加在大量环境噪声中直接从时域波形上看肉眼根本看不到明显的规律。DEMON谱做的事情就是把这个微弱的周期性幅度调制给“解调”出来再从频谱上找节奏。打个比方你听不清一个人说话的内容但只要他每秒钟固定敲两下桌子你总能从背景噪音里把这个节拍数出来——DEMON谱就是那个“数节拍”的动作。1.2 DEMON名称拆解与输出解读DEMON是Detection of Envelope Modulation On Noise的缩写翻译过来就是对噪声进行包络调制检测。它处理完后输出的不是常规功率谱而是一个“调制谱”横轴是调制频率纵轴是调制强度。谱线上第一个明显的峰对应轴频f_s后面间隔等距的峰对应叶频f_b n × f_s × Z这里的Z就是螺旋桨叶片数。因此根据谱峰位置可以直接估计螺旋桨轴频由相邻谱线间隔推算叶片数联合多个时刻的DEMON谱还能跟踪目标转速变化趋势。这些参数对目标分类识别意义很大。比如商船通常轴频低、叶片数4到6片有些高速巡逻艇轴频高、叶片数少这些差异会在DEMON谱上体现得非常清楚。1.3 这方法的工程价值在哪很多人第一次接触DEMON谱是从教材上觉得就是“带通滤波-检波-FFT”三件套简单得像课后习题。但真正放到工程里它有几个別的工具替代不了的优势计算量低一次处理几秒数据只需要做一次滤波、一次平方运算、一次FFT实时性极好稳健性强不依赖目标的线谱是否稳定对宽带空化噪声的检测效果尤其好可解释性强输出直接对应轴频、叶频等物理量不像某些端到端特征那样不可解释。所以纵使LOFAR谱低频线谱分析广泛应用于低频段目标检测DEMON谱在中高频段的螺旋桨参数估计上依然是标配两者通常配合使用。2. 仿真信号先造一个“标准答案”再动手处理2.1 三层结构的信号模型在实际声纳数据上调试算法前我强烈建议先做仿真验证。仿真最大的好处是“标准答案已知”——你知道真实的轴频、叶频是多少处理完直接对比就能定位算法每个环节是否正确。我们用的仿真模型可以分成三层背景噪声层模拟海洋环境噪声通常用高斯白噪声近似载波层模拟螺旋桨空化噪声本质是宽带随机信号可以看成白噪声经过带通滤波调制层模拟螺旋桨旋转对宽带载波的幅度调制由轴频基波、叶频倍频和一定的相位偏移叠加构成。观测信号可写成x(t) [n_b(t) × (1 m × e(t))] n_a(t)其中n_b(t)是宽带空化噪声e(t)是归一化的调制包络m是调制深度n_a(t)是环境噪声。调制深度m越大包络周期性越强DEMON谱上的线谱就越突出。2.2 参数设置与Matlab参考代码我习惯设置采样率20kHz时长30秒轴频2.5Hz叶片数5片。这样叶频基波是12.5Hz倍频在25Hz、37.5Hz处肉眼验证很方便。代码如下clear; close all; clc; rng(2024); fs 20000; % 采样率 20 kHz dur 30; % 信号时长 30 秒 N fs * dur; t (0:N-1) / fs; % ---- 载波层模拟空化噪声的宽带随机信号 ---- [b, a] butter(6, [300 8000]/(fs/2), bandpass); carrier filtfilt(b, a, randn(1, N)); % ---- 调制层轴频 叶频倍频叠加 ---- shaft 2.5; % 轴频 2.5 Hz Z 5; % 叶片数 5 env 0.6 0.8 * sin(2*pi*shaft*t); for k 1:Z env env 0.4 * cos(2*pi*shaft*Z*t k*0.5); end env env / max(abs(env)); % ---- 构造观测信号调制载波 环境噪声 ---- SNR -6; % 信噪比 -6 dB sig carrier .* (1 0.6 * env); sig sig / sqrt(mean(sig.^2)); noise randn(1, N); noise noise / sqrt(mean(noise.^2)); x sig * 10^(SNR/20) noise;这段代码生成的信噪比是-6dB意思是目标强度比背景噪声低大概一半。这个条件不算苛刻但也足够让你先体会到“从噪声里抠特征”的感觉。2.3 为什么仿真必须先于实测调试我见过不少同学拿真实录音直接跑DEMON跑出来一堆毛刺谱线第一反应是怀疑算法写错了其实是数据本身信噪比低、目标变速、多目标干涉等原因造成的。先跑仿真你手里有了“标准答案”算法任何问题都会立刻暴露然后用同样的处理流程去套实测数据这时再出问题就可以放心地往数据特性方向排查。另外仿真还有一个作用调节参数。你可以自由控制轴频、叶片数、信噪比、调制深度系统性观察它们对DEMON输出谱线的影响。这个手感积累对后续分析真实数据非常宝贵。3. 核心处理流程与Matlab源码实现3.1 算法链路全景标准的DEMON谱分析链路看起来简单但每个环节都有讲究。完整流程是带通滤波选取目标信号能量集中的频段避开低频强烈干扰和高频过度衰减区间包络提取采用平方或绝对值检波把幅度调制信息搬移到低频段低通滤波滤除载波残余分量保留调制谱频段通常只关心几百赫兹以内去直流与去趋势去掉包络平均值避免零频分量淹没低频线谱加窗与FFT对包络信号做功率谱估计峰值提取在调制谱中搜索等间距的谱峰估计轴频、叶频。这里我先做一个细节点说明包络提取为什么常用平方检波因为平方运算在数学上等效于“信号乘以自身”它会将载波频率翻倍同时把包络调制搬到低频。而包络的周期成分恰恰很低频后续低通滤波可以很干净地把载波翻倍项丢弃只留调制成分。绝对值检波也有类似效果但会引入更多高次谐波相比之下平方检波更干净。3.2 带通滤波器怎么取舍带通范围没有唯一标准答案取决于目标频率特性和环境干扰分布。以螺旋桨空化噪声为例它在中频段谱级相对平坦所以我通常选300Hz到8000Hz这一段。下限不能太低否则会混入螺旋桨线谱以外的低频机械噪声这些成分本身可能带有很强的周期性会和真正的螺旋桨调制混淆上限也不能太高高频段传播损失大远程目标的高频分量衰减严重留那么高只是徒增计算量。3.3 完整源码DEMON谱分析函数下面是完整可复用的MATLAB函数输入观测信号后直接返回调制谱的频率轴和功率谱密度function [f, pxx] demon_spectrum(x, fs, bpf, lpf_fc, ... win_len, overlap_ratio, nfft) % DEMON谱分析 % 输入 % x - 观测信号序列 % fs - 采样率 Hz % bpf - 带通滤波范围 [fL fH] Hz % lpf_fc - 低通滤波截止频率 Hz % win_len - 单帧长度单位秒 % overlap - 帧重叠比例 0~1 % nfft - FFT 点数建议为 2 的幂 % 输出 % f - 频率轴 % pxx - 平均调制功率谱 % 1. 带通滤波 [b1, a1] butter(4, bpf/(fs/2), bandpass); x filtfilt(b1, a1, x); % 2. 平方包络检波 env x .^ 2; % 3. 低通滤波 [b2, a2] butter(4, lpf_fc/(fs/2), low); env filtfilt(b2, a2, env); % 4. 分帧、加窗、FFT与平均 frame_len round(fs * win_len); hop round(frame_len * (1 - overlap_ratio)); win hann(frame_len, periodic); num_frames floor((length(env) - frame_len) / hop) 1; pxx_sum zeros(1, nfft/21); for idx 1:num_frames start (idx - 1) * hop 1; seg env(start : start frame_len - 1); seg seg .* win; seg seg - mean(seg); % 去直流 seg seg / sqrt(sum(win.^2)); % 能量归一化 spec abs(fft(seg, nfft)).^2; pxx_sum pxx_sum spec(1:nfft/21); end pxx pxx_sum / num_frames; f (0:nfft/2) * (fs / nfft); end这个函数把流程封得很干净可以直接复制到一个.m文件里使用。调用时大致的参数是bpf [300 8000]; % 带通范围 lpf_fc 1000; % 低通截止 win_len 4; % 4秒一帧 overlap 0.5; % 50%重叠 nfft 8192; % FFT点数 [f, pxx] demon_spectrum(x, fs, bpf, lpf_fc, win_len, overlap, nfft); figure; plot(f, 10*log10(pxx eps)); xlabel(调制频率 (Hz)); ylabel(调制谱功率谱密度 (dB)); xlim([0 30]); grid on;如果仿真参数是前面那些轴频2.5Hz、5叶片你会在2.5Hz附近看到明显峰12.5Hz、25Hz、37.5Hz处看到间隔相等的倍频峰。5叶片对应的叶频基波是12.5Hz叶片数越多倍频间隔越大。3.4 关于帧长和FFT点数的一个提醒有读者可能疑惑低通截止1kHz为什么nfft要取8192这段时间采样率20kHz8192点对应频率分辨率是20000/8192 ≈ 2.44Hz。这其实不够好——2.44Hz的分辨率意味着两个相距小于2.44Hz的谱峰没法区分而真实轴频可能是1.8Hz或3.2Hz2.44Hz根本无法稳定分辨。所以我实际使用时更喜欢让nfft自动取决于帧长而不是固定值。帧长4秒时nfft取2^nextpow2(fs * win_len)也就是8192还是2秒对应的4096如果是8秒帧长就取65536等一下还是算一下。fs*win_len80000点nextpow2是131072频率分辨率约0.15Hz那就很舒服了。代码里统一写成nfft 2^nextpow2(round(fs * win_len));这样频率分辨率严格由帧长倒数决定不存在“选错nfft”的问题。4. 关键参数选择四个直接影响结果的“旋钮”4.1 分析带宽的取舍原则带通范围的选择是在“有用信号能量”和“干扰抑制”之间做平衡。我前文推荐300~8000Hz但不代表这是万能组合。具体调整时可以参考三条经验如果目标在远距离优先压缩高频端比如取300~3000Hz避免高频噪声拉低整体信噪比如果近岸或港口环境低频干扰严重把低端抬到500Hz甚至800Hz以上如果同时想保留丰富的调制信息分析带宽最好跨一个倍频程以上比如2kHz到8kHz这样调制包络形态才能完整保留下来。一句话带通范围没有标准答案它是第一个需要根据数据反复试验的参数。4.2 帧长与频率分辨率的权衡帧长越长频率分辨率越高但同时时间分辨率变差目标的变速跟随能力变弱。这里的“时间分辨率”指的是你能把目标转速变化精确定位到什么时刻。我在工程上习惯按目标的机动性来定场景推荐帧长频率分辨率理由远距离稳定航行的商船8~10s0.1~0.125Hz转速稳定适合高分辨中近距离巡逻快艇2~4s0.25~0.5Hz转速变化较快需时间分辨强非平稳/机动目标1~2s0.5~1Hz只求趋势不求精细谱线帧长太短还有一个弊端低频谱线会“糊”掉。假如轴频是2.5Hz帧长只有0.5秒频率分辨率2Hz2.5Hz的峰几乎看不见。所以遇到高速目标时先算一下频率分辨率是否小于目标轴频的一半如果临界就需要延长帧长或放弃精细分辨。4.3 重叠率与平均次数重叠率的作用有两个一是提高帧利用率让短数据也能得到足够多的平均次数二是让帧与帧之间平滑连接减小频谱估计方差。经验值50%到75%都可以我常用75%因为计算量增加不多但谱线平滑度提升肉眼可见。平均次数呢理论上Welch法平均次数越多方差越小但代价是时间分辨率下降。工程上平均次数达到10帧以上后方差下降速度明显变慢。比如4秒帧长、75%重叠30秒数据大约能切出28帧平均次数非常充裕如果数据只有10秒4秒帧长就只能切出5~6帧这时可以降到2秒帧长换取更多平均次数。4.4 去直流、窗函数和归一化的细节这三个细节经常被初学者忽略但直接影响谱线形态。去直流必须在加窗之后、FFT之前做否则窗函数本身会带来一个常数的谱泄漏。更好的做法是直接从均值里减去而不是用高通滤波去直流因为高通滤波在低频端有可能把轴频附近的频谱也压下去。我在函数里用的是seg seg - mean(seg)很简单但很可靠。窗函数上Hann窗是DEMON谱的通用选择主瓣略宽但旁瓣衰减快适合线谱检测。矩形窗旁瓣太凶除非你清楚自己在做什么否则别用。FFT前的能量归一化seg / sqrt(sum(win.^2))是为了保证不同帧长、不同窗函数的功率谱幅度具有可比性。你要是只关注相对谱峰位置这个归一化可做可不做但如果你要比较不同航次、不同目标之间的调制强度就必须归一化否则谱级完全不可比。5. 实测数据中的常见问题与排查5.1 参考线谱混入翻倍的“假轴频”实测数据里最常见的问题是DEMON谱里出现一个很亮的窄峰但它并不是轴频。例如机器转速线谱通过带通滤波器泄漏进来或者变频器相关电磁干扰源直接耦合进采集系统都会产生一条窄带线谱而在平方检波后这个线谱会在它的两倍频处形成一个新的峰很容易被误判成轴频。排查方法很简单对比多个带通范围的DEMON输出。如果这个峰的位置跟着带通范围变化说明它是滤波边缘的边界效应如果峰的位置固定不变那大概率是外部干扰线谱不是螺旋桨调制产生的。再不行就对原始数据单独做一遍窄带频谱看这条线在原始谱上的确切频率。5.2 轴频半倍频程模糊叶频还是轴频还有一个常见的坑就是谱峰间隔刚好是轴频的两倍导致把“间隔的一半”误判为轴频。这种现象通常出现在调制包络以谐波为主、基波很弱的时候。比如实际轴频5Hz但4叶片螺旋桨在包络中以20Hz分量为最强5Hz的基波被淹没如果不了解先验信息很容易把20Hz当成轴频。这种情况下我会先计算相邻谱峰间隔再把间隔依次除以1、2、3看哪一组的最小公约数能同时解释所有峰的位置。叶片数通常是4到7的整数联合这个约束多数模糊可以消除。5.3 变速目标导致谱线展宽目标加速或减速时轴频不是常数长时间的FFT会把一个单一谱峰“涂抹”成一条宽带。这时固定帧长的DEMON谱就不够用了。我的处理习惯是先用短帧1~2秒快速扫一遍观察谱峰位置随时间的变化趋势再看是否做线性调频校正或分段估计。如果只是需要转速趋势短帧逐帧描点就够了不需要做特别的时频分析工具。5.4 低信噪比下的三个对策当目标距离远、信噪比低到调制峰几乎看不见时三条路可以走延长分析带宽收集更多空化噪声能量提高平方检波后的调制信号强度拉长帧长并用高重叠率做更多平均代价是时间分辨变差用多频段并联处理把多个子带的包络信号加权相加再统一做DEMON分析相当于空间分集的频域版本。前两条没什么技术门槛第三条是一个更进阶的技巧。具体实现时把原始信号按频段拆成4~8个子带每个子带各自平方检波、低通滤波再按信噪比加权叠加成一个总包络最后对这个总包络做FFT。这个方法在宽带空化噪声占优时效果很显著我实测过在-12dB左右仍能稳定提取轴频。6. 从DEMON延伸后续可深挖的几个方向6.1 多帧DEMON瀑布图单帧DEMON谱看的是静态特征多帧堆叠成时间-频率平面图就是DEMON瀑布图。它能直观展示轴频、叶频随时间的变化曲线特别适合做目标机动跟踪。在Matlab里用循环调用demon_spectrum函数把每一帧的结果按时间轴叠起来再用imagesc画二维图即可代码量很小视觉效果和数据信息量都提升一个档次。6.2 循环平稳谱分析DEMON谱是按“平方解调”的思路做的本质上是利用了信号的循环平稳特性。更严格的做法是计算循环谱密度在循环频率轴上直接检测轴频。循环谱密度能区分同一频率处的不同调制来源抗干扰能力更强但计算量和理解门槛都高不少。对于只想快速提取螺旋桨特征的项目DEMON谱够用做研究发论文可以考虑循环平稳方向。6.3 与LOFAR谱配合我一直主张DEMON和LOFAR搭配使用而不是二选一。LOFAR谱在低频段抓线谱、判断目标存在性DEMON谱在中高频段抓包络调制、估计螺旋桨参数。两边信息互相印证识别置信度比单看任何一方都高很多。我做目标识别的时候会把两类特征拼成一个特征向量再送进分类器。特征维度不高但物理意义明确分类效果比纯数据驱动方法更稳健尤其在小样本场合这种“物理特征轻量分类器”的组合非常能打。最后分享一点个人体会DEMON谱入门不难但真正做到“准”靠的是对数据的耐心和对参数的反复调试。仿真阶段建议把每个环节都可视化逐步确认带通、检波、低通、FFT每一步的中间结果处理实测数据时再多想想哪些峰是目标、哪些峰是干扰。拿一套已知轴频和叶片数的船载数据练手把流程跑通再换几组不同目标的数据试这套方法的“手感”很快就有了。