简介:针对水声目标识别中的噪声干扰与特征弱化问题,这份资源提供了一种基于自适应高斯平滑算法的MATLAB实现。算法根据声学信号的局部梯度或像素差异动态调整高斯核的尺寸与形状,以抑制水温、盐度、压力等环境因素引入的复杂噪声,同时保留边缘与目标细节,为后续短时傅里叶变换或小波变换等特征提取提供更干净的输入,从而提升目标识别的准确性和可靠性。压缩包内仅含1个m脚本,文件体积约1KB,代码精简却涵盖采样率转换、去直流偏置、梯度计算、自适应核构建及滤波操作等核心环节,便于读者从底层理解自适应平滑的完整流程,也可直接嵌入水声目标识别系统进行二次开发。目前已有423人学习,适合从事水下信号处理、目标识别相关研究的学生、工程师与算法开发者使用。
1. 水声目标识别里的预处理矛盾:降噪和保线谱为什么必须二选一
被动声呐的瀑布图(LOFAR图)上,目标船的辐射噪声在低频段拖出几条细亮的水平线谱,处理器想用高斯平滑压掉背景噪声,结果线谱也被抹成粗带;不压噪,模型又会被满屏的环境噪声和混响旁瓣带偏。这是水声目标识别落地时最常见的预处理矛盾,也是自适应高斯平滑算法的主场:让σ随局部梯度强度逐点变化,平坦噪声区开大σ猛降噪,线谱边缘区收小σ保特征。P37这条任务线里,我把这套预处理放在频谱图生成之后、模型输入之前,跑通的成本只有几行numpy计算。适合正在做被动声呐目标分类、想把频谱图增强写进识别pipeline的工程师参考,新手能照代码跑通,熟手可以拿着参数边界去调自己的数据。
2. 自适应高斯平滑原理:σ不再全局统一,让梯度给每个时频点发一个值
2.1 固定高斯核的短板:降噪和保线谱只能二选一
高斯平滑的本质是加权邻域平均,核函数里唯一能调的就是σ。σ决定邻域半径:σ越大,参与平均的时频点越多,随机噪声被摊得越薄,这是它作为降噪工具的全部底气。但目标线谱在频谱图上是沿时间轴延伸、沿频率轴只有一到两个bin宽的细结构,邻域一扩大,线谱就从"细线"摊成"粗带",相邻谐波彼此粘连,模型拿到手的判别特征已经变形。
水下目标辐射噪声的判别信息恰恰大多藏在这组线谱里。大型船舶的轴频通常在3到30赫兹量级,叶频是轴频乘以螺旋桨叶数,两者再拖出一串等间隔谐波。不同船型、不同推进工况的组合方式不同,这组线谱相当于频谱图上的"正文内容",而海洋环境噪声、混响旁瓣只是"纸张纹理"。固定σ的高斯平滑分不清正文和纹理:σ调大,正文被抹平;σ调小,纹理依旧满屏。一句话概括这个矛盾——降噪量和结构保持度由同一个参数控制,怎么调都是按下一个葫芦浮起一个瓢。
2.2 梯度当开关:σ值随边缘强度自适应缩放
自适应高斯平滑的出发点很朴素:先给每个时频点算一个"边缘强度",再按这个强度分配σ。边缘强度用梯度幅度近似,频谱图上线谱边界幅度跳变剧烈,梯度大;平静噪声区幅度起伏平缓,梯度小。于是σ的计算公式可以写成:
σ(p,q) = σ_min + (σ_max − σ_min) × (1 − g_norm(p,q))^α
g_norm是归一化到0到1的梯度幅度。g_norm接近1的强边缘位置,σ直接落到σ_min,几乎不动;g_norm接近0的平坦噪声区,σ顶到σ_max,全力降噪。α控制这条曲线的陡峭程度,α取1时线性过渡,取2以上时σ只在强边缘附近才快速收窄,中间地带保持较大平滑量。
这个公式解决的是加性噪声模型下的折中问题。P37数据上我用的参数范围是:σ_min取0.3到0.6个频率bin,σ_max取2.0到3.0个bin,α取2,梯度归一化截断在95分位。注意这里的最小单位是bin,不是赫兹——不同n_fft下每个bin的带宽不同,σ要跟着bin宽度换算,写死在赫兹上换一组参数就得重调。
提到等效实现,直接对每个像素做变核卷积代价太高,常见做法是分层高斯:把σ从σ_min到σ_max分成5到6档,对原图分别用固定σ平滑,每个像素按自己的σ_map给这几张结果分配加权系数再混合。6次固定卷积加一次加权求和,批量处理上百小时声呐数据的耗时完全可以接受,这也是后面代码里采用的做法。
2.3 和双边滤波、各向异性扩散怎么选
同样主打边缘保持的还有两位常客。双边滤波在距离权重之外乘了幅度差权重,实现简单,但在低信噪比频谱图上容易把孤立强噪点误判成边缘,平滑后留下盐粒状残留,需要再做一次中值收尾。各向异性扩散的边缘锐度最好,可Perona-Malik要调的迭代次数和梯度阈值两个超参对每段海况都得重来一遍,玄学成分偏高,放进批量流水线里难以把控。
相比之下自适应高斯是三者里最"工程友好"的:无迭代、输出确定可复现、超参就两个主值加一个梯度分位,每个参数都有明确物理含义。它的缺点是遇到弯曲粗边缘会出现轻微阶梯感,但水声频谱图以水平直线谱为主,这个缺点被数据形态天然绕开了。P37这条线里我把自适应高斯作为默认预处理,双边滤波只留作调试时的对照项。
3. 从WAV到自适应增强频谱图:可复现的P37预处理流水线
3.1 读进水声信号:采样率、分帧与STFT参数怎么定
水声原始数据通常是一段几十分钟的WAV,采样率从8k到96k都有,来自不同采集卡。处理前先统一重采样,线谱识别任务里16kHz足够覆盖轴频谐波所在的低频带,还能把后续FFT尺寸控制在合理范围。先读文件、做单通道和幅度归一化:
import numpy as np import scipy.io.wavfile as wavfile from scipy.signal import stft, resample_poly fs_in, data = wavfile.read("hydrophone_001.wav") if data.ndim > 1: # 多通道水听器,取第一个通道 data = data[:, 0] data = data.astype(np.float32) if data.dtype != np.float32: data = data / 32768.0 # int16 PCM 归一化到 [-1, 1] fs = 16000 data = resample_poly(data, fs, fs_in) # 重采样到16kHzresample_poly把数据按fs/fs_in的比值重采样,内部带抗混叠滤波,比直接线性插值可靠。注意wavfile.read返回的dtype:int16文件要除以32768,float32文件直接可用,写代码时加个判断比假设输入格式要省心得多。
STFT参数是这条流水线里第一次分叉。快速原型和线谱优先两套配置对比如下:
| 参数 | 快速原型 | 线谱优先 |
|---|---|---|
| fs | 48000 | 16000 |
| n_fft | 2048 | 32768 |
| hop | 512 | 4096 |
| 频率分辨率 Δf | 23.4Hz | 0.49Hz |
| 帧率 | 93.75帧/秒 | 3.9帧/秒 |
| 128帧patch对应时长 | 1.37s | 32.8s |
快速原型适合先验证流程,线谱优先才是真正能分辨轴频谐波间距的配置。轴频谐波间隔常在几十赫兹以内,Δf=23.4Hz时两根谐波粘在一起,自适应高斯平滑再怎么保边缘也无从保起。用线谱优先配置做STFT:
n_fft = 32768 hop = 4096 f, t, Zxx = stft(data, fs=fs, nperseg=n_fft, noverlap=n_fft - hop, window="hann", boundary=None) spec_db = 20.0 * np.log10(np.abs(Zxx) + 1e-8)noverlap=n_fft-hop让相邻帧有87.5%重叠,帧率3.9帧/秒对线谱的时间连续性足够。加窗选hann,旁瓣抑制好,避免强线谱能量泄漏到相邻bin,给后续梯度计算添乱。boundary=None去掉边缘补零引入的伪帧。spec_db转到dB域,这一步很关键,后面避坑章节会专门展开。
提示:如果测试数据里有带宽低于0.5Hz的极窄线谱,把n_fft提到65536,帧长4秒,代价是单帧计算量翻倍。P37数据上32768基本够用。
3.2 核心实现:梯度驱动σ_map,分层高斯加权混合
拿到dB频谱图后,自适应高斯平滑主体函数如下:
from scipy.ndimage import gaussian_filter def adaptive_gaussian_smooth(spec, sigma_min=0.4, sigma_max=2.5, grad_percentile=95.0, alpha=2.0, n_levels=6, blend_width=0.3): """对dB频谱图做自适应高斯平滑。 spec: (频率bin, 时间帧) 的二维数组 sigma 单位是bin,不是Hz;换n_fft时要按比例换算 """ spec = spec.astype(np.float32) gy, gx = np.gradient(spec) # 频率轴/时间轴的梯度 grad = np.sqrt(gy ** 2 + gx ** 2) # 梯度幅度 thr = np.percentile(grad, grad_percentile) # 分位数截断,防强干扰绑架 grad_n = np.clip(grad / (thr + 1e-8), 0.0, 1.0) sigma_map = sigma_min + (sigma_max - sigma_min) * \ (1.0 - grad_n) ** alpha # 分层高斯:每个固定σ平滑一次,按σ_map加权混合 levels = np.linspace(sigma_min, sigma_max, n_levels) out = np.zeros_like(spec) wsum = np.zeros_like(spec) for s in levels: w = np.exp(-0.5 * ((sigma_map - s) / blend_width) ** 2) out += w * gaussian_filter(spec, sigma=s, mode="nearest") wsum += w return out / (wsum + 1e-12) spec_enhanced = adaptive_gaussian_smooth(spec_db, sigma_min=0.4, sigma_max=2.5, alpha=2.0)逻辑分四步。第一步np.gradient同时算频率轴和时间轴的梯度,合成梯度幅度;第二步用分位数归一化把梯度压到0到1,取95分位的意思是只让5%的强边缘享受σ_min待遇,其余区域平滑量受控;第三步用上节公式生成σ_map;第四步分层混合,每个像素按自己对应的σ取邻近几档高斯结果的加权平均,比逐像素变核卷积快一个数量级,效果非常接近。
参数说明:sigma_min设为0.4个bin,强线谱边缘几乎不动;sigma_max设为2.5个bin,平坦区压噪明显。blend_width=0.3控制混合过渡带宽度,太小会在σ_map突变处留下条状伪影,太大则分级感消失,整个图糊成一团。mode="nearest"处理频谱图边缘,避免默认反射模式在高频截止处造出假梯度。整函数对4000×5000的频谱图耗时约两三秒,可以放心逐段批量跑。
3.3 平滑之后:频带裁剪、滑窗切帧与归一化
平滑完的频谱图还不能直接进模型,三件事要接着做。第一件是频带裁剪,目标线谱集中在低频,把0到2kHz之外的高频直接砍掉,既降计算量又避免模型被高频鱼群噪声带偏:
lim = int(2000 / (fs / n_fft)) # 2000Hz对应的bin数 ≈ 4096 spec_band = spec_enhanced[:lim, :] # 保留0~2kHz第二件是滑窗切帧。识别网络吃的是固定尺寸patch,线谱随时间漂移,窗口取128帧约33秒,既能覆盖足够多的谐波周期,又不会把航速变化拉进同一个patch:
def crop_patches(spec, crop_len=128, step=64): n_freq, n_time = spec.shape idx = range(0, n_time - crop_len + 1, step) return np.stack([spec[:, c:c + crop_len] for c in idx]) patches = crop_patches(spec_band) # (N, 4096, 128)第三件是逐patch归一化。每段录音的增益不同,直接喂CNN会让模型把绝对幅度当特征,归一化后只保留相对结构:
mean = patches.mean(axis=(1, 2), keepdims=True) std = patches.std(axis=(1, 2), keepdims=True) patches_norm = (patches - mean) / (std + 1e-6) patches_norm = patches_norm[:, np.newaxis, :, :] # (N, 1, 4096, 128)到这里预处理流水线闭环了:WAV进,形状为(N,1,4096,128)的张量出。顺便说一句,数据增强(SpecAugment的时频掩码、随机幅度增益)放在平滑之后做,不要放前面——先增强再平滑等于把增强造出的噪声又平滑掉了,增强白做。
4. 平滑频谱图怎么喂进识别模型:选型、训练与消融结果
4.1 模型选型:浅层CNN就够,ResNet不是必须
频谱图patch是单通道灰度图,判别信息集中在低频窄带上,结构比自然图像简单得多。P37任务里我默认先上一个四层卷积的浅网,而不是直接搬ResNet18。原因有三:patch尺寸4096×128比例细长,ResNet的224×224输入要求先缩放,缩放把频率方向的信息损失不少;浅网参数量小,万级patch的数据量下ResNet容易过拟合;浅网在CPU上也能跑推理,上船部署时不用换骨干。
import torch.nn as nn class LineSpecNet(nn.Module): def __init__(self, n_classes=5): super().__init__() self.features = nn.Sequential( nn.Conv2d(1, 32, 3, padding=1), nn.BatchNorm2d(32), nn.ReLU(), nn.MaxPool2d(2), nn.Conv2d(32, 64, 3, padding=1), nn.BatchNorm2d(64), nn.ReLU(), nn.MaxPool2d(2), nn.Conv2d(64, 128, 3, padding=1), nn.BatchNorm2d(128), nn.ReLU(), nn.AdaptiveAvgPool2d((1, 1)), ) self.head = nn.Linear(128, n_classes) def forward(self, x): x = self.features(x) return self.head(x.view(x.size(0), -1))两轮池化把4096×128压到1024×32,最后一层用AdaptiveAvgPool2d把空间维度聚合成1×1再进分类头,patch尺寸怎么换都能接。卷积核全部3×3、padding=1,感受野逐层扩大,第一层学局部线谱片段,第三层已经能看到谐波间隔这类全局结构。如果换数据集后patch高度低于512,记得同步减少池化层数,否则频率方向会被压过头。
4.2 训练配置与按航次划分数据
训练配置我固定用下面这套,换数据先不动它,优先调预处理参数:
| 配置项 | 取值 | 说明 |
|---|---|---|
| 优化器 | AdamW,lr=3e-4,wd=1e-4 | 小样本比SGD收敛快,wd压低过拟合 |
| 损失函数 | CrossEntropyLoss | 类别均衡时够用 |
| batch size | 64 | 单卡训练 |
| 轮数 | 60,早停在12轮无提升时 | 防止后段过拟合 |
| 学习率 | CosineAnnealing到1e-5 | 后期精细收敛 |
数据划分是这条流水线里最容易翻车的一步。水声数据天然按航次成组,同一航次里目标船、海况、噪声底高度一致,随机打散划分会让模型学到"认航次"而不是"认目标",验证分数虚高一截。正确做法是按voyage_id分组切割:
voyages = df["voyage_id"].unique() n = len(voyages) train_v = voyages[:int(n * 0.70)] val_v = voyages[int(n * 0.70):int(n * 0.85)] test_v = voyages[int(n * 0.85):] train_df = df[df["voyage_id"].isin(train_v)] val_df = df[df["voyage_id"].isin(val_v)] test_df = df[df["voyage_id"].isin(test_v)]注意切片前把voyages排序或固定随机种子,否则每次跑实验划分都不一样,前后指标没法对比。如果航次数太少(少于20个),可以放宽到按"天"分组,但绝不能按patch随机分。
注意:跨航次测试的F1一般比随机划分低5到15个点,这个差值不是模型不行,恰恰说明之前的划分在自欺欺人。
4.3 消融实验:不平滑、固定σ、自适应σ差多少
消融实验回答"平滑到底值不值得做"。P37类数据集上,我复现过的典型量级如下,具体数值随数据源浮动,但相对关系基本不变:
| 预处理方式 | 随机划分F1 | 跨航次F1 | 线谱保真度 |
|---|---|---|---|
| 无平滑 | 0.84~0.87 | 0.72~0.75 | 细线完整 |
| 固定σ=1.5 | 0.86~0.88 | 0.76~0.78 | 粗带粘连 |
| 自适应σ 0.4~2.5 | 0.88~0.90 | 0.80~0.83 | 细线完整,噪底下降 |
三个值得注意的点。第一,固定σ在随机划分下似乎还行,但跨航次F1比自适应低4到5个点,因为固定σ把线谱抹得粗细不一,模型学到的边缘形态过拟合了当前数据。第二,自适应的优势在跨航次测试里更明显,说明保留下来的线谱锐度是泛化特征而不是记忆特征。第三,消融时要保证除预处理外所有环节完全一致,我被自己坑过一次——换了预处理顺手改了随机种子,最后分不清提升来自平滑还是来自数据运气。后悔药是有的:所有消融共用一个固定种子列表,每组实验把五个种子的均值±方差报出来,比单次最高分可信得多。
5. 自适应高斯平滑常见问题与避坑记录:五个亲身踩过的坑
5.1 对dB频谱还是线性幅度谱做平滑
现象:用线性幅度谱做自适应平滑,线谱被抹得比背景还平,识别率没升反降。 原因:线性幅度下目标线谱只比噪声底高一点几倍,梯度幅度被宽带噪声的随机起伏盖过,σ_map把线谱区域当成平坦区,分配了σ_max;dB压缩把动态范围从几万压到几十,线谱相对噪声的凸起变得明显,梯度才真正反映结构边缘。 解决:一律在dB域做梯度计算、σ_map生成和平滑,整条链路别混域。如果最终要存成图片喂预训练模型,dB域数据线性映射到0到255即可,别再做指数还原。
5.2 时间轴和频率轴的σ共用一个值
现象:各向同性σ平滑后,相邻三根谐波粘成一条粗带,谐波间隔信息直接消失。 原因:频谱图两个轴物理意义不同。频率方向上线谱是窄脉冲,跨bin平滑等于直接加宽线谱;时间方向上线谱是持续横线,平滑只影响时间连续性,不破坏频率信息。各向同性处理把两者混为一谈。 解决:σ按轴分开设定。沿时间轴σ_t取2到4帧保持线谱时间连续,沿频率轴σ_f取0.4到1.0个bin。scipy的gaussian_filter(spec, sigma=(σ_f, σ_t))让两个轴独立设值,自适应版本里对频率梯度和时间梯度分别归一化、分别出σ_map,再合成二维σ_map。
5.3 梯度归一化被强干扰线谱"绑架"
现象:某段数据里混入一条其他船只的强窄带干扰,整张σ_map变小,自适应平滑退化成几乎不平滑。 原因:做max归一化时,单点极值把其余梯度全部压到0.01以下,σ_map在0到1范围内失去区分度,平坦区也拿不到大σ。 解决:用分位数归一化替代max归一化。grad_percentile取95,意思是最强的5%梯度才配享受σ_min,剩余95%梯度在正常区间内分布。如果数据里干扰稀疏且极强,分位数可以再降到90,让更多区域恢复平滑量。
5.4 随机划分数据集让模型学会了记住航次
现象:验证集F1到0.9以上,换一段新航次录音直接掉到0.6,让人怀疑模型是不是在背答案。 原因:声呐数据一天内采集的样本满足强相关性,同一航次内海况和目标工况几乎不变,随机划分等于把同一航次的patch一半放训练一半放测试,模型的捷径是记住航次级的声学指纹,而不是目标类级特征。 解决:按voyage_id分组做留出法,测试集必须来自模型从未见过的航次。4.2小节的具体代码可以直接抄,这里强调一个底线:报告指标时同时给随机划分和跨航次划分两个数字,后者才是能对外交差的。跨航次的模型,线上才不容易翻车。
5.5 平滑参数固定,模型对噪声风格过拟合
现象:A海区数据训的模型到B海区,F1掉10个点,σ调成B海区的又丢了A海区的成绩。 原因:σ_max固定等于在告诉模型"频谱图必然是某种平滑风格",换海区后噪声底宽度和线谱强度变化,同样的σ_map产生不同残余噪声形态,模型没见过这种风格就慌。 解决:把σ_max当成随机变量做预处理增强。每个batch从[2.0, 3.0]均匀采样一个σ_max,σ_min保持0.4不变,线谱的锐利度始终有保底,平滑力度却有变化。这个trick比加高斯噪声有效,因为它在预处理强度维度上做了数据增强,而高斯噪声只是往图上撒了一模一样的噪点。
6. 验证与进阶:证明平滑有效,并把σ_map变成模型的第二通道
验证一个预处理是否值得保留,我习惯三步走。第一步是看图说话,把原始频谱图、固定σ结果、自适应结果的LOFAR图并排打印,看线谱是否保持细锐、谐波间隔是否清晰、噪底是否下降,视觉上过不了关的预处理在特征上必然有损失。第二步是Grad-CAM热度图,把测试patch输入模型,看注意力是否落在谐波附近而不是噪声团块上;平滑如果有效,热度图会从散点状聚拢成沿谐波排列的条带。第三步是鲁棒性测试,往测试音频里混入-10dB白噪声重测F1,自适应平滑的模型通常比固定σ模型下降幅度小3到5个点,这是它能扛住新海况的最直接证据。
进阶一点的做法是把σ_map本身作为第二通道拼进输入。之前平滑是预处理,模型只能看到结果;把σ_map和频谱图堆叠成双通道输入,等于告诉模型"哪些区域是被强平滑过的、哪些区域保留了原始锐度",模型可以自行权衡两个来源的证据。我在P37任务线里的实验结果是双通道比单通道跨航次F1再高2个点左右,代价是输入通道从1变2,网络头两层的卷积参数略微增加。注意σ_map在训练和测试时必须用同一套归一化逻辑,否则这个通道的分布偏移会让模型产生新的过拟合,白忙一场。
我这几年做水声预处理的一个习惯是:每换一条数据源,先把不平滑的baseline跑出来,再往上加自适应平滑,两组指标差出来的部分才是这个预处理真正贡献的价值。否则哪天模型涨点,你都不知道该谢平滑算法还是谢数据集运气。希望帮到你。
本文还有配套的精品资源,点击获取