在水声信号处理这个圈子里DEMON谱分析几乎是每个做被动声呐目标识别的人绕不开的基本功。拿到一段水下噪声别人先听我们先看谱——尤其先看DEMON谱因为螺旋桨一转宽带噪声的包络就会被周期性调制这个调制频率就是目标的“节奏”相当于舰船的指纹。这篇文章我想从头到尾讲清楚DEMON谱分析怎么做包含物理原理、关键参数设计、Matlab完整可运行源码以及我实际跑数据时踩过的几个坑。无论你是刚接触水声信号的学生还是需要快速实现目标特征提取的工程师照着这篇文章的思路走一遍基本就能把DEMON谱这件事做扎实。1. DEMON谱分析的物理基础为什么螺旋桨噪声里藏着“指纹”1.1 舰船辐射噪声的三类来源水下目标辐射噪声通常可以分成三类机械噪声、螺旋桨噪声和水动力噪声。机械噪声来自主机、辅机、轴系等运动部件的振动通过船体耦合到水中主要表现为低频离散线谱。水动力噪声是水流经船体和附体时产生的湍流噪声频谱较宽但调制特征弱。真正让DEMON谱分析“管用”的是螺旋桨噪声。螺旋桨噪声的核心是空化噪声。简单说桨叶在高速旋转时叶片表面压力降低到一定程度后水会汽化形成气泡这些气泡到高压区又迅速溃灭产生一串宽带脉冲。这个过程随机性强能量分布宽直观听起来就是“沙沙”的宽带噪声。但关键点在于空化并不是均匀发生的它跟桨叶的旋转位置有关。桨叶转到某个角度时空化更剧烈噪声能量更大转到另一个角度时空化相对减弱。整个宽带噪声的包络就跟着螺旋桨的转角被周期性调制。这就是DEMON谱能被用于目标识别的最根本原因。如果我们把宽带噪声看成载波螺旋桨调制就是信号那么包络里就藏着轴频和叶频的信息。只要提取出这个调制周期就相当于拿到了目标的螺旋桨工作状态。1.2 轴频与叶频从调制周期到目标“心跳”轴频是螺旋桨每秒的转数单位Hz。比如螺旋桨转速150转/分钟轴频就是2.5Hz。叶频等于轴频乘以桨叶数比如4叶桨叶频就是10Hz。两者的物理意义要分清楚轴频对应的是“轴每转一圈”的调制周期叶频对应的是“每个桨叶叶片经过同一位置”时产生的空化脉冲重复率。实际DEMON谱上轴频位置不一定只有一个峰经常出现的是轴频的基波、二次谐波甚至更高次谐波。原因也不难理解螺旋桨的调制波形不是干净的正弦波而是类似脉冲串的形状脉冲串的傅里叶展开天然就包含丰富的谐波。叶频也有类似情况但幅值结构因船而异。所以在做目标识别时我一般不看单根谱线而是看整组谐波结构。如果你的谱图上只有孤零零一根峰先别急着下结论很可能是分析参数没调好或者数据长度不够。谐波族才是更可靠的识别特征。1.3 DEMON谱与LOFAR谱的分工一个看节奏一个看频率骨架刚接触水声信号的人经常把DEMON谱和LOFAR谱混在一起。两者的区别说穿了很简单LOFAR谱关注的是原始信号本身的低频线谱结构用来识别机械振动特征比如主机转速、轴系不平衡等。它是“静态”的频谱骨架。DEMON谱关注的是宽带噪声包络的调制频率用来识别螺旋桨工作状态。它提取的是“动态”的调制节奏。我见过不少项目把两者放在一起用LOFAR谱找机械线谱和轴频的倍频关系DEMON谱找轴频和叶频。两者互相印证目标识别的置信度才会高。单独拿DEMON谱去识别目标当然也可以但信息维度上会少一层遇到低速目标或调制深度小的目标时单靠DEMON谱容易漏判。2. 处理链路总览与关键参数设计2.1 从宽带噪声到调制频率完整的DEMON处理链路DEMON谱分析的标准链路可以概括为四步带通滤波、包络检波、低通滤波、谱分析。别小看这四步每一步都有讲究。带通滤波是为了从全频段信号里圈出空化噪声主导的频段把机械线谱和环境干扰挡在外面。包络检波是把调幅信息从宽带噪声里“解调”出来。低通滤波是为了只保留我们关心的低频调制分量去掉检波后产生的高频噪声。最后一步对包络信号做FFT或者功率谱估计在低频段找峰值得到轴频和叶频。这个链路看起来简单实际参数敏感度很高。我调试代码时早期版本直接把“带通滤波”这步省了结果包络谱里混入大量未滤除的线谱互调分量峰的位置乱七八糟。后来老老实实按链路走每个环节各司其职谱图才干净。2.2 带通频带怎么选避开线谱、圈住空化噪声带通频带选择是整个DEMON谱分析里最影响结果的一步。如果带通范围选得不对后面的分析再精致也白搭。选择的基本原则是尽量避开强机械线谱所在的低频段选空化噪声能量占主导的中高频段。工程上常见的做法是选5kHz到20kHz或者10kHz到30kHz。低频段比如几十Hz到几百Hz多数是机械线谱的天下如果把这些线谱也纳入检波强线谱和调制信号之间会产生大量互调分量谱面上会出现莫名其妙的假峰。带通滤波器的阶数也要注意。阶数太低过渡带太宽相邻频段的干扰滤不干净阶数太高滤波器群延迟会变大对短数据分析有影响。我自己常用的折中是四阶Butterworth带通滤波器幅频特性平坦群延迟可接受Matlab里一条命令就能生成。2.3 数据长度、帧长与频率分辨率的关系这是DEMON谱分析的硬约束。轴频通常在几赫兹到十几赫兹要想分辨出2.5Hz和2.6Hz这样的差别频率分辨率必须足够细。频率分辨率的计算公式是Δffs/N也就是1除以信号时长。如果你只拿到1秒钟数据频率分辨率是1Hz轴频2.5Hz和3Hz根本分不开。要分辨0.1Hz的频率间隔至少需要10秒数据。实测中我一般要求单段连续数据不低于15秒才能稳定看到轴频基波和二次谐波。拿到更长数据后是不是直接整段做一次FFT就行也不是。海洋环境噪声是非平稳的一段30秒的数据里某艘船可能突然转向、变速整段处理会把不同状态混在一起。更稳妥的做法是分帧加窗多帧平均既保留频率分辨率又平滑掉随机波动。帧长取10到20秒、重叠50%是比较实用的配置。3. 包络检波的核心原理与三种Matlab实现3.1 包络检波的数学原理宽带噪声怎么“显形”出调制信息包络检波是DEMON谱分析最核心的一步。假设带通滤波后的信号是x(t) m(t) · n(t)其中n(t)是宽带空化噪声m(t)是螺旋桨引起的慢变调制函数。我们的目标是从x(t)里把m(t)的频率成分提取出来。平方检波的思路很直接x²(t) m²(t) · n²(t)。这里n²(t)包含一个直流分量和快速起伏分量而m²(t)的频率成分恰好包含调制频率及其谐波。对x²(t)做低通和谱分析后调制频率就露出来了。这是经典的功率检波路线。为什么平方检波能保留调制频率你可以把m(t)展开成直流加各次余弦分量平方后这些分量之间会产生和差频项其中包含2倍轴频、轴频与叶频的差频等。虽然交叉项会带来一些额外谱线但主峰位置仍然对应轴频和叶频所以平方检波在工程上一直被大规模使用。3.2 平方检波、绝对值检波与Hilbert复包络的对比Matlab里实现包络检波有三种常见方式% 方式一平方检波 env_sq xb .* xb; % 方式二绝对值检波 env_abs abs(xb); % 方式三Hilbert变换求解析包络 env_hilbert abs(hilbert(xb));三种方式我都跑过对比说下实际感受平方检波计算量最小实时性最好对正弦型调制的提取效果稳定是传统DEMON谱分析的首选。它的缺点是动态范围被平方放大信号强的地方被抬得更高但这对后续谱分析影响不大因为我们关心的是峰值位置。绝对值检波相当于线性检波输出包络近似正比于m(t)本身动态范围比平方检波温和。它的问题是绝对值运算在零点附近有一个“折角”会产生额外的谐波和直流偏置谱图上的杂散成分会多一些。Hilbert变换求包络最“精准”得到的env_hilbert就是信号的真实幅度包络。但计算量明显高于前两者边界效应也需要处理。在做精细的调制深度分析时我会用Hilbert法常规识别用平方检波就够了。传统DEMON谱分析里多用平方检波还有一个原因平方检波后包络的直流分量与调制分量的比值理论上只跟调制深度有关方便做归一化对比。实测数据里各目标的绝对幅度差异很大归一化后更容易统一阈值。3.3 检波后的低通滤波处理为什么要做包络检波输出里除了低频调制分量还有大量高频成分。这些高频成分来源有两个一个是宽带噪声本身的快速起伏另一个是平方检波产生的倍频项。如果不对包络做低通滤波就做FFT高频成分会通过混叠和谱泄漏干扰低频段的判断。低通截止频率按你关心的最大调制频率来定。轴频和叶频一般不超过几十赫兹所以截止频率设在50Hz到100Hz是合理的。我习惯设50Hz保留足够的余量同时把无关高频压掉。滤波后再去直流均值谱图会干净很多。4. 完整仿真源码从生成水声信号到DEMON谱输出4.1 仿真信号设计思路用调制宽带噪声模拟真实螺旋桨辐射噪声在跑实海数据之前我强烈建议先用仿真把原理和代码流程理顺。仿真信号的设计思路是先产生宽带空化噪声再用螺旋桨调制函数去调制它。调制函数里同时放轴频基波、轴频二次谐波和叶频这样可以模拟出真实目标谱图上的典型结构。仿真参数设固定值采样率50kHz、信号时长30秒、轴频2.5Hz、桨叶数4片。这样叶频是10Hz。调制函数故意做得不完全正弦加入二次谐波和相位偏置更接近真实螺旋桨的脉冲型调制。环境噪声用低电平白噪声模拟让信噪比不至于高到失真。4.2 完整可运行的Matlab源码%% 水声信号DEMON谱分析完整示例 % 功能生成带螺旋桨调制特征的仿真水声信号并完成DEMON谱分析 % 环境MATLAB R2016b及以上版本需信号处理工具箱 clc; clear; close all; rng(2024); %% 1. 仿真信号生成 fs 50000; % 采样率 50 kHz T 30; % 信号时长 30 s t (0:1/fs:T-1/fs).; N length(t); % 螺旋桨参数 f_shaft 2.5; % 轴频 2.5 Hz对应150 RPM Z 4; % 螺旋桨叶片数 f_blade f_shaft * Z; % 叶频 10 Hz % 调制函数轴频基波、二次谐波、叶频共同作用 m 1 ... 0.60 * cos(2*pi*f_shaft*t) ... 0.35 * cos(2*pi*2*f_shaft*t) ... 0.45 * cos(2*pi*f_blade*t pi/6); % 生成宽带空化噪声白噪声经带通滤波 [b_noise, a_noise] butter(4, [5000 20000]/(fs/2), bandpass); band_noise filter(b_noise, a_noise, randn(N, 1)); % 合成接收信号低电平白噪声模拟环境背景 x band_noise .* m 0.05 * randn(N, 1); %% 2. 处理链路第一步带通滤波 % 仿真生成时已经带通一次这里再带通一次是为了模拟实际接收数据的完整处理流程 xb filter(b_noise, a_noise, x); %% 3. 包络检波平方检波 envelope xb .* xb; %% 4. 低通滤波保留50 Hz以下的调制信息 f_max_demon 50; [b_lp, a_lp] butter(4, f_max_demon/(fs/2), low); env_low filter(b_lp, a_lp, envelope); env_low env_low - mean(env_low); % 去直流避免零频分量掩盖真实峰 %% 5. 分段FFT与平均 seg_len 15 * fs; % 帧长15s频率分辨率0.0667Hz overlap 0.5; hop round(seg_len * (1 - overlap)); win hanning(seg_len, periodic); win win / sum(win); % 归一化避免平均谱幅度随窗函数变化 S_demon zeros(seg_len/2 1, 1); cnt 0; start_idx 1; while start_idx seg_len - 1 N idx start_idx : start_idx seg_len - 1; x_seg env_low(idx); X fft(x_seg .* win, seg_len); S_demon S_demon abs(X(1:seg_len/21)).^2; cnt cnt 1; start_idx start_idx hop; end S_demon S_demon / cnt; % 频率轴 f_demon (0:seg_len/2). * (fs / seg_len); %% 6. 绘制DEMON谱 figure(Color, w, Position, [100 100 900 500]); plot(f_demon, pow2db(S_demon), LineWidth, 1.2); xlim([0 30]); grid on; xlabel(调制频率 / Hz); ylabel(DEMON谱幅度 / dB); title(仿真水声信号的DEMON谱); %% 7. 自动标注主要峰值 db_demon pow2db(S_demon); [pks, locs] findpeaks(db_demon, MinPeakHeight, 0.3 * max(db_demon), ... MinPeakDistance, 1); hold on; plot(f_demon(locs), pks, ro, MarkerSize, 8, LineWidth, 1.5); for k 1:length(locs) text(f_demon(locs(k)) 0.3, pks(k), sprintf(%.2f Hz, f_demon(locs(k)))); end hold off;这段代码复制到Matlab里可以直接运行。运行后在0到30Hz范围内正常能看到三个突出的峰分别位于2.5Hz、5Hz和10Hz附近对应轴频基波、轴频二次谐波和叶频。第二个峰的高度不一定比第一个低因为调制函数里的系数是人为设的实际目标也经常出现二次谐波比基波强的现象。4.3 运行结果解读怎么从谱图上读出轴频和叶频拿到谱图后第一步不是看幅度而是看峰值的位置是否构成谐波族。2.5Hz、5Hz、10Hz这三个峰的关系很清晰2.5Hz是基频5Hz是它的2倍10Hz是它的4倍。其中2.5Hz就是轴频10Hz是叶频。判断轴频和叶频时要注意轴频是谐波族中的最小间隔不一定是最高峰。如果调制波形里二次谐波能量远大于基波最高峰可能会落在5Hz而不是2.5Hz。这时候如果直接取最高峰当轴频就错了。正确做法是找整组峰的“共同差频”也就是相邻峰的最小频率间隔。叶频的判断要结合目标信息。桨叶数是多少叶频和轴频的比值就该是多少。这个比值是整数通常在3到7之间。如果谱图上出现10Hz的峰而轴频是2.5Hz比值是4合理如果比值是非整数那这个峰大概率不是叶频而是别的干扰成分。5. 实操翻车点排查从谱图异常回到参数调优5.1 零频那座“大山”和它带来的假象我做DEMON谱最早的一个版本跑出来的谱图在0Hz处有一个巨大无比的峰轴频峰完全被它压住整个图谱就是一条从0Hz开始快速衰减的曲线什么信息都读不出来。原因后来定位得很清楚包络检波后的直流分量没有去掉。平方检波的输出里包含一个很大的直流成分它本质上是m²(t)的直流项加上n²(t)的直流项。如果不减均值直接做FFT这个直流分量会以零频线的形式铺满谱图甚至通过窗口旁瓣泄漏到相邻几个Hz把轴频基波淹没。解法也简单env_low - mean(env_low)。这一步相当于对包络信号做去直流处理等价于在频域把零频分量清零。别小看这一行代码少了它基本看不到任何峰。5.2 窗函数、谱泄漏与弱峰被淹没分段FFT时窗函数选择很关键。如果用矩形窗频率分辨率最理想但旁瓣衰减很差第一旁瓣只有-13dB左右。轴频峰如果幅度差距大旁边的弱峰很容易被矩形窗的旁瓣“盖住”。我在调试时试过一组对比同一份数据矩形窗和汉宁窗跑出来矩形窗在2.5Hz峰旁边多了一堆“锯齿”汉宁窗则干净得多。原因就是汉宁窗旁瓣衰减快能把泄漏压到-31dB以下。代价是主瓣变宽了约一倍频率分辨率略降。在这个应用里轴频和叶频的频率间隔通常足够大至少2到3Hz以上主瓣变宽一点完全能接受。所以我的建议是DEMON谱分析直接用汉宁窗除非你关心的两个峰频率间隔极小才考虑用主瓣更窄的窗函数。5.3 谐波族识别为什么不能只盯“最高峰”这个坑在仿真里不突出但到了实测数据里非常常见。很多目标的DEMON谱上二次谐波甚至三次谐波的幅度可能超过轴频基波。如果你程序里写的是“找最高峰当作轴频”结果就会偏到两倍频上去。我做过一次实测数据处理某目标轴频1.8Hz但谱图上最高峰在5.4Hz。后来把峰值列表全部打出来发现1.8Hz、3.6Hz、5.4Hz、7.2Hz一列等差峰等差间隔1.8Hz这才确定真实轴频是1.8Hz。单独看5.4Hz那个峰完全会误判。所以做自动识别时别用“全局最大峰”找轴频要用“谐波族匹配”找一组等间隔峰的最小间隔或者结合LOFAR谱上轴系的机械线谱去互相验证。5.4 常见DEMON谱异常现象的排查顺序我整理了一个排查表格按出现频率排序实操时可以按这个顺序查现象可能原因解决方案0Hz处巨大谱峰包络检波后未去直流env_low - mean(env_low)轴频峰被旁瓣淹没矩形窗泄漏严重改用汉宁窗峰位置偏移、频率不准帧长太短导致分辨率不足加长帧长或拼接连续数据高频段出现周期性假峰平方检波产生的倍频交叉项低通截止频率调低多个谐波族交错多目标信号叠加单独搜索每个谐波族结合目标轨迹峰极不稳定每次跑结果不同数据段选择不当、非平稳多帧平均或改用更平稳的时段遇到谱图异常时我通常是先回到数据本身看原始波形和LOFAR谱确认该频段内是不是真的只有宽带噪声。如果LOFAR谱上有明显的强线谱落进了带通范围第一步就该调带通而不是在后端折腾窗函数和平均次数。链路顺序决定问题排查顺序数据入口不对后面再怎么优化都是白搭。6. 从仿真到实海把DEMON谱做成能用的识别工具6.1 实测数据处理中必须补的功课仿真代码跑通了不等于实海数据就能直接套用。实测数据有几个特点需要额外处理。第一是信噪比明显降低。环境噪声、航船干扰、生物噪声都会叠加在目标信号上。仿真里我用的环境噪声只有0.05幅度实海环境经常是目标信号幅度跟噪声幅度在同一量级甚至更低。这时候需要增加多帧平均次数、延长观测时间把随机噪声压下去。第二是背景不平稳。海面波浪、风速变化会让噪声包络本身出现随机涨落这些涨落的谱可能会落在轴频附近。处理的办法是在选数据段时避开明显突变段比如突然出现大海浪或者近距离目标经过的时刻。如果一个数据段里包络整体幅度忽高忽低先做归一化处理再进DEMON分析。第三是目标调制深度可能很小。空化噪声调制深度跟航速、吃水深度都有关系低速目标调制深度可能只有0.1甚至更低。这时候谱峰会被噪声淹没。我的经验是尝试多个带通频段因为不同目标的最强调制频段可能不同。比如A目标在8kHz到12kHz调制最明显B目标在15kHz到20kHz最明显。只跑一个频段很容易漏掉。6.2 轴频自动估计与谐波族匹配流程把DEMON谱做成自动化工具时峰值搜索逻辑要完整。我的实现思路是第一步在0.5Hz到15Hz范围内搜索所有局部峰记录幅度和频率。第二步对每个候选峰计算它与其余候选峰的频率差看是否存在近似等间隔的谐波族。第三步取谐波族中所有相邻间隔的最小公约数作为轴频估计。第四步用估计的轴频反推叶频范围如果目标桨叶数未知搜索3到7倍轴频处是否有峰如果已知桨叶数直接检查对应倍频处是否有峰。这套逻辑麻烦一点但比“找最高峰”可靠得多。实测项目里它成功处理过多个轴频相近目标同时存在的情况——两个目标的轴频差只有0.3Hz靠单纯找峰完全分不开靠谐波族匹配才能各自锁定。6.3 DEMON谱的进一步扩展思路DEMON谱做到一定程度还可以往两个方向扩展。一是双谱分析。双谱可以检测频率成分之间的相位耦合关系。螺旋桨调制产生的谐波之间通常存在确定的相位关系双谱图上会出现明显的峰而随机噪声的相位关系混乱不会形成双谱峰。这相当于给DEMON谱加了一个“确认环节”提高抗干扰能力。双谱计算量大但只算切片的话运算量可控。二是DEMON谱与LOFAR谱联合识别。LOFAR谱上轴频的机械线谱可能和DEMON谱的调制轴频存在倍频关系。两个特征互相验证能大幅降低单谱分析的误判率。我在实际项目里目标识别输出的格式通常是一个特征向量轴频、叶频、桨叶数估计、LOFAR线谱频率组、谐波谱峰幅度比。这个向量后续可以接到分类器里做目标类型识别。最后分享一点个人心得。刚接触DEMON谱时我也走过“拿到数据直接跑FFT”的弯路谱图一团糟还以为是算法有问题。后来把每个环节拆开逐段验证带通后看时域包络是否明显有节奏、检波后看低通输出是否干净、最后才做谱分析。这种“逐级验证”的习惯帮我省了大量排查时间。你现在把这套流程捋清楚了后续接任何实测数据心里都有底。