企业官网建设流程全解析

简介：基于MATLAB的FK变换（傅里叶-基尔霍夫变换）完整资源包，面向光学成像、遥感图像处理、医学成像等方向的学习者与开发者，旨在解决复杂光学系统成像特性的分析与仿真问题。资源共10个文件，总大小仅926KB，类型覆盖4个docx文档、3个mat数据文件、1个m脚本、1个txt说明及1张运行结果图片。docx文档系统梳理FKT概述、程序手册与使用说明，m脚本fktran.m提供核心算法实现，mat文件保存处理所需数据，txt与jpg则可快速核对代码运行流程与结果。目前已有1779人学习。借助这套资源，读者可快速理解FK变换中二维傅里叶变换、光瞳函数卷积与逆变换的完整链路，直接运行脚本复现成像效果，并可根据文档示例将方法迁移到遥感图像理解、光学系统设计或OCT图像重建等实际任务中，兼具理论参考与工程实践价值，也可结合自身光学系统参数进行二次开发。 我们直接进入正题。

1. 为什么要从FK变换开始聊

做地震数据处理、探地雷达分析甚至实验力学振动信号整理的朋友，大概率都绕不过FK变换这个坎。我第一次接触FK变换是在处理二维地震记录时，目标很单纯：把斜向传播的干扰波（比如面波、声波）干掉，把同相轴更平缓的反射波保留下来。当时我用的就是matlab，一条fish命令来回试，最后跑通了，也踩了不少坑。

先说清楚FK变换到底处理的是什么问题。常规的FFT处理一维信号，只有时间轴t，得到的是频率f；FK变换处理的是二维信号，除了记录时间t，还有个空间轴x（比如检波点道号、测线位置），两个方向都做傅里叶变换后，得到的就是频率f和波数k（空间频率，代表波数，即每公里/每米有多少个波长）构成的谱。所以FK变换本质上是二维傅里叶变换，帮助我们把信号放到频率-波数平面上去区分波场。

为什么要在频率-波数域做？因为地震记录里的不同波场，在f-k平面里的分布区域有明显差异。反射波视速度大（同相轴平缓），能量集中在靠近k轴的小斜率区域；面波视速度低（同相轴陡峭），能量集中在斜率较大的锥形区域。两者在时间-空间域可能重叠得很厉害，但在f-k域可以分开，于是就可以用滤波算子在f-k平面上做切除，然后再逆变换回t-x域。这就是FK滤波“治“波场分离的最常见用法。

如果你正在做数据预处理、去噪、插值，或者只是想搞清楚二维阵列数据的频率成分，这篇内容应该能帮你少走不少弯路。我们以下都用matlab为操作环境，从物理原理讲到实际代码，再到排错经验，一步步拆开说。

2. FK变换的物理意义与matlab实现前的准备

2.1 直观理解f-k域里的“一个点代表一种波”

先做个简化比喻。一维信号FFT后，某个频率f上的能量大小说明信号里这个快慢变化的成分有多强；二维数据FK变换后，某个（f，k）坐标上的能量说明的是：这个数据里有多少能量是以“传播速度为f/k”这种特定的时空模式在移动。

因为速度v = 频率f / 波数k。于是，在f-k平面上，从原点出发的任意一条射线的斜率，就代表了波场的视速度。斜率越陡，速度越低；斜率越缓，速度越快。例如水平反射波同相轴在时空域是近水平的，FFT之后能量集中在k≈0附近，对应视速度接近无穷大，在f-k图上就是贴着f轴的一条窄带。而倾斜干扰波同相轴，例如线性面波，视速度是固定的几百米每秒，它的能量会沿一条射线分布。有了这个图像，我们就能设计扇形滤波器，把不需要的射线区域切掉。

这里必须提一个容易混淆的点：不要以为f-k域和t-x域是割裂的两个空间。事实上它们只是同一份波场数据的两种正交展开方式。你做的任何f-k域操作都会影响到时空域的结果，所以理解“域里面切了什么=域外面改了什么”很重要。

2.2 matlab实现前必备的三个基础操作

在写任何代码前，有准备工作要做。

第一步，数据摆放格式。FK变换要求输入是二维矩阵，而且建议是空间轴按行排列、时间轴按列排列，或者反过来，但必须保持恒定，因为它直接影响fft2输出的含义。习惯上地震行业用nt（时间采样点数）×nx（空间道数）的矩阵，行是固定一道的时间序列，列是固定时刻的各道振幅。绝大多数Segy读取工具导出的就是这种格式。

第二步，了解fft2的输出布局。matlab的fft2对矩阵做二维FFT，默认输出的第一个维度对应矩阵的行方向（即时间方向），第二个维度对应列方向（即空间方向）。很多坑都是这里产生的。如果你习惯用imagesc直接看abs(fft2(data))，会发现低频在四角，零频在左上角，视觉上反直觉。建议用fftshift做一次移位，把零频挪到矩阵中心，再显示，否则后面画扇形滤波边界的时候，坐标关系很容易搞乱。

第三步，搞清坐标轴。显示f-k谱时，横轴是波数k，纵轴是频率f（或者采样点数）。需要根据采样间隔dt、道间距dx，把数组坐标换算成物理坐标：

• 频率轴：f = (0 : nt-1) / (nt * dt)，单位Hz；
• 波数轴：k = (0 : nx-1) / (nx * dx)，单位1/m，或者换算成cycles/m；
• 如果用fftshift，需要生成从负到正的物理坐标范围：f = (-nt/2 : nt/2-1) / (nt*dt)，k同理。

这一步很多人偷懒不做，直接拿索引坐标去设计滤波器，结果滤波边界完全对不上波场位置，导致该去的没去掉，不该削的反射波反而被削了。

2.3 为什么有人用f-k插值而不是单纯去噪

除了压制干扰波，FK变换还有一个高频用法：数据规则化与插值。规则化是要把不规则采样的道间距变成等间距，插值是要把缺道补出来。时间域直接插值容易破坏波场的动力学特征，而f-k域可以利用波场的可预测性，在k方向做带宽限制或者反泄露处理后再反变换。

简单说原理：如果原始空间采样不足，存在空间假频（见后面第5节），插值算法在f-k域能识别哪些是假频能量、哪些是真实信号能量，通过迭代或者抗泄露算子把落在假频区的信号重建出来。这种用法在处理地面地震数据、GPR（探地雷达）数据甚至麦克风阵列信号时都非常常见。

3. matlab中FK变换的核心代码与参数设计

这部分我们写一个可以直接运行的示例流程。假设数据data变量是一个nt×nx的矩阵，代表了某个时刻采集的多道信号。

% 参数 dt = 0.002; % 时间采样间隔，单位秒 dx = 5; % 道间距，单位米 [nt, nx] = size(data); % 1. 二维傅里叶正变换 spec = fft2(data); % 2. 将零频移到中心，便于显示和设计滤波算子 spec_shift = fftshift(spec); % 3. 构建物理坐标轴 f_axis = (-nt/2 : nt/2-1) / (nt*dt); k_axis = (-nx/2 : nx/2-1) / (nx*dx); % 4. 显示f-k谱，注意转置为了符合imagesc显示习惯 figure; imagesc(k_axis, f_axis, abs(spec_shift)); xlabel('波数 k (1/m)'); ylabel('频率 f (Hz)'); axis xy; % 让y轴从小到大显示

这段代码是基础中的基础。显示之后，你能直观看到面波的能量团，通常是一个以原点为顶点、向高低频两侧扩展的扇形或V形区域。找到它的边界斜率，就是计算视速度范围的开始。

3.1 扇形滤波器设计：参数选择是关键

面波压制最常见的是扇形切除：把f-k平面中超过某个速度界限的能量置零。实现方式有两种思路。

第一种，直接在f-k域生成一个二维掩膜（mask）。比如我要保留视速度大于等于1500 m/s的成分，那么在f-k平面上，每个点对应的视速度v = f / k，我把|v| < 1500的点设为0，其余保留。这里注意k的正负都要考虑，因为波可以向左传也可以向右传，所以滤波器必须关于k=0对称（如果只是去除某个方向传播的波，则只保留一翼）。

% 用meshgrid生成二维网格 [K, F] = meshgrid(k_axis, f_axis); % 视速度计算，注意避免除以零 V = zeros(size(K)); idx = abs(K) > 1e-10; V(idx) = F(idx) ./ K(idx); V(~idx) = sign(F(~idx)) * 1e10; % 近似无穷大速度 % 设计低切滤波器：只保留速度绝对值大于等于vmin的信号 vmin = 1500; % 米/秒 mask = ones(size(K)); mask(abs(V) < vmin) = 0; % 保留近零波数区间（即垂直入射波附近），这在很多情况是必要的 k_zero_band = abs(K) < 0.01 / dx; % 归一化阈值，根据数据调整 mask(k_zero_band) = 1; % 加一点平滑过渡，避免吉布斯效应（见第5节） mask = imgaussfilt(mask, 1.5); % 5. 应用滤波 spec_filtered = spec_shift .* mask; % 6. 反变换 data_filtered = real(ifft2(ifftshift(spec_filtered)));

这里有几个坑要单独说明。

第一，直接硬截断会产生严重的吉布斯效应，表现为时间域信号上出现振荡拖尾、空间域出现条带状伪影。所以上面用imgaussfilt对mask做了一次轻度高斯平滑，让过渡带不要过于陡峭。平滑的程度要控制：太缓会把有效波也削掉一点，太陡又振铃，这个需要根据实际数据摸索。

第二，速度下限vmin的选取不能拍脑袋，要根据面波最大视速度和有效波最小视速度来定。如果目标反射波速度本来就低（例如浅层低速层里的折射波），vmin太高会伤害有效信号。我习惯先把f-k谱显示出来，用探针工具量出干扰波的视速度边界，再去填vmin，不要凭感觉。

第三，k_zero_band的处理。很多人在低频段保留全部信号，原因是浅层多次波、直达波通常都在低频小波数区，切太狠会影响后续反演。这个带保留的宽度（我这里用的是0.01/dx）属于经验值，如果你只关心深部中高频反射，可以把带设得更窄。

3.2 非扇形陷波：处理已知速度的规则干扰

有时候干扰波不是扇形铺满，而是事先知道它是某个固定视速度的声波（比如空气中传播的声波在地震记录上表现为直线同相轴），这时我倾向于做窄带陷波而不是整个扇形切除。做法是构造一个带阻掩膜，把落在v≈v_noise附近的窄带区域置零，同时保留其余区域。

具体操作：计算v_fk = F ./ K的矩阵（注意K=0的置大值），然后设定一个速度带宽v_tol，比如噪声速度是340 m/s，带宽设为±30 m/s，把|V - 340| < 30的点置零。同样建议加高斯过渡。这种做法比扇形切除更精细，副作用也更小。

但要注意，如果干扰波速度随时间变化（比如地形起伏导致视速度变化），这种固定速度陷波就不太行了，应回到扇形切除路线或用自适应方法。

3.3 空间抽稀问题：不充分空间采样该怎么办

无论如何，有一点必须在matlab代码里主动检查：空间方向的Nyquist波数。空间采样率dx决定了最大不混叠波数：

k_Nyquist = 1 / (2 * dx)

如果某频率f下存在波数绝对值大于k_Nyquist的成分，那这部分就是空间假频。假频在f-k谱上表现很具迷惑性：它不会出现在正确的速度射线位置，而是由于采样不够，被“折叠”到低波数区域，也就是看起来像是低速度能量。这与时间域的频率混叠是同一个道理。如果你用滤波直接切掉假频区，会丢失真实信息；如果不切，它又会污染有效信号。最彻底的办法是空间重采样或道距加密采集，但数据已经到手了，所以现实做法是：如果假频不严重，且目标频段不在假频区，可以考虑直接用带通滤波把假频频段的能量一并压制；如果假频严重，那常规FK滤波解决不了，要改用f-k反假频插值（比如基于抗泄露的迭代插值），但那是另一个话题，这里不展开。

4. 一个从数据构建到滤波的完整案例

为了让你直观看到整个过程，我构造一个合成的二维地震记录，包含一条水平反射同相轴和一组线性面波干扰。这样你可以自己跑一遍，看看f-k谱长什么样，滤波前后变化多明显。

% 参数 dt = 0.002; dx = 5; nt = 512; nx = 64; t = (0:nt-1)*dt; x = (0:nx-1)*dx; % 构建水平反射：时间方向子波 + 空间方向水平展布 w = ricker(30, dt, 25); % 自己实现的Ricker子波，峰值频率25Hz ref = zeros(nt, nx); t_ref = 0.25; % 反射波到达时刻 i_ref = round(t_ref/dt); ref(i_ref:i_ref+length(w)-1, :) = repmat(w(:), 1, nx); % 构建线性面波：斜率对应速度500 m/s slope = 1/500; % 每米延迟秒数 noise = zeros(nt, nx); for ix = 1:nx t_noise = 0.05 + x(ix) * slope; % 道间延迟 i_noise = round(t_noise/dt) + 1; if i_noise + length(w) - 1 <= nt noise(i_noise:i_noise+length(w)-1, ix) = 0.8 * w(:); end end % 混合 data = ref + noise;

这里我用了自定义的ricker子波，你可以写一个简单函数：

function w = ricker(fdom, dt, T) t = -T/2 : dt : T/2; w = (1 - 2*pi^2*fdom^2*t.^2) .* exp(-pi^2*fdom^2*t.^2); w = w / max(w); end

接下来就是显式地进行FK变换、滤波、反变换。把这套流程跑下来，你会在f-k谱上清楚地看到：反射波能量贴f轴且一般在较低频率；面波能量沿600 m/s速度射线分布，是一个V形或窄带区域。设计mask时保留|V|>800 m/s的成分，面波被有效压制，反射波基本保留。

这个案例虽然合成数据很干净，但跑通它以后，换成真实Segy数据无非就是多一步读数据、多一步道编辑。逻辑和代码大部分可以复用。

5. 常见问题与排查技巧实录

FK变换本身不复杂，但实际用起来几乎每个环节都有坑。这里把我在matlab里踩过的典型问题整理成一张速查表，再单独挑几个细说。

细说三个高频问题。

第一个就是吉布斯效应。很多新手觉得掩膜是逻辑判断，直接mask=zeros/ones然后一顿切就完事了。这样做在f-k谱上确实能看到干扰波被切掉了，但时间域里全是振铃。我建议无论如何都要加一点过渡带，哪怕只有几个像元的平滑。如果用imgaussfilt（比如σ取1.5～2个网格），对于道间距5米的勘探数据，实测下来既保住了主频能量，拖尾也明显减少。如果你想要更精确控制，也可以自己做斜坡：在掩膜边界附近定义一个过渡带宽Δv，用线性或余弦过渡把0和1连接起来。

第二个容易忽略的是f-k谱中的直流分量。matlab的fft2会把零频和零波数处的绝对值放到（1，1）位置，如果你在显示时没处理好，可能看到一片巨大的亮斑。滤波时如果对零频附近做硬置零，反变换出来信号整体均值会偏移，造成每道出现常数偏差。我的习惯是，如果不需要研究低频背景，保持零频附近一个很小的矩形区域不变，而不是置零。

第三个高频问题：真实数据往往不是纯二维均匀采样，比如弯线采集、变观测系统，它们的道间距不规则。直接塞进fft2会引入人为的采样不均匀噪声。我处理这种数据的做法是，先用matlab自带的scatteredInterpolant插值到均匀网格上，再做FK变换。但要注意插值会改变波场的频率成分，最好是先做带限插值，插值后做FK变换只用于分析，不在这个域直接做滤波，而是把滤波设计信息拿回时间域应用。

再说一个经验性技巧：f-k谱的线性干扰边界经常不是一条干净直线，因为振幅随频率有衰减，谱上能量看起来是“一头粗一头细”的形状。这时不要只用一个vmin去切，我常用两段式滤波：先做f-k扇形初滤，把强能量去掉，然后对剩余数据做一次自适应噪声压制（比如在时间域做预测反卷积），效果往往比单靠FK一步到位好。原因在于FK滤波是全局算子，无法很好处理时空变化的噪声背景，而组合策略可以弥补这个盲区。

另一个实战提示：在做FK滤波前，不要忘记先对数据做常规预处理，比如去均值、去线性趋势、带通滤波。否则低频漂移会严重影响f-k谱的显示比例，真实反射波能量在谱图里会被“淹没”，你根本看不清有效的速度区间。这个过程我每次都会做，算是固定前置流程。

6. 参数敏感性分析与个人调试心得

回到文章开头的那句话，FK变换的关键在于f-k谱上能量团的位置和边界。所以调试流程我一般这样走：

1. 先显示原始数据的f-k谱，把坐标轴物理单位标好，标注干扰波边界。
2. 选一个保守的vmin，先做一次滤波，对比时间域的残差（原数据减滤波数据），看残差里是否还有明显相干能量。
3. 逐步调整vmin和过渡带宽度，直到残差里大部分是非相干的随机噪声，而有效波的同相轴连续不破碎。
4. 如果同相轴出现“搓板状”振幅抖动，多半是滤波器过渡带过宽或者vmin太接近有效波速度下限。这时要回看f-k谱，重新确定两者的分离点。

关于vmin的具体取值，经验上我习惯取有效波最小视速度的0.7～0.9倍。为什么不是直接取有效波速度？因为f-k谱里速度是一个斜线，能量分布有一定宽度，如果直接卡在有效波边界附近，滤波器会切掉一部分有效波的低频或高频边缘。留出余量，配合过渡带的平滑，能大大减小损伤。这个系数看起来不起眼，但值得专门记录。

另外，我强烈建议把参数值固定为脚本开头的变量，不要散落在代码里乱改。因为滤波参数往往要做敏感性测试场景（比如vmin=1200/1500/2000），用变量名定义好，改起来安全很多。我早期经常在mask代码里直接写数字，后来一个数字填错，整个结果就废了，排查半天才发现是滤波器边界变量被覆盖了。

再补充一个道上平均法（trace averaging）和f-k域做互补的小技巧：如果面波速度不太低，可以先在时间域做相邻道相减（相当于空间高通滤波），把低速的水平相关噪声压掉，剩余数据再做f-k域低切，这种组合能大幅提升高速度弱信号的可见度。类似思路可以衍生出很多变体，核心原则是：在处理流程中把线性算子和非线性算子（比如中值滤波）搭配起来，往往能比单一域处理更稳健。

还有就是如果数据量很大，比如上千道、上万时间采样点，fft2本身并不慢，但显示f-k谱的大矩阵会占很多内存。我一般先decimate或者分频段处理，把数据降到可交互的规模，设计好滤波参数后，再对原始数据整体执行一次滤波。注意分频段处理时要保证各频段滤波器过渡带连续，否则整体反变换后可能出现频段接缝处的不自然起伏。

最后再分享一次真实数据调试的体会。当时的数据里，面波速度大概在400～600 m/s，反射波有效速度在2000 m/s以上，两者在f-k谱上离得很开，按理说非常好切。但实际滤波后反射波同相轴出现了明显的横向振幅波动，检查发现是道间存在振幅不一致，导致f-k滤波后横向道间差异被放大。解决方法是滤波前先做一道振幅均衡（AGC或道间能量归一化），处理后再恢复原始振幅关系。这个坑不在FK变换本身，却在工程链路里非常现实。

本文还有配套的精品资源，点击获取