简介F-K滤波是地震数据处理中依据频率与传播方向分离信号和噪声的常用手段尤其适合压制地滚波这类低频干扰。面向地震数据处理初学者这份MATLAB实现专注于压制记录中的低频地滚波噪声代码覆盖数据预处理、二维傅里叶变换、滤波窗口定义、频-空域滤波以及逆变换的完整流程并配有滤波过程示意图可对照理解每一步的操作目的。压缩包共2个文件包含1个m脚本和1张jpg效果图大小仅17KB结构简洁便于直接修改调试已有662人学习下载。通过这段代码读者能直观掌握F-K滤波在频-空间域按频率和方位角选择性衰减噪声的思路既能直接套用在实际地震记录上做降噪也能借此熟悉MATLAB在时频分析与信号处理中的典型操作。对于希望动手验证理论、快速入门的初学者来说是一份很有参考价值的实战样本。 干过地震数据处理的人应该都对面波压制这档子事不陌生。面波能量强、频率低、视速度慢混在有效反射波里特别烦人。以前处理这种问题要么在 τ-p 域切要么做 F-K 滤波——说白了就是把数据从时间-空间域搬到频率-波数域把那些“看着慢”的波直接切掉。今天我把在 MATLAB 里实现 f-k 滤波的完整思路、代码和踩坑记录整理出来争取让第一次接触的人也能顺利跑通。这套方法不挑行业地震勘探、探地雷达、声波阵列处理里都一样用。下面进入正题。1. 核心思路为什么要把数据搬到 f-k 域里做文章1.1 f-k 域就是“频率-视速度”联合空间先想一个问题一段地震记录横轴是道位置纵轴是时间你看到一个同相轴从左下角斜到右上角它到底是有效波还是面波肉眼能看出快慢但机器很难直接“画一条线”把慢的切掉。这时候就把它做二维傅里叶变换搬到 f-k 域去。f-k 域有两个轴频率 f单位 Hz和波数 k单位 1/m。一个线性同相轴 t t0 x/v做完二维傅里叶变换后能量会集中在满足 f v·k 的直线上。这里 v 就是视速度。面波视速度低比如 300~700 m/s对应的直线斜率小有效反射波视速度高比如 2000 m/s 以上对应的直线斜率大。所以到 f-k 域这些波直接就分开了。你可以把这个过程想成开运动会入场式所有人本来混在一起走但按“身高”视速度排完队高个儿和矮个儿自然分成两拨。f-k 域就是这个“排队后的操场”。1.2 扇形滤波器是怎么把低速波“切”掉的因为 f v·k波数 k 有正有负分别代表沿测线正向和反向传播所以低速干扰在 f-k 平面上分布在靠近 k 轴、远离 f 轴的两块三角区域里。要压掉速度低于某个阈值的波就要在 f-k 平面上构造一个掩膜mask把满足 |f/k| V_min 的点置零剩下的保留。因为这个掩膜画出来像个扇形低速三角区被挖掉高速区域保留所以叫扇形滤波器。有一点要特别注意如果掩膜边界是一条直线直接硬切逆变换回来会在时间-空间域产生振铃Gibbs 现象表现为同相轴旁边出现一系列“梳子齿”。解决办法是给掩膜边界加过渡带让强度从 0 到 1 平滑过渡。1.3 几个必须记住的硬性参数做 f-k 滤波有几个参数决定了滤波器的适用范围和效果奈奎斯特频率 f_nyq 1/(2dt)。高于这个频率的数据在时间域已经欠采样f-k 域里可能混叠。奈奎斯特波数 k_nyq 1/(2dx)。道间距 dx 越大k_nyq 越小大波数区域越容易混叠。换句话说道距太粗就别指望 f-k 滤波能救回来。视速度边界 V_min 和 V_max。V_min 是你要压制的最大速度V_max 是过渡带结束的速度。有效波速度如果和面波拉开一个量级这个参数特别好定。我在实际处理中一般是先看 f-k 频谱图画出几条等速度线f v·k 的直线看看有效波和干扰分别落在哪些区域再填参数。盲填 V_ min 容易把低速有效波也切掉。2. 动手实现从合成记录到完整滤波流程2.1 数据组织方式和矩阵维度约定MATLAB 里做二维 FFT要先把数据组织成矩阵。我最常用的约定是矩阵行代表时间采样点矩阵列代表地震道。data(nt, nx)nt时间采样点数nx空间道数这一步看起来简单但特别容易搞错。如果从 SEG-Y 文件读数据很多读取函数返回的矩阵是道数 × 采样点数即 data(nx, nt)和上面约定反着来。处理前必须先转置data data否则 FFT 之后 f 轴和 k 轴会对调滤波器方向完全反了。我因为这个吃过一次大亏后面会详细说。2.2 频率轴和波数轴到底怎么算做 f-k 滤波光有 FFT 结果不够还得有对应的 f 轴和 k 轴坐标不然你根本不知道频谱图哪个点对应什么频率、什么波数。时间采样率 dt单位秒道间距 dx单位米。fftshift 之后频率轴和波数轴的计算方式如下f (-nt/2 : nt/2-1) / (nt * dt); % 单位 Hz长度 nt k (-nx/2 : nx/2-1) / (nx * dx); % 单位 1/m长度 nx注意f 轴长度必须等于矩阵行数 ntk 轴长度必须等于矩阵列数 nx。对不上后面 meshgrid 一旦用错滤波器直接歪的。2.3 生成可验证的合成地震记录为了验证滤波器是否工作我习惯先生成合成数据一个视速度 2000 m/s 的有效波叠加上一个视速度 300 m/s 的面波。这样滤波之后理论上应该只剩高速成分。% 参数设置 dt 0.002; % 时间采样间隔 2 ms dx 10; % 道间距 10 m nt 512; % 时间采样点数 nx 128; % 道数 t (0:nt-1) * dt; x (0:nx-1) * dx; % 初始化数据 data zeros(nt, nx); % 有效反射波视速度 2000 m/s主频 30 Hzt0 0.2s for ix 1:nx t0 0.2 x(ix) / 2000; tau t - t0; idx abs(tau) 0.06; data(:, ix) data(:, ix) ... cos(2*pi*30*tau(idx)) .* exp(-(tau(idx)/0.03).^2); end % 面波视速度 300 m/s主频 10 Hzt0 0.1s for ix 1:nx t0 0.1 x(ix) / 300; tau t - t0; idx abs(tau) 0.12; data(:, ix) data(:, ix) ... cos(2*pi*10*tau(idx)) .* exp(-(tau(idx)/0.06).^2); end生成完可以先画一张 wiggle 图确认数据和实际地震剖面形态一致既有高频往右下方倾斜的有效波又有低频更“躺”的面波。3. 核心代码可以直接抄的 MATLAB 实现3.1 完整滤波函数 fk_filter下面是封装好的滤波函数。输入原始数据、时间采样间隔、道间距、速度边界和频率范围输出滤波后的数据。参数设计上保留了过渡带控制项具体怎么选后面会说。function data_filt fk_filter(data, dt, dx, vmin, vmax, fmin, fmax, taper) % FK_FILTER F-K 域扇形滤波 % 输入 % data - 输入数据矩阵维度 (nt, nx)行为时间列为空间 % dt - 时间采样间隔秒 % dx - 道间距米 % vmin - 完全压制速度阈值米/秒。视速度低于此值的成分被切除 % vmax - 完全保留速度阈值米/秒。视速度高于此值的成分完整保留 % fmin - 带通低截止频率赫兹 % fmax - 带通高截止频率赫兹 % taper - 过渡带控制量0~1。越大过渡带越宽振铃越轻 % 输出 % data_filt - 滤波后数据矩阵维度同 data [nt, nx] size(data); % 1. 二维 FFT并中心化 F fftshift(fft2(data)); % 2. 建立频率轴、波数轴网格 f (-nt/2 : nt/2-1) / (nt * dt); k (-nx/2 : nx/2-1) / (nx * dx); [f_grid, k_grid] meshgrid(k, f); % 3. 构造速度掩膜扇形区域 k_safe k_grid; k_safe(abs(k_safe) 1e-10) 1e-10; % 防止除零 ratio abs(f_grid ./ k_safe); % 视速度 mask ones(nt, nx); % 低速小斜率区域完全压掉 idx_cut ratio vmin; mask(idx_cut) 0; % 过渡带区域从 0 平滑过渡到 1余弦过渡 delta_v max(vmax - vmin, eps); idx_trans (ratio vmin) (ratio vmax); mask(idx_trans) 0.5 * (1 cos(pi * (ratio(idx_trans) - vmin) / delta_v)); % 4. 附加带通限制可选 idx_f (abs(f_grid) fmin) | (abs(f_grid) fmax); mask(idx_f) 0; % 5. 保护 k0 轴视速度为无穷大的垂直同相轴一般要保留 idx_axis abs(k_grid) 1e-6; mask(idx_axis) 1; % 6. 应用滤波并逆变换 F_filt F .* mask; data_filt real(ifft2(ifftshift(F_filt))); end这里有两个细节要解释一下。第一过渡带用的是余弦形式比线性过渡更平滑在我自己的测试里振铃更小。taper 参数可以控制过渡带宽度吗上面代码里过渡带宽度直接由 vmin 和 vmax 差值控制taper 参数其实是用在另一个实现版本里的这里为了不误导函数输入端 ta per 其实没直接用。我实际用的时候会直接调 vmin/vmax 的差值来控制过渡带宽度差值越大过渡越缓。你如果希望保留 taper 参数可以在构造过渡带积分的时候再乘上它但就目前这个函数来说vmax-vmin 已经足够你调了。第二保护 k0 轴这一点非常关键。在 f-k 平面上k0 表示视速度无穷大也就是所有道同时到达的那种水平同相轴。这种成分通常是折射波或直达水平能量一般不该被低速滤波器干掉。加一行保护能避免滤波把这类信号误伤。3.2 主程序调用流程调用方式很简单。接着上面生成的合成数据% 执行 f-k 滤波 data_filt fk_filter(data, dt, dx, 500, 700, 5, 80); % 显示三维对比或直接画 wiggle figure; subplot(1,2,1); imagesc(x, t, data); axis xy; title(滤波前); subplot(1,2,2); imagesc(x, t, data_filt); axis xy; title(滤波后);滤波前后的时间-空间剖面一对比效果非常直观低频、低视速度的面波被压掉了高速的反射波同相轴还在。如果你要观察 f-k 谱可以在函数外部单独做amp abs(fftshift(fft2(data))); amp_log 20*log10(amp eps); figure; imagesc(k, f, amp_log); axis xy; xlabel(波数 k (1/m)); ylabel(频率 f (Hz));在图上你加一条速度线确认滤波边界hold on; k_line linspace(0.5e-3, 20e-3, 100); f_line 500 * k_line; plot(k_line, f_line, r--, LineWidth, 1.5);这条红色虚线就是“视速度等于 500 m/s”的分界线。你会在图上看到面波能量集中在虚线以下有效波能量在虚线以上。这样调参数心里就有底了。4. 实测翻车记录这 5 个坑我全踩过理论归理论实际跑代码时问题特别多。我把踩过的坑整理成一份排查表希望你在复现的时候少走弯路。4.1 fftshift/ifftshift 用混恢复出来的数据是乱的这是我犯过最蠢的错误。二维 FFT 后我用 fftshift 把零频移到中心滤波完之后直接用了 fftshift 恢复原点没注意 MATLAB 里偶数长度下 fftshift 和 ifftshift 恰好相同一旦数据长度变成奇数比如有些数据切了某段道数两者就不一样了恢复出来的时间序列整体错位同相轴连续性完全破坏。这里的建议是固定写成 ifftshift(F_filt)不要嫌麻烦。习惯成自然任何中心化操作结束后都用对应的逆操作。4.2 滤波器边缘太陡结果全是振铃第一次写掩膜时我没有过度带直接在 ratio vmin 处硬切到 0。滤波结果一出来面波是没了但每个保留同相轴旁边多了一排平行“鬼影”画出来锯齿感严重。这是因为掩膜在频域产生一个陡峭台阶等效到时域是一个无限长的 sin c 函数和信号卷积自然就产生振铃。后来加上过渡带就正常了。过渡带宽度不要太小vmin 和 vmax 相差 100~300 m/s 都可以。滤波效果要是仍然不满意优先加大过渡带别一味去把 vmax 设成很大。4.3 道间距不均匀f-k 谱一塌糊涂野外采集的 SEG-Y 数据不是每道都能对准坏道、空道、地形原因造成的偏移都会让空间方向出现不规则采样。直接把这种数据做二维 FFTf-k 谱上到处都是“假能量”低速和高速波混在一起滤波器没法干净地切割。我的经验是先做道编辑或道插值规则化再进 f-k 滤波。现实处理中确实有一类流程是把 f-k 滤波放在规则化之后道理就在这。4.4 k0 轴上直流分量被误杀输出出现整体水平偏移扇形滤波器用 abs(f/k) 判断速度当 k0 时速度是无穷大理论上应该被保留。但代码里一旦直接做除法k0 要么产生 NaN要么被 ratio 判断成很大的数被保留要么因为 f0 和 k0 同时出现变成 NaN然后 mask 判断出错把这条轴上的能量全清零。结果就是滤波后的剖面每个道都少了一个常数分量表现为整道水平抬升或降低叠到叠加速度谱上会是一个高频假象。解决方式就是前面代码里的 k_safe 保护以及最后强制 idx_axis 置 1。这两个步骤缺一不可。4.5 忘记取实部输出结果带着小虚部逆 FFT 之后由于浮点误差结果矩阵多多少少有点虚部。直接拿去画图MATLAB 会给你报警告或者画出来全是噪点。所以在函数最后一行用 real() 包一下。取实部是标准操作别觉得多余加上它后续处理全程清爽。为了方便你排查我把这些问题的现象和解决办法汇总成一张速查表症状可能原因解法滤波后数据整体错位、同相轴跳变ifftshift 用错统一用 ifftshift 做逆中心化同相轴旁出现一排“鬼影”滤波器边缘太陡vmin/vmax 拉大加过渡带f-k 谱上高低速混叠严重有坏道或空间采样不均匀先做道插值/规则化再滤波滤波后整道水平偏移k0 轴能量被误删对 k0 轴单独置 1 保护剖面有噪点/警告虚部未取实部输出前 real(data_filt)5. 关于参数调试和项目扩展的几点个人经验参数怎么定这是大家问得最多的问题。我给的参考思路是先在 f-k 谱上拉三条速度线分别对应面波区、过渡区、有效波区。比如从谱上看到面波能量分布在 200~600 m/s 之间有效波从 1200 m/s 开始那 vmin 定 800 比较安全vmax 定 1200。这样面波被干净切除有效波完整保留过渡带刚好覆盖中间模糊区。另外这个滤波器函数稍微改一改还能做别的用。比如把 mask 反过来保留低速区、压掉高速区就成了“低速波分离器”可以用来提取面波之后做面波衰减或面波成像。再把掩膜换成只保留正 k 或负 k 区域就能分离正向和反向传播的波场这个在地震槽波处理和上下行波分离里也经常用到。最后分享一个小技巧真实数据处理前先用合成记录把参数摸透再上实际数据。合成记录的好处是答案已知你一眼就能看出滤波结果哪里对哪里不对比直接拿真实数据瞎调高效得多。F-K 滤波不是万能的它对空间采样均匀性要求高对高低速差异不够大时也很吃力。但只要参数设对、保护加到位它就是手里最趁手的那把刀。本文还有配套的精品资源点击获取