可用于单幅闭合干涉图相位恢复的正则化相位跟随技术

王贤敏，刘 东，臧仲明，吴 兰，严天亮，周宇豪，张与鹏

(浙江大学 光电科学与工程学院 现代光学仪器国家重点实验室，浙江 杭州 310027)

1 引 言

由于干涉仪具有响应速度快、测量精度高等优势，已被广泛应用于光学精密检测等行业[1-3]。但实际应用时，为了减少低频噪声对测量精度的影响，通常需要通过一定的调制手段，如移相干涉、载波干涉等方法[4]，将待测信息加载到干涉条纹图样中。移相干涉法可以获得较高的检测精度，但需要采集多幅间具有一定相位差的干涉图。这就要求外界环境必须非常稳定，环境中存在的机械振动、空气扰动等均会影响检测精度，因此一般不适合现场实时检测。而将载波干涉法与傅立叶分析技术相结合，仅需一幅干涉图即可获得被测信息，对噪声的抑制效果较好，比较适合对时间响应要求较高的应用场合。然而载波干涉需要在两路干涉波前之间加入倾斜或者离焦，当被测波前畸变较大时，会引入较大像差，从而影响检测精度。近年来，人们也提出了一些新颖的从干涉图中准确恢复相位的方法，如利用未知相移[5-6]的两幅干涉图反复迭代，也有基于单幅干涉图的多项式拟合[7]、基于黄氏希尔伯特变换[8]以及窗口傅立叶分析[9]等方法。这些方法有各自的优势，甚至只需要一幅干涉图就能恢复出相位，而且比之前提到的传统方法局限性更少，但仍然需要进一步的解包裹处理，或者只能检测小于一个波长的相位变化量。

在恢复单幅闭合条纹干涉图的相位信息这个问题上，Servin等人提出的二维正则化相位跟随技术(Regularized Phase Tracking，RPT)[10]是一种非常有效的方法，具有结果连续无跳变、自动去除高频噪声等优点，但RPT技术在很多方面都有改进的空间，如需要预处理干涉图从而去除噪声并归一化调制度、处理时间长、恢复精度不高、在处理鞍点和极值点等关键点时易出错、恢复效果尚依赖于扫描路径等。后期，研究者们从不同方面对RPT技术进行了改进，使得RPT技术在单幅、闭合条纹干涉图相位恢复方面得到了极大应用和发展。刘东等人[11]通过在第一次求解结果的基础上构造新的优化函数进行最优解的二次搜索，显著地提高了RPT技术的相位恢复精度。Tian等人[12]则在邻域内的相位拟合方式上进行改进，如用二次曲面拟合替代平面拟合，摆脱了对扫描路径的依赖性。Li等人[13-14]提出通用RPT技术，将背景项、调制度项以及相位均看成是邻域内的平面可拟合项，全部纳入到优化函数的优化范围内，省去了预处理中的背景去除和归一化调制度过程，可以直接处理实际干涉图。Deepan等人[15-17]通过提前求得的相位变化率信息，减少了优化函数中的未知数个数，使求解时间大幅减少。

本文介绍了RPT技术用于单幅干涉图相位恢复的基本原理，总结了近年来RPT技术的相关改进与发展，例举了采用RPT技术进行相位恢复的应用场合，并适当推测了RPT技术的未来发展方向。

2 正则化相位跟随技术

干涉仪得到的干涉图可以用式(1)来表示，

Ix,y=ax,y+bx,ycosφx,y+nx,y,

(1)

其中，I表示记录的干涉图的强度，a表示背景光，b表示调制度，φ表示相位信息，n表示干涉图中的噪声，下标x,y表示像素坐标。

在对单幅干涉图进行相位恢复前，需要对其进行去背景、降噪、归一化调制度等操作，得到式(2)。

(2)

目前已有的研究[18-23]已经可以很好地实现从式(1)到式(2)的过程，利用RPT技术恢复相位将在式(2)的基础上进行。一般来说，由于反三角函数的周期性，很难直接从式(2)恢复出连续无需解包裹的相位信息，而RPT技术可以实现。

2.1 基本原理

利用RPT技术恢复单幅干涉图的相位基于两点先验假设：

(1)恢复后的相位对应的强度信息与干涉图一致；

(2)待恢复的相位在空间上连续变化。

利用这两点先验假设，可以构造一个优化函数式(3)，表示恢复出的相位与先验假设之间的差异，优化过程则是不断缩小这个差异的过程。

λ[φ0(ε,η)-φe(ε,η)]2m(ε,η)} ,

(3)

式(3)中，(ε,η)表示点(x,y)周围邻域Nx,y内的像素坐标；φe(ε,η)为恢复该点相位过程中邻域内各点的相位计算值；φ0(ε,η)为该点恢复出来的相位；m(ε,η)为该点相位信息是否已被恢复的标志，为1时表示已恢复，为0时表示未恢复；λ为正则化系数。每一点(x,y)上的待恢复相位信息需在其邻域Nx,y内使式(3)中的代价函数Ux,y达到全局最小。在邻域Nx,y内，每一点(ε,η)的相位计算值φe(ε,η)是通过线性相位模型式(4)进行计算得到的。

φe(ε,η)=φ0+ωx(ε-x)+ωy(η-y) ,

(4)

从式(3)和式(4)来看，点(x,y)在其邻域内的优化函数Ux,y是关于φ0、ωx和ωy的函数，φ0、ωx和ωy则可以表示对邻域Nx,y中相位的估计情况。Ux,y与φ0、ωx和ωy的关系，可以用式(5)来表示。

Ux,y=Ux,y(φ0,ωx,ωy) .

(5)

梯度下降算法是求解凹函数极小值的有效搜索算法，针对式(5)，可以采用式(6)的梯度下降算法迭代求解使Ux,y为极小值时的p=[φ0ωxωy]T，式中τ为迭代步长，k表示第k次迭代。

(6)

利用RPT技术对一幅预处理后的干涉图的相位进行恢复时，一般遵循下面的过程：

步骤2：选择点(x,y)，首先检查该点的相位信息是否已经恢复，若未恢复，进入步骤3；若已恢复，重新选点；

步骤3：搜索邻域Nx,y是否有相位信息已经恢复的点，即找到m(ε,η;x,y)=1的点，若有，进入步骤4；若无，跳回步骤2；

步骤4：根据找到的m(ε,η;x,y)=1位置的相位信息[φ0ωxωy]T，将点(x,y)处的相位初始值设为φe(ε,η)=φ0+ωx(ε-x)+ωy(η-y)，利用式(6)迭代求解完成后，设置m(x,y)=1，进入步骤2，直至干涉图上所有点的相位都已被恢复。

2.2 技术分析

Servin等人认为利用RPT技术对单幅干涉图进行相位恢复时需要先对其进行归一化处理，如图1所示，将条纹强度数据归一化到-1～1之间。由于受到环境噪声、背景不均匀、调制度变化等因素的影响，探测器实际得到的干涉图，其图像

图1 干涉图灰度归一化处理 Fig.1 Interferogram grayscale normalization

强度数据往往复杂不规范，需要通过预处理将图像强度数据归一化。目前已有多种方法可以很好地归一化实际干涉图[18-23]，便于利用RPT技术恢复单幅干涉图的相位。

实际干涉图归一化后，在求解起始点的相位时，其邻域内并不存在已经恢复相位信息的点。此时的优化函数并不是一个严格的凹函数，利用梯度下降法搜索最优解时会陷入局部最小值。因此，一般在恢复起始点的相位信息时采用全局搜索的方式。而在后续点的相位恢复中，邻域内相位信息已经恢复的点的约束作用使优化函数转变为严格的凹函数，梯度下降法搜索到的极值点即为全局最优。由于斜率相反的两个相位平面可以有相同的强度分布，所以在对起始点求解时往往会找到两个全局最小值点，这两个点对应的相位平面具有相反的斜率，如果一开始选择了错误的斜率方向，那么整个干涉图恢复出来的相位结果就会和实际结果相反。应用线性拟合方式拟合邻域内各点相位时，是无法避免这个问题的。由于相同的原因，在后续遇到相位变化率为零处的鞍点时，很容易出现判断错误，此时后果更加严重。起始点最小值点选错，会使整个相位颠倒，但不影响整体判断，而鞍点处判断错方向会则使相位整体变形。在恢复鞍点处的相位时，如果邻域内存在较多已被正确恢复相位信息的点，那么由于这些点的约束作用，优化过程会自动往正确方向进行。预先设计好扫描路径，将鞍点有意识地放在最后并被大部分已经恢复相位的点包围，可以有效避免错误的发生。

利用RPT技术从单幅干涉图中恢复相位的过程本质上是从φ0、ωx和ωy3个参数出发，利用优化算法搜索使Ux,y达到全局最小的点的过程。应用梯度下降法搜索时，较多的未知数会使搜索过程曲折，达到最优解位置的时间变长。从2.1部分可以看出，式(6)中的步长τ会影响RPT技术从单幅干涉图中恢复出相位的速度。步长τ设置得大，能较快到达最优解的位置，但可能会在最优解附近陷入振荡；而步长τ设置得小，则需要经历更多次迭代才能到达最优解位置。程序中往往比较两次连续迭代结果的差值，如果差值小于阈值tol，则可以认为此时的φ0、ωx和ωy为最优解的位置。阈值tol的大小限制了理论精度，如果需要较高的精度，必须设置较小的步长τ和阈值tol，这意味着计算时间的增加。而RPT技术的精度有限，过小的步长τ和阈值tol除了增加计算时间外并不能显著地提高精度，因此没有必要选择非常小的参数。影响RPT技术求解相位快慢的因素除了步长τ和阈值tol，还有邻域Nx,y的大小，较大的Nx,y也会使计算时间增加，但适当增大邻域范围可以加强相位恢复过程中的约束作用，使后续恢复过程更稳定。

3 正则化相位跟随技术的发展

利用RPT技术从单幅干涉图中恢复相位信息的主要局限如2.2部分所述。针对这些局限，学者们主要从预处理和优化函数模型两个方面对RPT技术进行改进。

3.1 预处理与扫描路径

预处理的目的是从式(1)获得式(2)，实验中CCD采集到的干涉图中往往具有各种噪声，强度数据如同图1(a)所示。而RPT技术的实现是假设去除了噪声以及背景调制度的影响。因此为了对实际干涉图进行RPT处理，需要先进行一些预处理步骤，一般包括简单去噪、背景去除。其中简单去噪指中值滤波、均值滤波、高斯滤波等等；背景去除比较复杂，需要在干涉图对应的频域空间进行一系列操作。对干涉图的预处理过程一直也是研究热点，人们一直致力于对其进行高效降噪，现有的方法已经足以满足RPT技术需求，故其不作为本文重点，因此不做过多介绍。Quiroga等人在2001年提出的采用两个正交的带通滤波器共同作用的方法取得了很好的效果[19]。除此以外，Bernini等人将经验模式分解与Hilbert-Huang变换应用于条纹预处理中[20]，也得到了很好的滤波效果。

如2.2部分所述，合适的路径引导可以降低鞍点、极值点等关键点的选择错误概率，Servin等人提出利用条纹跟随(Fringe-follower RPT,FFRPT)[24]的方式规划扫描路径，引导相位恢复过程。FFRPT并没有在相位恢复过程中修改算法，只是指出有效的扫描路径规划，从而避免了相位恢复过程中的错误，并提出了基于条纹分布特点规划扫描路径的思想。FFRPT应用合理的扫描路径规划方式，有效地减小了关键点处的错误概率，可以从一幅复杂的干涉条纹中较好地恢复出相位信息。如图2所示，2(a)为计算机仿真生成的256 pixel×256 pixel的复杂条纹分布，2(b)为根据条纹分割出来的分布结果，2(c)为利用FFRPT技术从干涉图中恢复出来的相位分布。

图2 FFRPT解调结果，单位：波长 Fig.2 Results with FFRPT, Unit: Wavelength

3.2 优化函数模型

RPT技术的最主要特征就是在当前点某一邻域内通过构建优化函数，搜索当前点的相位信息，故优化函数模型的差异将直接影响其搜索效果。于是，很多研究人员都针对优化函数模型开展深入研究。

3.2.1 增加相移项

Servin等人提出，在经典RPT技术优化模型中，增加一个相移项[10](相移α一般在0.1π到0.3π之间)，可以加速优化函数收敛，如式(7)所示。应用式(7)进行求解时，一般是在第一次相位求解时设置一个不等于零的α，得到比较接近于真实情况的相位信息，再利用这些相位信息作为初始条件，并将α设置为零后进行第二次求解。

λ[φ0(ε,η)-φe(ε,η)]2m(ε,η)} .

(7)

图3 干涉图与利用恢复相位反求的干涉图 Fig.3 Interferograms before and after reversed phase recovery

干涉图与利用恢复相位反求的干涉图如图3所示，图3(a)为条纹图，图3(b)、3(c)和3(d)是α分别为0、0.1π和0.2π时，利用从干涉图3(a)中初次恢复出来的相位反求的条纹图，图3(e)为它们在同一位置(直线标记)上的截面视角。

由图3可以看出，α=0时，恢复出来的相位对应的强度分布与原图高度重合，而α=0.1π和α=0.2π时则出现了明显的不一致，这也是再次设置α=0进行二次求解的原因。

3.2.2 增加调制度项

经典RPT技术的一个局限是需要对条纹强度进行归一化处理，即要应用RPT技术求解相位必须先得到式(2)，而此过程所引入的误差必定会影响后续的相位恢复。Ricardo等人[25-26]将干涉图的调制度bx,y也作为一项未知数，并利用式(8)作为调制度拟合方程，将优化函数修正为式(9)。式中IB表示实际强度，Be表示迭代过程中的调制度值，B0、βx和βy则分别表示某点的调制度大小、调制度沿x方向的变化速率和沿y方向的变化速率，μ与λ表示正则化系数，为简洁起见，部分坐标(ε,η)已省略。

Be=B0+βx(x-ε)+βy(y-η) ,

(8)

{IB-Becos[φe+α]}2+

λ[φ0-φe]2m(ε,η)+

μ[B0-Be]2m(ε,η)} .

(9)

在这个模型中，Ricardo等人将调制度纳入求解过程，假设调制度变化连续且在一定区域内可看成平面。这样做不需要在预处理过程中归一化调制度，但缺点也很明显，即令求解更加复杂，因为解向量变为p=[φ0ωxωyB0βxβy]T，未知数个数增加了1倍，计算一幅干涉图的相位需要更长的时间。Ricardo等人为了解决计算时间过长的问题，将优化函数进行了修正，减少了未知数的个数，并利用高斯-赛德尔迭代方法进行最优值搜索[25]。

图4 对未归一化调制度干涉图的相位恢复(已包裹) Fig.4 Phase recovery of unnormalized interferogram(rewrapped)

如图4所示，干涉图整体的调制度(图4(b))呈现不均匀分布，通过该方法恢复出来的相位结果如图4(c)所示，处理时间只需原来的五分之一。

3.2.3 路径无关正则化相位跟随

在处理类似于图5中B点那样的极值点时，利用传统RPT技术求解相位信息时会发生错误，原因如2.2所述，即利用RPT技术进行单幅干涉图的相位恢复时，对于鞍点、极值点等关键点，线性拟合不具有判断相位变化方向的能力。路径无关正则化相位跟随(Path-Independent Regularized Phase Tracker[12],PIRPT)技术中将式(4)的线性拟合修改为式(10)所示的二阶拟合，即利用小曲面代替平面，并修正了相应的优化函数式(11)，式中的ωxx(·)、ωxy(·)和ωyy(·)表示曲面方程的二次项。利用曲面的二阶拟合，在鞍点、极值点处也可以判断出相位的变化方向，即可不需要单独设计扫描路线将关键点放到最后去处理，因此其相位恢复效果与扫描路径无关。图5(a)中的相位极值点处，二阶拟合方式比一阶拟合适应性更好。由于其能够判断相位变化的方向，PIRPT技术应用于单幅干涉图相位恢复中具有处理复杂干涉条纹的能力，算法稳定性较经典RPT技术也有了较大提高。但缺点是在计算过程中涉及到6个未知数的迭代优化，计算变得更加复杂，相应地，处理一幅干涉图所需的时间也大大增加。

图5 路径无关原理示意图 Fig.5 Path independent principle diagram

(10)

(11)

3.2.4 广义正则化相位跟随

Li等人[13-14]提出广义正则化相位跟随(Generalized Regularized Phase Tracker,GRPT)，即将式(1)中的背景a和调制度b均看成是邻域内的线性平面调制，将相位φ看成是二次曲面调制，如式(12)所示。优化函数则调整为式(13)，其中w(ε,η;x,y)是标准差为σ的高斯权重窗口，Li将窗口大小即邻域设为W×W，此时，标准差σ=W/4，式中已省略坐标(ε,η)。

(12)

λ[φ0-φe]2m} ,

(13)

(14)

这样做，可以直接对CCD记录得到的干涉图进行相位恢复，省略了背景归一化过程，但这样做的缺点也非常明显，优化过程中未知数个数达到了12个，计算时间变得冗长。Li等人利用了Levenberg-Marquardt(LM)和Broyden-Fletcher-Goldfarb-Shanno(BFGS)方法优化了搜索最优解过程[14]。GRPT对干涉图的相位恢复结果如图6所示，图6(a)为实际的干涉条纹，图6(b)为利用GRPT恢复的相位，图6(c)为恢复过程中各点的迭代次数。虽然GRPT技术已成功应用于多种干涉图的相位恢复，但在面对较为稀松的条纹图时处理效果一般较差。Li等人后期又将高斯窗口的大小修正为自适应调整模式，并在优化函数中增加了背景正则项和调制度正则项[13]，从而解决了这个问题。

图6 GRPT对干涉图的相位恢复(已包裹) Fig.6 Phase retrieval of interferogram with GRPT(rewrapped)

3.2.5 正则化稳频器

Tian等人[27]将稳频器(Regularized Frequency Stablizer, RFS)的概念引入到RPT技术中，充分利用邻域内恢复区域的信息，结果如图7所示。

图7 RFS与RPT对同一幅干涉图解调结果 Fig.7 Demodulation results of RFS and RPT in the same interferogram

RFS的恢复过程可以分为粗恢复和细化两个过程，粗恢复的结果比较粗糙，在此基础上细化可进一步去除噪声，达到较好的结果。

(15)

(16)

(17)

(18)

从式(15)～(18)可以看出，优化函数只有一个未知数，计算时间缩短。优化过程中减少的未知数即降低了计算复杂度，又使计算更稳定。

3.2.6 简易正则化相位跟随

Deepan等人[16-17]提出一种简易正则化相位跟随(Simplified Regularized Phase Tracker,SRPT)方法，其通过图像处理手段得到干涉图中的条纹密度和方向信息，从条纹密度和条纹方向出发推导对应相位在x和y方向的变化频率，利用推导出来的频率信息代替式(4)中的ωx和ωy有，

(19)

式中，N为条纹密度信息，θ(·)为条纹的方向。

图8 SRPT对单幅实际干涉图的相位恢复效果 Fig.8 Phase retrieval results of interferogram with SRPT

3.2.7 二次优化模型(RPT&GS)

RPT技术的一个局限在于它所恢复的结果精度不高，无法与移相法相比拟，这是由其基于数值分析的本质决定的。刘东等人[11]利用第一次恢复的结果作为初始条件，重新构造优化函数进行优化，将相位恢复精度提高到可与移相干涉仪相当的水平。刘东等人提出，在获取初次优化结果后，可以再构造一优化函数，如式(20)所示。 用其进行二次优化[11]。

(20)

在优化函数式(3)中，邻域内点的约束作用可以使相位朝着连续变化的方向进行求解，但也影响了求解精度。二次优化法的基本思想即，在第一次求解的基础上，利用优化函数式(20)，通过二次优化过程搜索出最优解。RPT&GS与RPT的效果对比如图9所示。可见RPT&GS的二次优化结果误差要比RPT低一个数量级。

图9 RPT&GS与RPT精度对比 Fig.9 Precision comparison of retrieval results by RPT&GS and RPT

4 正则化相位跟随技术的应用

RPT技术主要是针对无法进行调制的干涉图的相位提取问题提出的，故其精度始终无法与移相法等相媲美。但由于其仅需要一幅干涉条纹图即可恢复其中的相位信息，而无需任何附加调制，使得其应用也十分广泛。

4.1 应用于条纹偏折图的分析

莫尔条纹常用于分析光线偏折，分析的对象往往是偏折的带深度信息的条纹图，一般利用移相算法或傅立叶分析法来提取莫尔条纹偏折图的信息。但上述两种方法抗噪性均较差，从恢复出来的相位结果中确定梯度区域非常困难。Quiroga等人[28]将修正的RPT技术应用于莫尔条纹图分析过程中，以确定眼科镜片表面能量分布，具有较好的抗噪性。Villa等[29]在测量物体3-D形状轮廓时，提出了相位跟踪轮廓仪的概念，并应用了RPT相位恢复技术。

4.2 应用于Hartmann检验

Hartmann检验常用于检验大型望远镜，通过在瞳面放置均匀排布的阵列传感器对波前进行采样分析得到其像差分布。由于阵列各个面心与理想无像差成像时的偏离量理论上跟各面心所在位置的非球面波前曲率成正比，故通过分析记录得到的条纹分布情况与面心阵列理论情况相比较，即可得到光学系统引入的像差。Servin等人[30]将RPT相位恢复技术应用于Hartmann阵列图，来计算横向像差，即将得到的Hartmann阵列图作为干涉图进行相位恢复，推测横向像差分布情况。Servin等人将这种新的方法与人工计算以及傅立叶分析法进行比较，认为RPT技术在该场合中的应用具有优越性。

4.3 应用于方形光栅条纹图分析

方形光栅可以把两个正交方向的表面偏差信息相乘于一幅图像中，从而不需要旋转光栅。RPT技术由于良好的抗噪特性被应用到恢复方形光栅条纹图的相位分布，Villa等人[31]在分析方形光栅条纹图时将优化算法修改为拟牛顿优化算法，并在扫描路径上利用基于条纹质量引导的算法进行有效规划，得到了较好的结果。

4.4 应用于光测弹性仪分析

光测弹性应力分析是分析透明物体应力与应变的重要实验技术，具有悠长的发展历史[32-33]，常被用于大型高架桥梁等大型复杂结构的自重效应分析中。应力与应变条件下透明材料内部的双折射率发生变化，在光照下有不同的条纹分布。通过分析条纹的形状，间接可以得到物体应力与应变信息。光测弹性分析有等倾线与等色线之分，前者指试样上主应力方向相同的点组成的轨迹，后者则是指试样上两个主应力差值相同的点的轨迹即最大切应力相同的点的轨迹。RPT技术应用于分析光弹性条纹[34]，尤其是应用于分析等倾线与等色线处[35]的相位信息，在光测弹性分析中发挥了重要的作用。

5 结束语

将RPT技术应用于单幅干涉图的相位恢复，可以得到连续的相位分布，省去解包裹处理过程，潜在的滤波作用使其对干涉图的噪声也有很强的适应能力。从其首次被成功应用于单幅干涉图相位恢复以来，研究者们在此基础上做出了许多改进，这些改进使RPT技术在单幅干涉图相位恢复上保存自己原有优势的同时，也渐渐摆脱了一些弊端，适用性越来越强，甚至可以利用RPT技术实现动态观测。RPT技术可以从单幅干涉图中恢复出连续相位的特有优势一直吸引着研究者们进行更深入的研究。其大致可以分为两种：一种是构建更加准确同时也更加复杂的相位优化函数，在优化过程中进行修正，这是早期大多数研究努力的方向；另一种则是放弃了复杂化的相位优化函数，尽量利用图像处理方法得到更多的相位衍生信息，求解相位时则充分利用这些信息简化计算过程，使整体过程趋于简短高效，是近些年来趋于主流的方法。从现有的基础及趋势来看，将单幅干涉图相位恢复中的RPT技术纳入到图像处理以至图像分析甚至图像理解中更为合适。最开始，RPT技术更多被看做是一种数学方法，利用建模和优化的思想找到最优解，从而恢复相位。而在现在，我们更多地会去考虑条纹分布与相位分布合理性的问题，这已经从单纯的数学约束问题跨入了广义上图像理解的领域。在图像理解上，机器学习方法正被视作解决问题的好方法，并成为研究的热门方向，相关的理论与硬件基础都已达到一个历史高点。借助机器学习的理论与工具，处理图像问题会更加得心应手，将它们应用于干涉条纹相位恢复的应用研究中，对于扫描路径规划、相位变化率及方向判断、相位恢复结果评价等方面都有很大的意义，是目前尚未有太多人涉足但非常具有潜力的研究方向。