企业官网建设流程全解析

去年处理某工区一批多源同时激发采集的混合地震记录时，我被一个看似简单的问题卡了很久：面波、浅层折射和深层反射在时间-频率平面里基本重叠，常规带通滤波一压，面波是没了，但浅层反射的弱同相轴也跟着断了。后来我尝试把同步压缩小波变换（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 整体数据流的设计逻辑

我最后定的流程分五步，每一步都有明确的职责边界：

1. 读取混合剖面，组织成nt×nx的矩阵，逐道计算SST，得到三维时频体。
2. 在SST三维体上检测能量脊线，为面波和反射波分别构建时频掩模，逆变换得到两个初步分离分量。
3. 初步反射分量里通常还残留一些方向性强的面波尾巴，把它和初步面波分量一起送入曲波域。
4. 在曲波域里按尺度、角度做阈值收缩，把两类分量里的方向性伪影和噪声进一步剥离。
5. 输出最终面波分量、反射分量以及残差，残差再反馈回反射分量，保证弱信号能量不丢失。

这个顺序是反复试出来的。先做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谱在主导频段的相关系数，相关系数高说明弱信号细节保得好。

合成数据上的对比结果如下表：

可以看到，单用曲波效果甚至不如单用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以下。原因很简单，弱反射在曲波域里幅度本身就很低，全局阈值一刀切，它和噪声一起被归零。

后来我把策略改成“大尺度高阈值、小尺度低阈值”，并且分两轮迭代。第一轮用较高的c0=2.5，直接压面波；第二轮把c0降到0.8，只会处理第一轮剩下的残差。这样弱反射的损伤大幅下降，合成数据上SNR又多提升了2dB左右。

还有一个我吃过亏的细节：阈值操作不要用abs(coeff) > thr后取abs(coeff)重建，必须保留原始复数系数。我第一次图省事，把系数模值置零再逆变换，结果重建波形相位全乱，剖面上出现大量随机抖动。

6.3 生产级执行的效率与测试建议

整炮数据逐道算SST非常慢。我处理一炮360道、每道6000样点的数据，单核跑要将近一小时，生产上完全没法接受。后来改成parfor并行循环，12核的机器跑十到十五分钟能出结果，提升明显。如果你有GPU，wsst本身没有很好的GPU加速路径，但曲波变换的正反变换可以考虑用GPU实现，能省下第二级的大部分时间。

内存也要精打细算。三维SST体如果用double精度存储，大规模工区的数据量可以达到几十GB。建议只保存掩模和低频段投影，或者把整炮按滑窗切成若干重叠子块分别处理，处理完再拼回去。滑窗接缝处用余弦权重叠加，避免出现横条伪影。

对刚接触这套方法的读者，我的建议是先不要一上来就跑全炮。抽20到40道、时间窗截取前一半，快速试参数，看时频谱和分离剖面是否合理，参数稳定后再放全炮。这一条能帮你省下大量无效计算时间。

7. 一点个人体会

折腾了这么一大圈，我最深的体会是：单一工具在混合数据分离问题上都有自己的天花板，SST擅长时频，曲波擅长空间结构，两者结合不是炫技，而是信息维度互补的必然选择。现在深度学习方法在很多数据分离任务里确实很火，我也试过一些，但标签短缺、混合模式变化时，确定性方法反而更稳，参数含义也更好解释。

最后再分享一个小技巧：调参时不要只盯着分离剖面看，把“反射分量+残差”叠回去和原始混合剖面比，如果残差频谱和原始面波频谱高度吻合，说明SST掩模还有泄漏，阈值要加；如果残差是白噪形态，说明分离已经到底了。这样判断收敛比肉眼找同相轴快得多，也更客观。