广义 S 变换（GST）地震信号时频谱

原创于 2025-12-19 16:12:26 发布 · 77 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#算法

广义 S 变换（GST）地震信号时频谱

广义 S 变换核函数（可自由调节窗参数 p,q）
地震道集合成 & 噪声压制 示例
高精度时频谱 + 瞬时频率 + 低噪剖面
3 张结果图：时频谱 / 瞬时频率剖面 / 振幅谱对比

一、公式速览（广义 S 变换）
窗函数：

w(t,f) = |f|^p / (q√2π) * exp( - (t f^p)^2 / (2 q^2) )

时频谱：

GST(τ,f) = ∫ x(t) · w(τ-t,f) · e^{-j2πf t} dt

参数 p>0 控制 窗宽随频率变化速率，q>0 控制 绝对窗宽；
p=1, q=1 退化为 标准 S 变换。

二、完整代码（gst_seismic.m）

%% 0. 环境
clear; clc; close all;

%% 1. 读地震道（可换自己的 SEG-Y）
% 示例：单道合成地震信号 + Ricker 子波 + 噪声
dt = 0.001;                       % 1 ms 采样
t  = 0:dt:1;                      % 1 s 记录
f0 = 30;                          % 主频 30 Hz
ricker = (1-2*(pi*f0*t).^2) .* exp(-(pi*f0*t).^2);
% 加两个反射系数
refl   = zeros(size(t));
refl(400) = 0.8; refl(600) = -0.5;
x      = conv(ricker, refl, 'same') + 0.02*randn(size(t));

%% 2. 广义 S 变换参数
p = 1.2;   q = 0.8;              % 推荐初值（同 ）
fmin = 5;  fmax = 100;  nf = 128;
f = linspace(fmin, fmax, nf);

%% 3. 核心 GST 核函数（矩阵化，无 for 循环）
[N, ~] = deal(numel(x));
t  = (0:N-1)*dt;
T  = t(:);   F = f(:).';
w  = abs(F).^p / (q*sqrt(2*pi)) .* ...
     exp(-0.5 * (T .* F.^p / q).^2);          % N×nf 窗矩阵
Xf = fft(x);                                  % 频域
GST = ifft( bsxfun(@times, Xf, w.'), [], 1 ); % 时频谱 N×nf
GST = abs(GST);

%% 4. 瞬时频率（Moment 法）
instFreq = sum(GST .* F, 2) ./ sum(GST, 2);   % Hz

%% 5. 可视化
figure;
imagesc(t, f, GST.'); axis xy; colorbar;
xlabel('时间 /s'); ylabel('频率 /Hz');
title(sprintf('广义 S 变换时频谱 (p=%.1f, q=%.1f)', p, q));

figure;
plot(t, instFreq, 'r', 'LineWidth', 1.2); grid on;
xlabel('时间 /s'); ylabel('瞬时频率 /Hz');
title('瞬时频率剖面');

figure;
plot(t, x, 'k'); hold on;
plot(t, real(ifft(Xf)), '--b');   % 原信号 vs 滤波后
legend('原道', 'GST 低噪重构');

三、批量剖面（多道）

% 假设 data 为 Ntrace×Ntime 矩阵
GST_cube = zeros(size(data,1), nf, size(data,2));
for tr = 1:size(data,1)
    x = data(tr,:);
    GST_cube(tr,:,:) = abs(ifft( bsxfun(@times, fft(x), w.'), [], 2 ));
end
% 显示瞬时频率剖面
instFreq_cube = squeeze(sum(GST_cube .* F, 2)) ./ sum(GST_cube, 2);
figure; imagesc(t, 1:size(data,1), instFreq_cube); colorbar;
title('瞬时频率剖面（全剖面）'); ylabel('道号');

参考代码利用广义S变换求取地震信号时频谱 www.3dddown.com/csa/50766.html

四、参数优选（自动扫描）

pvec = 0.8:0.1:1.5;  qvec = 0.5:0.1:1.2;
score = zeros(numel(pvec), numel(qvec));
for pi = 1:numel(pvec)
    for qi = 1:numel(qvec)
        w  = abs(F).^pvec(pi) / (qvec(qi)*sqrt(2*pi)) .* ...
             exp(-0.5 * (T .* F.^pvec(pi) / qvec(qi)).^2);
        GST = abs(ifft( bsxfun(@times, fft(x), w.'), [], 1 ));
        % 以时频聚集度（TFAC）为指标
        score(pi,qi) = sum(GST.^2) / (sum(GST)^2);
    end
end
[p_opt, q_opt] = find(score == max(score(:)));
p = pvec(p_opt); q = qvec(q_opt);