简介面向地震勘探、信号处理与时频分析研究者的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变换存在精确的逆变换我做完时频域滤波之后能无损地回到时间域。下面是三种方法的直观对比。方法窗函数时间-频率分辨率频率轴物理含义逆变换短时傅里叶固定宽度窗全局固定需人工折中明确Hz有但窗的选取影响重构质量连续小波随尺度伸缩低频宽、高频窄需要额外换算有但基函数选择影响结果S变换高斯窗随频率反向变化低频宽、高频窄明确Hz有且与傅里叶谱直接对应我个人的体会是如果只是要看定性时频特征三个方法都能用但如果要做定量分析比如从时频谱里提取频散曲线、做时频滤波之后再进行后续反演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) Σ_{k0}^{N-1} X(kn) · e^(-2π²k²/n²) · e^(i2πkm/N)这里X(k)是信号x(t)的离散傅里叶变换DFTn是频率索引m是时间索引。关键点在于对每一个固定的频率n先对X做循环移位把第kn个频率分量挪到第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变换的高斯窗严格跟随频率倒数变化这是其理论上的优势但在实际地震数据处理里有时你会觉得低频端的时间分辨率不够或者高频端的频率分辨率不够。这时可以对高斯窗引入一个调节参数gammagauss 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完整S矩阵内存只用前N/2频率的内存102416 MB8 MB4096256 MB128 MB81921 GB512 MB163844 GB2 GB实际地震记录动辄几万采样点如果直接做完整的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实现直接、结果可量化验证非常适合作为自己深入研究时频分析的第一个完整工具。本文还有配套的精品资源点击获取 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变换存在精确的逆变换我做完时频域滤波之后能无损地回到时间域。下面是三种方法的直观对比。方法窗函数时间-频率分辨率频率轴物理含义逆变换短时傅里叶固定宽度窗全局固定需人工折中明确Hz有但窗的选取影响重构质量连续小波随尺度伸缩低频宽、高频窄需要额外换算有但基函数选择影响结果S变换高斯窗随频率反向变化低频宽、高频窄明确Hz有且与傅里叶谱直接对应我个人的体会是如果只是要看定性时频特征三个方法都能用但如果要做定量分析比如从时频谱里提取频散曲线、做时频滤波之后再进行后续反演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) Σ_{k0}^{N-1} X(kn) · e^(-2π²k²/n²) · e^(i2πkm/N)这里X(k)是信号x(t)的离散傅里叶变换DFTn是频率索引m是时间索引。关键点在于对每一个固定的频率n先对X做循环移位把第kn个频率分量挪到第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变换的高斯窗严格跟随频率倒数变化这是其理论上的优势但在实际地震数据处理里有时你会觉得低频端的时间分辨率不够或者高频端的频率分辨率不够。这时可以对高斯窗引入一个调节参数gammagauss 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完整S矩阵内存只用前N/2频率的内存102416 MB8 MB4096256 MB128 MB81921 GB512 MB163844 GB2 GB实际地震记录动辄几万采样点如果直接做完整的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实现直接、结果可量化验证非常适合作为自己深入研究时频分析的第一个完整工具。本文还有配套的精品资源点击获取