MATLAB高分辨率时频分析：同步压缩小波变换与曲波变换分离混合地震数据

MATLAB环境中应用高分辨率二维时频分析方法——同步压缩小波变换与曲波变换在混合地震数据分离中的实践

做地震数据处理的人，大概都躲不过一件事：手里的资料不“干净”。野外采集回来的数据常常是多种成分搅在一起的混合地震数据——有效反射信号里叠着面波、随机噪声，甚至多源同时激发时不同震源的能量互相串扰。要把这些成分逐个拆开，时频分析工具是绕不开的基本功。可传统手段在分辨率上总差口气，直到我接触了同步压缩小波变换和曲波变换的组合方案，才算是找到了既能看清细节、又能稳定重构的分离路径。

这套方法本质上解决的是两件事：一是把频带上靠得极近的成分区分开，二是把空间上方向各异的波场分开。同步压缩小波变换负责前者，它在连续小波变换的基础上做了一步“重新分配”，把模糊的时频能量压回真正的瞬时频率位置，分辨率直接拉开一个档次；曲波变换负责后者，它有天然的强方向选择性，能像“滤波器组”一样按角度和尺度切分二维图像。两者接在同一个MATLAB流程里，对付混合地震数据分离效果很稳。这篇文章我会把整个思路、原理、代码实现和踩过的坑一起倒出来，做勘探处理、探地雷达或者任何涉及波动信号分离的朋友，都应该能用得上。

1. 问题从哪来：混合地震数据为什么要做高分辨率分离

1.1 混合地震数据的几种典型成因

我平时遇到的数据混合问题，至少有三个来源，处理思路侧重点各不一样，先说清楚才能对症下药。

第一种是环境噪声与有效信号的混合。野外采集时，风吹草动、工业干扰、车辆振动都会变成随机噪声，这些噪声在时频谱上表现为“撒了一地”的能量，低频段还有面波这种强线性干扰，能量大、频带宽，把浅层反射信号完全压在底下。

第二种是多源同时激发的混叠。现在高效采集经常做同时激发，两个震源间隔几秒轮流干活，记录上是两套波场叠在一起，要在时间域层面把它们分离开，靠的是两种波场之间的统计差异和某些可分离特征。

第三种是地质体的散射波与反射波混合。复杂构造区里头，绕射波、散射波混在正常反射里，常规处理把它们吃掉就丢失了裂缝和断点信息，硬留又会干扰速度分析。这类问题更像“信号里找影子”，对分辨率要求最高。

传统频域滤波遇到这三种情况都很尴尬。低截滤波压不住与有效波频带重叠的噪声，常规频带分离又把同相轴切成畸变的“面条”，更别说多源混叠时双方频谱几乎完全重叠，分无可分。

1.2 传统时频方法卡在哪

短时傅里叶变换估计大家都用过，它的缺点一句话概括：窗口一固定，时间和频率的分辨率就锁死了。想要时间定位准，频率必然变宽；想要频率分得细，时间窗就得拉长，瞬时事件就看不清。地震信号本身就是非平稳的，用固定窗去匹配变化多端的波场，先天吃亏。

Wigner-Ville分布分辨率确实高，但双线性变换会产生交叉项干扰。两个频率成分中间的“幽灵能量”搞得图面一团糟，后续做提取、分离时一并引入噪声，越处理越花。S变换把短时傅里叶和连续小波拧在一起，低频窄窗让频率分辨率提升不少，但代价是低频段时间分辨率很差，对地震道的深层走时定位不友好。

常规小波变换解决了窗自适应的问题，可它输出的也不是纯频率域——是尺度域。尺度与频率之间存在映射关系，却始终隔着一层，而且小波系数在时频平面上的能量“弥散带”比较宽，两个相近频率的高斯包络一旦靠近就糊成一个峰，这为后续做门槛分离埋下了隐患。

1.3 为什么这套组合能补齐短板

同步压缩小波变换的核心改进，是把连续小波系数沿尺度方向做一次“重排”。它先对每个时间-尺度格点估算瞬时频率，然后把该点所有能量往瞬时频率对应的位置集中，结果就是一张能量高度聚焦的二维时频谱。频率靠得再近，只要瞬时频率有差异，理论上就能被拉开；而且它是可逆变换，分离完成后原信号能近乎无损地重建出来，就冲这一点，已经比一堆“分析专用”的时频方法强了一大截。

曲波变换又是另一维度的补充。它是面向二维图像的稀疏表示，基函数呈细长的“针状”，在尺度上和方向上都能自适应匹配弯曲同相轴。对地震剖面而言，反射同相轴在局部表现为近似直线的纹理，曲波系数把它表示成极少数几个大系数，而面波、随机噪声和混叠波场往往落在不同的尺度和方向通道里，于是分离操作变成了“按方向通道丢系数”的干净活。

这两种变换一个管频率、一个管方向，组合起来正好是“高分辨率二维时频分析”这个题目的完整落点。我在MATLAB里把两者做成流水线，先后处理混合数据，分别解决不同维度的分离需求，整个过程参数少、逻辑直，改起来也快。

2. 两个核心算法在MATLAB里的落地

2.1 同步压缩小波变换的原理和代码实现

同步压缩小波变换我习惯先备一句话的公式记忆法：小波变换算出能量分布，再求相位把能量“拧”到真实频率上去。数学上，先用连续小波变换得到不同尺度下的系数，然后对时间变量做局部回归估计瞬时频率，最后把相同瞬时频率对应的系数加起来，得到同步压缩后的时频系数。整个过程等于是给模糊的小波谱做了一次“锐化”。

MATLAB里这个功能其实有内置函数可以直接用，别自己从底层造轮子。wsst就是同步压缩小波变换的主函数，输入一维信号和采样率，输出为高分辨率的时频稀疏表示。重构时对应使用iwsst。实测下来它的门槛参数和尺度设置都能直接传，不需要手动做复杂调参。

一个典型的调用片段像这样：

% 读取单道地震数据，x 为列向量，fs 为采样率 fs = 1000; % 采样率，单位 Hz x = data(:, 1); % 取第一道 % 执行同步压缩小波变换 [sst, fr] = wsst(x, fs); % sst 为时频系数矩阵，fr 为频率向量 % 查看时频图 figure; imagesc(t, fr, abs(sst)); axis xy; xlabel('时间 (s)'); ylabel('频率 (Hz)');

这里要提醒一句：输出给wsst的数据最好先做首尾端的去均值处理，否则边界上的微小直流偏置在压缩过程里被无限放大，会让最左边和最右边的时间边缘出现一整条高能量伪影。我一开始没注意，切了段带常量的道集进去，结果低频带全是红的，排查了半天才发现只是均值没去干净。

具体的门槛参数wsst可以直接指定，我用的是默认值'default'。如果后续要控制压缩程度，可以手动传阈值向量。该值控制的是有效瞬时频率估计的可靠程度，数值越小条件越苛刻，压缩就越激进，但对噪声也越敏感。经验上，SNR 低于 5 dB 的数据我会把阈值放宽，保住信号完整性再说。

2.2 曲波变换的原理和代码实现

曲波变换的直观理解，可以类比成“很多不同长度、不同角度的刷子从各个方向扫过图片”。它把二维图像分解成多个尺度和多个方向上的系数，每个系数代表一个局部区域在这个特定方向和尺度上的强度。地震剖面上的反射轴是平滑曲线，用少数几个方向的基函数就能拟得非常好；随机噪声则散布在所有方向所有尺度上，系数值小而均匀；两种波若方向不同，在曲波域天然分居不同的“抽屉”里。

MATLAB不默认内置曲波工具箱，但业内开源的曲波变换实现并不少，基本都是拆成两步：分解得到系数、重构合成图像。方向数、尺度层数、窗口大小三个参数，是控制整个分离效果的关键。

以下是我常用的一段剖面级曲波分解与重构示例：

% 读取二维剖面 M，假设已经完成道间归一化 [M, dt] = read_seismic_section('synthetic_mixed.dat'); % 正向曲波变换（伪代码，实践中请使用正确的工具包函数） [C, S] = curvelet_decomp(M); % C 为系数元胞数组，S 为结构描述 % 按方向通道做系数筛选 C_thr = C; for scale = 2 : length(S) % 第一层是粗尺度，一般不动 for ori = 1 : size(C{scale}, 3) % 保留该方向上能量较大的系数，清零其余 C_thr{scale}(:, :, ori) = hard_threshold(... C{scale}(:, :, ori), 0.4 * max(abs(C{scale}(:, :, ori)))); end end % 逆变换重建，得到分离后的重构剖面 M_rec = curvelet_reconstruct(C_thr, S);

这里有个值得花精力调的地方：方向通道对应的是同相轴的倾角区间。要想让曲波变换真正替你干活，分离之前得搞清楚目标信号的倾角范围。以反射波和面波为例，反射波通常近水平或稍有倾角，面波尤其是浅层面波往往形成大倾角、低频的强线性条纹，两者在曲波域中落在方向通道的极端与中间位置上，做系数筛选时不要全局统一阈值，应该按方向分组分别设定保留水平。

同时也要注意，曲波变换本身的基函数不具备频率紧致性，不同尺度的系数之间会有重叠。如果只做一次硬阈值就把尺度1和尺度2之间的连通性切断，重构时会出现“麻点状”纹理伪影。我后来在筛选之后补了一次小窗口的平滑滤波，效果明显改善，重构剖面的空间连续性回到可用状态。

3. 分离策略设计：两段式流程与参数配合

3.1 分离流程的整体架构

这套方法最具实战价值的一点，是把两个变换的优势拼接成了一条流水线。我的核心设计思想是先频率后方向：越靠前的处理越“粗”，只负责把主要能量成分拉开，避免杂成分干扰后续精细处理。

流程第一步是用同步压缩小波变换处理单道或小窗口内的信号，把数据拆成带限成分。因为同步压缩小波变换提供的是逐道、逐时刻的高分辨率时频系数，能够把频带高度重合但瞬时频率不同的信号成分分到不同频率通道。第二步是对分离后的目标频段剖面施加曲波域分解，按方向筛选系数，这一步负责把前面混在过去频段里的面波模态和随机干扰在空间条纹上进一步切开。

具体来说，我在模拟混合数据上做过如下拆分：

• 用wsst将含噪声混合地震道逐道分解，取有效频带范围（比如 20~80 Hz）内的系数做重构，得到第一步去噪后的道集。
• 在第一步重构剖面上取 f-k 谱确认剩余的同相轴方向分布，确定反射波与侧面波在频率-波数域中各自占据的扇区。
• 用曲波变换对剖面分解，根据上一步确认的方向信息保留反射波对应方向通道，对面波、散射干扰通道系数做出衰减。
• 逆曲波变换重构剖面，再与第一次结果做差比照，观察泄露量。

实际处理过程里，我几乎总是先跑一版默认参数作为基线，再针对数据特征迭代。这个思路建议做成按“分析-决策-处理-再验证”循环，能在很短时间内对齐数据本身的内容差异。

3.2 同步压缩小波变换的关键参数

同步压缩小波变换有四个核心参数最影响效果：VoicesPerOctave（每倍频程内的尺度数）、NumOctaves（倍频程数量）、频率轴采样密度、阈值参数。前面两个参数决定时频谱分辨率；分辨率越高，压缩越精细，但是计算量会指数级上涨。

实测经验是：VoicesPerOctave设到16以上就能得到很紧致的谱图，设置成32时边际收益已经很小；NumOctaves则依据目标频带跨度来定，地震数据有效频带低频到高频跨度最多两个数量级，设成8~10足够覆盖。

阈值参数是干分离时我最常微调的点。硬阈值会把系数中小于阈值的能量整个清零，适合分离强相干成分；软阈值连续萎缩，适合噪声压制的场景。做混合地震数据分离时我建议先用硬阈值把能量主干找出来，对剩下的低幅值成分再单独做一遍软阈值，两级组合能兼顾分离彻底性和细节保持。

这里提供一个参数配置的对照表，方便大家照着自己的数据情况做首轮设置：

3.3 曲波变换的关键参数

曲波变换的三个主要参数分别是：分解尺度层数、每层方向数、窗函数尺寸。分解尺度层数决定多分辨率的粒度，方向数决定倾角分辨能力，窗类型影响交叉泄露。设置原则是让基函数尺寸尽量匹配同相轴在横纵方向的长度展布比。

比如处理 512×512 的合成剖面，我用 4 层尺度，每一层的方向数设置为 16，即能够在 0~180 度范围内区分 11.25 度间隔的波场方向，足够分离水平反射波和倾角约 30~45 度的面波条纹。窗函数则选用各向异性较强的类型，它对方向的选择性更好，代价是各方向系数间有轻微旁瓣泄露。

与曲波变换系数筛选直接相关的还有模态幅值比问题。混合数据里有效信号与噪声幅值比往往在 1:5 到 5:1 之间波动，统一用百分比阈值容易照顾不到弱信号。我在筛选时采用了“方向分组自适应阈值”的做法：按目标方向通道的最大系数乘一个系数权重，非目标通道则执行均值减若干倍标准差的下限截断。这样既保证提取，又避免误伤。

3.4 质量验证指标怎么选择

分离做得好不好，不能光盯着剖面眼睛看，得有量化指标。我常用三个指标判断：

一是重构残差率。把分离后的各成分叠加起来与原始数据求差异，残差能量占比在 5% 以下算合格；若超过 10%，通常说明筛选太狠把有效信号削掉了。这个指标不仅是衡量恢复程度的标尺，也能反向指导门槛的设定。

二是时频谱聚集度。取同步压缩后目标成分的主脊，计算能量波动带的宽度，越窄说明时频聚集越好，分离越干净。具体计算可用主脊附近能量占整个成分能量的比例替代，省事且直观。

三是同相轴连续性指数。在分离后的剖面上沿着目标同相轴追踪，计算振幅跳变次数和相位翻转次数。指数高说明重构后的同相轴完整，反之说明损坏了数据连续性。这个指标虽然有点“软”，但对评价后续成像表现非常有效，尤其在波动方程正演数据上，它的区分度比残差率更高。

4. 实施中遇到的坑和排查日志

4.1 边界“滤波器裙边”效应

用wsst处理第一道和最后一道数据时，谱图边缘常出现强伪影。第一次跑合成记录时，边缘上突然多出低频横纹，我还以为是数据里藏了面波，排查半天才锁定是边界效应。

连续小波变换在信号两端的数据不足，瞬时频率估计会严重偏移，压缩后这些偏移能量被扔到低频段，形成裙边。解决方式就两步：一是先用编辑函数对称延拓数据，处理完再裁掉延拓段；二是对wsst输入前强制detrend去除趋势项，避免零频成分放大。

对称延拓的长度在边界长度的 20% 附近即可，太长计算浪费，太短延拓泄漏还会渗进有效区域。我最后固定用 64 个采样点延拓，512 点地震道数据分离开来边缘几乎看不到痕迹。

4.2 同步压缩的“虚假频率”问题

如果数据本身噪声强或者存在轻微的非线性调频，同步压缩过程可能把能量推到不存在的瞬时频率处，产生虚假频脊。排查日志里最经典的状况是：时频谱上多出一条不随时间变化、强度忽大忽小的恒定谱线。

这个问题源自瞬时频率估计时用了局部相位差分，相位差分对噪声敏感，一旦有局部幅值扰动，相位导数就可能出现在“真正的脊”两侧对称的旁瓣。解决的办法是用高斯窗平滑瞬时频率场，再做压缩。在wsst里可以通过改变尺度方向上的插值方式降低敏感度，但更直接的方法是先对原始道做一次轻度低通，削掉高频毛刺，再进同步压缩。

4.3 曲波域的系数交叉泄露

曲波域筛选时，目标方向通道与非目标方向通道之间会因窗函数的旁瓣发生交叉泄露。这个现象在强面波与弱反射同向重叠时尤其明显。我在某次处理时发现分离后的剖面上，目标反射附近有浅浅的斜纹残余，正是不该出现的面波旁瓣。

修复办法是：在系数筛选之前先对整张曲波系数图做一次“方向性中值滤波”，沿着垂直于目标方向的方向滑窗取中值，能有效衰减低幅值的旁瓣能量。滑窗宽度选成目标同相轴宽度的 3~5 倍，效果最好。

4.4 常见问题速查表

4.5 计算量的控制经验

同步压缩小波变换配合曲波变换整套跑下来，计算负担不小。512×512 剖面、32 voices、4层曲波尺度，单次完整分离在普通台式机上大约需要 15~20 秒。如果要处理几十个 shot gather，累积时间就很可观了。

我优化过两轮。第一轮是数据层面：每个 shot gather 单独分离，只分离含目标事件的时窗段，事后把分离结果拼接回原长度；第二轮是算法层面：曲波变换每层的方向数按倾角范围动态削减，水平构造为主的数据不需要每个尺度都配满 16 个方向，能减少约 35% 的时间开销。

还有个小技巧值得分享：如果多道数据的噪声背景相近，可以抽 20% 的道数做参考道集，用参考道集的曲波域系数分布统计出统一阈值，再应用到全部道集，能省下一次全道集统计的计算，且阈值空间更平滑，剖面的空间一致性反而更好。

5. 个人实操体会与下一步可以玩的扩展

这一部分没有太多高深理论，就说点实在的操作心得。

第一，这套组合方法的真正价值不在某一个变换有多强，而在“级联”的先后顺序。我试过先做曲波域分离再做同步压缩小波重构，结果显示面波压制效果差不多，但弱反射的保幅性明显变差。原因是曲波对幅度敏感的系数压缩在前，后续小波重构又把那些微弱的信号脊能量吸收走了一部分。顺序换回来之后，残差率降低了 2 个百分点，这道坎是几次试验才试出来的。

第二，门槛参数宁可从宽到紧，不要一开始就追求“干干净净”。单位数据处理时，我常常先用宽松阈值跑一遍看主成分，再逐步缩小阈值观察成分变化，直到找到一个“再收就伤信号”的临界点。用这个临界点作为产量参数，在整个数据集上得到的结果是最稳定的。这个调参习惯比任何默认参数都可靠。

第三，去多次试验替代参数寻优时，先固定同步压缩小波的尺度数不变，只动曲波的阈值。两个变换共同作用时，如果同时变动多个参数，很难判定哪个参数在影响结果。固定一个变量，逐个扫另一个，效率更高。

做这套流程至今，我还留了两个扩展方向没做完，但已经验证可行。一个是把同步压缩得到的瞬时频率曲面当成先验信息，反馈给曲波变换的方向筛选，让曲波域筛序具备频率引导能力；另一个是把全流程包成批处理函数后接入并行池，把大测线的分离时间从小时级别压到分钟级别。这两个方向都不算特别难，但这个领域数据形态千变万化，想真正吃透高分辨率时频方法，还是得一条测线一条测线慢慢磨，经验就藏在每次阈值调整和每张残余剖面里。