从原始与平滑图像恢复高斯模糊Sigma：卷积核计算异常排查

问题背景 #

我需要恢复高斯模糊的sigma值，手头有原始图像（vp_true）和平滑后的图像（vp_smooth），尝试用反卷积方法提取高斯核，但得到的结果完全不具备高斯形态。使用的代码如下：

import numpy as np


nz = 141
nx = 681
h = 25

x = [i*h / 1000 for i in range(nx)]
z = [i*h / 1000 for i in range(nz)]

vp_true = np.fromfile("vp_true", np.float32).astype(np.float64).reshape(nz, nx, order='F')
vp_smooth = np.fromfile("vp_smooth", np.float32).astype(np.float64).reshape(nz, nx, order='F')

fft_vp_true = np.fft.fft2(vp_true)
fft_vp_smooth = np.fft.fft2(vp_smooth)

fft_g = np.fft.fftshift(fft_vp_smooth) / np.fft.fftshift(fft_vp_true)
g = np.real(np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(fft_g))))

fig = plt.figure()
ax = plt.gca()

im = ax.pcolorfast(x, z, g, cmap="jet")
ax.invert_yaxis()
ax.set_aspect('equal')

plt.tight_layout()
plt.savefig("test.pdf")
plt.show()

错误原因分析 #

1. 频域除法的数值稳定性缺失
   原始图像的频域中存在大量接近0的分量，直接做除法会导致噪声被极端放大，完全扭曲卷积核的形态，这是朴素反卷积的典型问题。

2. 未处理核的归一化与移位
   高斯核是总和为1的对称核，直接反卷积得到的结果既没有归一化，也可能因FFT循环卷积的特性出现中心偏移，导致形态异常。

3. 不必要的频域移位操作
   虽然fftshift本身逻辑没问题，但在频域相除前后的移位操作没有实际意义，反而增加了不必要的复杂度，且未解决核心的数值不稳定问题。

修正方案 #

步骤1：加入正则化的反卷积代码 #

使用Tikhonov正则化抑制噪声，同时修正核的移位与归一化：

import numpy as np
import matplotlib.pyplot as plt

nz = 141
nx = 681
h = 25

x = [i*h / 1000 for i in range(nx)]
z = [i*h / 1000 for i in range(nz)]

# 读取原始与平滑图像
vp_true = np.fromfile("vp_true", np.float32).astype(np.float64).reshape(nz, nx, order='F')
vp_smooth = np.fromfile("vp_smooth", np.float32).astype(np.float64).reshape(nz, nx, order='F')

# 计算频域变换
fft_true = np.fft.fft2(vp_true)
fft_smooth = np.fft.fft2(vp_smooth)

# Tikhonov正则化：避免除零与噪声放大，lambda_reg需根据数据调整
lambda_reg = 1e-3
fft_g = fft_smooth / (fft_true + lambda_reg * np.max(np.abs(fft_true)))

# 逆傅里叶变换得到实值核
g = np.real(np.fft.ifft2(fft_g))

# 将核的中心移到图像中心（修正循环卷积的移位）
g = np.fft.fftshift(g)

# 归一化核，确保总和为1（符合高斯核的属性）
g = g / np.sum(g)

# 可视化结果
fig = plt.figure()
ax = plt.gca()
im = ax.pcolorfast(x, z, g, cmap="jet")
ax.invert_yaxis()
ax.set_aspect('equal')
plt.colorbar(im)
plt.tight_layout()
plt.savefig("test.pdf")
plt.show()

步骤2：直接拟合高斯核获取sigma值 #

如果已知核是高斯形态，直接拟合比反卷积更稳定：

from scipy.optimize import curve_fit

# 提取核的中心区域（根据实际核大小调整窗口，示例用31x31）
kernel_size = 31
center_z, center_x = nz // 2, nx // 2
g_crop = g[center_z - kernel_size//2 : center_z + kernel_size//2 + 1,
           center_x - kernel_size//2 : center_x + kernel_size//2 + 1]

# 定义二维高斯函数模型
def gaussian_2d(xy, sigma_z, sigma_x, amplitude, offset):
    x, y = xy
    return offset + amplitude * np.exp(-(x**2/(2*sigma_x**2) + y**2/(2*sigma_z**2)))

# 创建拟合用的网格坐标（以核中心为原点）
x_grid, z_grid = np.meshgrid(np.arange(kernel_size), np.arange(kernel_size))
x_grid -= kernel_size // 2
z_grid -= kernel_size // 2
xy_data = np.vstack((x_grid.ravel(), z_grid.ravel()))

# 拟合高斯参数
popt, _ = curve_fit(gaussian_2d, xy_data, g_crop.ravel(), 
                    p0=[5, 5, np.max(g_crop), np.min(g_crop)])

# 转换为物理单位（h=25单位是毫米，除以1000转为千米）
sigma_z_pix, sigma_x_pix = popt[0], popt[1]
sigma_z = sigma_z_pix * h / 1000
sigma_x = sigma_x_pix * h / 1000

print(f"恢复的z方向sigma值：{sigma_z:.4f} km")
print(f"恢复的x方向sigma值：{sigma_x:.4f} km")

内容的提问来源于stack exchange，提问作者Alexandre Hoffmann