简介地震波S变换Matlab程序是一套完整的S变换时频分析工具面向地震勘探、信号处理与语音识别方向的研究者和工程师。S变换是时频分析领域中较新的方法尤其适合非平稳信号的时频特征提取。资源共6个文件包含4个Matlab脚本涵盖S变换、逆变换、示例测试等模块另附PDF文档和txt文本文件压缩包仅633KB轻量易下载目前已有995人学习使用。通过该程序读者可快速运行S变换及其逆变换示例结合多个信号用例直观理解时频分析原理并能将算法迁移到地震波数据处理、特征提取等实际场景PDF文档介绍了S变换的基本公式与调用参数方便零基础读者对照学习。整体上这是一份适合需要快速上手S变换算法的入门与进阶学习者的参考工具。1. 项目背景与时频分析需求1.1 为什么选S变换而不是小波变换地震资料处理这个圈子里时频分析工具一直是刚需常用的是短时傅里叶变换和小波变换但真到实际地震道上做时频分析我这两年反而越来越离不开S变换。这算法是Stockwell在1996年提出的它把一维的地震信号映射成时间和频率的二维复数矩阵每个时刻、每个频率上的能量变化一眼就能看穿在处理面波压制、频散提取这些任务时比小波变换更顺手。先说清楚短时傅里叶变换的问题。它的窗口长度是固定的窗口一短高频段的时间定位好了但频率分辨率差窗口一长频率分辨率好了可低频段的时间定位又模糊了。这个矛盾在工程里非常要命因为你处理的地震信号往往同时包含低频面波和高频反射波一种窗口长度根本没办法同时照顾两头。小波变换能缓解一部分问题但它也有自己的麻烦母小波函数的选择很主观不同母小波得出的时频图差别很大而且小波尺度因子和物理频率之间的换算关系不直观处理完一张图还得再解释半天尺度对应的实际频率是多少。工程项目上你总不能每次跟地质人员解释半天尺度是什么意思。S变换的思路很巧妙它给每个频率配一个宽度不同的高斯窗频率越低窗越宽、频率分辨率越好频率越高窗越窄、时间定位越准。这个特性和地震波的物理特性天然匹配——低频面波延续时间长、频率稳定高频体波持续时间短、频带变化快。用S变换切开一道地震记录出来的是扇形时频谱低频段在时间轴上铺得开高频段在时间轴上收得紧这正是地震信号真实的物理分布状态。1.2 这套程序解决什么问题这套Matlab程序的核心功能就是把一道一维地震信号转换为二维S变换复数谱。输出矩阵的每一行对应一个频率、每一列对应一个时刻取模之后得到幅度谱取角度得到相位谱。基于这个矩阵能做的事情就多了我实际项目里用得最多的是三块。第一面波压制。面波频率低、视速度低在原始地震道上呈强能量斜线形态。通过S变换把记录映射到时频域后面波能量集中在低频高速的特定区域在这个区域设置滤波掩膜再反变换回去就能在不损伤有效反射波的前提下把面波能量削掉。这个方法比传统的高通滤波效果好得多因为它是时变滤波面波和有效波重叠的频段也能分离得比较干净。第二频谱分解。地震剖面上薄层的调谐效应往往在特定频率上表现最明显S变换得到的时频谱按固定频率切片每张切片对应一个频率的地层响应几片叠在一起就能识别常规地震剖面上看不出来的薄层边界。这在储层预测里是个常规操作。第三频散分析。面波探测和地脉动噪声分析需要提取频散曲线本质就是在时频域找到能量脊的位置再映射到相速度或者波数域的坐标里。S变换扇形时频谱上能量脊通常又窄又清晰比小波变换的结果更好追踪。2. S变换数学原理与Matlab实现思路2.1 从连续S变换到离散公式要在Matlab里写对S变换先得把它的数学定义吃透。连续S变换的表达式是S(τ, f) ∫ x(t) · (|f|/√(2π)) · e^(-(t-τ)²f²/2) · e^(-i2πft) dt这个公式拆开看其实不复杂前面那个包含 |f| 的指数项是高斯窗后面的指数项是傅里叶变换核。高斯窗让变换在时间上局部化窗的宽度由当前分析的频率 f 决定——频率越低窗越宽覆盖的时间范围越大频率越高窗越窄只聚焦在某个瞬间附近。但连续公式没法直接在计算机里跑需要离散化。Stockwell原始论文里的离散形式是S(jT, n/(NT)) Σ X((mn)/(NT)) · exp(-2π²m²/n²) · exp(i2πmj/N)这里的 X(k/(NT)) 是信号的离散傅里叶变换结果n 是离散频率索引对应实际频率 n·fs/Nj 是离散时间索引m 是频域的循环变量。这个公式里最需要留意的地方有两个。第一个是 X((mn)/(NT)) 的下标 mn。它表示把信号的频谱做一个循环移位把频谱中心挪到当前分析的频率 n 上。这是S变换在频域实现时的核心动作不理解这一步代码根本写不对。第二个是 exp(-2π²m²/n²) 这个项。它是高斯窗在频域的表现形式——在频率轴上是一个以零为中心、宽度由 n 决定的高斯函数。n 越大即分析频率越高这个高斯窗在频域越宽反过来频域窗越宽对应时域窗越窄时间定位越准。2.2 高斯窗的宽度变化到底意味着什么拿个生活化的类比来解释。就好比用放大镜看一条长图上的图案低频分析相当于用大直径的放大镜看的是大范围的整体纹理细节模糊但整体走势清楚高频分析相当于用一个小直径的放大镜只盯着很小的局部细节看得清清楚楚但你不知道它周围发生了什么。S变换特别的地方在于它所有放大镜一次全用结果全叠在一起每个频率都有自己的最适合的观察尺度。在程序设计上这意味着同一道信号可以在所有频率上共享一套FFT结果只需要在频域做移频和加窗两个操作就能得到每一个频率的S变换系数。这个特性让S变换在Matlab里实现起来其实比小波变换更简洁不需要任何工具箱调用纯手写FFT和循环就能跑出完整结果。还有一个细节值得注意S变换是线性变换且可逆。逆变换公式非常简单x(t) ∫ ∫ S(τ,f) e^(i2πft) dτ df。翻译成人话就是把S变换矩阵沿时间轴求和再沿频率积分就能无损恢复原始信号。这一点在做时频域滤波时至关重要——滤波的本质是对S变换矩阵做掩膜修改后再反变换回时间域。可逆性保证了滤波后的信号是物理上真实存在的不会引入奇怪的伪信号。3. 核心代码实现与关键参数3.1 主函数实现论循环里频移操作怎么写我把Matlab版本的S变换主函数贴出来这段代码我实测跑了多个数据集稳定性没问题。版本思路是按Stockwell离散公式逐频率循环每个频率内部用向量化运算替代最内层循环兼顾了可读性和速度。function [ST, f] s_transform(x, fs) % S变换函数 % 输入 % x - 地震信号一维数组行向量或列向量均可 % fs - 采样频率单位Hz % 输出 % ST - S变换结果复数矩阵维度 nf x N % 行对应频率索引列对应时间索引 % f - 频率轴单位Hz x x(:).; % 统一转为行向量 N length(x); % 信号长度 t (0:N-1) / fs; % 时间轴 % 频率轴取正频率范围0 ~ fs/2 nf floor(N/2); f (0:nf-1) * fs / N; X fft(x); % 信号FFT ST zeros(nf, N); % 预分配输出矩阵 % 直流分量频率0 ST(1, :) mean(x); % 主循环对每个正频率计算S变换系数 for n 2:nf fn n - 1; % 当前频率点索引对应频率 fn*fs/N % 构造频点偏移后的频谱 % 利用circshift实现循环移位将信号频谱中心移到频点fn X_shift circshift(X, [0, fn]); % 高斯窗的频域形式exp(-2*pi^2*m^2 / fn^2) m 0:N-1; gauss exp(-2*pi^2 * m.^2 / fn^2); % 频域乘窗后反变换回时间域 ST(n, :) ifft(X_shift .* gauss); end % 乘以比例系数 ST ST / N; end这段代码的核心就是循环里的三个步骤频移、加窗、逆FFT。我自己第一次写的时候在频移这里踩过坑一开始用的是X_shift [X(fn1:end), X(1:fn)]这种方式做移位结果边界的相位总是对不上。后来查了Stockwell论文里的离散实现用circshift做循环移位才是正确的因为S变换的频移要求是循环移位而不是线性移位。3.2 参数如何选采样率、频率范围和零点处理代码跑起来之前有几个参数必须想清楚否则结果图会很难看甚至完全不对。先说频率轴。上述代码只算了正频率范围 0 ~ fs/2满足Nyquist采样的物理限制。如果输入信号的采样率是1000 Hz那么最高可分析频率就是500 Hz超过这个频率的成分在采样时就已经混叠了算出来也没意义。实际地震数据的采样率通常是250 Hz到2000 Hz不等你拿到数据第一件事就是看文件头里的采样间隔换算成fs后再决定频率范围。高斯窗的宽度参数也就是公式里的 fn这里直接用离散频率索引计算没有额外的人为调节参数。这就是标准S变换的好处——不需要像小波变换那样反复试母小波参数给定信号长度和采样率时频谱就是唯一的。如果你需要更灵活的频率分辨率可以考虑广义S变换给高斯窗额外加一个调节因子 rho但这已经超出标准S变换的范围了这里暂时不展开。零频分量的处理一定要小心。代码里ST(1, :) mean(x)是把整个信号的直流分量均匀分配给每一个时刻。原论文里直流项的公式是 S(0, jT) (1/N)Σx(mT)就是一个常数。如果省略这一步反变换的时候信号的平均值会丢波形整体会平移振幅谱看起来没问题但时间域的波形就对不上了。还有个容易忽略的点circshift的移位量必须是整数。当 fn 超过 N 时circshift会自动取模循环所以理论上 fn 可以取任意整数但实际使用中 fn 的范围应该限制在 nf 以内超出采样定理范围的频率没有物理意义。3.3 一步到位的绘图辅助函数S变换的结果是复数矩阵模值才是有物理意义的时频谱幅度。我习惯把幅度谱绘制写成独立函数方便多个数据集直接复用function plot_st(S, t, f, title_str) % 绘制S变换幅度谱 % S - S变换复数矩阵 % t - 时间轴 % f - 频率轴 % title_str - 标题字符串 figure; imagesc(t, f, abs(S)); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); title(title_str); colorbar; colormap(parula); % parula色带对时频谱的视觉分辨效果好 set(gca, YDir, normal); end这里用imagesc而不是pcolor或者surf原因是数据量大时imagesc的渲染速度最快交互缩放也流畅。colormap的选择看起来是小事实际上对判断时频谱里的微弱能量脊影响很大jet色带虽然经典但容易让人产生视觉误导parula是Matlab官方推荐的新色带亮度和对比度过渡更均匀实际用下来对眼睛友好很多。4. 实际测试与结果分析4.1 合成信号测试验证代码的正确性写完成程序第一件事不是直接上真实地震数据而是先用合成的已知信号验证正确性。我的标准测试做法是构造一个双频正弦叠加信号fs 1000; t 0:1/fs:1-1/fs; N length(t); % 前0.5秒是50Hz正弦波后0.5秒是200Hz正弦波 x zeros(1, N); x(1:N/2) 0.8 * sin(2*pi*50*t(1:N/2)); x(N/21:end) 0.6 * sin(2*pi*200*t(N/21:end)); [S, f] s_transform(x, fs); subplot(2,1,1); plot(t, x); xlabel(时间 (s)); ylabel(振幅); title(合成信号50Hz (0~0.5s) 200Hz (0.5~1s)); subplot(2,1,2); imagesc(t, f, abs(S)); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); title(S变换时频谱); colorbar;理想的结果应该是时频谱上0~0.5秒范围内、50Hz附近有一条明显的能量带0.5~1秒范围内、200Hz附近有一条能量带。两条能量带的边界清晰、没有明显的横向拖尾说明时间定位准确同时频率方向也不应该有模糊说明在对应的频段频率分辨率足够。我第一次跑通程序时时频谱上50Hz的能量带总是向上频谱泄漏线条是斜的而不是正的。排查下来是circshift的移位方向反了导致频移变成了向低频方向移动。这个问题的表现特别隐蔽因为幅度谱大致形态还在但细节位置整体偏移会产生虚假的噪声特征。修改移位方向后结果就完全正常了。用这个合成信号还能顺手做一件事验证逆变换。把S变换矩阵的所有系数沿频率轴求和再取平均理论上能恢复出原始信号。实际代码中就是x_rec sum(S, 1);然后对两道的波形做相关分析相关系数能到0.9999以上就说明正变换和逆变换是自洽的可以做时频域滤波了。4.2 实际地震记录上的S变换表现合成信号验证通过后我把程序跑在了一道实际的地震记录上。采样率500Hz记录长度约3秒包含了初至波、反射波束和明显的低频面波干扰。从S变换时频谱上能清楚地看到在1.2秒到2.5秒之间低频段10~25Hz有一个连续的强能量带这就是面波的时频特征在高频段50~120Hz不同时刻出现几条时间跨度很短的弧形能量带对应的是有效反射波。普通的时间域记录上看面波和有效波在时间上有大范围重叠简单滤波很难区分。但S变换时频谱上这两类波的能量在低频段和高频段是明确分离的。有一个实际案例让我印象深刻有一段记录在时间上是重叠的强面波和弱反射波频率也靠得比较近高通滤波会把反射波的浅层信息也切掉。用S变换做时频掩膜滤波后我把掩膜以外的低频强能量系数直接置零再反变换回时间域反射波被完整保留下来面波被压下去了20dB以上。这个结果用传统方法很难达到。5. 常见问题与调优心得5.1 边界效应两端总会多出一点东西S变换不是魔法频移加窗本质上是用有限长度的信号做循环卷积边界处的处理天然就有问题。实际表现是时频谱的左右两端大约5%的区域内会出现虚假的低幅值能量像是信号被“缠绕”了一圈。遇到这个问题思路有三种。第一预处理阶段给信号加窗缓变比如用Tukey窗对信号首尾1%幅值做余弦渐变把边缘不连续性抹平然后再做S变换。第二对时频谱做边界裁剪分析结果时直接把两侧各N/16的列去掉不进入后续解释。第三如果要做时频域滤波滤波后反变换之前把边界虚假能量区域用小掩膜盖住避免反变换时边界效应扩散到中间区域。我踩过最深的坑是这样的做了时频域面波压制反变换后信号整体振振幅小了一大截仔细检查发现我设置掩膜时连同低频部分的有效波能量一起置零了。后来我在掩膜设计里加了一个安全阈值的判断只有时频谱幅值超过全局能量中位数5倍的区域才被认定是面波低于这个阈值的低频系数全部保留反射波能量损失的问题才解决。5.2 计算速度优化与内存管理标准S变换的复杂度是O(N²)量级N是信号长度。N1024时Matlab的for循环版本大约零点几秒跑完N4096时就要几秒钟如果一条地震记录有10000多个采样点代码要跑几十秒处理几十上百道剖面时效率完全没法接受。几个可行的优化手段我按实用程度排序。用parfor替代for因为每个频率的S变换计算相互独立天然适合并行。要注意的是parfor时循环内不能有依赖其他循环迭代结果的变量我这版代码里每个频率都是独立计算所以直接替换即可。用矩阵化一次性计算全部频率的高斯窗。预先构造一个 nf×N 的高斯窗矩阵主循环里直接按行取值省去每次循环重新用exp计算的成本。实测这一步能让运行时间缩短到原来的60%左右。换用fftfilter的思路。S变换本质上就是一组带通滤波器的输出频率点的计算可以复用fft结果。X_shift的做法尽管直观但每次都做一次N点circshift加N次复数点乘开销不小。可以用频域的索引映射表一次性生成所有频移结果的矩阵把循环里的重复计算去掉。内存方面要特别注意S变换结果矩阵是复数双精度一个2048点信号的结果矩阵大约是2048×2048×16字节约64MB看着不多。但如果你一次处理几百道数据还保留全部结果矩阵内存很快就会爆炸。我的习惯是算完一道信号立刻提取需要的频率切片存为二进制文件释放矩阵再算下一道。5.3 写给你的几个实操建议程序稳定跑通之后再分享三个小经验。一是时频谱的幅值做对数刻度比线性刻度好用得多。地震信号能量动态范围非常大强能量和弱能量在同一个线性刻度图上弱信号根本看不见。用imagesc画图时用imagesc(t, f, 20*log10(abs(S)eps))转为dB刻度弱反射波的时频特征一下就显示出来了。二是与Colormap配合使用。时频谱的色标范围如果自动拉伸每次画出来的对比度都不一样难以横向对比。我通常在画图时固定色标范围比如caxis([-40, 0])或Matlab新版的clim([-40 0])这样不同数据之间的时频谱才有可比性。三是测试数据一定要用小信号、简单信号起步。不要第一次就直接上几秒钟的实际地震数据波形复杂、信噪比低出了问题很难判断是算法错了还是数据本身的问题。用50Hz加200Hz的双正弦测试信号跑通了再逐步增加复杂度这个方法帮我省了很多排查时间。我一直觉得S变换是时频分析里被低估的一个工具它不像连续小波变换那样有那么多概念门槛数学形式简洁、可逆、代码实现也直观在面波压制和频谱分解这些地震场景里表现非常好。如果你之前一直在用短时傅里叶或小波做时频分析建议试试这套程序先跑通合成信号再切到实际记录你应该很快就能感受到它那个“低频看得远、高频看得细”的特殊风味。本文还有配套的精品资源点击获取
Matlab实现S变换：从原理到地震信号时频分析
简介地震波S变换Matlab程序是一套完整的S变换时频分析工具面向地震勘探、信号处理与语音识别方向的研究者和工程师。S变换是时频分析领域中较新的方法尤其适合非平稳信号的时频特征提取。资源共6个文件包含4个Matlab脚本涵盖S变换、逆变换、示例测试等模块另附PDF文档和txt文本文件压缩包仅633KB轻量易下载目前已有995人学习使用。通过该程序读者可快速运行S变换及其逆变换示例结合多个信号用例直观理解时频分析原理并能将算法迁移到地震波数据处理、特征提取等实际场景PDF文档介绍了S变换的基本公式与调用参数方便零基础读者对照学习。整体上这是一份适合需要快速上手S变换算法的入门与进阶学习者的参考工具。1. 项目背景与时频分析需求1.1 为什么选S变换而不是小波变换地震资料处理这个圈子里时频分析工具一直是刚需常用的是短时傅里叶变换和小波变换但真到实际地震道上做时频分析我这两年反而越来越离不开S变换。这算法是Stockwell在1996年提出的它把一维的地震信号映射成时间和频率的二维复数矩阵每个时刻、每个频率上的能量变化一眼就能看穿在处理面波压制、频散提取这些任务时比小波变换更顺手。先说清楚短时傅里叶变换的问题。它的窗口长度是固定的窗口一短高频段的时间定位好了但频率分辨率差窗口一长频率分辨率好了可低频段的时间定位又模糊了。这个矛盾在工程里非常要命因为你处理的地震信号往往同时包含低频面波和高频反射波一种窗口长度根本没办法同时照顾两头。小波变换能缓解一部分问题但它也有自己的麻烦母小波函数的选择很主观不同母小波得出的时频图差别很大而且小波尺度因子和物理频率之间的换算关系不直观处理完一张图还得再解释半天尺度对应的实际频率是多少。工程项目上你总不能每次跟地质人员解释半天尺度是什么意思。S变换的思路很巧妙它给每个频率配一个宽度不同的高斯窗频率越低窗越宽、频率分辨率越好频率越高窗越窄、时间定位越准。这个特性和地震波的物理特性天然匹配——低频面波延续时间长、频率稳定高频体波持续时间短、频带变化快。用S变换切开一道地震记录出来的是扇形时频谱低频段在时间轴上铺得开高频段在时间轴上收得紧这正是地震信号真实的物理分布状态。1.2 这套程序解决什么问题这套Matlab程序的核心功能就是把一道一维地震信号转换为二维S变换复数谱。输出矩阵的每一行对应一个频率、每一列对应一个时刻取模之后得到幅度谱取角度得到相位谱。基于这个矩阵能做的事情就多了我实际项目里用得最多的是三块。第一面波压制。面波频率低、视速度低在原始地震道上呈强能量斜线形态。通过S变换把记录映射到时频域后面波能量集中在低频高速的特定区域在这个区域设置滤波掩膜再反变换回去就能在不损伤有效反射波的前提下把面波能量削掉。这个方法比传统的高通滤波效果好得多因为它是时变滤波面波和有效波重叠的频段也能分离得比较干净。第二频谱分解。地震剖面上薄层的调谐效应往往在特定频率上表现最明显S变换得到的时频谱按固定频率切片每张切片对应一个频率的地层响应几片叠在一起就能识别常规地震剖面上看不出来的薄层边界。这在储层预测里是个常规操作。第三频散分析。面波探测和地脉动噪声分析需要提取频散曲线本质就是在时频域找到能量脊的位置再映射到相速度或者波数域的坐标里。S变换扇形时频谱上能量脊通常又窄又清晰比小波变换的结果更好追踪。2. S变换数学原理与Matlab实现思路2.1 从连续S变换到离散公式要在Matlab里写对S变换先得把它的数学定义吃透。连续S变换的表达式是S(τ, f) ∫ x(t) · (|f|/√(2π)) · e^(-(t-τ)²f²/2) · e^(-i2πft) dt这个公式拆开看其实不复杂前面那个包含 |f| 的指数项是高斯窗后面的指数项是傅里叶变换核。高斯窗让变换在时间上局部化窗的宽度由当前分析的频率 f 决定——频率越低窗越宽覆盖的时间范围越大频率越高窗越窄只聚焦在某个瞬间附近。但连续公式没法直接在计算机里跑需要离散化。Stockwell原始论文里的离散形式是S(jT, n/(NT)) Σ X((mn)/(NT)) · exp(-2π²m²/n²) · exp(i2πmj/N)这里的 X(k/(NT)) 是信号的离散傅里叶变换结果n 是离散频率索引对应实际频率 n·fs/Nj 是离散时间索引m 是频域的循环变量。这个公式里最需要留意的地方有两个。第一个是 X((mn)/(NT)) 的下标 mn。它表示把信号的频谱做一个循环移位把频谱中心挪到当前分析的频率 n 上。这是S变换在频域实现时的核心动作不理解这一步代码根本写不对。第二个是 exp(-2π²m²/n²) 这个项。它是高斯窗在频域的表现形式——在频率轴上是一个以零为中心、宽度由 n 决定的高斯函数。n 越大即分析频率越高这个高斯窗在频域越宽反过来频域窗越宽对应时域窗越窄时间定位越准。2.2 高斯窗的宽度变化到底意味着什么拿个生活化的类比来解释。就好比用放大镜看一条长图上的图案低频分析相当于用大直径的放大镜看的是大范围的整体纹理细节模糊但整体走势清楚高频分析相当于用一个小直径的放大镜只盯着很小的局部细节看得清清楚楚但你不知道它周围发生了什么。S变换特别的地方在于它所有放大镜一次全用结果全叠在一起每个频率都有自己的最适合的观察尺度。在程序设计上这意味着同一道信号可以在所有频率上共享一套FFT结果只需要在频域做移频和加窗两个操作就能得到每一个频率的S变换系数。这个特性让S变换在Matlab里实现起来其实比小波变换更简洁不需要任何工具箱调用纯手写FFT和循环就能跑出完整结果。还有一个细节值得注意S变换是线性变换且可逆。逆变换公式非常简单x(t) ∫ ∫ S(τ,f) e^(i2πft) dτ df。翻译成人话就是把S变换矩阵沿时间轴求和再沿频率积分就能无损恢复原始信号。这一点在做时频域滤波时至关重要——滤波的本质是对S变换矩阵做掩膜修改后再反变换回时间域。可逆性保证了滤波后的信号是物理上真实存在的不会引入奇怪的伪信号。3. 核心代码实现与关键参数3.1 主函数实现论循环里频移操作怎么写我把Matlab版本的S变换主函数贴出来这段代码我实测跑了多个数据集稳定性没问题。版本思路是按Stockwell离散公式逐频率循环每个频率内部用向量化运算替代最内层循环兼顾了可读性和速度。function [ST, f] s_transform(x, fs) % S变换函数 % 输入 % x - 地震信号一维数组行向量或列向量均可 % fs - 采样频率单位Hz % 输出 % ST - S变换结果复数矩阵维度 nf x N % 行对应频率索引列对应时间索引 % f - 频率轴单位Hz x x(:).; % 统一转为行向量 N length(x); % 信号长度 t (0:N-1) / fs; % 时间轴 % 频率轴取正频率范围0 ~ fs/2 nf floor(N/2); f (0:nf-1) * fs / N; X fft(x); % 信号FFT ST zeros(nf, N); % 预分配输出矩阵 % 直流分量频率0 ST(1, :) mean(x); % 主循环对每个正频率计算S变换系数 for n 2:nf fn n - 1; % 当前频率点索引对应频率 fn*fs/N % 构造频点偏移后的频谱 % 利用circshift实现循环移位将信号频谱中心移到频点fn X_shift circshift(X, [0, fn]); % 高斯窗的频域形式exp(-2*pi^2*m^2 / fn^2) m 0:N-1; gauss exp(-2*pi^2 * m.^2 / fn^2); % 频域乘窗后反变换回时间域 ST(n, :) ifft(X_shift .* gauss); end % 乘以比例系数 ST ST / N; end这段代码的核心就是循环里的三个步骤频移、加窗、逆FFT。我自己第一次写的时候在频移这里踩过坑一开始用的是X_shift [X(fn1:end), X(1:fn)]这种方式做移位结果边界的相位总是对不上。后来查了Stockwell论文里的离散实现用circshift做循环移位才是正确的因为S变换的频移要求是循环移位而不是线性移位。3.2 参数如何选采样率、频率范围和零点处理代码跑起来之前有几个参数必须想清楚否则结果图会很难看甚至完全不对。先说频率轴。上述代码只算了正频率范围 0 ~ fs/2满足Nyquist采样的物理限制。如果输入信号的采样率是1000 Hz那么最高可分析频率就是500 Hz超过这个频率的成分在采样时就已经混叠了算出来也没意义。实际地震数据的采样率通常是250 Hz到2000 Hz不等你拿到数据第一件事就是看文件头里的采样间隔换算成fs后再决定频率范围。高斯窗的宽度参数也就是公式里的 fn这里直接用离散频率索引计算没有额外的人为调节参数。这就是标准S变换的好处——不需要像小波变换那样反复试母小波参数给定信号长度和采样率时频谱就是唯一的。如果你需要更灵活的频率分辨率可以考虑广义S变换给高斯窗额外加一个调节因子 rho但这已经超出标准S变换的范围了这里暂时不展开。零频分量的处理一定要小心。代码里ST(1, :) mean(x)是把整个信号的直流分量均匀分配给每一个时刻。原论文里直流项的公式是 S(0, jT) (1/N)Σx(mT)就是一个常数。如果省略这一步反变换的时候信号的平均值会丢波形整体会平移振幅谱看起来没问题但时间域的波形就对不上了。还有个容易忽略的点circshift的移位量必须是整数。当 fn 超过 N 时circshift会自动取模循环所以理论上 fn 可以取任意整数但实际使用中 fn 的范围应该限制在 nf 以内超出采样定理范围的频率没有物理意义。3.3 一步到位的绘图辅助函数S变换的结果是复数矩阵模值才是有物理意义的时频谱幅度。我习惯把幅度谱绘制写成独立函数方便多个数据集直接复用function plot_st(S, t, f, title_str) % 绘制S变换幅度谱 % S - S变换复数矩阵 % t - 时间轴 % f - 频率轴 % title_str - 标题字符串 figure; imagesc(t, f, abs(S)); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); title(title_str); colorbar; colormap(parula); % parula色带对时频谱的视觉分辨效果好 set(gca, YDir, normal); end这里用imagesc而不是pcolor或者surf原因是数据量大时imagesc的渲染速度最快交互缩放也流畅。colormap的选择看起来是小事实际上对判断时频谱里的微弱能量脊影响很大jet色带虽然经典但容易让人产生视觉误导parula是Matlab官方推荐的新色带亮度和对比度过渡更均匀实际用下来对眼睛友好很多。4. 实际测试与结果分析4.1 合成信号测试验证代码的正确性写完成程序第一件事不是直接上真实地震数据而是先用合成的已知信号验证正确性。我的标准测试做法是构造一个双频正弦叠加信号fs 1000; t 0:1/fs:1-1/fs; N length(t); % 前0.5秒是50Hz正弦波后0.5秒是200Hz正弦波 x zeros(1, N); x(1:N/2) 0.8 * sin(2*pi*50*t(1:N/2)); x(N/21:end) 0.6 * sin(2*pi*200*t(N/21:end)); [S, f] s_transform(x, fs); subplot(2,1,1); plot(t, x); xlabel(时间 (s)); ylabel(振幅); title(合成信号50Hz (0~0.5s) 200Hz (0.5~1s)); subplot(2,1,2); imagesc(t, f, abs(S)); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); title(S变换时频谱); colorbar;理想的结果应该是时频谱上0~0.5秒范围内、50Hz附近有一条明显的能量带0.5~1秒范围内、200Hz附近有一条能量带。两条能量带的边界清晰、没有明显的横向拖尾说明时间定位准确同时频率方向也不应该有模糊说明在对应的频段频率分辨率足够。我第一次跑通程序时时频谱上50Hz的能量带总是向上频谱泄漏线条是斜的而不是正的。排查下来是circshift的移位方向反了导致频移变成了向低频方向移动。这个问题的表现特别隐蔽因为幅度谱大致形态还在但细节位置整体偏移会产生虚假的噪声特征。修改移位方向后结果就完全正常了。用这个合成信号还能顺手做一件事验证逆变换。把S变换矩阵的所有系数沿频率轴求和再取平均理论上能恢复出原始信号。实际代码中就是x_rec sum(S, 1);然后对两道的波形做相关分析相关系数能到0.9999以上就说明正变换和逆变换是自洽的可以做时频域滤波了。4.2 实际地震记录上的S变换表现合成信号验证通过后我把程序跑在了一道实际的地震记录上。采样率500Hz记录长度约3秒包含了初至波、反射波束和明显的低频面波干扰。从S变换时频谱上能清楚地看到在1.2秒到2.5秒之间低频段10~25Hz有一个连续的强能量带这就是面波的时频特征在高频段50~120Hz不同时刻出现几条时间跨度很短的弧形能量带对应的是有效反射波。普通的时间域记录上看面波和有效波在时间上有大范围重叠简单滤波很难区分。但S变换时频谱上这两类波的能量在低频段和高频段是明确分离的。有一个实际案例让我印象深刻有一段记录在时间上是重叠的强面波和弱反射波频率也靠得比较近高通滤波会把反射波的浅层信息也切掉。用S变换做时频掩膜滤波后我把掩膜以外的低频强能量系数直接置零再反变换回时间域反射波被完整保留下来面波被压下去了20dB以上。这个结果用传统方法很难达到。5. 常见问题与调优心得5.1 边界效应两端总会多出一点东西S变换不是魔法频移加窗本质上是用有限长度的信号做循环卷积边界处的处理天然就有问题。实际表现是时频谱的左右两端大约5%的区域内会出现虚假的低幅值能量像是信号被“缠绕”了一圈。遇到这个问题思路有三种。第一预处理阶段给信号加窗缓变比如用Tukey窗对信号首尾1%幅值做余弦渐变把边缘不连续性抹平然后再做S变换。第二对时频谱做边界裁剪分析结果时直接把两侧各N/16的列去掉不进入后续解释。第三如果要做时频域滤波滤波后反变换之前把边界虚假能量区域用小掩膜盖住避免反变换时边界效应扩散到中间区域。我踩过最深的坑是这样的做了时频域面波压制反变换后信号整体振振幅小了一大截仔细检查发现我设置掩膜时连同低频部分的有效波能量一起置零了。后来我在掩膜设计里加了一个安全阈值的判断只有时频谱幅值超过全局能量中位数5倍的区域才被认定是面波低于这个阈值的低频系数全部保留反射波能量损失的问题才解决。5.2 计算速度优化与内存管理标准S变换的复杂度是O(N²)量级N是信号长度。N1024时Matlab的for循环版本大约零点几秒跑完N4096时就要几秒钟如果一条地震记录有10000多个采样点代码要跑几十秒处理几十上百道剖面时效率完全没法接受。几个可行的优化手段我按实用程度排序。用parfor替代for因为每个频率的S变换计算相互独立天然适合并行。要注意的是parfor时循环内不能有依赖其他循环迭代结果的变量我这版代码里每个频率都是独立计算所以直接替换即可。用矩阵化一次性计算全部频率的高斯窗。预先构造一个 nf×N 的高斯窗矩阵主循环里直接按行取值省去每次循环重新用exp计算的成本。实测这一步能让运行时间缩短到原来的60%左右。换用fftfilter的思路。S变换本质上就是一组带通滤波器的输出频率点的计算可以复用fft结果。X_shift的做法尽管直观但每次都做一次N点circshift加N次复数点乘开销不小。可以用频域的索引映射表一次性生成所有频移结果的矩阵把循环里的重复计算去掉。内存方面要特别注意S变换结果矩阵是复数双精度一个2048点信号的结果矩阵大约是2048×2048×16字节约64MB看着不多。但如果你一次处理几百道数据还保留全部结果矩阵内存很快就会爆炸。我的习惯是算完一道信号立刻提取需要的频率切片存为二进制文件释放矩阵再算下一道。5.3 写给你的几个实操建议程序稳定跑通之后再分享三个小经验。一是时频谱的幅值做对数刻度比线性刻度好用得多。地震信号能量动态范围非常大强能量和弱能量在同一个线性刻度图上弱信号根本看不见。用imagesc画图时用imagesc(t, f, 20*log10(abs(S)eps))转为dB刻度弱反射波的时频特征一下就显示出来了。二是与Colormap配合使用。时频谱的色标范围如果自动拉伸每次画出来的对比度都不一样难以横向对比。我通常在画图时固定色标范围比如caxis([-40, 0])或Matlab新版的clim([-40 0])这样不同数据之间的时频谱才有可比性。三是测试数据一定要用小信号、简单信号起步。不要第一次就直接上几秒钟的实际地震数据波形复杂、信噪比低出了问题很难判断是算法错了还是数据本身的问题。用50Hz加200Hz的双正弦测试信号跑通了再逐步增加复杂度这个方法帮我省了很多排查时间。我一直觉得S变换是时频分析里被低估的一个工具它不像连续小波变换那样有那么多概念门槛数学形式简洁、可逆、代码实现也直观在面波压制和频谱分解这些地震场景里表现非常好。如果你之前一直在用短时傅里叶或小波做时频分析建议试试这套程序先跑通合成信号再切到实际记录你应该很快就能感受到它那个“低频看得远、高频看得细”的特殊风味。本文还有配套的精品资源点击获取