企业官网建设流程全解析

做地震数据处理的人，几乎没有人能绕开频率波数域变换（f-k域变换）这个工具。我最早接触它是在某区块三维地震资料的线性干扰处理上，当时面波和声波混在炮集里，常规带通滤波完全压不住，后来把数据变换到f-k域，用了一个简单的扇形滤波器，效果立竿见影。从那以后，f-k域去噪就成了我处理流程里的常客。这篇文章把整套方法从原理到MATLAB实现完整讲一遍，包括我踩过的坐标陷阱、振铃问题，以及最终沉淀下来的可复现代码，给你一个能直接套用到自己数据上的参考。

1. 从x-t域到f-k域：为什么倾斜同相轴是一根"斜线"

1.1 两次傅里叶变换叠出来的二维频谱

频率波数域变换本质上就是二维傅里叶变换。对地震道集这个二维矩阵（第一维是时间t，第二维是空间道号x）做二维傅里叶变换，我们就得到了f-k域振幅谱。一个最简单直观的理解方式是把过程拆成两步：先对每一道做时间方向的一维傅里叶变换，得到频率域的分量；再对每一个频率分量沿空间方向做一维傅里叶变换，得到波数域的分量。两步叠在一起，就得到了 $U(f,k)$。

数学表达式是：

$$U(f,k)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} u(t,x) e^{-i2\pi(ft+kx)} dt dx$$

这里面f是时间频率，单位Hz；k是空间波数，单位是1/m，表示能量沿空间方向的周期性变化快慢。很多初学者会把波数想象成"空间频率"，这个类比很准确：就像时间频率描述信号随时间震荡的快慢，波数描述地震波场沿测线方向震荡的快慢。

f-k域的关键价值在于：它把x-t域里"看起来都差不多的波形"按视速度重新排列了。一个以视速度v传播的线性同相轴，比如面波、直达波，在x-t域里是一组彼此平行、有一定倾斜角度的波形；变换到f-k域后，它的能量会高度集中到一条穿过原点的直线上，这条直线的斜率正好等于视速度的倒数。斜率越缓，视速度越低；斜率越陡，视速度越高。

这个特性是f-k域去噪的根本依据。不同视速度的波在二维谱上占据不同的方向区域，因此我们可以用带方向性的滤波器把它们分开。这是时间域滤波和单道频率域滤波做不到的。

1.2 为什么去噪要选f-k域而不是单独的时间域滤波

单道滤波的一个天然局限是"只看一个维度"。假设有效反射和线性干扰恰好都分布在30Hz附近，频率域带通滤波只能把这两个频率成分一起滤掉或一起保留，完全没有区分能力。但它们的视速度可能是完全不同的：一次反射波视速度通常在1500m/s以上，面波视速度只有300到1000m/s，声波大约340m/s，折射波可能有一千到两千。只要视速度不同，它们在f-k域里就落在不同的方向区域。

所以f-k滤波本质上是"速度域滤波"。它利用的是地震信号与噪声在视速度上的差异，而不是仅频率上的差异。对于线性相干干扰这种"频率重叠但速度不重叠"的噪声，f-k滤波几乎是首选。

还有一个实际原因：f-k域实现极其简单。二维快速傅里叶变换在MATLAB里就是一个fft2调用，滤波也只是频谱点乘掩膜，不存在复杂的最优化求解，参数物理意义非常明确。相比之下，f-x域预测滤波、小波变换去噪、曲波变换等方法各有优势，但一旦需要快速迭代、直观调整，f-k域仍然是最顺手的工具之一。

2. MATLAB里做f-k正反变换：先避开三个坐标陷阱

2.1 最简可运行的变换代码

在动手设计滤波器之前，先把正变换、显示、反变换这套基础流程跑通。下面这段代码我用得最多，它完成了完整的"x-t域 -> f-k域 -> 反变换回x-t域"流程：

% data: 二维地震道集, 行为时间, 列为道 % dt: 时间采样间隔(秒), dx: 道间距(米) nt = size(data, 1); nx = size(data, 2); % 正变换 spec = fft2(data); spec_shift = fftshift(spec); % 把零频与零波数移到矩阵中心 % 构建频率轴与波数轴 freq = (-nt/2 : nt/2-1) / (nt * dt); % Hz k = (-nx/2 : nx/2-1) / (nx * dx); % 1/m % 显示振幅谱(对数幅度) figure; [KK, FF] = meshgrid(k, freq); amp = 20 * log10(abs(spec_shift) + 1e-8); imagesc(k, freq, amp); xlabel('波数 k (1/m)'); ylabel('频率 f (Hz)'); axis xy; % 让y轴频率从低到高 colorbar;

反变换同样简单：

% 假设spec_filtered是我们在f-k域做过滤波后的频谱 data_recon = real(ifft2(ifftshift(spec_filtered)));

这里面有三个容易出错的坐标细节，我逐个说。

2.2 陷阱一：fft2的维度顺序到底谁是第一维

fft2是按矩阵维度做变换的，第一维对应行、也就是时间方向，第二维对应列、也就是空间道方向。所以变换结果spec里，行方向是频率f，列方向是波数k。这个顺序和后面meshgrid构造坐标矩阵时必须严格一致。

很多人在这一步栽跟头，因为习惯性把"道"放在前面。如果数据矩阵是nx行、nt列存储，那么fft2之后第一维反而成了波数、第二维成了频率，滤波掩膜也要跟着转置，这是最常见的坐标错乱来源。我的建议是：从一开始就统一用"nt行、nx列"的存储方式，并在代码里加一行注释写清楚。

构造坐标网格时，不要直接写[F, K] = meshgrid(freq, k)，因为meshgrid的输出行数等于第一个输入长度、列数等于第二个输入长度，而spec是nt行、nx列。正确写法是：

[KK, FF] = meshgrid(k, freq);

这样FF和KK都是nt行、nx列，FF(i,j)对应spec中的第i个频率、KK(i,j)对应第j个波数。

2.3 陷阱二：freq和k轴的原点、方向与奈奎斯特边界

频率轴和波数轴都有一个奈奎斯特极限。时间采样间隔dt决定了最大有效频率是1/(2dt)，比如dt=2ms时最高频率是250Hz；道间距dx决定了最大有效波数是1/(2dx)，比如dx=10m时最高波数是0.05 1/m。超过这个范围的频率或波数都会发生折叠，也就是混叠。

MATLAB的fftshift把零频挪到矩阵中心，所以shifted之后频谱的坐标范围是从负的奈奎斯特值到正的奈奎斯特值附近。用(-nt/2 : nt/2-1)这样构造轴是常规写法，注意当nt为偶数时没有正的中心点；如果nt是奇数，则需要用(-(nt-1)/2 : (nt-1)/2)。为避免这种边界差异，我通常建议把道集长度凑成偶数，处理流程会省心很多。

振幅谱显示时还有个细节是取对数后加一个小常数，否则零值处的log(0)会显示成黑块，影响观察。

2.4 陷阱三：fftshift与ifftshift不是什么时候都能互换

在MATLAB里，对于偶数长度数组fftshift和ifftshift效果一样，但奇数长度时不一样。规范做法是正变换使用fftshift，反变换使用ifftshift，不管长度奇偶都按这个对应关系写。虽然地震数据长度几乎都是偶数，但既然养成正确习惯不费事，就坚持用ifftshift。

验证坐标是否正确的保险方法是做一次"圆整测试"：

spec_test = fft2(data); data_back = real(ifft2(ifftshift(fftshift(spec_test)))); max(abs(data_back(:) - data(:)))

这个数值应该接近机器精度。如果差得很大，说明fftshift和ifftshift的使用有问题，或者正反变换之间多了一重/少了一重shift。

3. 各类噪声在f-k域的长相：识别比滤波更重要

f-k域滤波效果好不好，很大程度上取决于你能不能一眼看懂频谱上的能量分布。拿到一张f-k谱，首先要问的并不是"该用什么滤波器"，而是"这些能量团分别代表什么波"。

3.1 线性相干干扰：面波、声波、浅层折射

线性干扰的共同特征是x-t域里同相轴近似直线，视速度恒定。这种波在f-k谱上是穿过原点的一条能量带，斜率为视速度的倒数。

面波视速度低，一般300到1000m/s，在谱上表现为靠近水平方向、斜率很小的窄带能量，通常集中在低频段。声波速度快一些，空气中约340m/s，水中约1500m/s，能量带比面波陡。浅层折射视速度变化大，但同一批次炮记录里往往也有固定的优势视速度。

识别这些干扰时，我习惯先看原始炮集的x-t显示，估计干扰的视速度范围，再到f-k谱上对照找对应的能量条带。如果在x-t图上看到一组斜率一致、频率上明显低于有效反射的强线性波，那么在f-k谱上它的能量条带一定落在某个扇形区域内，这样后续滤波器参数就有了依据。

3.2 随机噪声与异常振幅脉冲：背景地毯和十字线

随机噪声在f-k域的特征是"铺满全谱"。如果噪声是白噪声，它的二维谱近似均匀；如果是有色噪声，能量会集中在某个频率范围，但方向上没有偏好。这种背景能量会抬高整个谱的底噪，但不会形成明显的条带或团块。

异常振幅脉冲噪声则完全不同。如果一个高频脉冲只是出现在单道或者极少数道上，它在空间方向上是窄的，所以会在f-k谱上沿波数轴展宽成一条"竖向"能量线；如果脉冲在时间上也短，还会沿频率轴展宽，最后形成十字形或放射状能量线。这种形态很容易识别，处理时可以用矩形切除或者中值滤波配合，而不是用扇形滤波器。

3.3 多次波与规则噪声的谱形态差异

多次波的问题要复杂一些。多次波的时距曲线也近似双曲线，经过NMO校正后会变得接近水平，因此如果直接在原始炮集或CMP道集上做f-k滤波，多次波和一次波在谱上往往是重叠的，单靠f-k切除很难干净分离。

实际项目中处理多次波，我更多是在NMO校正之后做f-k域滤波：一次波NMO后近似水平，视速度接近无穷大，能量集中在k=0附近；剩余动校正量较大的多次波仍然倾斜，能量偏到波数轴两侧。这样就能用扇形滤波器保留低速部分、切除中高速倾斜能量。当然这种做法要配合反NMO使用，流程上多两步，但效果往往比直接切除好很多。

4. f-k域去噪滤波器设计：从扇形切除到局部陷波

4.1 扇形滤波器：按视速度通放带切出有效波

扇形滤波器是f-k去噪里最常用的工具。它的思想是：给定视速度阈值vmin和vmax，保留所有视速度在区间内的能量。因为f和k之间的比值就是视速度，所以边界在谱上是两条穿过原点的直线。

下面这段代码构造了一个完整的扇形通放掩膜：

function mask = fan_filter_mask(k, freq, vmin, vmax, taper_frac) [KK, FF] = meshgrid(k, freq); % 避免除零 KK_safe = KK; KK_safe(KK == 0) = 1e-12; ratio = FF ./ KK_safe; % ratio = f/k, 即视速度 mask = zeros(size(FF)); % 核心通放带: vmin < ratio < vmax mask(ratio >= vmin & ratio <= vmax) = 1; % 在边界加渐变过渡, 减少振铃 if taper_frac > 0 taper_low = max(vmin*(1 - taper_frac), vmin*0.5); taper_high = vmax*(1 + taper_frac); % 下边界渐变 transition = (ratio - taper_low) / (vmin - taper_low); transition(transition < 0) = 0; transition(transition > 1) = 1; transition = 0.5 - 0.5*cos(pi * transition); mask(ratio > taper_low & ratio < vmin) = transition(ratio > taper_low & ratio < vmin); % 上边界渐变 transition_high = (taper_high - ratio) / (taper_high - vmax); transition_high(transition_high < 0) = 0; transition_high(transition_high > 1) = 1; transition_high = 0.5 - 0.5*cos(pi * transition_high); mask(ratio > vmax & ratio < taper_high) = transition_high(ratio > vmax & ratio < taper_high); end end

使用方式很简单，正变换得到spec_shift后点乘掩膜：

mask = fan_filter_mask(k, freq, 1500, 8000, 0.15); spec_filtered = spec_shift .* mask; data_clean = real(ifft2(ifftshift(spec_filtered)));

注意mask的维度必须和spec_shift完全一致，这也是前面坐标陷阱提到的原因。

4.2 矩形切除与定向陷波：把特定噪声"抠掉"

扇形滤波器解决的是"有效信号视速度区间"问题，但实际数据里往往还有局部强噪声。比如一组强线性干扰的频谱能量条带虽然落在扇形通放带内，但强度远高于周围信号，这时候用扇形滤波器反而会把它保留下来。

我的处理习惯是先用扇形滤波器做粗去噪，拿到结果后对比剩余频谱。如果还能看到明显的局部能量团，就用矩形切除或任意多边形掩膜做局部陷波：

% 例如切除波数0.005~0.02、频率5~20Hz的一块区域 mask_local = ones(size(spec_shift)); % 找到对应索引区间 band = abs(KK) >= 0.005 & abs(KK) <= 0.02 & ... FF >= 5 & FF <= 20; % 在区域边缘加简单的余弦过渡 mask_local(band) = 0.3; spec_filtered2 = spec_shift .* mask_local;

局部切除比扇形滤波更灵活，但也更容易引入局部振铃。我强烈建议不要直接硬置零，而是把切除区域内幅值衰减到某个比例，或者做渐变过渡。硬置零会导致频谱突变，反变换后噪声道周围会出现周期性假象。

4.3 切除之后的重建：硬边界振铃的抑制

频谱域的硬边界切除，在反变换回x-t域后必然会产生吉布斯振铃。表现就是信号边缘出现等间隔的伪波形，尤其在强反射同相轴两端。这个问题是所有频率域滤波的共性难题，f-k域也同样面对。

抑制振铃我常用的手段有三个：

第一是掩膜渐变。前面代码里的taper_frac参数就是干这个的，过渡带越宽，振铃越弱，但滤波选择性会下降。一般取0.1到0.2比较均衡。

第二是时窗衰减。滤波前对道集两端做时间方向的余弦斜坡衰减，把数据先变到零附近再变换，可以显著降低频谱泄漏引起的振铃。

第三是空域混合。把滤波后的结果和原始道集做加权混合，在信号强、噪声弱的地方多用原始数据，在噪声强的地方用滤波数据。这种做法保幅性更好，但需要额外的质量控制步骤。

5. 合成道集上的完整测试：滤波器参数如何定

5.1 构造带噪声的合成记录

在实际数据上调试参数之前，先用合成记录验证滤波器逻辑是效率最高的方式。下面这段代码生成一个包含双曲线有效反射、线性干扰、随机噪声的合成道集：

nt = 512; nx = 64; dt = 0.002; dx = 10; t = (0:nt-1)*dt; x = (0:nx-1)*dx; data = zeros(nt, nx); % 有效反射: 双曲线同相轴, 速度2500m/s, 零偏移距时间0.2s v_ref = 2500; t0 = 0.2; for ix = 1:nx tr = sqrt(t0^2 + (x(ix)/v_ref)^2); data(:,ix) = data(:,ix) + ... exp(-((t-tr)*50).^2) .* sin(2*pi*30*(t-tr)); end % 线性干扰: 视速度800m/s, 主频10Hz v_noise = 800; for ix = 1:nx tn = x(ix)/v_noise; data(:,ix) = data(:,ix) + ... 0.8 * exp(-((t-tn)*30).^2) .* sin(2*pi*10*(t-tn)); end % 随机噪声 data = data + 0.1*randn(nt, nx);

这个模型贴近真实炮集的基本组成：有效反射是双曲线、线性干扰是恒定视速度斜直线、随机噪声让频谱背景抬高。

5.2 f-k滤波后信号恢复效果评估

用前面写的fan_filter_mask处理，选择保留视速度1500m/s以上能量：

[KK, FF] = meshgrid(k, freq); spec = fftshift(fft2(data)); mask = fan_filter_mask(k, freq, 1500, 8000, 0.15); data_clean = real(ifft2(ifftshift(spec .* mask)));

如果合成记录的真实无噪信号已知，可以直接计算滤波前后的信噪比：

% 构造无噪信号 clean = zeros(nt, nx); for ix = 1:nx tr = sqrt(t0^2 + (x(ix)/v_ref)^2); clean(:,ix) = clean(:,ix) + ... exp(-((t-tr)*50).^2) .* sin(2*pi*30*(t-tr)); end snr_before = 10*log10(sum(clean(:).^2) / sum((data-clean).^2)); snr_after = 10*log10(sum(clean(:).^2) / sum((data_clean-clean).^2)); fprintf('滤波前SNR: %.2f dB, 滤波后SNR: %.2f dB\n', snr_before, snr_after);

我实测这样的合成例子，滤波前信噪比大约在一两个dB量级，滤波后可以提升到10dB以上。要注意的是，滤波器如果过度切除，会把有效反射的高视速度成分也去掉，信噪比反而下降。所以在合成记录上多做几组参数试验，能帮你形成一个直观感受：到底多大的扇形角度是安全的。

5.3 保幅特性与边界影响的平衡

地震资料处理的后期要做AVO分析或反演，所以前面滤波步骤尽量保幅。f-k滤波是一个线性算子，对振幅的影响是确定的、可预测的，但也存在边界效应。数据边缘道的能量在空间方向傅里叶变换时会被摊薄，造成边缘道滤波前后振幅失真。

在实际处理中，我对保幅性的控制思路是：先做一次f-k滤波得到噪声模型，再用原始数据减去噪声模型得到剩余信号。也就是把"从数据中切掉噪声"变成"估计噪声并从数据中减去"。这样做的好处是滤波器的切除尺度不需要设计得很极端，只要能把大部分噪声特征估计出来就行，剩余能量里保留下来的有效信号更完整。

具体实现就是：

mask_noise = 1 - mask; % 只保留噪声区 noise_model = real(ifft2(ifftshift(spec .* mask_noise))); data_clean2 = data - noise_model;

这个思路和个性化去噪的习惯很接近，能保留更多有效信号细节。

6. 实际地震资料中绕不开的坑与我的处理习惯

6.1 道间距不均匀和空道：先规则化再变到f-k域

f-k变换隐含假设是空间方向均匀采样。实际采集数据很少完美满足这个条件：坏道、空道、野外变观造成的道间距抖动，都会让f-k谱出现虚假能量。

我的经验是：处理前先做道编辑和规则化。空道用相邻道插值补齐，或者至少充零并在频谱中做相应处理；道间距不均匀的，先内插到统一网格。否则你在f-k域看到一个强能量条带，很可能不是地质信号也不是真实噪声，而是采样不均匀造成的假频。

6.2 空间混叠：高频远偏移道折叠到低视速度区

空间混叠是f-k去噪最容易踩的坑。当地震信号频率较高、道间距较大时，高波数成分会折叠回低波数区域，在谱上表现为看起来像低速噪声的假能量。比如本来真实的视速度是2000m/s的反射，折叠后可能出现在600m/s的扇形区里，这时候盲目用扇形滤波器切除低速区，会把有效信号误伤。

判断是否混叠的一个实用方法是看频率和道间距是否满足 $\Delta x < v_{app} / (2 f_{max})$。比如目标视速度1500m/s、最高频率60Hz，那么道间距不能大于12.5m。如果实际道间距大于这个值，要么先做空间插值加密道距，要么把处理频率上限限制到安全范围。

6.3 大数据分块处理：重叠与拼接

三维数据体做f-k滤波时，一次性fft2一整块数据往往内存不够，或者计算时间过长。我的做法是按炮集按偏移距分块处理。分块时有两个细节：

一是块与块之间要有重叠。因为f-k滤波在块边界会产生振铃，两个相邻块分别滤波后如果直接拼接，边界上会出现明显的接缝。我通常让相邻块在空间上重叠10到20道，滤波后只取中间部分，弃掉靠近块边缘的若干道，再用线性斜坡混合重叠区域。

二是块大小要方便FFT。MATLAB的fft在块大小为2的幂时效率最高，所以把每个块的道数和时间采样点凑成512、1024这类长度可以显著加速。

6.4 参数记忆与作业标准化

最后分享一个我觉得很重要的习惯：把每次处理的参数记下来。哪个工区、哪一批炮集，用了什么样的扇形速度区间、过渡带宽、是否做规则化，这些信息都整理成表格。因为f-k滤波的参数挑选高度依赖数据质量，一批数据的最佳参数换到另一批数据上可能完全不可用。有记录才能回溯、对比、复用。

我现在做实际项目，f-k域去噪的流程基本固定成四步：规则化与道编辑，频谱检查与参数标定，扇形滤波或陷波切除，质量控制与噪声模型差减。每一步都有对应的质量图件。这套流程用下来，处理效率和数据可靠性都比以前随手调参数高得多。如果你正打算在自己的数据集上尝试f-k去噪，建议也从这四步开始，而不是直接跳到最后一步。