简介广义S变换GST是一种面向非平稳信号的时频分析方法相比短时傅立叶变换可在时间与频率分辨率之间更灵活地权衡因此常用于地震事件识别、地质勘探和地球内部结构分析等场景。这份压缩包提供的是GST核心算法的C语言实现面向具备C编程基础和信号处理理论的研究人员、工程师及高年级学生既可用于实际数据分析也可作为算法学习的参考代码。资源包仅含1个C源文件大小约2KB代码非常精简其中通常包括信号预处理、广义S变换积分路径与权重计算、尺度参数调节、结果后处理等关键模块方便用户根据自身数据和需求修改参数、移植到其他工程或扩展功能。目前已有224人学习下载。借助该程序用户可以快速获得信号的时频表示深入理解非平稳信号的局部特征尤其对地震数据中的异常识别与地下结构分析具有直接的参考价值。1. 为什么地震数据要放弃短时傅里叶改用广义S变换地震道和VSP记录这类信号最麻烦的不是幅值小而是频率成分随时间剧烈变化。短时傅里叶变换STFT一旦把窗长定死高频段的时间分辨率和低频段的频率分辨率总有一个要牺牲。广义S变换Generalized S Transform, GST通过一个随频率变化的窗函数和额外尺度参数在同样的数据上做到“低频看谱、高频看到达时”这让它在地震薄层检测、衰减属性提取里比STFT更常用。gst.zip里这份gst.c就是该算法的C实现适合那些想脱离MATLAB、把时频分析直接嵌进C/C处理流程的工程师和研究者。如果你是做地震资料处理、或是在搞非平稳信号特征提取读这篇能把它编译起来、调好参数、看懂输出并知道结果里的每个数值对应什么物理量。2. GST数学原理与gst.c的算法骨架2.1 从S变换到广义S变换改动到底在哪传统S变换可以写成S(τ,f) ∫ x(t) · (|f| / √(2π)) · exp( - (t-τ)² f² / 2 ) · exp( -i2πft ) dt这里的高斯窗宽度固定为 1/|f|也就是说频率越高时窗越窄。这个特性让S变换在低频时频率分辨率好在高频时时间分辨率好但两个方向的调整都是被动的。广义S变换把窗宽分母里的 |f| 改成 |f|^p其中 p 就是尺度参数S_G(τ,f) ∫ x(t) · (|f|^p / √(2π)) · exp( - (t-τ)² f^(2p) / 2 ) · exp( -i2πft ) dtp1 时它退化成标准S变换p1 时窗随频率变快的速度加大时间分辨率更强p1 时则提高频率分辨率、削弱时间分辨率。gst.c 里这个 p 通常以命令行参数或宏定义的方式暴露出来也正是摘要中所说的“尺度参数”。为什么称“广义”因为标准S变换是 p 固定为 1 的特例。改成 p 带来的直接效果是我们可以根据地层Q值或信号衰减特性主动调节窗宽而不是让窗宽跟着频率粗细随意变化。在地震数据中深层信号的频率普遍偏低此时让时间窗加宽一些反而能稳定提取低频段瞬时属性。对浅层高分辨率目标增大 p 则更容易看清高频到达时。2.2 gst.c 里应当出现的核心模块一个可用的 gst.c 至少包含五个模块输入读取、参数解析、高斯窗构造、GST核心循环、输出转储。下面是一个典型的窗构造函数double *make_gauss_window(int n, double f, double df, double p) { double *w (double *)calloc(n, sizeof(double)); double sigma (f 0) ? 1.0 / pow(fabs(f), p) : 1.0 / pow(df, p); for (int t 0; t n; t) { double tau (t - n / 2) * dt; w[t] (double)(fabs(f) / sqrt(2.0 * PI)) * exp(-0.5 * tau * tau / (sigma * sigma)); } return w; }这里 sigma 是窗宽的时间量纲表示。注意 f0 时要单独处理否则 pow(0, p) 会出问题一般直接用 df 作为低频保护常数。dt 应当作为全局变量传入而不是像网上下到的某些版本那样写成常量。C语言里这类计算密集代码最容易出错的就是浮点边界条件比如 f 为负、n 为奇数、p 为小数时pow的底数不能为负。2.3 双循环完整计算离散化与 FFT按定义直接离散化时GST 的计算可以看作两层外层遍历每个频率点内层遍历每个时间点做卷积。朴素写法是for (int k 0; k n; k) { double f (k n/2) ? k * df : (k - n) * df; if (fabs(f) 1e-10) continue; for (int t 0; t n; t) { double sum 0; for (int m 0; m n; m) { double tau (m - t) * dt; sum x[m] * gauss(tau, f, p) * cexp(-2 * PI * I * f * t * dt); } out[k * n t] sum * dt; } }这个三重循环是教学原型复杂度 O(N³)不用于实际数据。gst.c 如果追求效率会采用“时域乘窗 FFT”的等价方案把信号与高斯窗相乘再做 FFT取对应频点。空间换时间而且能复用现成的 FFT 库。数据规模在几千采样点时朴素写法也能接受但地震道长度通常上万还是建议用 FFT 方案。实际工程里FFT 库有 FFTW、KissFFT 或者直接用 fftwf 的单精度版本gst.c 一般自带一个简单的 radix-2 实现。2.4 输出矩阵的组织方式GST 结果可以按“行对应时间、列对应频率”或反过来存。假设输出为一个 n×n 复数矩阵文件里常见保存格式是每行先输出时间索引然后是全部频率点的模值或实部。读取这种矩阵时最干净的方式是把它当作普通二维数组但关键是要知道保存的是模值还是复数值。如果后续要提取瞬时相位就不能只存模值必须存实部和虚部两列。常见保存格式含义适合用途每行时间索引 所有频率模值振幅谱时频矩阵画等值线图、峰值频率追踪每行时间索引 频率索引 实部虚部复数矩阵瞬时相位、滤波重构按二进制 double 顺序排列FFTW 直接输出大数据量高密度存储gst.c 多数情况下会提供一个简单文本输出方便验证。如果要做生产级处理可以在此基础上加一个-binary选项。若你拿到的版本里输出是“列对应频率”读数据时把矩阵转置一下即可。3. 编译 gst.c 并跑通第一个时频分析3.1 编译前先看代码结构拿到 gst.zip 后先解压unzip gst.zip cd gst ls -la如果只有一个 gst.c没有 Makefile不要慌。用文本编辑器打开 gst.c重点看 main 函数开头的注释或 printf通常写着Usage: gst input_file output_file dt p这比读完整代码快得多。有些版本会把 dt 固定为 1.0并通过编译期宏去改在文件头部会看到类似#define DT 0.004的定义。C语言文件读写操作最常见的就是 fopen/fscanfgst.c 里这两块一般不会写得多复杂但要注意它是否做了文件长度检查。如果输入的采样点数和程序内部预设的 N 不一致绝大多数实现会直接崩溃或产生越界访问这是第一个需要排查的坑。3.2 GCC 编译的命令与依赖最基础的编译命令是gcc gst.c -o gst -lm -O2其中-lm链接数学库因为pow和sqrt都在 libm 里。-O2对浮点密集代码提升明显。如果 gst.c 里用了 FFTW则需要提前安装sudo apt install libfftw3-dev gcc gst.c -o gst -lfftw3 -lm -O2Windows 上则是在 VSCode 里配好 C 语言环境后用gcc gst.c -o gst.exe -lm。要特别注意gst.c 如果在代码开头声明了complex类型而你没有使用 C99 标准编译会报错。在 gcc 后面加-stdc99或-stdgnu11通常能解决问题gcc gst.c -o gst -stdc99 -lm -O2如果报错“cexpundeclared”说明编译器版本默认没有开启 C99 的复数支持同样用-stdc11解决。3.3 生成测试信号并运行我们用 10Hz 正弦波做冒烟测试。生成 512 点数据采样率 200Hzawk BEGIN { for (i0; i512; i) print sin(2*3.14159265359*10*i/200.0); } sine.dat然后运行./gst sine.dat st_out.txt 0.005 0.8其中 0.005 是采样间隔秒0.8 是尺度参数 p。程序若无报错st_out.txt会生成一个矩阵。用head查看前几行head -5 st_out.txt输出第 0 行通常是直流分量第 1 行到第 10 行会有明显峰值。如果第 0 行直接是 0说明程序已经做了去均值处理。如果第 0 行是很大的常数那要注意后续属性提取时把直流列删掉。3.4 验证结果数出峰值位置找每个时间点的最大模值用 Python 最方便import numpy as np # 假设矩阵行时间、列频率 mat np.loadtxt(st_out.txt) max_idx np.argmax(mat[:, 1:], axis1) # 跳过直流列 freqs max_idx * (200.0 / 512) # 索引转频率 print(freqs[:20])期望输出在 10Hz 附近。如果峰值出现在 0Hz多半是直流分量没移除需要把输入信号减掉均值。如果峰值出现在 25Hz说明程序输出矩阵的排列是“行频率、列时间”需要调整索引方向。这一步能快速确认程序的数据布局比直接读全部代码有效。4. 参数调优与常见坑从合成信号到地震道4.1 尺度参数 p 的影响对比p 的常用范围是 0.5 到 1.5超出这个范围时高斯窗不是过宽就是过窄。用 chirp 信号频率从 5Hz 线性扫到 60Hz跑不同 p 值可以看到时频脊的形态变化p 值时窗宽度随频率变化趋势实际效果适合场景0.5窗宽随频率缓慢缩小频率分辨率好时间分辨差深层长时间衰减分析1.0标准 S 变换平衡但高频时间分辨率仍一般默认试跑1.2窗宽快速缩小高频到达时清晰低频谱模糊浅层高分辨率层序解释1.5极度依赖频率只有单频附近少数点有可靠结果特殊 Q 补偿分析建议先以 p1.0 跑一遍再根据目标层主频调整。地震信号主频越低p 应该越接近 1.0 甚至低于 1.0否则高频段的时间窗太窄导致该频段仅有几个有效采样点振幅畸变。经验是主频低于 20Hz 时 p 取 0.7~0.9主频在 30~60Hz 时 p 取 1.0~1.2。4.2 直流分量与负频率的处理很多 gst.c 实现输出矩阵包含全部 n 列其中第 0 列是直流第 1~n/2 列是正频率第 n/21~n-1 列是负频率。如果开发者不处理负频率FFT 得到的复数谱会关于中心对称时频谱上出现“镜像双峰”。处理负频率的正确方式是只计算正频率的 GST 系数对负频率做共轭对称补充或者干脆只输出正频率段。检查方法用单频信号跑出来的时频谱如果对称出现两个峰值说明程序把正负频率都原样输出了。对于地震数据这不算致命但会让后续属性提取混乱因为峰值追踪时可能选到镜像频率。修正时把频率索引大于 n/2 的区域清零即可。4.3 归一化为什么换个采样率幅值就变GST 定义中连续积分对应离散化后必须乘上采样间隔 dt。下面这段代码展示改正out[k * n t] sum * dt; // 乘 dt 而不是直接赋 sum如果不乘 dt当采样率从 100Hz 变成 1000Hz 时输出幅值会差 10 倍。很多论文代码不写这一步因为它们只关心相对幅值。但在地震衰减属性里要比较不同井的频谱差异必须保证绝对幅值正确。验证手段是对 GST 结果做逆变换如果能近似恢复原信号说明归一化正确。实际 gst.c 里我一般会加一个全局scale变量把 dt、窗系数、输出缩放分开管理。4.4 边界效应与延展策略当高斯窗覆盖范围超出信号长度时边界之外按 0 处理会让时频谱在首尾时间点出现非真实能量。解决办法是“两侧补零后再变换变换后裁剪回到原长度”。补零的量至少为窗宽的两倍。更保守的做法是用镜像延展信号这样边界处的相位变化更连续。我的处理流程是输入信号 → 去均值 → 左右各延展100点(镜像) → GST变换 → 裁剪边界 → 输出如果 gst.c 不支持延展可以先在外面做延展再喂给程序。对地震道而言边界效应会影响浅层和深层各几十毫秒的时频属性尤其在做薄层分析时不能忽略。如果输出矩阵的首尾时间点出现明显的强振幅竖条就是边界处理的典型症状。4.5 常见编译与运行错误排查现象原因解决输出全是 nan未初始化数组或 pow 函数域错误用 calloc 分配检查 p 和 f 范围程序直接崩溃输入文件点数与内部 N 不一致读文件后检查长度或动态分配内存峰值频率偏一格频率轴从 1 开始索引导致检查公式里的频偏修正时间轴整体偏移窗函数中心取错应取 (n-1)/2 而非 n/2这些坑在把 MATLAB 算法移植到 C 时几乎都会踩一遍最常见就是“数组维度1”和“FFT 结果顺序”。特别是 FFTW 的输出顺序是正频率从 0 到 Nyquist然后是负频率如果不做 fftshift直接按自然顺序读频率轴就是错的。拿到新 gst.c 不要急着上真实地震数据先用合成信号做回归测试。5. 用 GST 输出做地震同相轴时频属性提取5.1 峰值频率扫描的快速实现峰值频率Peak Frequency是最容易从 GST 时频矩阵中提取的属性。对每个时间索引 t在正频率范围内找到最大模值对应的频率值。下面这段 Python 脚本可以直接读 gst.c 的输出import numpy as np # 假设矩阵已经按行时间、列频率保存 spectrum np.loadtxt(st_out.txt) t_end, n_freq spectrum.shape df 50.0 / (n_freq - 1) # 奈奎斯特频率50Hz仅示例 peak_freq np.zeros(t_end) for t in range(t_end): mag np.abs(spectrum[t, :]) peak_freq[t] np.argmax(mag[1:]) * df这里[1:]跳过直流。如果峰值频率出现大的跳变可以加一个 3 点中值滤波或移动平均抑制由噪声引起的单点抖动。在实际 C 程序里也可以直接在主循环里记录每个时间点的最大模值位置省掉二次扫描。5.2 瞬时带宽估计与时频聚焦度瞬时带宽反映信号的衰减快慢可以用 GST 结果计算每个时间点的二阶矩B(t) sqrt( Σ ( f - fp(t) )² |G(t,f)|² / Σ |G(t,f)|² )从振幅谱矩阵出发一行 numpy 代码就能完成。当目标层下方出现强衰减时带宽通常会抬升因为高频部分被吸收。此时把 p 值调大一些可以让高频段的时间分辨率更好但注意不要引入窗截断伪影。更实用的做法是同时输出“最大振幅频率”和“质心频率”如果两条曲线在某个时间窗出现分叉常被解释为含气层响应这对油气检测很有帮助。5.3 与 CWT 结果交叉验证 GST 参数最后一个常用技巧取同一地震道用 CWT 得到主频曲线再与 GST 峰值频率曲线作差。两条曲线差异应该控制在 5Hz 以内否则说明 p 值选择或边界处理有问题。原因是 CWT 的尺度轴与 GST 的频率轴存在映射关系当 GST 的窗函数相对频带有偏置时峰值频率会系统性偏移。校正方法很简单计算两条曲线的均值差将 GST 频率轴乘一个接近 1 的校正因子重新运行一遍即可。这样能把 gst.c 调到你所在工区的主频附近后续提取的时频属性才有横向可比性。本文还有配套的精品资源点击获取
