1 引言

在当前许多高度、面形计量技术中,相位信息通常作为直接测量参数,经相位高度映射后,可得出被测物体的高度分布。由于相位值大多是通过反正切函数取得的,会被截断在(-π,π]之间,因此为了得到正确的结果,需在进行相位高度映射之前将截断相位(或称为包裹相位)展开成为连续分布的真实相位。二维相位展开就是从二维的截断相位中恢复物体真实相位的过程,因此,二维截断相位展开(或称为相位解包裹、相位解缠)是决定测量结果成败及精度的关键环节。

在无噪声等干扰信号的理想情况下,相位展开只需要沿着截断相位的行和列逐点进行展开即可。实现原理是:从某一点作为相位展开起始点沿着截断相位的行或列方向,比较相邻两个像素间的相位值,若相邻像素点之间的相位差小于-π,则后一点相位值加2π;若相邻像素点之间的相位差大于π,则后一点相位值减2π。逐行(列)展开后,再逐列(行)进行同样的展开操作。Itoh条件指出:当相邻像素间真实相位的相位差绝对值小于π时,真实相位的差分等于截断相位的差分再截断,而在实际测量过程中,由于噪声、欠采样等的影响,该条件并不能总是得到满足,最终将由于截断级次的判断错误而造成误差传播,展开为错误的结果^[1] 。

为了解决上述的问题,学者们已经提出了很多相位展开的算法。其中空间相位展开算法大致分为与路径无关的最小化目标函数算法和与路径有关的路径跟踪算法^[2]。与路径无关的相位展开典型算法包括最小二乘法^[3-5]、细胞机自动算法^[6-7]、最小不连续(MD)法^[8-9]等。大部分的空间相位展开算法均为与路径有关的算法,如枝切(BC)法^[10-12]、质量图导向法^[13-15]、掩膜切割法^[16]、快速相位展开(FPU)算法^[17-18]等。

本文对连续相位受噪声干扰的情形进行仿真模拟,分别利用空间相位展开算法中的枝切法、最小不连续法以及快速相位展开算法将所得截断相位进行展开。通过比较三种算法各自的性能,得出了最小不连续法的抗噪性能最强、但耗时也最长的结论。随后,提出将枝切法和快速相位展开算法分别与最小不连续法融合,既可以利用枝切法与快速相位展开算法展开速度快的优点,又可以利用最小不连续法展开准确性高的性能,仿真模拟与实际实验的相位展开结果证实了融合后的两种算法展开速度快且具有较高的准确度,能够有效运用到二维空间相位展开中。

2 三种基本相位展开算法原理

空间相位展开算法是基于Itoh条件提出的,理想情况下可根据该条件直接展开成为正确的连续相位。然而在实际测量中,由于噪声、欠采样以及相位信息本身不连续等原因的影响,获得的相位信息无法在所有点处都满足Itoh条件,将不满足Itoh条件的点称为残点,在残点处将会出现展开错误。由于空间相位展开实质上是一个积分过程,任一点的展开错误都将传播到后续展开过程中,因此若展开路径经过残点,就会引起误差传播,造成相位展开失败。为了在相位展开过程中减少残点的影响,路径跟踪算法通过选取合适的路径展开相位,避开不可靠区域,从而减小展开相位的误差,本研究所介绍的枝切法和快速相位展开算法即基于此思想。最小化目标函数算法通过建立目标函数,使目标函数获得最小值来展开相位。本研究中最小不连续法即基于此思想。

在相位展开过程中,由于残点的存在,当展开路径穿过残点时会导致展开错误,并在之后的展开路径中造成误差传播。枝切法通过识别残点,连接残点形成枝切线,利用枝切线引导相位展开路径从而避开不可靠区域。具体实现过程为:

1)在2×2闭环中,根据相位展开结果与路径无关的条件识别残点^[6],并将残点进一步分为正残点与负残点^[19],如图1(a)所示,图中A、B、D、I和K表示正残点位置,F、H和N表示负残点位置,O、E为图像边界。

2)通过平衡残点形成枝切线,使得所形成的枝切线上正残点和负残点数目相等,或者通过将残点与边界相连接形成枝切线。如图1(b)所示,图中枝切线EFGCABD是通过将残点与边界相连接形成的,枝切线HI与枝切线IJLKMN通过平衡正负残点数目所形成,两枝切线相交于I点形成一条枝切线。

3)在枝切线外找一个点作为相位展开起始点,将其四周非枝切线上的点进行相位展开,直到所有非枝切线上的点全部展开后,再展开相邻的枝切线上的点。

图 1. 枝切法中的(a)残点和(b)枝切线分布

Fig. 1. Distribution of (a) residues and (b) branch cuttings in BC algorithm

下载图片 查看所有图片

将相邻像素间的相位跳跃级次称为跳跃数,在理想情况下,连续相位全场跳跃数总和为0,但在实际测量中由于噪声、欠采样以及相位本身不连续等原因使得真实连续相位存在一定的跳跃数。根据最小范数法则^[20-21],当目标函数—全场跳跃数总和达到最小时,正是最小范数P=1的特殊情况。最小不连续法通过一系列迭代来完成截断相位的展开过程,最终找到一个满足全场跳跃数总和(即不连续性)达到最小的相位分布,作为与真实连续相位最为接近的展开相位。迭代过程中使用“结点”“边”“环”来辅助得到最终全场的跳跃数,结点为围绕一个像素的4个角点,边由相邻两个结点连接形成,环为由多条边围成的回路,增长环的正边值数目大于负边值。如图2所示,像素点用圆圈表示,结点用小方块表示,边表示为连接两个相邻结点的箭头。实现过程为:

1)根据截断相位分布,初始化水平跳跃数组、垂直跳跃数组以及结点数组;

2)由三个数组共同决定是否可以添加边,当添加的边生成一个增长环时,删除形成该环的边,并更新三个数组,由于图2中所形成的环为逆时针方向,因此环内所有像素点的相位值均增加2π,以此减少相位的不连续性;

3)重复过程2),直到没有新边生成;

4)利用最终的水平跳跃数组、垂直跳跃数组以及相位展开起始点计算得到展开相位。

图 2. 最小不连续法中的结点、边、环示意图

Fig. 2. Schematic of node, edge and loop in MD algorithm

下载图片 查看所有图片

相位值的二阶导数反映了对应点的凹凸度,可以较好地表征相位的不连续性。对于每一个像素点,分别在其8邻域的水平、垂直以及两个斜方向上计算该像素点的二阶导数,将计算得到的4个二阶导数求和,用该值表征像素点可靠度,图3(a)所示数值为各个像素点对应的可靠度,值越大,可靠度越小。在快速相位展开算法中,相邻像素点在水平方向和垂直方向形成边,边的可靠度定义为形成边的两个像素点可靠度之和,如图3(b)中相邻两个像素点之间的值代表各边的可靠度。该算法利用边的可靠度规划展开路径,可以限制低可靠度区域的误差传播。将以同一个像素点作为相位展开起始点展开的像素集合定义为组,图3(c)~(e)中相同颜色代表同一组。组的可靠度定义为组内所有像素可靠度之和,初始化时每个截断相位属于一个独立的组。实现过程为:

1)将每条边按照可靠度大小进行排序。

2)选取可靠度最大的边,判断构成该边的两个像素点是否属于同组。若两个像素点不属于同组,则认为该边有效。

3)比较两个像素点所在组的可靠度大小,根据所在组可靠度较大的像素点相位展开所在组可靠度较小的像素点,如图3(c)中以F点展开E点,以A点展开B点,以G点展开H点。若处在可靠度较小组中像素点除被展开像素点外还有其他像素点,则将处在可靠度较小组中的像素点展开前后的相位差加在所在组中其他像素点上,如图3(e)中以D点展开G点后,将G点展开前后相位差加在H点上。展完之后,合并两组。

4)判断下一条可靠度较高的边,重复过程2)、3),直到所有的边全部判断结束,所有像素点均展开完毕。

图 3. 快速相位展开算法的(a)像素点可靠度, (b)边可靠度定义和(c)~(e)相位展开过程

Fig. 3. Definition of (a) reliability of pixel points, (b) reliability of edges, and (c)~(e) process of phase unwrapping of FPU algorithm

下载图片 查看所有图片

快速相位展开算法实质上将截断相位分布按照边可靠度排序结果逐步展开,最终实现所有像素点属于同一个组。

枝切法的优点是展开速度快,但该方法在残点分布密集的地方枝切线容易闭合,形成孤立区域,造成区域性展开误差。最小不连续法虽然可以得到较平滑的展开相位,抗局部噪声能力较强,但在高噪声区域会严重扭曲相位,且大量迭代过程会耗费较长时间。快速相位展开算法速度快,但可能会在展开过程中穿过孤立区域,在一定区域内造成误差传播。通过对比研究,可以了解各个算法之间的优点。融合各个算法的优点形成新的算法即为本研究的思想。

3 计算机仿真模拟抗噪实验

用计算机模拟生成标准连续的相位图像,其函数关系式为:

f(x,y)=20exp-x2+y24+2x+y,(-3<x<3,-3<y<3),(1)

在X和Y方向上分别生成512个抽样点,即图像大小为512 pixel×512 pixel,满足抽样定理。其连续相位与截断相位如图4所示。

图 4. 仿真模拟的(a)连续相位及其(b)截断相位

Fig. 4. Simulation of (a) continuous phase and (b) its corresponding wrapped phase

下载图片 查看所有图片

在图4(a)所示的连续相位上增加随机噪声A,其参数如表1所示,模拟了叠加噪声逐渐增强的6种情况(分别为表1中噪声A的1,1.33,1.67,2,2.33,2.67倍)。

使用三种算法对分别受6种噪声干扰的截断相位进行展开,对应展开相位如图5所示,从左到右第1列代表带有6种噪声的截断相位,后三列分别代表利用枝切法、最小不连续法、快速相位展开算法的相应展开结果。利用三种算法展开6种截断相位的结果误差和消耗时间分布如图6所示。

图 5. 利用三种算法对6种噪声的展开结果。噪声分别为(a) A、(b) 1.33A、(c) 1.67A、(d) 2A、(e) 2.33A和(f) 2.67A

Fig. 5. Unwrapped results of six noises by using three algorithms. The noises are (a) A, (b) 1.33A, (c) 1.67A, (d) 2A, (e) 2.33A, (f) 2.67A, respectively

下载图片 查看所有图片

图 6. 三种算法展开6种噪声干扰截断相位的(a)误差变化与(b)耗时变化

Fig. 6. Changing of (a) error and (b) time consuming that using three algorithms to unwrap six kinds of wrapped phases with noise

下载图片 查看所有图片

由图6可知,当所加噪声比较小时,三种算法均可以正确展开截断相位。随着噪声的增大,枝切法与快速相位展开算法的误差均快速增大,但最小不连续法的误差变化不大。枝切法由于连接枝切线在高噪声区域自我闭合形成孤立区域,从而导致有部分像素点展开错误。快速相位展开算法则因为在噪声区域无法选择正确的展开路径,且在噪声分布密集的区域发生误差叠加,从而展开为错误的结果。最小不连续法在高噪声区域由于平滑展开相位而导致部分相位展开错误。最小不连续法由于不断迭代,其所消耗时间远大于其他两种算法。

从抗噪性能模拟实验中可以看出,三种相位展开算法均有一定的抗噪性能,其中,最小不连续法是最后在增强噪声中展开失败的,相比其他两种算法,其抗噪性能最强。但最小不连续法也有其缺点,即展开相位所用时间长于其他两种算法。

4 两种融合的相位展开算法原理以及性能

4.1 融合算法原理

使用枝切法对相位进行展开时,由于在高噪声区域连接残点形成枝切线时会自我闭合形成孤立区域。在展开相位的过程中,若展开路径穿过由枝切线构成的孤立区域时,可能将残点引起的展开错误向其他像素点传播开来,但它具有展开速度快的优点。最小不连续法的抗噪性能较强,但其缺点是在算法实现过程中由于大量使用递归,内存消耗大,运行时间长。故将枝切法与最小不连续法结合,利用枝切法将孤立区域之外的点全部展开,之后再使用最小不连续法展开孤立区域内的点并进行优化,此时大部分截断已使用枝切法展开,只需要少量迭代即可展开整幅截断相位图。因此会达到速度较快且精度较高的效果。

快速相位展开算法中,由于没有枝切线的引导,在展开过程中遇到高噪声区域时无法选择正确的路径,因此会出现展开路径穿过孤立区域引起误差传播的情况。而最小不连续法由于大量迭代使得相位展开耗时长,故将两种算法融合。对整幅相位图设定一个合适的阈值,大于阈值的区域认为是高质量区域,小于阈值的区域认为是低质量区域。先使用快速相位展开算法展开高质量区域,之后再使用最小不连续法展开低质量区域并进行全局优化,这样既避免了快速相位展开由于高噪声区域无法选择正确的路径而引起的相位展开错误,又解决了最小不连续法耗时长的问题。

为了表达的方便,将枝切法与最小不连续法融合后的算法记为BC-MD,将快速相位展开算法与最小不连续法融合后的算法记为FPU-MD。

4.2 计算机仿真模拟5种相位展开算法性能比较

计算机模拟生成一个大小为512 pixel×512 pixel、相位值为5.6倍peaks函数的连续相位,并在其上叠加了大小为2A的噪声,所得到的连续相位和截断相位如图7所示。随后用5种相位展开算法分别对相位进行展开,比较各算法的性能。

利用三种基本算法和两种融合算法分别展开上述截断相位,其误差标准差和其他算法相对于所耗时间最长的算法—最小不连续算法的比率如表2所示。5种算法的展开结果及其误差如图8所示,图8右边的色度图用于表征相位值。

图 7. 模拟加噪声的(a)连续相位及其(b)~(c)截断相位

Fig. 7. Simulation of (a) continuous phase and its (b)~(c) corresponding wrapped phase with noise

下载图片 查看所有图片

图 8. 分别使用(a) BC、(b) MD、(c) FPU、(d) BC-MD和(e) FPU-MD算法得到的相位展开结果(第1行)及其误差(第2行)

Fig. 8. Phase unwrapped results (first row) and their errors (second row) using (a) BC, (b) MD, (c) FPU, (d) BC-MD, and (e) FPU-MD algorithms, respectively

下载图片 查看所有图片

由表2及展开结果可知,相比于最小不连续法,融合后的两种算法展开相位时所消耗的时间均已减少,且展开的相位误差相对于枝切法和快速相位展开算法均较小。从最小不连续法特点来分析,可以看出融合后的两种算法在迭代过程中所形成的边数均小于使用最小不连续法所形成的边数。融合后的算法展开结果仍有误差,其原因是最小不连续法的目的是使得目标函数最小,因此在高噪声区域具有平滑相位的特性,将展开相位进行了平滑,与原相位相比产生了误差。

4.3 实物实验

本研究还使用两个实物实验数据对融合后算法的有效性进行证明。和相位测量轮廓术^[22-23]等其他三维传感技术相比,傅里叶变换轮廓术^[24-25] 具有只需要一帧变形条纹图即可恢复出物体三维形貌的优点。因此本研究使用傅里叶变换轮廓术获得物体变形条纹。实物图像大小为1024 pixel×1280 pixel,分别投影一幅正弦条纹图到两个被测物体表面,变形条纹为图9(a)所示的佛像和图9(c)所示的面具,其截断相位通过傅里叶条纹分析方法计算得到,

图 9. 佛像的(a)变形条纹及其(b)截断相位;人脸面具的(c)变形条纹及其(d)截断相位

Fig. 9. (a) Deformed fringes and (b) wrapped phase of a Buddha sculpture; (c) deformed fringes and (d) wrapped phase of a face mask

下载图片 查看所有图片

如图9(b)、(d)所示。实验结果中以最小不连续法展开相位的消耗时间为参考量进行比较。

图 10. 分别使用(a) BC、(b) MD和(c) BC-MD算法对佛像截断相位的展开结果

Fig. 10. Unwrapped results of the wrapped phase of Buddha sculpture using (a) BC, (b) MD and (c) BC-MD algorithms, respectively

下载图片 查看所有图片

分别利用枝切法、最小不连续法以及两种算法融合后的算法展开佛像相位,其结果如图10所示。利用枝切法和融合算法展开相位所消耗的时间分别是最小不连续法所消耗时间的1.64%和13.82%。最小不连续法和融合后的算法在实现过程中所形成的边数分别为84795704和10618675,环数为778和247。由于融合算法中的枝切法已经将孤立区域之外的像素点全部展开,此时全场的跳跃数已大大减少,因此使用融合后的算法进行相位展开时所形成的边数少于最小不连续法所形成的边数。而由于枝切法的展开路径是连续的,因此其形成的环数也相应减少。

图 11. 分别使用(a) FPU、(b) MD和(c) FPU-MD算法对人脸截断相位的展开结果

Fig. 11. Unwrapped results of the wrapped phase of face mask using (a) FPU, (b) MD and (c) FPU-MD algorithms, respectively

下载图片 查看所有图片

分别利用快速相位展开算法、最小不连续法以及两种算法融合后的算法展开人脸相位,其结果如图11所示。利用快速相位展开算法和融合算法展开相位所消耗的时间分别是最小不连续法所消耗时间的1.11%和51.84%。最小不连续法和融合后的算法在实现过程中所形成的边数分别为284445573和149575906,环数分别为2097和17865。使用融合后的算法所形成环数多于最小不连续法的原因是使用快速相位展开算法展开高质量区域后,由于快速相位展开算法的展开路径是离散的,展完之后整幅图像相位被分成多个离散小区域,因此会将本来在最小不连续法中形成的一个大环在融合后的算法中分散形成多个小环。而因为快速相位展开算法已经将高质量的区域展开,此时全场的跳跃数已减少,因此使用融合后的算法进行相位展开时所形成的边数少于最小不连续法所形成的边数。

5 结论

本文介绍了枝切法、最小不连续法以及快速相位展开算法三种空间相位展开算法以及两种融合后的算法。利用模拟实验分析讨论了三种基本算法的优点和缺点,结果表明枝切法与快速相位展开算法均具有展开速度快的特点,但其抗噪性能均较差。而最小不连续法相比于其他两种算法具有较强的抗噪性能,但其缺点是展开相位所消耗时间较长。

基于三种算法的性能比较,将枝切法和快速相位展开算法分别与最小不连续法进行融合,得到两种新的算法。新算法结合了原有算法的优点,避免了原有算法的缺陷,融合后的算法由于在使用最小不连续法之前,整场的跳跃总数已经减少,因此所形成的边数减少,进而节省了时间,因此使用新算法的展开结果既具有较快的速度,又具有较高的精度。计算机仿真模拟和实物实验结果证明了两种新算法的可行性。