从“锯齿”到“光滑”：相位解包裹 (Phase Unwrapping) 算法深度解析

在光学仿真（如超表面设计）或光学测量（如干涉仪、结构光 3D 扫描）中，我们经常遇到一个恼人的现象：

当我们计算波前的相位时，得到的不是一个连续光滑的曲面，而是一张布满锯齿状跳变的“条纹图”。

这是因为数学上的反三角函数（如 atan2）只能输出 $(-\pi, \pi]$ 之间的主值。这种被强行截断的相位，被称为包裹相位 (Wrapped Phase)。而将它恢复为真实的、连续的绝对相位 (Unwrapped Phase) 的过程，就是相位解包裹。

这不仅是图像处理问题，更是恢复物理真实的关键步骤。

---

01. 数学定义：迷失的 $2\pi$

假设真实的物理相位是 $\phi(x, y)$，它是一个连续变化的函数。
我们的探测器或算法计算出的包裹相位 $\psi(x, y)$ 与真实相位的关系为：

$$\psi(x, y) = \mathcal{W}[\phi(x, y)] = \phi(x, y) + 2\pi k(x, y)$$

其中：

• $\mathcal{W}[\cdot]$ 是包裹算子（取模运算）。
• $k(x, y)$ 是一个未知的整数（阶次）。

相位解包裹的核心任务：就是对每一个像素点 $(x, y)$，找到正确的整数 $k(x, y)$，使得恢复后的 $\phi(x, y)$ 满足空间连续性。

---

02. 一维解包裹：Itoh 条件

在 1D 情况下，解包裹非常简单，主要依赖 Itoh 条件。

核心假设：如果采样率足够高，相邻两个采样点的真实相位差绝对值小于 $\pi$。
$$|\Delta \phi| < \pi$$

算法步骤：

1. 计算相邻像素的差分：$\Delta \psi_i = \psi_i - \psi_{i-1}$
2. 如果 $\Delta \psi_i > \pi$，说明发生了下跳变，真实相位应该减去 $2\pi$。
3. 如果 $\Delta \psi_i < -\pi$，说明发生了上跳变，真实相位应该加上 $2\pi$。
4. 对修正后的差分进行累加（积分），即可恢复 $\phi$。

代码实现：

import numpy as np

def unwrap_1d(wrapped_phase):
    diff = np.diff(wrapped_phase)
    # 修正跳变：将差分值限制在 (-pi, pi] 之间
    diff_corrected = np.arctan2(np.sin(diff), np.cos(diff))
    # 累加恢复
    unwrapped = np.concatenate(([wrapped_phase[0]], wrapped_phase[0] + np.cumsum(diff_corrected)))
    return unwrapped

---

03. 二维解包裹：噩梦的开始

到了 2D 图像，问题变得复杂得多。我们不能简单地逐行进行 1D 解包裹，因为噪声和断点会导致误差传播，形成贯穿整张图的“拉丝”错误。

这就引入了 2D 解包裹的核心概念：残点 (Residues)。

什么是残点？

在一个 $2 \times 2$ 的像素小闭环中，计算包裹相位的梯度积分。理论上，在一个闭合路径上积分一圈，结果应该是 0。
但如果该区域存在噪声或相位断裂，积分结果可能是 $\pm 2\pi$。

• $+2\pi$：正残点（源）。
• $-2\pi$：负残点（汇）。

残点的存在意味着该点的解包裹结果与积分路径有关 (Path Dependent)。 如果你的积分路径穿过了一对残点之间，你的解包裹结果就会出错。

---

04. 主流算法流派

为了处理残点，业界发展出了两大流派：

1. 路径跟踪类 (Path Following) —— 代表：枝切法 (Goldstein Branch Cut)

思路：既然穿过残点连线会出错，那我就把这些残点连起来，设为“禁区”（枝切线，Branch Cut），积分路径绕着走，绝对不穿过它。

• 优点：如果不穿过切线，解是精确的。
• 缺点：如果残点太多（噪声大），切线会连成网，把图像分割成孤岛，导致无法解包裹。

2. 最小范数类 (Minimum Norm) —— 代表：最小二乘法 (DCT/FFT based Least Squares)

思路：不纠结于局部路径，而是寻找一个全局最优解 $\phi$，使得其梯度 $\nabla \phi$ 与包裹相位梯度 $\nabla \psi$ 的差的平方和最小。
这就转化为求解泊松方程 (Poisson Equation)：
$$\nabla^2 \phi = \rho$$
其中 $\rho$ 是由包裹相位计算出的拉普拉斯算子。利用离散余弦变换 (DCT) 或 FFT 可以快速求解。

• 优点：速度极快，对高斯噪声鲁棒，得到的表面非常光滑。
• 缺点：会抹平真实的相位陡峭变化，低估斜率。

---

05. Python 实战：质量导向路径跟踪算法

在实际工程（如 scikit-image）中，最常用的是一种折中方案：质量导向 (Quality Guided)。
它优先处理“质量好”（梯度平缓、信噪比高）的像素，最后处理“质量差”的像素，从而阻断误差传播。

我们可以直接使用 scikit-image 库：

import numpy as np
import matplotlib.pyplot as plt
from skimage.restoration import unwrap_phase

# 1. 构造一个模拟的 2D 真实相位 (高斯山峰)
N = 256
x = np.linspace(-3, 3, N)
y = np.linspace(-3, 3, N)
X, Y = np.meshgrid(x, y)
true_phase = 5 * np.exp(-(X**2 + Y**2)) * 5  # 高度为 25 的山峰

# 2. 生成包裹相位 (Wrapped Phase)
# 加上一点噪声
noise = np.random.normal(0, 0.5, (N, N))
wrapped_phase = np.angle(np.exp(1j * (true_phase + noise)))

# 3. 执行解包裹
unwrapped_phase = unwrap_phase(wrapped_phase)

# 4. 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

axes[0].imshow(true_phase, cmap='jet')
axes[0].set_title('True Phase (Ground Truth)')

axes[1].imshow(wrapped_phase, cmap='hsv') # 使用循环色图展示条纹
axes[1].set_title('Wrapped Phase (-pi to pi)')

axes[2].imshow(unwrapped_phase, cmap='jet')
axes[2].set_title('Unwrapped Result')

plt.tight_layout()
plt.show()

代码解读

• np.angle(exp(1j * phi)) 是生成包裹相位的标准数学方法。
• unwrap_phase 函数内部通常使用了基于网络流 (Network Flow) 或质量图导向的算法，比单纯的 Itoh 积分要健壮得多。

---

06. 总结

算法类型	核心思想	适用场景	Python 库
Itoh (1D)	简单积分	简单的无噪声一维信号	numpy.unwrap
枝切法 (Branch Cut)	避开残点连线	噪声较少、对精度要求极高的干涉图	自写或科研专用库
最小二乘 (Least Squares)	求解泊松方程	噪声较大、追求全局平滑、速度要求高	需基于 DCT/FFT 实现
质量导向 (Quality Guided)	好的区域先走	通用性最强，图像处理首选	skimage.restoration

掌握相位解包裹，是你从“光学仿真新手”迈向“计算成像专家”的重要一步。下次看到黑白条纹图，别忘了它背后藏着一个连续的光滑世界。

---