4. 频率域向上延拓算法实现

频率域向上延拓的核心思想，是把空间域中“位场满足拉普拉斯方程”这一物理约束，转化为波数域中一个极其简洁的乘法算子。空间域向上延拓本是一个积分运算，计算量大且难以直接离散；而进入频率域后，延拓就退化为“对每个波数分量乘以一个衰减因子”。理解这一点，是写出正确程序的前提。本节围绕四个环节展开：二维离散傅里叶变换的落地、延拓算子的构造、离散化误差的来源与控制、以及可运行的编程实现步骤。

一、从空间域到波数域：变换与算子

设观测面为 z=0，位场为 f(x,y)，向上延拓高度为 h（h>0，向上为正）。对 f 做二维离散傅里叶变换得到频谱 F(kx,ky)，则延拓后频谱为：

F_up(kx, ky) = F(kx, ky) · H(kx, ky)
H(kx, ky) = exp( -h · sqrt(kx² + ky²) ) = exp(-h·k)

其中 k=sqrt(kx²+ky²) 为径向波数。该算子是一个纯实数、恒小于等于 1 的低通衰减器：波数越大（对应越短波长、越浅源信息），衰减越剧烈；k=0 的直流分量衰减因子为 1，保持不变。这正是向上延拓“压制浅部高频、突出深部低频”的数学本质。逆变换回空间域即得延拓结果。

需要特别强调符号约定：若采用向下延拓，则指数取正号 exp(+h·k)，算子大于 1，会放大高频噪声，属于不适定问题。本节只讨论向上延拓，全文与图示均采用“向上为正、算子取负指数”的约定。

二、离散化误差的来源与控制

公式简单，但直接照抄往往得到错误结果，问题几乎都出在离散化细节上。主要有四类误差源，需要在编程时逐一处理。

第一，波数坐标的排列顺序。FFT 输出的频谱零频位于数组左上角（索引 0,0），而正负波数分布在四个角。若直接用索引计算 k，会把高频误当低频。标准做法是先构造与频谱同形的波数网格，再用 fftshift 把零频移到中心，或直接在未移位的网格上按“前半为正、后半为负”的规则生成 kx、ky。以 N 点采样、采样间隔 dx 为例，角波数增量为 dk=2π/(N·dx)，第 i 个索引对应的波数为：i≤N/2 时取 i·dk，i>N/2 时取 (i−N)·dk。

第二，采样与延拓高度的匹配。延拓算子对高频的衰减是指数级的，当 h 远大于网格间距时，高频分量被压到数值噪声以下，此时结果主要由低频决定，这是物理上合理的；但若网格太稀（dx 过大），Nyquist 波数偏低，短波长信息在采样阶段就已丢失，延拓后无法恢复。经验上应保证 dx 明显小于目标探测深度对应的波长尺度，且延拓高度 h 不宜超过测区尺度的数倍，否则边缘效应会显著增强。

第三，边界效应与周期延拓假设。FFT 隐含把数据视为周期信号，实测数据在边界处不连续，会在延拓后产生沿边界的条带状假异常。常用控制手段是：先对数据做扩边（余弦 taper、镜像对称或最小曲率外推），扩边宽度一般取延拓高度的 2～3 倍，延拓后再裁回原测区。这是实际生产中最容易被忽略、却对成果质量影响最大的一步。

第四，数值精度与直流分量。建议全程使用双精度浮点；对含明显区域背景场的数据，可先去掉均值（直流分量），延拓后再加回，避免大数相减带来的相对误差。此外，若数据存在坏点或空洞，应先插值填补，否则单个异常值经 FFT 会污染整个频谱。

三、编程实现步骤

下面给出以 Python/NumPy 实现的完整流程，参数为教学示例（网格 256×256、dx=dy=50 m、延拓高度 h=200 m），实际使用时按数据替换即可。

import numpy as np

def up_continue(grid, dx, dy, h, pad_ratio=2.0):
    # 1) 去均值，减小直流分量误差
    mean_val = grid.mean()
    g = grid - mean_val

    # 2) 余弦扩边，抑制周期假设引起的边界假异常
    pad = int(max(h / dx, h / dy) * pad_ratio)
    g = np.pad(g, pad, mode='reflect')
    ny, nx = g.shape

    # 3) 构造波数网格（未移位，按 FFT 索引规则）
    kx = 2*np.pi*np.fft.fftfreq(nx, d=dx)
    ky = 2*np.pi*np.fft.fftfreq(ny, d=dy)
    KX, KY = np.meshgrid(kx, ky)
    K = np.sqrt(KX**2 + KY**2)

    # 4) 延拓算子（向上为正，取负指数）
    H = np.exp(-h * K)

    # 5) 正变换 → 乘算子 → 逆变换
    F = np.fft.fft2(g)
    F_up = F * H
    g_up = np.real(np.fft.ifft2(F_up))

    # 6) 裁回原测区并恢复均值
    g_up = g_up[pad:pad+grid.shape[0], pad:pad+grid.shape[1]]
    return g_up + mean_val

实现中有三个易错点值得单独说明。其一，np.fft.fftfreq 返回的正是按 FFT 索引排列的波数，直接用它构造网格即可，无需再手动 fftshift，这样可避免符号混乱。其二，逆变换后必须取实部，理论上结果应为实数，但浮点运算会留下极小虚部。其三，扩边与裁剪的索引必须严格对应，否则结果会整体错位。

运行后如何判断结果合理？可从三方面观察：延拓后异常幅值整体减小、等值线更平滑、局部小尺度跳跃被抹平；原图中的浅部孤立高值（如浅埋磁性体引起的尖峰）明显减弱，而区域性宽缓异常形态基本保留；边界处不出现明显条带。若幅值反而增大或出现高频振荡，应首先检查算子指数符号与波数网格是否正确。

把上述流程封装成函数后，改变 h 即可快速获得不同高度的延拓结果，为后续边界增强与场源分离提供稳定的输入。需要牢记的是：算法本身只有一行乘法，真正决定成果质量的是波数网格的正确性、边界处理是否充分以及参数 h 与网格尺度的匹配关系。