自适应高斯平滑：水声目标识别中降噪与保线谱兼得

简介针对水声目标识别中的噪声干扰与特征弱化问题这份资源提供了一种基于自适应高斯平滑算法的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返回的dtypeint16文件要除以32768float32文件直接可用写代码时加个判断比假设输入格式要省心得多。STFT参数是这条流水线里第一次分叉。快速原型和线谱优先两套配置对比如下参数快速原型线谱优先fs4800016000n_fft204832768hop5124096频率分辨率 Δf23.4Hz0.49Hz帧率93.75帧/秒3.9帧/秒128帧patch对应时长1.37s32.8s快速原型适合先验证流程线谱优先才是真正能分辨轴频谐波间距的配置。轴频谐波间隔常在几十赫兹以内Δf23.4Hz时两根谐波粘在一起自适应高斯平滑再怎么保边缘也无从保起。用线谱优先配置做STFTn_fft 32768 hop 4096 f, t, Zxx stft(data, fsfs, npersegn_fft, noverlapn_fft - hop, windowhann, boundaryNone) spec_db 20.0 * np.log10(np.abs(Zxx) 1e-8)noverlapn_fft-hop让相邻帧有87.5%重叠帧率3.9帧/秒对线谱的时间连续性足够。加窗选hann旁瓣抑制好避免强线谱能量泄漏到相邻bin给后续梯度计算添乱。boundaryNone去掉边缘补零引入的伪帧。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_min0.4, sigma_max2.5, grad_percentile95.0, alpha2.0, n_levels6, blend_width0.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, sigmas, modenearest) wsum w return out / (wsum 1e-12) spec_enhanced adaptive_gaussian_smooth(spec_db, sigma_min0.4, sigma_max2.5, alpha2.0)逻辑分四步。第一步np.gradient同时算频率轴和时间轴的梯度合成梯度幅度第二步用分位数归一化把梯度压到0到1取95分位的意思是只让5%的强边缘享受σ_min待遇其余区域平滑量受控第三步用上节公式生成σ_map第四步分层混合每个像素按自己对应的σ取邻近几档高斯结果的加权平均比逐像素变核卷积快一个数量级效果非常接近。参数说明sigma_min设为0.4个bin强线谱边缘几乎不动sigma_max设为2.5个bin平坦区压噪明显。blend_width0.3控制混合过渡带宽度太小会在σ_map突变处留下条状伪影太大则分级感消失整个图糊成一团。modenearest处理频谱图边缘避免默认反射模式在高频截止处造出假梯度。整函数对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秒既能覆盖足够多的谐波周期又不会把航速变化拉进同一个patchdef crop_patches(spec, crop_len128, step64): 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), keepdimsTrue) std patches.std(axis(1, 2), keepdimsTrue) 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_classes5): super().__init__() self.features nn.Sequential( nn.Conv2d(1, 32, 3, padding1), nn.BatchNorm2d(32), nn.ReLU(), nn.MaxPool2d(2), nn.Conv2d(32, 64, 3, padding1), nn.BatchNorm2d(64), nn.ReLU(), nn.MaxPool2d(2), nn.Conv2d(64, 128, 3, padding1), 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、padding1感受野逐层扩大第一层学局部线谱片段第三层已经能看到谐波间隔这类全局结构。如果换数据集后patch高度低于512记得同步减少池化层数否则频率方向会被压过头。4.2 训练配置与按航次划分数据训练配置我固定用下面这套换数据先不动它优先调预处理参数配置项取值说明优化器AdamWlr3e-4wd1e-4小样本比SGD收敛快wd压低过拟合损失函数CrossEntropyLoss类别均衡时够用batch size64单卡训练轮数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.870.72~0.75细线完整固定σ1.50.86~0.880.76~0.78粗带粘连自适应σ 0.4~2.50.88~0.900.80~0.83细线完整噪底下降三个值得注意的点。第一固定σ在随机划分下似乎还行但跨航次F1比自适应低4到5个点因为固定σ把线谱抹得粗细不一模型学到的边缘形态过拟合了当前数据。第二自适应的优势在跨航次测试里更明显说明保留下来的线谱锐度是泛化特征而不是记忆特征。第三消融时要保证除预处理外所有环节完全一致我被自己坑过一次——换了预处理顺手改了随机种子最后分不清提升来自平滑还是来自数据运气。后悔药是有的所有消融共用一个固定种子列表每组实验把五个种子的均值±方差报出来比单次最高分可信得多。5. 自适应高斯平滑常见问题与避坑记录五个亲身踩过的坑5.1 对dB频谱还是线性幅度谱做平滑现象用线性幅度谱做自适应平滑线谱被抹得比背景还平识别率没升反降。 原因线性幅度下目标线谱只比噪声底高一点几倍梯度幅度被宽带噪声的随机起伏盖过σ_map把线谱区域当成平坦区分配了σ_maxdB压缩把动态范围从几万压到几十线谱相对噪声的凸起变得明显梯度才真正反映结构边缘。 解决一律在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跑出来再往上加自适应平滑两组指标差出来的部分才是这个预处理真正贡献的价值。否则哪天模型涨点你都不知道该谢平滑算法还是谢数据集运气。希望帮到你。本文还有配套的精品资源点击获取