通过解调使用正则化相位跟踪技术进行相位解包裹

1. 绪论

光学计量学通常使用光学干涉仪来测量各种物理量。1,2 根据应用的不同，可以使用多种类型的干涉仪，但它们的共同目标是产生一个由被测物理量调制的条纹图案。使用这种光束编码程序可以检测到的物理量范围非常广：深度测量、应变分析、温度梯度、表面变形等等。

干涉图可以用以下数学表达式表示：

$I(x,y)=a(x,y)+b(x,y)cos[phi(x,y)]+n(x,y)(1)$

其中 $a(x,y)$ 是缓慢变化的背景照明，$b(x,y)$ 是幅度调制，它也是一个低频信号，$\varphi (x,y)$ 是与被测物理量相关的相位项。项 $n(x,y)$ 是一个加性噪声，尽管它也可能是乘性的，就像电子散斑干涉测量法 (ESPI) 的情况一样。然后将连续干涉图在 CCD 摄像机上成像，并使用视频帧抓取器数字化，以便在数字计算机中进行进一步分析。计算机辅助条纹分析的目的是自动检测干涉图上发生的二维相位变化 $\varphi (x,y)$，这是由于被测物理量的空间变化引起的。

有许多技术可以用来测量有趣的空间相位变化 $phi(x,y)$。其中，我们可以提到相位步进干涉测量法 (PSI)，3,4 它至少需要三个相位步进干涉图。干涉图之间的相位差必须在整个干涉图中已知。在这种情况下，可以估计每个可分辨图像像素的调制相位。当干涉仪的大气湍流和机械条件在获得三个相位步进干涉图所需的时间内保持恒定，PSI 是首选技术。当上述要求无法满足时，如果将载波条纹引入条纹图案以获得载波频率干涉图，则可以仅分析一个干涉图。然后可以使用诸如傅里叶变换、5 同步、6,7 空间相位偏移、8-11 相位锁定环路 (PLL)、12,13 等等这些众所周知的技术来分析该干涉图。

除了 PLL 技术外，所有其他方法（傅里叶变换、同步、PSI 和空间相位偏移）都会给出检测到的相位包裹（真实相位的模 $2\pi$），因为相位估计过程中涉及反正切函数。包裹相位和非包裹相位之间的关系可以表示为：

$\varphi (x,y)={\varphi }_W(x,y)+2\pi k(x,y),(2)$

其中 $\varphi (x,y)$ 是非包裹相位，${\varphi }_W(x,y)$ 是包裹相位，$k(x,y)$ 是一个整数值校正场。对于从高质量条纹数据计算出的相位图，解包裹问题很简单；在这些相位图中，水平和垂直方向上相邻相位样本之间的绝对相位差小于 π，除了预期的 $2\pi$ 不连续点。因此，解包裹是一个简单的问题，即在遇到的每个不连续点上添加或减去 $2\pi$ 偏移量，14,15 或者对包裹相位差进行积分。16,17

当相邻像素之间的绝对相位差在反正切函数的不连续点以外的点大于 π 时，解包裹变得更加困难。这些不连续点可能是由高频、高幅度噪声、相位跳跃不连续和条纹图案中的区域欠采样引入的。Ghiglia 等人 18 考虑通过在开始解包裹过程之前隔离这些错误的不连续点来解包裹相位。当围绕大小为 L 的正方形路径的包裹相位差之和不等于零时，就会检测到错误的不连续点或相位不一致。不一致会产生相位误差（意外相位跳跃），这些误差会沿着解包裹方向（在顺序解包裹中）传播。因此，解包裹过程变得依赖于路径，即，根据所选的解包裹方向，可能会获得不同的解包裹相位场。

Ghiglia 等人 17 通过应用 Fried19 和 Hudgin20 的相位梯度最小二乘积分 21-23 的思想，为获得鲁棒的路径无关相位解包裹器迈出了重要的一步。参考文献 17 中所需的相位梯度是作为包裹估计相位的 x 和 y 方向上的包裹相位差获得的。然后对这个包裹的梯度场进行最小二乘积分，以获得所寻找的连续相位。最近，Marroquin 等人 24,25 通过在所寻找的解中添加一个正则化项，以势能范数的形式来扩展包裹相位梯度最小二乘积分技术。当使用这种技术时，可以滤除解包裹相位中的一些噪声，以及对无效相位数据（如孔洞）进行插值，这些数据在内部具有明确的行为。

最小二乘积分 17 或其正则化扩展 25 的一个缺点是假设相邻像素之间的相位差小于 π 的绝对值。也就是说，这些技术将包裹相位的包裹差视为真实梯度场。不幸的是，当解包裹严重噪声的相位图时，情况并非如此。这里获得的相位梯度实际上是在高噪声和高相位梯度区域中包裹的。这是在 ESPI 中遇到的非常嘈杂的情况，其中可能存在如此大量的噪声，以至于相邻像素之间的相位差可能大于 π 弧度，而不是预期的 $2\pi$ 跳跃。

这里提出的解包裹系统可以看作是先前报道的 PLL 系统的推广。12,13 正则化相位跟踪 (RPT) 系统在估计相位中具有明显更高的信噪比，相对于傅里叶变换、同步或 PSI 技术而言。原因是 RPT 系统的行为类似于一个极窄带通滤波器，跟踪宽带条纹图案信号。迄今为止，RPT 系统已在多个不同应用中得到报道，26-28 范围从单个闭合条纹干涉图的解调到 PSI。在所有这些应用中，RPT 技术已被证明比以前众所周知的干涉测量技术 3-18 对噪声更鲁棒，并且对条纹图案边界不敏感。

所提出的 RPT 方法的一个优点是它在干涉图边缘没有相位失真。当干涉图中存在无效条纹数据区域（如阴影和孔洞）时，这一点特别有用。当使用 ESPI 分析机械零件中的应变时，这是一个常见的情况。但是，这里提出的 RPT 技术仅限于解调平滑波前。也就是说，我们没有考虑分段连续函数的解调。

然而，Servin 等人 26 已经将类似的 RPT 技术应用于直接解包裹相位图（无需生成两个中间条纹图案）；这里介绍的 RPT 更易于理解，并且对噪声更鲁棒。

2. 使用正则化相位跟踪技术的解调解包裹

如上所述，使用这里提出的技术解包裹给定相位图的第一步是将包裹相位转换为两个相位偏移的条纹图案。我们可以通过使用被解包裹的相位图的余弦和正弦来获得这些条纹图案，即：

$\begin{array}{rl}{I}_C(x,y)& =cos[{\varphi }_W(x,y)],\\ I_S(x,y)& =sin[{\varphi }_W(x,y)],\end{array}$ (3)

其中 ${\varphi }_W(x,y)$ 是被解包裹的相位图。

现在，解包裹问题可以看作是使用 RPT 技术 28 对两个相位偏移的条纹图案进行解调。现在，我们概述 RPT 技术的基本原理，该技术用于解调-解包裹这些相位偏移的条纹图案。

鉴于所提出的技术使用正则化相位跟踪器，我们简要解释对相位估计问题进行正则化的动机。为了在条纹分析中正则化相位检测问题，有必要找到一个合适的能量或成本函数，该函数至少使用两个有助于约束估计相位场的项。这些项与以下内容相关：

(1) 估计函数与观测值之间的保真度。
(2) 调制相位场的平滑度。

然后，我们假设所寻找的相位函数是使成本函数最小化的函数。

特别是在这里提出的解包裹 RPT 技术中，我们假设局部条纹图案 [公式 (3)] 可以被认为是空间单色的，因此它们的辐照度可以被建模为一个在相位空间中由平面调制的余弦函数。该余弦函数的幅度必须接近观测到的辐照度 [陈述 (1)]。这样的相位平面必须适应条纹图案中的每个区域，因为它的局部频率在二维空间中连续变化。所提出的成本函数的第二项指的是预期解包裹检测相位的平滑度和连续性 [陈述 (2)]。具体来说，所提出的成本函数，由解包裹相位 ${\varphi }_0(x,y)$ 在每个位置 $(x,y)$ 处最小化，是：

其中 L 是一个具有有效条纹数据（良好的幅度调制）的二维格点，$N_{x,y}$ 是围绕坐标 $(x,y)$ 的邻域区域，其中正在解包裹相位，$m(x,y)$ 是一个指示场，如果位置 $(x,y)$ 已经被解包裹，则等于 1，否则等于 0。函数 ${\omega }_x(x,y)$ 和 ${\omega }_y(x,y)$ 是沿 x 和 y 方向估计的局部频率。最后，λ 是一个正则化参数，它与 $N_{x,y}$ 的大小一起控制检测到的解包裹相位的平滑度。

公式 (4) 中的前两项试图使局部条纹模型在最小二乘意义上尽可能接近邻域 $N_{x,y}$ 中的观测到的辐照度 [陈述 (1)]。第三项在仅使用由 $m(x,y)$ 标记的先前检测到的像素时，强制执行平滑度和连续性的假设 [陈述 (2)]。换句话说，第三项强化了这样一个事实，即解包裹的相位值 ${\varphi }_0(x,y)$ 必须在最小二乘意义上尽可能接近邻域 $N_{x,y}$ 中的适应相位平面 $$\begin{gathered} p(x,y,ϵ, \\ eta) \end{gathered}$$。还要注意，局部相位平面同时适应观测数据（通过其余弦和正弦模型）以及由 $m(x,y)$ 标记的已经估计的相位值 ${\varphi }_0(x,y)$。这样做是为了找到与观测到的条纹辐照度相兼容的最平滑相位。

为了解调-解包裹由公式 (3) 给出的条纹图案，我们必须找到成本函数 $U(x,y)$ 关于场 ${\varphi }_0(x,y)$、${\omega }_x(x,y)$ 和 ${\omega }_y(x,y)$ 的最小值。为此，我们提出了一种解调算法，该算法通过以下顺序方式优化每个 L 中的 $U(x,y)$，然后获得我们所谓的第一个全局解包裹相位估计。

对 L 上的第一个全局解包裹相位估计的执行方式如下：首先，将指示函数 $m(x,y)$ 设置为零 $[m(x,y)=0在L中]$。然后选择一个种子或起点 $(x_0,y_0)$ 在 L 内部，开始对条纹图案进行解调。然后，函数 $U(x_0,y_0)$ 关于 ${\varphi }_0(x_0,y_0)$、${\omega }_x(x_0,y_0)$、${\omega }_y(x_0,y_0)$ 进行优化。已访问的位置被标记为已解包裹；也就是说，我们将 $m(x_0,y_0)$ 设置为 1。一旦种子像素被解包裹，顺序相位解包裹将按以下步骤进行：

1. 选择 L 内部的 $(x,y)$ 像素（随机或以预定的扫描顺序）。
2. 如果 $m(x,y)=1$，则返回第一个语句。如果 $m(x,y)=0$，则测试以确定 $m(x^{\prime },y^{\prime })=1$ 是否适用于任何相邻像素 $(x^{\prime },y^{\prime })$。如果没有任何相邻像素已被估计，则返回第一个语句。如果 $m(x^{\prime },y^{\prime })=1$ 适用于相邻像素，则将 $[{\varphi }_0(x^{\prime },y^{\prime })$、${\omega }_x(x^{\prime },y^{\prime })$、${\omega }_y(x^{\prime },y^{\prime })]$ 作为初始条件来最小化 $U(x,y)$ 关于 $[{\varphi }_0(x,y)$、${\omega }_x(x,y)$、${\omega }_y(x,y)]$。
3. 设置 $m(x,y)=1$。
4. 返回第一个语句，直到 L 中的所有像素都被估计出来。

直观地理解第一个迭代的另一种方式是将其视为一个晶体生长 (CG) 过程，其中新的分子以最小化局部晶体能量的特定方向添加到主体中，该方向由相邻和先前定位的分子决定。在图 1 中，我们可以看到一个关于先前定义的二维区域的空间分布的特定示例，在这些区域中，当我们选择特定的顺序扫描时，条纹图案正在被解调。

为了优化位置 $(x,y)$ 处的 $U(x,y)$ 关于 $({\varphi }_0$、${\omega }_x$、${\omega }_y)$，我们使用了一个简单的梯度下降：

其中 $\tau =0.05$ 是步长，k 是迭代次数。通常只需要公式 (6) 中的两次迭代（除了对起始种子点的解调，这可能需要大约 10 次迭代）。这是因为公式 (6) 中的初始条件是从一个邻域解包裹像素中获取的，因此初始条件已经非常接近公式 (4) 的稳定点。在实践中，公式 (6) 中第一个关系中的 τ 参数可以乘以大约 10 来加速梯度搜索的收敛速度。

第一个全局相位估计 [当公式 (6) 与 CG 算法一起使用时] 通常非常接近实际的调制相位；如果需要，可以执行额外的全局迭代来改进解包裹过程。可以执行额外的迭代，再次使用公式 (6)，但现在将优化后的三元组 $({\varphi }_0$、${\omega }_x$、${\omega }_y)$ 作为同一位置 $(x,y)$ 的初始条件（而不是像在第一个全局 CG 迭代中那样在邻域位置）。请注意，对于额外的解包裹迭代，公式 (4) 中的指示函数 $m(x,y)$ 现在在 L 中等于 1。因此，只要 L 中的所有位置在每次全局迭代中都被访问，就可以以任何所需的顺序扫描格点。在实践中

仅需要三到四次额外的全局迭代才能在二维格点 L 中的每个位置 $(x,y)$ 处达到 $U(x,y)$ 的稳定最小值。

为了方便起见，我们可以评论一下成本函数 $U(x,y)$，它是在先前关于解包裹的研究中提出的，其中使用了 RPT 系统。26 我们使用了以下形式的成本函数：

其中 $V(x)=[x-2\pi int(x/\pi ){]}^2$，$int(y)$ 取 y 的整数部分。公式 (7) 中的所有其他变量和函数与为公式 (4) 定义的变量和函数相同。如参考文献 26 中所示，用于优化公式 (7) 的顺序技术与本文中使用的技术相同。公式 (7) 相对于公式 (4) 的两个主要缺点是，包装运算符 $V(\cdot )$ 不如公式 (4) 中的 sin(.) 和 cos(.) 函数那样容易相对于三元组 $({\varphi }_0$、${\omega }_x$、${\omega }_y)$ 推导出来，而且计算机模拟表明，从公式 (4) 推导出的优化系统（RPT 系统）比从公式 (7) 推导出的优化系统对噪声更鲁棒。关于从公式 (4) 推导出的 RPT 系统的噪声鲁棒性更高的一个假设可能是由于使用了两个数据场，$I_C(x,y)$ 和 $I_S(x,y)$。[鉴于公式 (6) 中给出的 RPT 动态系统的非线性，对这一点进行数学证明可能非常困难]。在公式 (6) 中给出的 RPT 动态系统中，我们保留了与 RPT 系统 26-28 通常相关的相同边缘鲁棒性。

3. 讨论

现在我们可以给出该 RPT 解包裹系统的频域解释。在 CG 算法中，由于从估计的三元组 $({\varphi }_0,{\omega }_x,{\omega }_y)$（来自一个相邻检测到的像素）获得梯度搜索系统的初始条件，使得 RPT 系统寻找与观测到的辐照度相兼容的最平滑调制相位。在频域中，RPT 解包裹技术可以被解释为一个窄带通自适应滤波器，跟踪相位偏移干涉图的瞬时相位和频率。被检测到的信号的频谱可能比跟踪 RPT 检测滤波器宽得多。该窄带通相位跟踪系统具有很大的惯性，因此它可以平滑地移动到条纹图案的二维频谱中。具有自适应窄带通滤波器跟踪更宽信号频谱的另一个优点是检测到的相位的信噪比有很大提高。换句话说，使用窄带通相位跟踪系统，在解调可能具有明显更宽频谱的信号时，该系统具有自适应窄带通滤波器的噪声抑制能力。解包裹相位中的信噪比增益是众所周知的，并且在现在广泛使用的电信中，相位跟踪系统通常被使用。

最后，正如我们在第 4 和 5 节中看到的那样，RPT 技术对边界边缘效应几乎不敏感，这在标准相位移技术中是一个严重的缺点，在这些技术中，在实践中必须对条纹图案使用卷积低通滤波器来滤除一些相位噪声。

如上所述，参数 λ 和邻域 $N_{x,y}$ 的大小与检测到的相位的带宽以及 RPT 算法的鲁棒性相关。例如，对于非常嘈杂的条纹图案，$N_{x,y}$ 的大小应该很大，以便 RPT 系统能够跟踪平滑的调制相位，并忽略快速变化的噪声。如果邻域 $N_{x,y}$ 的大小太小，系统可能会被相位噪声困住。换句话说，$N_{x,y}$ 的大小与要恢复的预期相位的尺度以及噪声量有关。另一方面，一旦 $N_{x,y}$ 的大小被选择，公式 (4) 中 λ 参数的值就不那么重要了。在整个过程中使用了 $\lambda =1$。这里给出了计算机模拟和实验结果。RPT 技术的计算速度与邻域 $N_{x,y}$ 的大小和格点 L 的大小成线性关系。在所呈现的实验中，$N_{x,y}$ 的大小范围从 9 × 9 到 17 × 17 像素。

就像在 CG 过程中一样，RPT 技术的第一次迭代非常关键。如果它成功，第一次 RPT 迭代将使整个解包裹系统移动到正确的吸引子，因此进一步的细化计算（计算额外的全局迭代）将使估计的相位仅略微移动。如果 RPT 系统到达错误的吸引子，RPT 系统将被困在该局部最小值中，并且进一步的迭代可能无法将系统拉出来。在这种情况下，可能需要尝试另一个邻域大小 $N_{x,y}$ 来达到所需的吸引子。

4. 计算机模拟

我们进行了一次计算机模拟，将这里提出的 RPT 解包裹器与众所周知的相位解包裹技术（即包裹相位差的最小二乘积分）进行比较。18 图 2 展示了一个无噪声的计算机生成的包裹相位 $phi(x,y)$。我们向该相位添加了范围在 (0.0, 1.55) 弧度内的白噪声和均匀分布的随机噪声。使用该噪声相位，我们生成了三个相位偏移的条纹图案，如下所示：

我们使用 128 × 128 个像素，每个像素有 256 个灰度级来计算这些条纹图案。对噪声条纹图案进行了 3 × 3 卷积平均器进行低通滤波。此外，对 3 × 3 平均器进行了三次卷积，并与相位偏移的

条纹图案 [公式 (8)] 进行卷积。进行低通滤波是为了从三个可用的条纹图案中获得最佳的相位估计，方法是使用最小二乘解包裹技术。得到的条纹图案如图 3 所示。然后，使用以下公式估计这些条纹图案的调制相位：

其中运算符 $LPF[\cdot ]$ 表示低通滤波过程。然后使用包裹相位差的最小二乘积分对该包裹相位进行解包裹，结果如图 4 所示。如上所述，对解包裹的相位进行了重新包裹，以便与图 2 中显示的真实调制相位进行比较。我们可以看到真实相位的动态范围显著地

下降了。这种下降是由于在严重噪声的情况下，包裹差运算符是真实相位梯度的糟糕估计器，因为高噪声会导致一些相位差出现在 $(-\pi ,\pi )$ 弧度范围之外，因此包裹运算符将给出错误的相位梯度。我们尝试了使用 3 × 3 平均器低通滤波器对公式 (8) 中的条纹图案进行不同次数的卷积，但发现估计的解包裹相位在所有情况下都更差。

另一方面，我们使用这里提出的 RPT 解包裹方法对由以下公式给出的噪声包裹相位进行了测试：

正如我们所看到的，这次我们没有使用低通滤波过程，因为我们的技术以这种方式效果更好。这里提出的 RPT 解包裹器得到的解包裹相位（再次重新包裹以用于比较目的）如图 5 所示。正如图 5 所示，这里提出的 RPT 解包裹器在很大程度上保留了真实相位的动态范围。对于第一个全局 CG 迭代，使用了 11 × 11 个像素的邻域 $N_{x,y}$。随后，需要三次全局迭代来改进估计的解包裹相位。如上所述，在应用 RPT 解包裹器之前，不需要对条纹图案进行低通滤波。这是一个优点，因为我们不必担心低通滤波器的最佳通带以平滑条纹图案。

第五章是关于实验结果的。

5. 实验结果

图 6 展示了一组三幅相位偏移的电子散斑干涉图。ESPI 图像采用减法 ESPI 模式从一块边缘夹紧的铝板上获得。该板用一个放置在它后部表面的烙铁头加热。这些干涉图之间的相位偏移是 $pi/2$ 弧度。我们使用以下三步公式来获得图 7 中所示的相位图：

${\varphi }_W(x,y)=arctan\left[\frac{I_1(x,y)-I_3(x,y)}{2I_2(x,y)-I_1(x,y)-I_3(x,y)}\right],(11)$

其中 ${\varphi }_W(x,y)$ 代表估计的包裹相位图。

图 8(a) 和 8(b) 展示了公式 (11) 中给出的包裹相位的正弦和余弦条纹图案。这些条纹图案是输入到公式 (6) 中所示的 RPT 解包裹系统中的，没有任何预处理。图 9 展示了解包裹的相位（重新包裹以与图 6 中的条纹图案进行比较），该相位是使用这里提出的 RPT 解包裹技术获得的。图 10 展示了检测到的解包裹相位的网格图。对于第一个全局 CG 迭代，使用了 5 × 5 个像素的邻域 $N_{x,y}$。随后，使用了四次全局迭代来改进估计的解包裹相位。获得图 9 中所示的解包裹相位所需的计算时间约为 1 分钟，使用的是 180 MHz 的奔腾-Pro 机器。

图 10. 图 9 中所示的解包裹相位的网格图。