图像复原

Illustrated by GAloDos（Pixiv），Source

图像增强以“更适合观察或处理”为目标，并不要求解释图像为什么变差；图像复原则先建立退化模型，再利用关于成像系统和噪声的知识估计原图。若原图为 $𝑓$，退化算子为 $\mathrm{\mathscr{H}︀}$，噪声为 $𝜂$，观测图像统一写成

复原算法接收 $𝑔$ 以及对 $\mathrm{\mathscr{H}︀}$、$𝜂$ 的估计，输出原图估计 $\widehat{𝑓}$。同一幅模糊图可能由运动、离焦或大气扰动造成；退化原因不同，逆运算也不同。复原本质是一个反问题：即使已知 $𝑔$ 和对 $\mathrm{\mathscr{H}︀}$、$𝜂$ 的估计，通常也不存在唯一且稳定的原图 $𝑓$，因为退化会抹平细节，而噪声又混入观测。因此复原不能从“选哪个锐化滤波器”开始，而应从噪声模型、点扩散函数与先验约束开始。后文的逆滤波、维纳滤波和约束最小二乘，正是对“如何约束这个反问题”给出的三种不同回答。[1]

图 1 退化过程从 $𝑓$ 经系统 $𝐻$ 并叠加噪声 $𝜂$ 得到 $𝑔$；复原算子 $𝑅$ 利用模型与观测生成 $\widehat{𝑓}$。

图中的 $𝑓\left(𝑥,𝑦\right)$ 是未知的清晰原图，$\mathrm{\mathscr{H}︀}$ 是造成模糊或几何失真的退化系统，$𝜂\left(𝑥,𝑦\right)$ 是成像过程中叠加的噪声，$𝑔\left(𝑥,𝑦\right)$ 是传感器实际记录的退化图像。复原阶段并不能直接取得 $𝑓$，只能依据对 $\mathrm{\mathscr{H}︀}$ 与 $𝜂$ 的认识设计复原算子 $\mathrm{ℛ︀}$，从 $𝑔$ 得到原图估计 $\widehat{𝑓}\left(𝑥,𝑦\right)$。因此，$\widehat{𝑓}$ 通常只是满足模型和先验约束的估计，并不必然与 $𝑓$ 完全相同。

一、噪声模型

噪声模型描述随机灰度扰动可能取什么值、各值出现的概率，以及它是否依赖位置或原图。概率密度不是噪声图像本身：同一分布可产生许多不同样本，参数则决定均值、方差、偏斜程度和尾部。选定某种模型是假设而非事实——它以成像机制为依据，但必须通过平坦区样本或残差直方图检验；若实际噪声是混合型或位置相关，单一密度模型只是近似。

1. 高斯、瑞利、伽马噪声

$𝜇=0$ 时噪声不系统改变平均亮度，但会增加局部方差；$𝜎$ 越大，灰度扰动分布越宽。高斯模型允许任意实数值，所以加入有限范围图像后不应立刻裁剪，否则两端样本被堆到 0 与 255，实际分布不再是原高斯分布。

瑞利分布只在一侧取值，适合描述由两个正交高斯分量合成的幅值型随机量。它与高斯噪声的对称钟形不同，不能仅用相同均值与方差就认为两者可互换。

当形状参数改变时，伽马分布可从强烈右偏逐渐接近钟形。不同资料也会使用“尺度参数” $𝜃=\frac{1}{𝑎}$，读取公式和代码时必须确认第二参数到底是率还是尺度。

图 2 对称、单侧偏斜和单调衰减的密度具有不同统计特征；仅观察一幅噪声样本可能难以直接分辨，需要结合直方图与成像机制。

2. 指数、均匀、椒盐噪声

椒盐噪声不是在原灰度上加一个小扰动，而是把少量像素直接替换成极端值，因此均值滤波往往会把极端值扩散到邻域；中值和自适应中值更符合这种离群点模型。只有盐噪声或只有椒噪声时，调和、逆调和均值与最大/最小值滤波可以有针对性地处理。

3. 周期噪声

周期噪声可能来自电气干扰、扫描机构振动或重复采样误差。它在空间域表现为规则条纹、波纹或网纹，在频域却集中到少数位置，所以频域陷波通常比大范围空间平滑更有选择性：空间低通会同时模糊真实边缘，陷波只压制已识别的噪声频率。

4. 噪声参数的估计

有原始无噪声图像时，可以直接用差值 $𝜂=𝑔-𝑓$ 估计噪声；实际更常见的是只有退化图像。此时应选取若干内容近似恒定的小区域，因为平坦区的灰度变化主要来自噪声，而不是物体纹理。

图 3 先选取缺少真实边缘与纹理的局部块，再从样本直方图、均值和方差估计噪声参数。

对包含 $𝑛$ 个样本 ${𝑧}_{𝑖}$ 的区域，样本均值和无偏方差为

若局部真实亮度未知，可以先减去局部均值；若怀疑含有脉冲离群值，可用中位数和中位绝对偏差

得到稳健尺度估计，高斯情形常用 $\widehat{𝜎}\approx 1.4826\text{MAD}$。

椒盐概率可由直方图在 0、$𝐿-1$ 的异常计数估计，但真实图像本身也可能含纯黑或纯白区域，所以应结合空间孤立性与多个平坦块。周期噪声则在去除平滑背景后检查频谱峰值，其频率由峰相对中心的偏移估计，振幅可由峰值复系数估计。

patch = gray[40:90, 120:180].astype(np.float64)
mean = patch.mean()
variance = patch.var(ddof=1)

median = np.median(patch)
mad = np.median(np.abs(patch - median))
robust_sigma = 1.4826 * mad

二、使用空间滤波复原图像

空间复原假设噪声可由局部邻域统计削弱。窗口越大，随机波动平均得越充分，但边缘和细节损失也越严重。选择滤波器时应先匹配噪声模型，再决定窗口大小。

1. 均值滤波器

算术平均适合零均值随机噪声，却会平滑边缘。几何平均对大值的影响弱于算术平均，通常保留细节稍好；实现时应对数求和，避免大量乘法溢出。调和平均能抑制盐噪声，但窗口内出现 0 时分母发散，不适合椒噪声。逆调和平均在 $𝑄>0$ 时消除椒噪声，在 $𝑄<0$ 时消除盐噪声；选错符号会强化相反类型的脉冲。

图 4 (a) 无噪声参考图；(b) 高斯噪声图像；(c) 算术平均结果；(d) 几何平均结果；(e) 盐噪声图像；(f) 调和平均结果；(g) $𝑄=-1.5$ 的逆调和平均结果。

def geometric_mean_filter(image: np.ndarray, size: int) -> np.ndarray:
    padded = np.pad(image, size // 2, mode="reflect")
    windows = np.lib.stride_tricks.sliding_window_view(
        padded,
        (size, size),
    )
    safe = np.maximum(windows, 1e-6)
    return np.exp(np.mean(np.log(safe), axis=(-2, -1)))

def contraharmonic_mean_filter(
    image: np.ndarray,
    size: int,
    order: float,
) -> np.ndarray:
    padded = np.pad(image, size // 2, mode="reflect")
    windows = np.lib.stride_tricks.sliding_window_view(
        padded,
        (size, size),
    )
    safe = np.maximum(windows, 1e-6)
    numerator = np.sum(safe ** (order + 1.0), axis=(-2, -1))
    denominator = np.sum(safe**order, axis=(-2, -1))
    return numerator / np.maximum(denominator, 1e-12)

2. 统计排序滤波器

中值滤波取排序后的中间值，适合双极椒盐噪声。最大值滤波扩张亮区域、消除孤立椒点；最小值滤波扩张暗区域、消除孤立盐点。中点滤波取窗口最大值与最小值的平均

更适合均匀或高斯型随机噪声，不适合大量脉冲，因为两个极端值恰好被保留。

修正阿尔法均值把 $𝑚𝑛$ 个样本排序，删除最低和最高各 $\frac{𝑑}{2}$ 个，再对剩余 $𝑚𝑛-𝑑$ 个求平均：

$𝑑=0$ 时退化为算术平均；$𝑑=𝑚𝑛-1$ 时接近中值。它在高斯噪声与少量脉冲同时存在时提供折中。

图 5 (a) 椒盐噪声图像；(b) 中值滤波结果；(c) 最小值滤波结果；(d) 最大值滤波结果；(e) 中点滤波结果；(f) 修正阿尔法均值滤波结果。

kernel = np.ones((3, 3), dtype=np.uint8)
median = cv2.medianBlur(noisy_u8, 3)
minimum = cv2.erode(noisy_u8, kernel)
maximum = cv2.dilate(noisy_u8, kernel)
midpoint = 0.5 * (
    minimum.astype(np.float32)
    + maximum.astype(np.float32)
)

3. 自适应滤波器

在平坦区，${𝜎}_{𝐿}^2$ 接近噪声方差，系数接近 1，输出趋向局部均值；在真实边缘处，局部方差远大于噪声方差，系数较小，输出更接近原像素。将比值限制到 1 可避免噪声方差估计略大于局部方差时产生反向增强。

自适应中值滤波逐步扩大奇数窗口。对窗口最小值 ${𝑧}_{\text{min}}$、最大值 ${𝑧}_{\text{max}}$、中值 ${𝑧}_{\text{med}}$ 和中心像素 ${𝑧}_{𝑥𝑦}$，算法分两层判断：

1. 若 ${𝑧}_{\text{min}}<{𝑧}_{\text{med}}<{𝑧}_{\text{max}}$，说明中值不是脉冲，进入下一步；否则扩大窗口，直到最大尺寸。
2. 若 ${𝑧}_{\text{min}}<{𝑧}_{𝑥𝑦}<{𝑧}_{\text{max}}$，中心不是脉冲，保留中心；否则以中值替换。

固定中值窗口在高密度脉冲下可能所有样本都被污染；自适应扩大窗口能找到未污染中值，又在低噪声区域保持较小窗口。

图 6 (a) 方差从左向右增大的高斯噪声；(b) 自适应局部降噪结果；(c) 高密度椒盐噪声；(d) 自适应中值滤波结果。

def adaptive_local_filter(
    image: np.ndarray,
    size: int,
    noise_variance: float,
) -> np.ndarray:
    source = image.astype(np.float32)
    mean = cv2.boxFilter(source, -1, (size, size))
    mean_square = cv2.boxFilter(source * source, -1, (size, size))
    variance = np.maximum(mean_square - mean * mean, 0.0)
    ratio = np.minimum(
        noise_variance / np.maximum(variance, 1e-6),
        1.0,
    )
    return source - ratio * (source - mean)


def adaptive_median_filter(
    image: np.ndarray,
    maximum_size: int = 7,
) -> np.ndarray:
    if maximum_size < 3 or maximum_size % 2 == 0:
        raise ValueError("maximum_size 必须是大于等于 3 的奇数")

    source = image.astype(np.uint8)
    output = source.copy()
    radius = maximum_size // 2
    padded = np.pad(source, radius, mode="reflect")

    for row in range(source.shape[0]):
        for col in range(source.shape[1]):
            center = source[row, col]

            for size in range(3, maximum_size + 1, 2):
                half = size // 2
                padded_row = row + radius
                padded_col = col + radius
                window = padded[
                    padded_row - half:padded_row + half + 1,
                    padded_col - half:padded_col + half + 1,
                ]
                minimum = int(window.min())
                maximum = int(window.max())
                median = int(np.median(window))

                # A 层：中值不是脉冲时，才检验中心像素。
                if minimum < median < maximum:
                    # B 层：中心不是脉冲则保留，否则换成中值。
                    output[row, col] = (
                        center
                        if minimum < center < maximum
                        else median
                    )
                    break

                # 到达最大窗口仍未通过 A 层，只能采用当前中值。
                if size == maximum_size:
                    output[row, col] = median

    return output

三、使用频域滤波复原图像

周期噪声的能量集中于少数离散频率，而真实图像的大部分频谱仍可保留，因此频域陷波比空间窗口平均更有选择性。空间滤波需要为了压制整幅条纹而使用较大窗口，通常会模糊全部边缘；陷波只在噪声峰附近衰减。

1. 陷波滤波

设频谱高为 $𝑃$、宽为 $𝑄$，则中心化后的零频位置为 $\left(\lfloor \frac{𝑄}{2}\rfloor ,\lfloor \frac{𝑃}{2}\rfloor \right)$。第 $𝑘$ 个噪声峰相对中心的偏移为 $\left({𝑢}_{𝑘},{𝑣}_{𝑘}\right)$，到正、负共轭峰的距离分别为

高斯陷波带阻的一般形式可写为

每个因子在对应峰中心为 0，远离峰时趋近 1；正负峰成对出现，保证对实图像的共轭对称性。陷波带通为 ${𝐻}_{\text{NP}}=1-{𝐻}_{\text{NR}}$，可先提取周期干扰，再从退化图像中减去。

下面的实际图像叠加两项周期噪声：一项产生偏移 $\left(0,\pm 48\right)$ 的峰，另一项产生 $\left(\pm 31,\pm 18\right)$ 的共轭峰。处理链完整展示了从空间诊断到复原的过程。

图 7 (a) 无噪声参考图；(b) 含两组周期噪声的退化图；(c) 中心化频谱及成对噪声峰；(d) 成对高斯陷波传递函数；(e) 陷波后的频谱；(f) 逆变换得到的复原图像。

def gaussian_notch_reject(
    shape: tuple[int, int],
    offsets: list[tuple[float, float]],
    sigma: float,
) -> np.ndarray:
    rows, cols = shape
    y, x = np.ogrid[:rows, :cols]
    center_y = rows // 2
    center_x = cols // 2
    transfer = np.ones(shape, dtype=np.float64)

    for offset_y, offset_x in offsets:
        for sign in (-1.0, 1.0):
            notch_y = center_y + sign * offset_y
            notch_x = center_x + sign * offset_x
            distance2 = (y - notch_y) ** 2 + (x - notch_x) ** 2
            transfer *= 1.0 - np.exp(
                -distance2 / (2.0 * sigma**2)
            )
    return transfer

spectrum = np.fft.fftshift(np.fft.fft2(noisy))
notch = gaussian_notch_reject(
    noisy.shape,
    offsets=[(0.0, 48.0), (31.0, 18.0)],
    sigma=2.6,
)
filtered = spectrum * notch
restored = np.fft.ifft2(np.fft.ifftshift(filtered)).real

2. 最优陷波滤波

普通陷波把指定频率固定衰减，却没有利用噪声在空间中的强弱变化。最优陷波先用陷波带通提取干扰估计 $\widehat{𝜂}$，再在每个局部区域选择权重 $𝑤\left(𝑥,𝑦\right)$：

令局部期望用窗口平均近似，使复原结果局部方差最小，可得

分子是观测与噪声估计的局部协方差，分母是噪声估计的局部方差。噪声在某区域较强时权重增大，较弱时减小；分母过小时应令权重为 0 或加正则项，避免数值爆炸。

def optimal_notch_subtraction(
    degraded: np.ndarray,
    noise_estimate: np.ndarray,
    window: int = 15,
) -> np.ndarray:
    g = degraded.astype(np.float32)
    n = noise_estimate.astype(np.float32)
    size = (window, window)

    mean_g = cv2.blur(g, size)
    mean_n = cv2.blur(n, size)
    covariance = cv2.blur(g * n, size) - mean_g * mean_n
    variance_n = cv2.blur(n * n, size) - mean_n * mean_n
    weight = covariance / np.maximum(variance_n, 1e-6)
    weight = np.clip(weight, 0.0, 1.5)
    return g - weight * n

四、线性位置不变退化与退化函数的估计

1. 线性位置不变退化

$ℎ$ 是一个理想点目标经过系统后形成的图样，所以称点扩散函数；$𝐻$ 是光学传递函数。位置不变意味着同一个点目标放到图像不同位置时，扩散形状只发生平移。实际镜头边缘像差、景深变化和非刚性运动可能破坏这一假设，此时单一全局 $𝐻$ 只能近似。

图 8 点扩散函数把每个输入点扩展并叠加，噪声随后进入观测；频域中卷积转为 $𝐻𝐹$。

2. 退化函数估计

a. 观察法

观察法在退化图像中寻找结构简单、原形状可推断的区域。例如孤立亮点的模糊形状近似 PSF，笔直边缘的边缘扩展函数求导后可估计线扩散函数。选区频谱 ${𝐺}_{𝑠}$ 与对理想选区的估计 ${𝐹}_{𝑠}$ 可给出

分母接近 0 的频率不能可靠相除，选区也必须足够小，使退化近似位置不变。

b. 试验法

若能控制成像系统，可输入已知点源、细线、棋盘格或正弦图样。点源响应直接测量 PSF；不同空间频率正弦图样的幅度衰减和相位移动可逐点测量 $𝐻$。试验法通常比仅观察未知图像可靠，但要求设备状态、焦距、曝光和运动条件与待复原图像一致。

c. 建模法

建模法从物理运动或光学过程推导 $𝐻$。以曝光时间 $𝑇$ 内的平面匀速运动为例，若相机与场景的相对位移为 $\left({𝑥}_0\left(𝑡\right),{𝑦}_0\left(𝑡\right)\right)$，时间平均观测为

利用傅里叶平移性质，传递函数为

若总位移为 $\left(𝑎,𝑏\right)$，即 ${𝑥}_0\left(𝑡\right)=𝑎\frac{𝑡}{𝑇}$、${𝑦}_0\left(𝑡\right)=𝑏\frac{𝑡}{𝑇}$，令 $𝑠=𝑢𝑎+𝑣𝑏$，则

其中 $\text{sinc}\left(𝑠\right)=\frac{sin\left(𝜋𝑠\right)}{𝜋𝑠}$，$𝑠=0$ 时极限为 1。幅度沿运动方向的频率坐标呈 sinc 形，并在非零整数 $𝑠$ 处出现零点；相位因子对应模糊核中心位置。

图 9 匀速直线路径产生线状点扩散函数，其频率响应沿相应方向出现周期零点；这些零点是逆滤波不稳定的根源。

下面给出一个实际数值例子。令总运动长度为 $23$ 像素、方向角为 $18°$，则水平与竖直位移近似为

在曝光期间沿这条线段均匀采样并归一化，便得到线状 PSF；将其与清晰图像卷积，得到具有明显斜向拖影的运动模糊图像。再叠加标准差为 $3$ 个灰度级的零均值高斯噪声，得到后续复原算法共同使用的观测图像。这样，物理假设、模型参数、PSF 与实际退化结果便一一对应，而不只是停留在传递函数公式上。

图 10 (a) 清晰参考图像；(b) 使用长度 $23$ 像素、方向角 $18°$ 的匀速运动 PSF 得到的无噪声模糊图像；(c) 在 (b) 上叠加标准差为 $3$ 的高斯噪声后得到的实际退化观测。

def uniform_motion_transfer(
    shape: tuple[int, int],
    displacement_x: float,
    displacement_y: float,
) -> np.ndarray:
    rows, cols = shape
    v = np.fft.fftfreq(rows)[:, None]
    u = np.fft.fftfreq(cols)[None, :]
    projection = u * displacement_x + v * displacement_y
    return np.sinc(projection) * np.exp(-1j * np.pi * projection)

motion_h = uniform_motion_transfer(
    gray.shape,
    displacement_x=22.0,
    displacement_y=7.0,
)

非匀速或曲线运动时，仍可对一般积分数值求和：在曝光区间采样轨迹点，计算每个点的相位因子并取平均。若轨迹未知，可以从陀螺仪、相邻视频帧、图像中直线拖影的长度与方向，或盲去卷积联合估计原图和 PSF。

五、逆滤波

即使准确知道 $𝐻$，也不能保证准确复原。由

直接相除得到

当 $𝐻$ 很小，噪声项被 $\frac{1}{𝐻}$ 巨大放大；当 $𝐻=0$，该频率的原图信息已完全丢失，不存在有限逆；此外，PSF 参数的微小误差也会在小 $𝐻$ 附近被放大。因此“知道退化函数”并不等于逆问题已经稳定。

截断逆滤波只在 $\left|𝐻\right|\ge 𝜀$ 的位置相除，其余频率置零或不作逆补偿：

$𝜀$ 越小，恢复的频率更多，噪声放大也更严重；越大则稳定但模糊残留更多。

def truncated_inverse(
    degraded: np.ndarray,
    transfer: np.ndarray,
    threshold: float,
) -> np.ndarray:
    spectrum = np.fft.fft2(degraded)
    inverse = np.zeros_like(transfer, dtype=np.complex128)
    stable = np.abs(transfer) >= threshold
    inverse[stable] = 1.0 / transfer[stable]
    return np.fft.ifft2(inverse * spectrum).real

六、维纳滤波

维纳滤波不要求把退化完全反转，而是在噪声放大与去模糊之间最小化均方误差。以下推导把原图与加性噪声视为零均值、二阶平稳且互不相关的随机场，并假设退化系统线性位置不变。设使用线性频域估计 $\widehat{𝐹}=𝐴𝐺$，目标为

代入 $𝐺=𝐻𝐹+𝑁$，展开期望并对 ${𝐴}^{\ast }$ 求极值，可得正规方程

其中 ${𝑆}_{𝑓}=𝐸\left[{\left|𝐹\right|}^2\right]$、${𝑆}_{𝜂}=𝐸\left[{\left|𝑁\right|}^2\right]$ 是信号和噪声功率谱。因此

在高信噪比频率，$\frac{{𝑆}_{𝜂}}{{𝑆}_{𝑓}}$ 很小，维纳滤波接近逆滤波；在低信噪比或 $𝐻$ 很小的位置，分母中的噪声项限制增益。未知功率谱时常用常数 $𝐾\approx \frac{{𝑆}_{𝜂}}{{𝑆}_{𝑓}}$，但它只是白噪声与平均信号功率的粗略近似。

def wiener_restore(
    degraded: np.ndarray,
    transfer: np.ndarray,
    noise_to_signal: float | np.ndarray,
) -> np.ndarray:
    spectrum = np.fft.fft2(degraded)
    numerator = np.conjugate(transfer)
    denominator = np.abs(transfer) ** 2 + noise_to_signal
    estimate = numerator / np.maximum(denominator, 1e-12) * spectrum
    return np.fft.ifft2(estimate).real

图 11 (a) 清晰参考图像；(b) 无噪声运动模糊；(c) 叠加高斯噪声后的观测；(d) 截断逆滤波结果；(e) 维纳滤波结果；(f) 约束最小二乘结果。

七、约束最小二乘滤波

维纳滤波需要信号与噪声功率谱。约束最小二乘改用平滑先验：在观测残差符合噪声能量的条件下，让估计图像经过高通算子 $𝑝$ 后的能量尽量小：

用拉格朗日乘子把约束并入目标，并在频域对 ${\widehat{𝐹}}^{\ast }$ 求导，得到

所以

$𝑃$ 常取离散拉普拉斯核的 DFT，使高频振荡受到更强惩罚。$𝛾$ 小时接近逆滤波，细节多但噪声敏感；$𝛾$ 大时更平滑、稳定，却可能过度模糊。$𝜀$ 若来自噪声方差估计，可用迭代调整 $𝛾$，使最终残差能量与预期噪声能量相符。

def constrained_least_squares(
    degraded: np.ndarray,
    transfer: np.ndarray,
    penalty: np.ndarray,
    gamma: float,
) -> np.ndarray:
    spectrum = np.fft.fft2(degraded)
    denominator = (
        np.abs(transfer) ** 2
        + gamma * np.abs(penalty) ** 2
    )
    estimate = (
        np.conjugate(transfer)
        / np.maximum(denominator, 1e-12)
        * spectrum
    )
    return np.fft.ifft2(estimate).real

八、几何均值滤波

几何均值复原滤波器在逆滤波与维纳滤波之间作乘法插值：

$𝛼=1$ 时得到逆滤波，$𝛼=0$、$𝛽=1$ 时得到维纳滤波；$𝛽$ 进一步调节噪声项强度。由于两个因子具有相同的 ${𝐻}^{\ast }$ 相位，更稳定的实现可合并为

def geometric_mean_restore(
    degraded: np.ndarray,
    transfer: np.ndarray,
    noise_to_signal: float | np.ndarray,
    alpha: float,
    beta: float = 1.0,
) -> np.ndarray:
    power = np.abs(transfer) ** 2
    denominator = (
        np.maximum(power, 1e-12) ** alpha
        * np.maximum(
            power + beta * noise_to_signal,
            1e-12,
        ) ** (1.0 - alpha)
    )
    restoration = np.conjugate(transfer) / denominator
    estimate = restoration * np.fft.fft2(degraded)
    return np.fft.ifft2(estimate).real

图 12 (a) 同时含运动模糊与高斯噪声的退化图像；(b) 取 $𝛼=0.45$ 的几何均值滤波结果。