去年处理某工区一批多源同时激发采集的混合地震记录时我被一个看似简单的问题卡了很久面波、浅层折射和深层反射在时间-频率平面里基本重叠常规带通滤波一压面波是没了但浅层反射的弱同相轴也跟着断了。后来我尝试把同步压缩小波变换SST和曲波变换Curvelet组合起来做两级分离效果才明显好转。这篇文章就围绕这套方法讲清楚它为什么有效、怎么在MATLAB里落地、以及实际工程里最容易踩的几个坑。适合正在做混合炮分离、面波压制或复杂波场成像预处理的地球物理从业者和相关专业的研究生参考。先说清楚标题里的“二维时频分析”到底指什么。在这套思路里SST负责把单道地震信号映射到一张高分辨率的“时间-频率平面”上这是一维信号在二维平面上的表示曲波变换负责对二维剖面本身进行多尺度、多方向分解相当于在“空间-方向平面”上做分析。两者合在一起一个管时间方向上的频率区分一个管空间方向上的几何形态区分恰好互补。文章不会堆公式重点放在原理直觉、参数设置和可复现的代码骨架希望你看完能直接在自己项目里试。1. 先从“混合地震数据”说起为什么常规方法被逼到墙角1.1 混合记录从哪来混合地震数据最常见的一种来源是多源同时激发。为了提升采集效率野外会把多组震源放在不同位置同时或准同时激发所有能量叠加到同一张记录里单炮之间互相“串扰”。这种混合炮在时间域、频率域、空间域上都叠在一起常规按炮分离的思路完全失效必须先做一次数据分离把所有被污染的炮集还原成干净的拟单炮记录。另一种更容易遇到的场景是单炮内部的强干扰与弱信号混叠。面波能量通常能超过反射波几十上百倍而且以低频、低速、带频散的方式出现在记录上浅层折射和高速噪声也常常横在有效波上面。这里头的难点是面波和反射波的主频范围往往是部分重叠的时间窗也重叠你要是直接高通滤波面波是压掉了浅层反射的主频也丢了一半。这两种情况本质是同一个问题不同波场在时间-频率-空间三个维度上耦合分离的目标是尽量保全弱信号、压掉强干扰而单纯在时间域或频率域取一刀没法同时满足这两个要求。1.2 传统STFT和小波为什么不够短时傅里叶变换用固定窗长去截信号窗一固定时间分辨率和频率分辨率就只能此消彼长。对地震信号这种时变非平稳信号来说浅层反射到深层反射的频谱变化很快固定窗根本追不上。更麻烦的是混合数据里两个分量的谱峰挨得很近时STFT的谱会把两个峰涂成一团你想在两个峰之间找一条干净的分割线基本做不到。连续小波变换比STFT好一些它的窗长会随着尺度自动伸缩低频看长窗、高频看短窗符合地震信号的物理习惯。但小波系数的能量仍然会在时间-频率平面上发散存在“涂抹”现象。两个波场如果频率差不过半个倍频程时频谱上的两个脊线就会贴在一起做掩模时很容易互相串扰。这个“差一点就能分开”的状态最让人难受。我当初实际处理时用STFT和CWT都试过分离结果要么是面波残留过多要么是弱反射被削得没法做后续速度分析。所以后来才转向同步压缩小波变换想看看能不能把时频谱的脊线“收紧”让两个分量真正在平面上分开。2. 同步压缩小波变换把时频模糊“挤”回脊线上2.1 SST到底在压缩什么同步压缩小波变换的全名听起来吓人核心思想其实很直接。连续小波变换算完之后每个时间-尺度点上的系数不仅幅度有意义它的相位沿时间方向求导还能得到一个局部瞬时频率估计。SST做的就是把那些瞬时频率估计接近的系数沿着频率轴重新“堆积”到它们各自的脊线位置上去。打个比方CWT像一台老式投影仪对焦不准光斑散得到处都是SST相当于根据每一束光的方向信息重新把光收回合适的针孔里。针孔越细时频谱上的两个邻近分量就越容易分开。对于地震记录里常见的调频波、衰减波、频散面波这个操作非常有效因为它们本身的瞬时频率变化是有规律的系数重排后能量高度集中。这个“把散掉的能量集中回来”的过程就是标题里“高分辨率二维时频分析”的核心含义。SST最后得到的矩阵仍然是时间×频率的二维谱只是这个谱的分辨率远高于普通CWT谱。2.2 MATLAB中wsst的用法边界与参数取舍MATLAB从R2016b开始提供了原生的同步压缩小波变换函数wsst不需要额外装工具箱调用方式很直接% 输入一段单道数据 x采样间隔 dt 秒 [sst, f] wsst(x, amor, VoicePerOctave, 12);返回的sst是频率点数×时间点数矩阵f是对应的频率向量。‘amor’代表复Morlet小波对地震这种带瞬态特征的信号比较友好因为它有较好的频率局部性。VoicePerOctave决定每个倍频程内插多少个“声音”说人话就是频率轴的细分程度。我在实际项目里的经验是VoicePerOctave取8到16之间比较稳妥。取太小时频率轴太粗面波和反射的脊线还是粘着取到32甚至64时时频谱上会出现很多细碎的伪峰尤其是低频端信号本身信息量不足纯粹是插值插出来的纹理。你在合成数据上看到细密纹理时先别高兴那多半不是高分辨率是把噪声也一并“锐化”了。wsst是单道处理函数整炮数据要逐道循环。所有道的f向量是一致的这给后面构建立体掩模提供了方便。需要注意逆变换函数isst的使用重建时小波类型和VoicePerOctave要和正变换保持一致否则重建能量对不上xrec isst(sst, amor, VoicePerOctave, 12);从实践看isst重建不是百分之百无损残余误差通常在几个百分点动态范围大的记录误差会更明显。所以分离流程里一定要有残差处理环节不能天真地以为逆变换出来就是干净分量。2.3 从单道时频到剖面二维时频体单道SST只能告诉你“这一道在什么时间有哪些频率成分”但混合数据分离需要知道整个剖面的规律。所以我会把每一道都算一遍SST堆成一个三维数据体三个维度分别是时间、频率、道号。拿到这个三维体后有两个重要的观察方式。第一个是固定某个时间窗看频率-道数切片能直接看到面波的能量在低频段横着铺开反射波的能量在中等频率段呈双曲状分布两者在切片上往往能看出不同的空间走势。第二个是固定某个频率切片看时间-道数切片等于做了一个时变带通能看到该频率分量在空间上的连续性强不强。分离的线索就藏在这个三维体里面波集中低频段反射波跨较宽频段面波空间连续性强但视速度低反射波在空间上也有连续性但斜率不同。先用SST谱找到两类波在时间和频率上的主导区域生成掩模再用曲波变换去处理空间方向上的残留形成两级分离。这个三级结构是整个流程的地基。3. 曲波变换从“逐道时频”到“空间几何”的补位3.1 面波与反射波在曲波域的差异化表征曲波变换是一种多尺度几何分析工具它和普通小波最大的区别是“各向异性”。小波用方形窗口去“看”图像看到的是一个个方块曲波用细长条形的基函数去匹配图像里的曲线边缘基函数本身带有方向。这套设计是为了能稀疏表示方向性强、弯曲但局部可近似为直线的结构。地震剖面其实是典型的“曲线纹理图像”反射同相轴是弯曲的曲线面波是低速、大倾角的斜条纹随机噪声则是没有方向规律的散点。把这些结构送到曲波域后面波的能量会集中在少数几个尺度、少数几个角度的子带里反射波因为倾角范围更宽散布在多个子带随机噪声对各个子带几乎均摊但幅度小。这种差异天然适合做方向滤波。你在曲波域里保留某些尺度角度、衰减另一些尺度角度就同时完成了带通滤波和方向滤波比单纯在f-k域做扇形滤波要细腻得多。f-k域滤波在强能量面波附近容易产生空间截断伪影曲波域由于是多尺度多方向分解伪影被分摊到各个子带视觉上干净不少。3.2 MATLAB环境下曲波实现的选用MATLAB本身没有内置曲波变换函数通常要借助开源的曲波变换工具箱。这类工具箱的核心函数一般分为正变换和逆变换两部分我习惯统一封装成curvelet_forward和curvelet_inverse两个函数方便替换不同实现版本。返回的系数是嵌套元胞结构大致是C{j}{a}表示第j个尺度、第a个方向角的系数矩阵。有一点必须提前说这类工具通常对图像尺寸有要求边长最好是2的幂或者至少是偶数且满足某些整除条件。地震剖面动辄上千道、上万采样点直接扔进去会有边界填充问题。我的做法是先对剖面做对称延拓到nextpow2尺寸变换完再裁剪回来。延拓宽度一般取30~50个样点既能压制边界效应又不会带来太多计算量。系数是复数里面有幅度也有相位。你在做阈值时如果只保留幅度相位信息会丢失重建出来的波形会“噼里啪啦”地乱跳。正确做法是保留原系数只把幅度小于阈值的那些位置置零。后面代码部分还会强调这一点。3.3 为什么单靠Curvelet不能一步到位那能不能跳过SST直接用曲波变换分离混合数据我试过效果不理想。原因是混合数据里两类波如果频率和时间窗重叠得很厉害单纯靠空间方向的差异很难把它们分开。比如多源同时激发产生的两炮记录波场形态几乎一样方向也相似唯一区别可能是激发时刻或频谱内容这些信息在曲波域里被摊到了无数个尺度角度子带中反而看不清楚。曲波变换擅长对付“空间形态不同的波场”但不擅长对付“时间线上挨得很近、频率上互相纠缠”的暂态信号。SST正好相反它沿着每一道做时频锐化能看出激发时刻、频率调制趋势这些时间维度的差异但对同一时刻来自不同空间方向的波场无能为力。一个管时间频率一个管空间方向两个合起来才有完整的信息维度。4. 两级级联分离流程SST掩模先行Curvelet兜底4.1 整体数据流的设计逻辑我最后定的流程分五步每一步都有明确的职责边界读取混合剖面组织成nt×nx的矩阵逐道计算SST得到三维时频体。在SST三维体上检测能量脊线为面波和反射波分别构建时频掩模逆变换得到两个初步分离分量。初步反射分量里通常还残留一些方向性强的面波尾巴把它和初步面波分量一起送入曲波域。在曲波域里按尺度、角度做阈值收缩把两类分量里的方向性伪影和噪声进一步剥离。输出最终面波分量、反射分量以及残差残差再反馈回反射分量保证弱信号能量不丢失。这个顺序是反复试出来的。先做SST是因为它能最有效地利用时频差异做粗分离粗分离之后信噪比已经上来了此时再做曲波方向滤波阈值参数更容易定不容易把弱反射误伤。反过来先做曲波强面波和弱反射在空间形态上如果相似曲波会把弱反射也当成面波一起压掉。4.2 SST时频掩模的构建与软/硬掩模选择掩模构建是整个第一级分离成败的关键。最简单的做法是画一条频率阈值线低频算面波、高频算反射但这在强噪声和频散面波面前太粗糙。我的做法是先做脊线检测对三维时频体中每一个时间点在频率维度上搜索局部极大值把幅值超过全局能量1%的极大值点连成脊线。脊线找到之后以每条脊线为中心开一个频带窗口。面波脊线的频带窗口用相对带宽10%~20%反射脊线窗口用相对带宽20%~30%因为反射波频带本身就比面波宽。窗口之间如果有重叠重叠区域的归一化幅度按距离加权谁离脊线近倾向谁。这样构建出来的掩模是软掩模边缘有渐变过渡逆变换重建时不会因为硬截断而产生类似吉布斯现象的震荡。硬掩模我也试过二值化的优点是处理速度快缺点是重建结果在掩模边界附近有轻微振铃在强反射区域尤其明显。软掩模多花一点点计算时间但剖面干净程度提升一个档次值得。4.3 Curvelet域阈值收缩的具体参数第二级曲波处理的核心是阈值收缩。我使用的阈值参考了经典的高斯噪声阈值公式thr c * sigma * sqrt(2 * log(N))其中N取该尺度子带系数总个数sigma用子带系数的中值绝对偏差估计sigma median(abs(coeff(:) - median(coeff(:)))); % 再除以0.6745归一化c是经验调节因子。第一轮处理时c取1.5~2.5主要目的是把面波强能量子带压下去第二轮迭代时降到0.5~1目的是保护弱反射。尺度之间也要做区别低尺度的低频大尺度子带面波能量集中阈值适当提高高尺度的高频细节子带有效信号强度弱阈值必须放低否则反射同相轴的细节纹理都被抹平。还有一个容易被忽略的参数是角度数。角度数太少方向分辨能力不够面波残留会比较严重角度数太多每个子带包含的系数太少统计阈值不稳定重建结果出现“钉板”纹理。我常用的配置是4到6个尺度16到24个角度先在一个小范围数据块上试再决定最终值。4.4 主流程MATLAB代码框架下面这段代码是浓缩后的骨架实际项目里建议拆成函数文件。用小波同步压缩函数时以MATLAB R2020b及以后版本为准不同版本函数签名略有差异。% 输入 data: nt x nx 混合单炮记录 % 采样间隔 dt 0.002 s nt size(data, 1); nx size(data, 2); % 第I级逐道SST构建三维时频体 vpO 12; % 先跑一道拿到频点数量 [~, f] wsst(data(:, 1), amor, VoicePerOctave, vpO); nf length(f); sstCube zeros(nt, nf, nx, single); % 单精度省内存 for ix 1:nx sst wsst(data(:, ix), amor, VoicePerOctave, vpO); sstCube(:, :, ix) single(sst); end % 注意 wsst 返回矩阵是 nf x nt预分配时按 nt x nf 存循环里转置实际写的时候要把维度对齐看清楚wsst返回的矩阵是nf×nt我预分配成nt×nf×nx所以存的时候要转置。每条道算完SST后f向量完全一致这可以放心。接下来构建脊线掩模和逆变换maskReflect zeros(nt, nf, nx, single); for ix 1:nx sst squeeze(sstCube(ix, :, :)); % 重新变回 nf x nt % 简单起见低频段置为面波掩模中高频段置为反射掩模 lowIdx f 15; % 面波主导区 midIdx f 15 f 80; % 反射主导区 mask zeros(size(sst)); mask(lowIdx, :) 0.3; mask(midIdx, :) 0.85; mask(f 80, :) 0.5; maskReflect(:, :, ix) single(mask); % 逆变换得到反射分量 dataReflectPre(:, ix) isst(sst .* mask, amor, VoicePerOctave, vpO); end dataSurfacePre data - dataReflectPre;这是一个高度简化的雏形真正的脊线检测和软掩模逻辑比这个复杂得多但方向就是这样。先用低频置弱、中高频置强的方式跑通全流程再去细化掩模生成算法。第二级曲波处理框架如下% 第II级曲波域方向阈值收缩 function C curvelet_shrink(dat, c0, numIter) datPad padarray(dat, [32 32], symmetric); C curvelet_forward(datPad); for iter 1:numIter for js 1:numel(C) for ja 1:numel(C{js}) coeff C{js}{ja}; med median(coeff(:)); sigma median(abs(coeff(:) - med)) / 0.6745; if isempty(sigma) || sigma eps continue; end N numel(coeff); % 大尺度低频子带阈值高高频细节子带阈值低 if js 2 c c0; else c c0 * 0.4; end thr c * sigma * sqrt(2 * log(N)); % 保留相位只用幅度判断置零 C{js}{ja} coeff .* (abs(coeff) thr); end end end datClean curvelet_inverse(C); datClean datClean(33:end-32, 33:end-32); end % 对初步反射分量再做一次曲波处理 dataReflect curvelet_shrink(dataReflectPre, 2.0, 2); % 残留低能量再补回反射分量 dataResidual data - dataReflect - dataSurfacePre; dataReflectFinal dataReflect 0.6 * dataResidual;这个流程跑通之后面波分量就是dataSurfacePre经过曲波处理后的结果反射分量就是dataReflectFinal。残差按比例补回的系数要小心补多了面波又回来补少了弱反射继续丢。我后面的参数敏感性部分会具体讲怎么调。5. 验证与评测合成数据和实测记录我们都试了什么5.1 合成混合记录怎么造为了验证方法有效性我合成了一张混合剖面参数尽量贴近陆上资料。反射波用主频30Hz的零相位雷克子波沿双曲线走时排列面波用主频10Hz的相速度频散信号速度900m/s在剖面上呈大角度斜纹最后加了高斯白噪信噪比定在6dB左右。生成方式不复杂关键是把两个波场的时频特征做成真实可区分的状态否则测试就没有说服力。反射波和面波在时间窗上大面积重叠传统带通滤波只能分掉一部分这正是我想模拟的困难场景。5.2 关键指标与评价结果评价分离质量我用两个指标。第一个是信噪比改进量定义成SNR 10 * log10( sum(ref^2) / sum((ref - sep)^2) )其中ref是干净反射波sep是分离出的反射分量单位dB。第二个是时频特征保持度我直接看分离后的SST谱和干净反射的SST谱在主导频段的相关系数相关系数高说明弱信号细节保得好。合成数据上的对比结果如下表方法分离后SNR/dB时频相关系数面波残留程度原始混合数据6.10.52严重单用带通滤波9.80.67中等单用SST掩模分离13.50.83少量单用曲波方向滤波10.20.71中等SST曲波两级分离17.90.92极低可以看到单用曲波效果甚至不如单用SST这印证了前面的判断混合数据里的关键差异在时频维度空间方向滤波只能作为补充。两级联合作业的收益不是简单相加而是先解决主要矛盾再处理残余矛盾。5.3 实测记录的情况合成数据过了关我又在某项目组的实际采集记录上做了验证。那段记录是多源同时激发方式采集的单炮表面上看满屏都是斜纹干扰浅层有效反射几乎被淹没。用这套流程处理后面波斜纹明显压平浅层反射同相轴能连续追踪了后续速度分析道集上低速能量的影响也小了很多。实测和合成的最大差异在于实测记录的波场形态更复杂面波不是一条干净的斜线而是多模式、变速度、带散射的宽带结构。SST掩模在第一级分离时偶尔会把强散射能量的尾巴误判成有效信号这个时候曲波的第二级“兜底”作用就显得特别重要它能把那些空间形态异常的碎片单独揪出来。6. 绕不开的坑与工程落地建议6.1 wsst参数的常见天花板VoicePerOctave并非越大越好这个前面提过我这里再强调一个细节当你把VoicePerOctave提到24以上低频端的频率轴加密会带来伪脊线原因是信号低频成分本身没有足够周期去支撑瞬时频率估计。如果检测到低频出现密密麻麻的细线先别急着调高参数而是把参数降回16试试。边界效应值得一提。wsst在信号两端各有一段时间长度与小波支撑有关的影响区这段区域里计算出的瞬时频率不可靠掩模打上去会失真。我在项目里的处理是先对每道做镜像延拓大约64个样点算完SST和逆变换后再裁掉效果比在原始数据上硬算好很多。还有一个容易出问题的地方是逆变换重建误差。isst重建出来的信号和原信号不是完全相等的误差在强能量段附近尤其明显。因此分离之后一定要算一下data - dataReflect - dataSurface把残差单独画出来看如果残差里还有规律的斜纹说明掩模泄漏了需要把掩模边界调得更紧。6.2 Curvelet阈值与迭代的联动风险曲波阈值参数里最坑的是“高阈值清得干净但伤弱信号”。我在合成数据上试过把c直接设成3.5面波确实完全消失但反射波的主频也被削掉1/4时频相关系数掉到0.7以下。原因很简单弱反射在曲波域里幅度本身就很低全局阈值一刀切它和噪声一起被归零。后来我把策略改成“大尺度高阈值、小尺度低阈值”并且分两轮迭代。第一轮用较高的c02.5直接压面波第二轮把c0降到0.8只会处理第一轮剩下的残差。这样弱反射的损伤大幅下降合成数据上SNR又多提升了2dB左右。还有一个我吃过亏的细节阈值操作不要用abs(coeff) thr后取abs(coeff)重建必须保留原始复数系数。我第一次图省事把系数模值置零再逆变换结果重建波形相位全乱剖面上出现大量随机抖动。6.3 生产级执行的效率与测试建议整炮数据逐道算SST非常慢。我处理一炮360道、每道6000样点的数据单核跑要将近一小时生产上完全没法接受。后来改成parfor并行循环12核的机器跑十到十五分钟能出结果提升明显。如果你有GPUwsst本身没有很好的GPU加速路径但曲波变换的正反变换可以考虑用GPU实现能省下第二级的大部分时间。内存也要精打细算。三维SST体如果用double精度存储大规模工区的数据量可以达到几十GB。建议只保存掩模和低频段投影或者把整炮按滑窗切成若干重叠子块分别处理处理完再拼回去。滑窗接缝处用余弦权重叠加避免出现横条伪影。对刚接触这套方法的读者我的建议是先不要一上来就跑全炮。抽20到40道、时间窗截取前一半快速试参数看时频谱和分离剖面是否合理参数稳定后再放全炮。这一条能帮你省下大量无效计算时间。7. 一点个人体会折腾了这么一大圈我最深的体会是单一工具在混合数据分离问题上都有自己的天花板SST擅长时频曲波擅长空间结构两者结合不是炫技而是信息维度互补的必然选择。现在深度学习方法在很多数据分离任务里确实很火我也试过一些但标签短缺、混合模式变化时确定性方法反而更稳参数含义也更好解释。最后再分享一个小技巧调参时不要只盯着分离剖面看把“反射分量残差”叠回去和原始混合剖面比如果残差频谱和原始面波频谱高度吻合说明SST掩模还有泄漏阈值要加如果残差是白噪形态说明分离已经到底了。这样判断收敛比肉眼找同相轴快得多也更客观。