频率-波数域F-K偏移成像单点定位仿真程序

资 源 简 介

本项目基于MATLAB环境完整实现了频率-波数（F-K）域偏移成像算法，专门用于单点目标的精确定位与成像仿真。程序首先包含一个正演模拟模块，利用波动方程或走时公式生成单点散射体在时空域（Time-Space Domain）的原始回波信号，该信号在B-scan图像中呈现典型的双曲线绕射特征。核心处理部分采用F-K偏移（亦称Stolt偏移）技术，具体步骤包括：对原始数据矩阵执行二维快速傅里叶变换（2D-FFT）将其转换至频率-波数域；根据标量波动方程的色散关系，计算垂直波数与频率及水平波数的非线性映射关系；利用插值算法（如最近邻或双线性插值）将数据从时间频率轴重新映射到垂直波数轴，以校正波场传播的几何畸变；最后执行二维傅里叶逆变换（2D-IFFT）将数据变换回空间-深度域。该项目能够演示能量从双曲线轨迹聚焦回真实目标空间位置的完整过程，有效验证F-K算法在提高横向分辨率和还原目标真实位置方面的能力，适用于地震勘探、探地雷达（GPR）或合成孔径雷达（SAR）的基础教学与算法验证。

详 情 说 明

FK_Frequency_Wavenumber_Migration_Single_Target

项目介绍

本项目基于 MATLAB 环境完整实现了频率-波数（Frequency-Wavenumber, F-K）域偏移成像算法，该算法也被称为 Stolt 偏移。项目旨在演示和验证如何利用波动方程的色散关系，将时空域（Time-Space Domain）中呈现双曲线特征的绕射波场，通过频域映射归位到真实的几何空间位置。

该程序专门针对单点目标（Singular Point Target）进行设计，包含从正演数据生成到最终成像可视化的全流程，适用于探地雷达（GPR）、地震勘探或合成孔径雷达（SAR）领域的算法教学与基础研究。

功能特性

• 正演模拟系统：内置基于运动学原理的回波生成器，利用双曲线走时方程和 Ricker 子波合成单点目标的 B-scan 原始数据。
• F-K Stolt 偏移核心算法：完整实现了基于 FFT 的偏移流程，包括频域变换、色散关系映射、插值校正及雅可比能量补偿。
• 高分辨率成像：通过波场延拓和重聚焦，显著提高了目标的横向分辨率，能够将双曲线能量收敛至真实的物理坐标。
• 多维数据可视化：提供原始时空域数据、F-K 域频谱分布以及最终深度-空间域成像结果的直观对比展示。

算法原理与实现细节

本项目中的 main.m 文件严格按照 F-K 偏移的数学原理编写，具体实现逻辑如下：

1. 仿真参数配置

程序首先建立物理模型和计算网格，设定介质波速为 1000 m/s，中心频率为 100 MHz。定义了横向（x）和深度（z）方向的网格大小，并根据奈奎斯特采样定理设定时间步长 dt。同时，精确定义了单点散射体在空间中的坐标 (target_x, target_z)。

2. 正演模拟（B-scan 数据生成）

在进行偏移之前，程序首先生成原始回波数据（raw_data）：

• Ricker 子波生成：构造一个中心频率为 100 MHz 的 Ricker 子波作为发射信号。
• 双曲线走时计算：对于每一个地面接收位置（道），根据几何关系计算从发射点到目标点再返回接收点的双程走时。
• 信号叠加：采用“脉冲响应叠加”的思路，将子波按照计算出的走时延迟累加到时间-空间矩阵中。这模拟了点目标在 B-scan 图像中典型的双曲线绕射特征。

3. F-K (Stolt) 偏移处理

这是本项目的核心模块，包含以下关键步骤：

• FFT 预处理：对原始数据矩阵的时间维和空间维进行补零操作（至 2 的幂次），以优化 FFT 计算效率。
• 2D-FFT 变换：利用二维快速傅里叶变换将数据从时空域转换至频率-波数域（$f-k_x$），并对频谱进行中心化平移。
• 网格构建：构建原始的频率轴（$omega$）和水平波数轴（$k_x$），以及用于重采样的目标垂直波数轴（$k_z$）。
• Stolt 映射关系计算：

根据标量波动方程的色散关系，计算从输出域 $(k_x, k_z)$ 到输入域频率 $omega$ 的逆映射关系。代码中明确实现了以下公式： $W_{mapped} = frac{v}{2} cdot text{sign}(k_z) cdot sqrt{k_x^2 + k_z^2}$ 这一步决定了如何重新排列频域能量以校正几何畸变。

• 插值重采样 (Interpolation)：

使用 interp2 函数，根据计算出的映射频率 $W_{mapped}$，将频谱数据从原始的 $(k_x, omega)$ 网格插值到新的 $(k_x, k_z)$ 网格上。

• 雅可比缩放 (Jacobian Scaling)：

为了在坐标变换过程中保持能量守恒，程序计算并应用了雅可比行列式比例因子。代码实现了如下缩放因子，以补偿从 $omega$ 域到 $k_z$ 域变换带来的振幅变化： Scale Factor $= frac{v^2 cdot |k_z|}{4 cdot |W_{mapped}|}$

• 2D-IFFT 变换：将插值并缩放后的频谱通过二维傅里叶逆变换还原回空间域，取实部得到最终的深度剖面图。

4. 结果可视化

程序最后通过一个包含三个子图的窗口展示处理结果：

• 左图：原始未偏移数据（Time-Space），展示目标的双曲线绕射形态。
• 中图：F-K 域频谱幅度（dB 刻度），展示信号在频率-波数域的能量分布。
• 右图：F-K 偏移后成像（Depth-Space），展示能量聚焦后的点目标，并用白色十字标记了真实的设定位置以验证精度。

关键代码分析

• fft2 / ifft2：执行二维傅里叶变换和逆变换，是连接时空域与频率波数域的桥梁。
• fftshift / ifftshift：将零频分量移至频谱中心，便于直观分析波数与频率的关系及后续的对称网格构建。
• meshgrid：用于生成二维网格矩阵，支撑后续的插值和并行矩阵运算。
• interp2：实现 Stolt 偏移中最关键的插值步骤，程序选用线性插值（linear）并将超出范围的消失波置零。
• 矩阵向量化操作：程序中大量的数学运算（如距离计算、映射关系推导）均采用 MATLAB 的矩阵化操作，避免了低效的循环，保证了计算速度。

系统要求

• MATLAB R2016b 或更高版本。
• 无需额外的工具箱（Toolbox），主要依赖 MATLAB 基础数值计算功能。

使用方法

1. 确保 MATLAB 的当前工作目录包含 main.m 文件。
2. 在 MATLAB 命令行窗口输入 main 并回车，或直接点击编辑器中的“运行”按钮。
3. 程序将自动执行仿真与处理，并弹出一个图形窗口展示处理前后的对比结果。