简介：面向瞬变电磁法（TEM）数据处理中的小波去噪需求，该MATLAB脚本适合地球物理勘探、信号处理方向的工程师与研究者，适用于矿产资源勘查、地下水探测及环境工程等场景。针对TEM原始信号易受环境随机干扰与系统误差影响的问题，脚本基于小波变换对观测数据进行去噪处理，旨在保留信号主要特征、提升深部结构成像精度。压缩包内含1个m文件，大小仅1KB，为纯源码形式，便于直接查看算法逻辑和参数设置。已有304人学习下载，适合需要快速上手TEM数据预处理或比较小波阈值去噪效果的读者。通过该脚本可了解从数据读取、小波分解到系数处理与重构的完整流程，并可根据实际数据特性调整软阈值、硬阈值等策略；清晰的结构也便于进一步扩展自适应阈值或贝叶斯去噪方法，是瞬变电磁信号降噪的实用参考。
1. 瞬变电磁信号里的小波去噪：从一条被噪声淹没的衰减曲线说起
瞬变电磁（TEM）测量里，信号最值钱的部分往往也最难看。发射电流关断后，接收到的二次场电压从几百毫伏一路衰减到几微伏，早期道在几微秒内完成主要形态变化，晚期道却要持续到几毫秒甚至更久，整条曲线的动态范围接近 100 dB。随机噪声一旦压过晚期信号，衰减曲线在做视电阻率换算时就会偏离实际地电断面。小波去噪恰好满足“不破坏早期突变、不抹平晚期趋势”的要求，它在时频平面上按尺度把信号和噪声分离，再按幅度阈值逐层裁剪，因此成为 TEM 数据预处理里的常用手段。mainTEM1_good 这类脚本名，前半段是处理入口标识，good 后缀基本对应“加了小波去噪后的处理版本”。下面这套流程把该环节最常见的做法拆开讲，包括小波基、分解层数、阈值规则和验证方法，可以直接套到自己的 TEM 剖面上。
2. 瞬变电磁信号为什么适合小波去噪：非平稳信号的时频分布
2.1 TEM 衰减曲线是天然的非平稳信号
TEM 观测到的二次场电动势 V(t) 不是稳态信号，而是发射电流关断后的暂态响应。一次场切断相当于给地下介质一个阶跃激励，接收线圈记录到的电动势随时间近似呈幂律衰减。早期道近地表信息集中在前沿附近，波形陡；晚期道深层信息体现在尾部，波形缓。这样一个从陡到缓、幅度跨几个数量级的序列，用固定窗傅里叶分析会顾此失彼。窗口短则频率分辨率不足，窗口长则早期和晚期的时域特征被强行平均。小波变换靠伸缩和平移两个参数扫描信号：分析高频分量时时间窗自动变窄，分析低频时自动展宽。这种自适应时频划分与 TEM 衰减曲线的形态相匹配，所以瞬变电磁数据处理里，小波去噪比普通数字滤波更常用。
2.2 傅里叶低通滤波的两个硬伤
许多处理流程早期都用傅里叶低通滤掉高频噪声，但它在 TEM 上有两个问题绕不开。第一是频带重叠。TEM 有用信号本身是宽频带的，晚期段的低频成分与低频人文电磁噪声频谱高度重叠，线性滤波器无论怎么切通带，都会把信号的低频尾巴连在一起削掉。第二是振铃。低通滤波器在信号陡变处会产生吉布斯振荡，TEM 早期前沿本身就是阶跃型陡变，滤波后的曲线在真实地层响应之前出现假起伏，这种假象在后期反演里很难识别。小波阈值去噪把“滤除”变成“按系数幅度剪裁”，不对整个频带一刀切，保留下的是超过阈值的时频原子，因此上面两个硬伤不会以同样形态出现。
表：傅里叶低通与小波阈值去噪的关键差异
| 对比维度 | 傅里叶低通 | 小波阈值去噪 |
|---|---|---|
| 频率分辨率 | 全局固定 | 高频时间分辨率高、低频频率分辨率高 |
| 时域定位 | 无 | 每个尺度对应特定时窗 |
| 信号前沿 | 被钝化，易产生振铃 | 边缘奇异点可以保留 |
| 噪声模型 | 隐含平稳假设 | 不要求平稳，按尺度分别估计 |
| 主要参数 | 截止频率、阶数 | 小波基、分解层数、阈值 |
2.3 分解、阈值、重构：小波去噪的固定骨架
所有常见小波去噪方法都沿用三步：先做小波分解，得到最后一层近似系数和各层细节系数；再对细节系数做阈值处理；最后逆变换重构信号。噪声接近白噪声假设时，它在各个细节层都有能量，而 TEM 信号的能量相对集中在幅度较大的少数系数上，因此把小于门限的系数置零或压小，重构后就是相对干净的信号。下面这段代码不完成去噪，只用来观察各层细节系数的幅度分布，这是决定小波参数前必经的一步。
% 观察分解后各层细节系数的分布 % v 为归一化后的 TEM 衰减电压序列 [C, L] = wavedec(v, 5, 'db4'); for j = 1:5 d = detcoef(C, L, j); fprintf('level %d: max=%.3f median=%.4f\n', ... j, max(abs(d)), median(abs(d))); end这里detcoef提取第 j 层细节系数，用median代替均值，是因为细节系数里可能混有真实信号的尖峰，后者会把均值抬高。如果level 1的max值比level 2高一个数量级，说明记录里还有明显高频脉冲，直接进阈值处理会留下脉冲残余。常见做法是先做一次中值滤波剔除孤立尖峰，再进入小波去噪。
3. 用 MATLAB 复现 mainTEM1 的小波去噪主流程
3.1 进入去噪前的两个数据准备动作
仪器输出的时间道和电压序列，通常需要先处理两个问题再进小波模块。一是剔除发射电流关断前的样本，只保留关断后的衰减段，否则一次场残留会生成很大的首点尖峰，小波分解后扩散到相邻时间道。二是处理 NaN 和溢出的负值。晚期道的负值可能是噪声，也可能是仪器动态范围溢出，直接保留会让细节系数出现假大值；常见做法是用两端有效点做线性插值，或把该道权重置零后做叠加。这两步做完再进入小波去噪，参数重复性会明显变好。
3.2 核心代码：分解、阈值、重构一个函数写完
下面给一个可直接调用的去噪函数，结构对应 mainTEM1_good 里最核心的那段小波处理。输入时间向量t和电压向量v，输出去噪后的电压vDen。
function vDen = temWaveletDenoise(t, v, varargin) % TEM 衰减曲线小波去噪 % t: 时间道向量 (s) % v: 电压向量 (V) % 可选参数: wname 小波基, nlevel 分解层数, thrRule 阈值规则 p = inputParser; addParameter(p, 'wname', 'db4', @ischar); addParameter(p, 'nlevel', 5, @isscalar); addParameter(p, 'thrRule', 'universal', @ischar); addParameter(p, 'thrFun', 'soft', @(x) any(strcmp(x, {'soft','hard'}))); parse(p, varargin{:}); opt = p.Results; v = v(:); t = t(:); v(isnan(v) | isinf(v)) = 0; % 补齐无效点 amp = max(abs(v)); % 记录幅度用于还原 s = v / amp; % 归一化到 [-1,1] N = length(s); [C, L] = wavedec(s, opt.nlevel, opt.wname); % 用第 1 层细节系数的 MAD 估计噪声标准差 cd1 = detcoef(C, L, 1); sigma = median(abs(cd1)) / 0.6745; % 阈值: 默认通用阈值 sigma*sqrt(2*lnN) lambda = sigma * sqrt(2 * log(N)); if strcmp(opt.thrRule, 'rigrsure') lambda = thselect(s, 'rigrsure'); end Cden = C; for j = 1:opt.nlevel cd = detcoef(C, L, j); cdD = wthresh(cd, opt.thrFun, lambda); % 软阈值收缩 Cden = wthcoef('d', Cden, L, j, cdD); % 写回该层细节 end vDen = waverec(Cden, L, opt.wname) * amp; end调用示例：vDen = temWaveletDenoise(t, v)，即使用默认的 db4、5 层、通用阈值、软阈值。代码中关键位置的作用如下表：
| 环节 | 代码位置 | 作用 |
|---|---|---|
| 归一化 | s = v/amp | 防止大动态范围影响小波变换的数值稳定性 |
| 噪声估计 | median(abs(cd1))/0.6745 | 鲁棒估计噪声标准差，不受少量信号尖峰干扰 |
| 阈值计算 | sigma*sqrt(2*log(N)) | 通用阈值，门限随点数增长而自适应 |
| 软阈值处理 | wthresh(cd,'soft',lambda) | 系数向零收缩，避免硬阈值带来的伪振荡 |
| 系数写回 | wthcoef('d',Cden,L,j,cdD) | 只更新第 j 层细节，近似层不受影响 |
| 重构还原 | waverec(...) * amp | 乘回幅度，恢复到原信号量纲 |
median(abs(cd1))/0.6745里的 0.6745 是正态分布四分位点倒数。用中位数绝对偏差（MAD）估计噪声标准差，比直接用std更稳，因为细节系数里会混有真实信号边缘，这些大系数会把方差估计显著拉高。
3.3 参数变化带来的可观察差异
nlevel取 4 和取 8 的差异主要在晚期段：层数少时晚期平滑不足，层数太多时近似层也被压缩，早期电压幅值会系统性下降。thrFun选'hard'时，wthresh只做置零不做收缩，去噪后的曲线在过渡处会出现小台阶；TEM 晚期信号幅值低，这种台阶会影响后续求导，所以默认用软阈值。如果记录里有明显工频干扰（50 Hz 及其谐波），需要先做窄带陷波再进小波模块，否则工频会集中泄漏到某几个细节层，让阈值估计偏离实际噪声水平。
4. 小波基、分解层数与阈值：TEM 去噪最值得调的三组参数
4.1 小波基：先看对称性和消失矩
MATLAB 小波工具箱里，正交小波族最常用的是dbN与symN。dbN的消失矩等于 N，symN的消失矩也等于 N，但对称性更好。消失矩越高，对光滑信号的压缩效率越高，代价是对跳变信号容易产生过冲；对称性越好，重构后的相位畸变越小。TEM 的早期是陡变，晚期是缓变，参数选两头都不如选中间档。
表：常用小波基在 TEM 去噪上的表现
| 小波基 | 消失矩 | 对称性 | 早期突变 | 晚期光滑 | 计算量 |
|---|---|---|---|---|---|
| db2 | 1 | 差 | 保留最锐 | 噪声残余多 | 最小 |
| db4 | 2 | 较差 | 好 | 尚可 | 小 |
| db8 | 4 | 差 | 轻微过冲 | 好 | 中 |
| sym8 | 4 | 近似对称 | 好 | 好 | 中 |
| coif3 | 6 | 近似对称 | 较好 | 很光滑 | 中 |
实际处理里，我一般从 db4 起步，因为它在早期和晚期之间比较折中。若晚期段去噪后仍有明显毛刺，换 sym8 往往比把分解层数加到 7、8 更有效。做多测点剖面时，统一用 sym8 全剖面处理，比逐条测线调小波基更容易保证成像结果的一致性。
4.2 分解层数：跟着采样点数走
分解层数决定最后一层近似系数的频带下限，也决定噪声被拆分到多少个细节层。TEM 单测点的时间道通常在 200～1024 之间，经验上取 4～6 层。N 少于 256 时取 3～4 层；N 在 1024 左右取 6 层就已经足够。层数过多并不会提高去噪能力，反而会把接近低频的信号也当成噪声压掉。判断方法很简单：用 5 层与 6 层各跑一遍，对比去噪后曲线的尾部，如果 6 层结果在尾部出现整体下压，说明已经压到信号本身，退回 5 层。
注意：层数上限受点数限制，
wavedec到达上限后会提示 level 超过wmaxlev，这时不是调层数，而是该考虑数据是否做过有效叠加。
4.3 阈值规则：压随机噪声选 universal，保晚期弱信号选 rigrsure
thselect提供universal、rigrsure、heursure、minimaxi四种规则。universal门限是sigma*sqrt(2*log(N))，对肉眼可见的毛刺最有效，但随着 N 增大门限抬高，晚期弱信号容易被压过头。rigrsure基于 Stein 无偏风险估计，门限相对保守，适合晚期信号幅值低、变化平缓的情况。heursure是两者结合，适合判断不了噪声强弱时先用一轮；minimaxi最保守，去噪幅度也最小。
更好的做法是分层用不同阈值。TEM 早期突变的信息集中在第 1、2 层，晚期噪声散布在更深层：
% 分层阈值: 前两层用保守阈值, 深层用通用阈值 for j = 1:opt.nlevel cd = detcoef(C, L, j); if j <= 2 lam = 0.7 * thselect(s, 'rigrsure'); % 保护早期突变 else lam = sigma * sqrt(2 * log(N)); % 压随机噪声 end cdD = wthresh(cd, 'soft', lam); Cden = wthcoef('d', Cden, L, j, cdD); end这里第 1、2 层的阈值缩小到 rigrsure 的 0.7 倍，主要是不让早期峰被削成圆顶；深层保持通用阈值，针对晚期毛刺做压缩。0.7是一个常用的缩放系数，范围 0.6～0.9 都可以试，最终以晚期段无毛刺、早期峰值不下降为准。
5. 去噪效果怎么验证：晚期段残差与分段去噪技巧
5.1 用晚期段残差量化压噪幅度
去噪效果的验收不能只看图像。TEM 曲线尾部通常已到噪声底，可以把这一段当作近似纯噪声来量化：
idx_late = t > 0.5 * max(t); % 后期噪声主导段 rms_before = std(v(idx_late)); rms_after = std(v(idx_late) - vDen(idx_late)); fprintf('晚期RMS: 原始=%.4e 残差=%.4e\n', rms_before, rms_after);std受个别尖峰影响大时，可以用 MAD 替代。残差标准差比原始标准差压掉一半以上，同时早期道峰值变化小于 5%，这套参数基本可以固定下来。只盯着平滑度调参数容易把真实地质响应一起抹平，所以“早期峰值”和“晚期残差”两个指标要同时看。
5.2 分段小波去噪：保护晚期道弱信号
整条曲线用一个全局阈值时，总是存在一个矛盾：阈值大了晚期毛刺压得干净但早期被削，阈值小了早期保持但晚期仍然毛茸茸。几个 TEM 流程里常见的做法是把曲线按时间道分成早、晚两段，各用不同参数处理后再拼接：
% 按时间道切分, 以 1 ms 为示例边界 idx_split = find(t >= 1e-3, 1, 'first'); d_early = temWaveletDenoise(t(1:idx_split), v(1:idx_split), ... 'nlevel', 4, 'thrRule', 'universal'); d_late = temWaveletDenoise(t(idx_split:end), v(idx_split:end), ... 'nlevel', 6, 'thrRule', 'rigrsure'); % 边界处取 6 个点做线性过渡拼接 w = linspace(0, 1, 6)'; v_final = [d_early(1:end-6); ... w .* d_late(1:6) + (1 - w) .* d_early(end-5:end); ... d_late(7:end)];切分点一般选在双对数曲线的拐弯处，而不是峰值附近，否则两侧斜率差会在拼接处被小波重构放大。1 ms 只是示例，实际要根据关断时间、发射频率和延迟时间道重新选。拼接后观察整条曲线是否存在台阶；若存在，把过渡区从 6 个点扩到 10～12 个点，同时让两段的软阈值参数保持一致，台阶通常就会消失。这个方法既能保住早期地质响应的锐度，也不亏晚期弱信号，是瞬变电磁数据日常处理里最常用的一套收尾方式。
本文还有配套的精品资源，点击获取
