简介:水声目标识别中的辐射噪声常掩盖目标信号,希尔伯特变换解调可提取信号包络与瞬时特征。这份由脚本和文本数据组成的微型资源包,聚焦三种解调方法的对比实验,面向水声信号处理学习者、算法验证人员以及涉及目标识别与噪声分析的科研场景,帮助理解从实值信号获取同相与正交分量并揭示隐藏特征的分析路径。压缩包共两个文件,分别是.m脚本和.txt记录文件,总大小仅2KB,代码精简,便于逐行研读与二次修改;其中脚本对应三种解调方式的对照实现,文本则可能保存了仿真信号与实际船舶辐射噪声信号的实验数据或结果记录。已有682人学习该资源,适合用于快速搭建仿真对比环境、评估不同解调方法在复杂水声环境下的性能差异,也可作为课程设计或论文实验的参考模板。
1. 辐射噪声包络里的周期量,是水声目标识别最该优先看的特征
被动声呐拿到的舰船辐射噪声,表面上看就是一段宽带随机信号。直接对原始时域信号做 FFT,得到的是连续谱叠加上少量线谱;而螺旋桨叶片周期切割水流引起的空化噪声调制,全部藏在宽带噪声的幅度包络里,在原始功率谱上几乎不可见。用希尔伯特变换解调把包络提取出来,再对包络做频谱分析,得到的就是 DEMON 谱,轴频和叶频会以清晰的低频峰值形式出现在谱图上。这篇文章把这条链路完整拆开:辐射噪声为什么自带周期调制、希尔伯特变换解调在数学上做了什么、DEMON 谱的参数怎么设才能把 1 Hz 附近的轴频找出来,以及这些特征如何落到水声目标识别上。适合刚接手被动声呐数据分析的工程师,也适合想把手里噪声采集程序升级成特征提取管线的开发者。
2. 希尔伯特变换解调的原理:从解析信号到辐射噪声包络谱
2.1 辐射噪声为什么自带周期调制:轴频与叶频的物理来源
舰船辐射噪声的主要成分是螺旋桨空化噪声,叶片旋转时,每个叶片周期性地经过同一空间位置,空化强度也随叶片方位角周期变化。于是接收到的噪声可以建模为宽带平稳随机信号与一个周期函数的乘积:宽带空化噪声充当载波,螺旋桨的旋转充当幅度调制源。调制基频是轴频,数值上等于螺旋桨转速乘以轴数;高次调制频率是叶频,等于轴频乘以叶片数。
这个物理机制直接决定了信号处理路径。空化噪声的能量集中在几百赫兹到几十千赫兹,而轴频通常在 1~10 Hz,叶频在 10~60 Hz,大型商船甚至更低。目标频率和载波频率相差三个数量级,直接在原始辐射噪声谱里找轴频几乎不可能。但包络信号是慢变的,频率范围正好落在目标特征所在的低频段。这就是为什么水声目标识别必须先做解调,再做谱分析。
2.2 解析信号与包络提取:希尔伯特变换解调这一步的数学含义
对带通滤波后的实信号 x(t),希尔伯特变换 H[x(t)] 可以构造解析信号:
z(t)=x(t)+j·H[x(t)]
z(t) 的模就是 x(t) 的包络:
a(t)=|z(t)|=sqrt(x²(t)+H²[x(t)])
这一步是严格意义上的正交解调,包络 a(t) 保留了调制信息,去掉了高频载波。得到包络后,去掉直流分量再做 FFT,就得到 DEMON 谱:
X_demon(f)=|FFT(a(t)-mean(a(t)))|
频率分辨率由帧长决定,Δf=1/T。要分辨 1 Hz 附近的轴频,单帧分析时长至少 1 秒,实际工程中通常取 4~10 秒。帧长取 4 秒时分辨率是 0.25 Hz,帧长取 10 秒时是 0.1 Hz,足以区分相邻目标的轴频差异。DEMON 这个缩写本身就是 Detection of Envelope Modulation on Noise,整条链路的每一环都对应名字里的单词:检测、包络、调制、噪声。
2.3 平方解调、绝对值解调与希尔伯特变换解调的差别
经典的模拟解调电路常用全波整流加低通滤波来提取包络,数字实现里也有人直接取 abs(x(t)) 或 x²(t) 再滤波。这两种做法在数学上并不精确。全波整流会引入载波频率二倍频的残余分量,平方解调则在包络频谱中额外产生谐波失真,对后续轴频、叶频的峰值检测会造成干扰。
希尔伯特变换解调在数学上是精确的包络提取,不产生额外谐波,也不改变包络频谱的相对幅度关系。代价是多一次 FFT 和一次逆 FFT 的计算。以 44.1 kHz 采样率、4 秒帧长为例,单帧希尔伯特变换在普通桌面 CPU 上耗时毫秒级,对 Desktop 端的实时或半实时处理完全不是瓶颈。这也是为什么这个标题把辐射噪声、希尔伯特变换解调、DEMON 三个词放在一起时,工程上最顺的路径就是先带通、再解调、后谱分析。
3. 用 Python 复现辐射噪声 DEMON 谱:解调参数与完整代码
3.1 输入数据与预滤波:希尔伯特变换解调前的带通选择
解调之前必须先做带通滤波。原因有二:一是辐射噪声的调制能量集中在特定频段,低频段流噪声和机械振动与螺旋桨调制无关,不滤掉会稀释包络中的调制成分;二是希尔伯特变换对全频带信号提取包络时,宽带噪声会把调制峰淹没在随机波动里。滤波的目标是只留下调制能量集中且调制深度比较大的频段。
import numpy as np from scipy.io import wavfile from scipy.signal import butter, sosfiltfilt, hilbert, windows fs, data = wavfile.read('radiated_noise.wav') data = data / np.max(np.abs(data)) # 归一化到 [-1, 1] sos = butter(4, [1000, 8000], btype='bandpass', fs=fs, output='sos') x = sosfiltfilt(sos, data)带通范围取 1~8 kHz 是水面舰辐射噪声的常见选择:低于 1 kHz 的部分包含流噪声和机械线谱,高于 8 kHz 的高频段水声传播衰减快,信噪比低。sosfiltfilt 做的是零相位滤波,正反各过一次,输出没有群延迟,避免包络时间对齐出问题。
3.2 希尔伯特变换提取包络并计算 DEMON 谱
滤波之后进入核心链路:分帧、降窗、希尔伯特变换取包络、去直流、FFT。
frame_len = int(4 * fs) # 单帧 4 秒,频率分辨率 0.25 Hz hop = int(2 * fs) # 50% 重叠 nfft = 4096 win = windows.hann(frame_len).astype(np.float64) demon_acc = np.zeros(nfft // 2 + 1) frame_count = 0 for start in range(0, len(x) - frame_len + 1, hop): seg = x[start:start + frame_len].astype(np.float64) * win env = np.abs(hilbert(seg)) # 解析信号的模 = 包络 env = env - np.mean(env) # 去掉直流分量 spectrum = np.abs(np.fft.rfft(env, nfft)) demon_acc += spectrum frame_count += 1 demon = demon_acc / frame_count freqs = np.fft.rfftfreq(nfft, 1 / fs)帧长 4 秒对应 0.25 Hz 的频率分辨率,能清晰分辨 1~10 Hz 的轴频。50% 重叠是为了提高弱调制瞬态的捕获概率,代价是计算量翻倍,桌面端依然可以接受。nfft 取 4096 而帧长是 176400 点,rfft 会自动对包络补零到 4096 点,等效于频谱插值。需要注意:补零只让谱线看起来更密,并不会提升真实分辨率,真实分辨率只由帧长决定。如果希望 DEMON 谱更平滑,优先增加帧长,而不是加大 nfft。
3.3 分段平均:处理辐射噪声非平稳性的关键
单帧 DEMON 谱的随机波动非常大。辐射噪声不是平稳信号:目标机动、航速变化、海况干扰都会让包络结构在秒级尺度上变化。直接拿一帧包络做谱分析,轴频峰可能被随机尖峰掩盖。整段数据分帧计算谱再做平均,本质上是 Welch 平均周期图法在包络域的应用,平均 N 次后随机波动幅度约降为原来的 1/sqrt(N)。
平均次数也不是越多越好。超过 20 帧平均会把目标机动引起的轴频漂移抹成宽带隆起,反而丢失识别线索。推荐的平衡点是 5~20 帧,具体取决于数据长度和目标航速稳定性。处理长记录时,建议先把数据按 30 秒一段切分,段内做 50% 重叠的谱平均,段与段之间保留独立谱输出,这样既压了随机波动,又保留了目标的时变信息。
| 参数 | 推荐值 | 作用与影响 |
|---|---|---|
| 采样率 | 不低于 44.1 kHz | 决定可分析的载波频带上限 |
| 带通范围 | 1~8 kHz | 低于 1 kHz 流噪声干扰大,高于 8 kHz 衰减快 |
| 帧长 | 2~10 秒 | 决定频率分辨率,轴频低时取长帧 |
| 重叠率 | 50%~75% | 提高弱调制捕获概率,增大计算量 |
| nfft | 大于等于帧长 | 只做插值,不提升真实分辨率 |
| 平均次数 | 5~20 帧 | 压随机波动,过多会抹平轴频漂移 |
注意:nfft 小于帧长时 rfft 会自动截断包络,导致频谱混叠。宁可让 nfft 大于帧长,也不要小于帧长。
4. 从 DEMON 谱到水声目标识别:轴频叶频特征的提取与判别
4.1 轴频、叶频与谐波结构:目标身份的指纹
DEMON 谱里最稳定的特征不是某个固定频率,而是谐波结构。轴频 f0 处出现基频峰,2f0、3f0 处出现谐波峰;叶频同样带谐波,且叶频与轴频之比正好是叶片数。这个比例关系由螺旋桨物理结构决定,航速变化只改变轴频的绝对值,不改变叶频与轴频之比。所以识别判据应当基于频率比例和等差谐波序列,而不是绝对频率匹配。
实际谱图上轴频的 2 次、3 次谐波往往比基频更突出,这是因为空化调制波形不是纯正弦,而是接近周期性脉冲串。找峰值时如果只挑最大峰,常会错把二次谐波当基频。正确做法是先找一组满足等差关系的谱线,频率间隔一致且间隔最小,那个最小间隔才是真正的轴频。
4.2 自动提取 DEMON 谱线峰值的代码
from scipy.signal import find_peaks # demon 为上一步得到的平均 DEMON 谱 max_val = demon.max() peaks, props = find_peaks( demon, height=0.3 * max_val, # 峰高阈值,低于 30% 最大峰的忽略 distance=int(0.5 / (freqs[1] - freqs[0])) # 最小间距 0.5 Hz,抑制旁瓣 ) peak_freqs = freqs[peaks] # 在 60 Hz 以下搜索轴频候选,做谐波匹配 cand = {} for fp in peak_freqs[peak_freqs <= 60]: harmonics = [k * fp for k in range(2, 7)] match = sum( 1 for h in harmonics if np.min(np.abs(peak_freqs - h)) < 0.5 # 容差 0.5 Hz ) cand[fp] = match shaft_freq = max(cand, key=cand.get) if cand else None blade_freq = None if shaft_freq: # 叶频候选:轴频整数倍且在 200 Hz 以内 blade_cands = peak_freqs[(peak_freqs > 3 * shaft_freq) & (peak_freqs < 200)] blade_freq = max(blade_cands, default=None)find_peaks 的 height 参数过滤低幅度杂峰,distance 参数保证同一个峰的旁瓣不会被当作独立峰。0.5 Hz 的容差覆盖了目标低速机动时的轴频漂移范围,如果数据来自航速稳定的目标,可以收窄到 0.1 Hz。叶片数估计直接取 blade_freq / shaft_freq 后四舍五入,商船常见 4~7 叶,这个数值本身就是一个强判别特征。
4.3 识别判据与混淆场景:轴频漂移和环境干扰
拿到轴频和叶频后,水声目标识别通常走一条层次化判据。首先确认谐波序列是否满足等差关系,这是区分螺旋桨调制与随机干扰的最有效手段。其次确认叶频与轴频的比值是否在合理叶片数范围,5~7 是商船典型值,鱼雷等高速小目标则明显偏高。最后做多段判决:同一目标连续多段 DEMON 谱的轴频应该平滑变化,环境噪声产生的随机谱峰各段位置不相关。
容易混淆的情况主要有三种。海浪拍击产生的包络调制频率在 0.1~0.5 Hz,低于轴频下限,加 1 Hz 高通即可排除。船上机械线谱是固定频率的窄带成分,不随转速平滑变化,与轴频的比例关系不稳定,可通过多航速数据验证。双目标同频是最难的情况,此时 DEMON 谱上会出现两组谐波序列相互交叠,需要结合波束形成分方位处理后再提取,这已经超出单一通道解调的范畴。
5. 解调频带选择、边界效应与链路验证的三个实战技巧
5.1 解调频带怎么定:调制密度优先原则
1~8 kHz 是通用起点,但不同目标的调制能量分布差异很大。低速水下目标可能低至 500 Hz 以下,高速水面艇的调制能量集中在 5~15 kHz。逐段试滤波频段效率太低,更快的做法是按 1/3 倍频程分带,每个子带分别做希尔伯特变换解调,计算包络的调制密度:
md = std(env) / mean(env)
取调制密度最大的子带作为正式解调频带。这个指标直接衡量包络中调制成分的相对强度,比主观试听和谱峰目测都客观。算出最优先导频带后,再以它为中心取相邻 1~2 个倍频程合并为最终带通范围。
5.2 把希尔伯特变换的边界效应压下去
hilbert 基于 FFT 实现,隐含着信号周期延拓的假设,帧的首尾样本在边界上不连续,会在包络两端产生振荡,幅度可达正常包络的十几倍,对 DEMON 谱的低频段污染尤其严重。压制手段有两个层次。第一是分帧时加汉宁窗,窗把帧两端压到接近零,边界不连续量随之消除。第二是在每帧解调后丢弃首尾各 5% 的包络样本再进入 FFT,代价是有效帧长缩减,频率分辨率略降,但对边界振荡的抑制非常彻底。两者同时使用时,丢弃比例控制在 5%~10% 即可,再多会浪费数据。
处理中场记数据时还要注意:如果用重叠分帧,窗函数叠加区的幅度补偿在包络域会引入额外的调制,建议先对去窗后的包络做同步补偿,或者直接取重叠帧中靠近中心的那一段包络参与平均,避免同一段数据被不同窗权重重复统计。
5.3 先仿真后真实数据的验证顺序
真实辐射噪声里干扰因素太多,链路调不通时很难判断是解调参数不对还是数据本身没有调制。标准做法是先用仿真信号验证整条链路,再上真实数据。构造一个调制深度 m=0.3、轴频 3.2 Hz、叶频 16 Hz(5 叶桨)的宽带噪声信号,带通到 1~8 kHz 后走一遍完整流程,DEMON 谱的 3.2 Hz 和 16 Hz 处必须出现尖峰,且 6.4 Hz、9.6 Hz 处有谐波;如果峰位不对,先查频率分辨率,再查边界效应。链路验证通过后,把同一套参数直接套到真实数据上,剩下的工作就只是调带通范围和平均次数了。这个顺序能把参数调试和代码错误分开,对 Desktop 端快速原型验证尤其有效。
本文还有配套的精品资源,点击获取