CN108170636B - 一种从线性载频移相干涉图中恢复出真实相位的方法 - Google Patents

Abstract

本发明涉及一种从线性载频移相干涉图中恢复出真实相位的方法。首先，选取三幅移相干涉图，通过三幅图像的像素尺寸确定分组；然后，利用评价函数计算出所有分组数据的全局移相量，再利用公式得到每幅图之间的全局移相量；最后，利用最小二乘原理得到压包相位分布，再经过展开操作、二维一次项拟合、消除线性载频，最终恢复出真实相位分布。通过本发明提出的算法可以准确快速的恢复出待求的真实相位分布，所需要的移相次数少，缩短了采集数据的时间，算法对线性载频的频率、干涉图的对比度以及背景光强的均匀性均没有限制，因此算法的鲁棒性较高，可以广泛地应用于干涉计量与测试、条纹投影、全息干涉等领域。

Description

一种从线性载频移相干涉图中恢复出真实相位的方法

技术领域

本发明涉及移相干涉测量技术邻域，具体涉及一种从线性载频移相干涉图中恢复出真实相位的方法。

背景技术

移相干涉术是一种应用于光学干涉测量中的技术。该技术的精度高，稳定性好，硬件实现方式多样，四十多年来，始终受到研究者们的广泛关注。在移相干涉术中，需要采集一幅以上的干涉图，一般来说，若要在存在背景光强、调制度不一致的情况下，能够准确恢复的真实相位，需要采集干涉图的数量大于等于三幅，并且在这些干涉图之间需要引入一定的全局移相量。经典的移相算法包括固定步长的移相算法、可变步长的移相算法、随机移相的迭代算法等。近年来，研究者们提出了多种全新的移相算法，包括2帧的正则化条纹的移相算法、未知或未校正移相量的移相算法以及通用移相算法等，这些算法能够仅采集1-2幅干涉图，并从中求解真实相位分布，但是对于背景光强以及调制度需要有一些额外的要求。产生全局移相量的方法，最常见的方法是通过压电陶瓷促动器推动参考镜面进行移动，除此之外，还有波长调谐、随机移相、倾斜移相等多种方法。这些产生全局移相量的方法各有特点，在光学测量、数字全息术、条纹投影技术等方面有着重要的作用。在一些情况下，全局移相量是一个已知量，通过先验信息给出预期值，但是这个值在实际实验中经常会产生误差，例如环境振动的影响。在另一些情况下，全局移相量本身即为未知量，需要计算一系列干涉图之间全局移相量。

对于上述全局移相量的求解，Farrell和Player提出了基于李萨如图进行拟合的方法，Brug提出了计算干涉信号的相关强度的算法，Goldberg和Boker等人提出了基于单点傅里叶变换的方法，该方法通过比较两幅干涉图之间载频频率的功率变化，计算全局移相量，但是所有的干涉信号都是有限长度的，载频频率不是单一频率的，+1级频谱会出现频谱泄露，仅通过计算单一频率的功率变化不能完全反映出全局移相量的变化，造成了计算精度下降；Guo和Rong等人提出了基于能量最小傅里叶变换的算法，然而，测量采集的干涉图受到切趾效应、背景光强不均匀、信号包络、随机噪声等因素的影响，+1级频谱的边带可能会与零级频谱发生混叠，这一问题显著降低了计算结果的准确性，尤其是载频频率不高的时候影响明显；Vishnyakov和Levin等人提出了将两幅干涉图相减再进行傅里叶变换的方法，可以有效避免混叠现象的出现，使得该算法在载频很低的情况下，也可以有较高的计算精度，然而该算法在计算过程中，需要使用三幅干涉图，从而限制了应用的范围；此外，Bai和Rao提出了一种无须傅里叶变换的全局移相量计算方法，对干涉图的载频没有要求，可以在背景光强不均匀以及存在包络调制的情况下计算全局移相量，但仅用于恢复全局移相量为nπ/2时的相位分布；因此，现有的全局移相量求解算法仍然存在一些局限。

发明内容

针对现有技术的不足，本发明提出了一种从线性载频移相干涉图中恢复出真实相位的方法，通过对线性载频干涉图进行计算，在不存在测量噪声的理想情况下，对于全局移相量未知的移相干涉图，能够从中完全准确的计算出全局移相量，并且对载频频率没有特别要求，能够提高现有方法的鲁棒性和求解的精确度。

本发明的技术方案如下：

一种从线性载频移相干涉图中恢复出真实相位的方法，包括如下步骤：

1)从干涉测量仪器的图像采集系统中获取至少三幅线性载频移相干涉图，且两幅图像的像素尺寸均为P行和Q列，对应像素点(x,y)的光强分布分别定义为M₁(x,y)、M₂(x,y)和M₃(x,y)，每两幅干涉图之间的全局移相量是未知的，根据线性载频干涉图的定义，上述干涉图光强分布用下式表示：

其中，a(x,y)表示背景不均匀分布，b(x,y)表示调制度，本方法中认为这两项分布在三幅移相干涉图中保持一致，

表示相位分布，k₁和k₂分别为x和y方向的角波数，代表了线性载频的大小，

为待求的相位分布，假设M₁(x,y)的全局移相量为0，则Δ₁表示M₂(x,y)相对于M₁(x,y)的全局移相量，Δ₂表示M₃(x,y)相对于M₁(x,y)的全局移相量；

2)选取其中的M₁(x,y)和M₂(x,y)，根据该两幅图像的像素尺寸合理地选择N行和T列的光强数据，每组光强数据分别记作

及

其中，上标n＝1,2…N、t＝1,2,…T分别表示行数据和列数据的分组序数，下标1和2表示分别取自第一和第二幅干涉图；

3)计算每一组行数据

和每一组列数据

对应的评价函数，利用评价函数求解每组行数据和列数据对应的全局移相量dⁿ和d^t；

4)根据下式计算出M₁(x,y)和M₂(x,y)之间全局移相量的平均值d₁，也就是公式(2)中Δ₁的计算值，

5)同理，采用步骤2)至步骤4)相同的流程，计算出M₁(x,y)和M₃(x,y)之间全局移相量的平均值d₂，第一幅移相干涉图的全局移相量设定为d₀＝0；

6)根据前述所得的d₀、d₁、d₂、M₁、M₂及M₃，利用最小二乘原则，并结合公式计算得到包含线性载频的压包相位分布r(x,y)；

7)对步骤6)中得到的r(x,y)进行相位展开操作，得到包含线性载频的解包相位分布r’(x,y)；

8)对步骤7)中得到的包含线性载频的解包相位进行二维的一次项拟合，并消除线性载频，从而得到不包含线性载频的真实相位分布

所述步骤3)中求解每组行数据和列数据对应的全局移相量dⁿ和d^t过程为：

①选择一组行数据

进行计算，初始计算时n＝1；

②根据步骤①计算出的行数据分别利用下列公式计算出I_m(x)和I_p(x)：

③针对步骤②中得到的I_m(x)使用一维希尔伯特变换计算得到I_mc(x)：

其中，H[.]表示取一维希尔伯特变换；

④依据步骤②中得到的I_p(x)和步骤③中得到的I_mc(x)，利用下式计算评价函数：

E(δ′)＝∫_l|I_p(x)-I_mc(x)·cot(δ′/2)|dx， (8)

其中，l为该组数据的采样点数，δ’＝δ₀+iε，i表示步数，移相量的初值δ₀和步长ε依经验确定；

⑤依据步骤④中得到的评价函数绘制出评价函数的曲线，并找到最小值点，该点代表的全局移相量即为从这组数据中求得的全局移相量dⁿ＝δ’；

⑥判断是否已经按照上述流程处理完所有选择出的行数据组，如果没有完成，则n＝n+1，回到步骤①，代入下一组数据，重复步骤①至步骤⑤的流程，直至确定出所有的dⁿ；

⑦同理，采用步骤①至步骤⑥相同的流程处理完所有选择出的列数据组，确定出所有列数据对应的全局移相量d^t。

所述公式(8)中步数i的取值范围为200～400。

所述步骤6)中压包相位分布r(x,y)的计算过程为：

①利用最小二乘原则，将前述所得的d₀、d₁、d₂、M₁、M₂及M₃代入最小二乘方程(9)，计算出D₁、D₂、D₃，

其中，i＝1,2,3；

②将步骤①计算得到的D₂、D₃，代入下式即可计算出包含线性载频的压包相位分布r(x,y)：

本发明与现有技术相比，有益效果在于：通过本发明提出的方法可以计算任意两幅线性载频移相干涉图之间的全局移相量，通过三幅移相干涉图，可以准确快速的恢复出待求的真实相位分布，所需要的移相次数少，缩短了采集数据的时间，方法对线性载频的频率、干涉图的对比度以及背景光强的均匀性均没有限制，因此方法的鲁棒性较高，可以广泛地应用于干涉计量与测试、条纹投影、全息干涉等领域，具有优异的社会效益和经济效益。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1是本发明中提及的方法的流程图。

图2是本发明实施例中所使用的模拟干涉图一。

图3是本发明实施例中所使用的模拟干涉图二。

图4是本发明实施例中所使用的模拟干涉图三。

图5是本发明实施例中模拟干涉图中所设定的真实相位分布。

图6是本发明实施例中计算出的三行数据的评价函数曲线，并标出了最小值的位置。

图7是本发明实施例中所以计算出的压包相位分布。

图8是本发明实施例中经过解包的相位分布，其中包含了线性载频。

图9是本发明实施例中经过平面拟合以后的相位分布，其中已经消除了线性载频的部分。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细说明。

实施例：

如图2、3和4所示，分别为设定的待模拟采集的三幅线性载频移相干涉图，已知：图像像素尺寸为256x256，干涉图中的背景强度a(x,y)＝0.5(x²+y²)，调制度b(x,y)＝-0.1(x²-y²)+0.9，x和y两个方向上的角波数分别为k_x＝24，k_y＝20，其中待求解的第二和第三幅干涉图相对于第一幅图的移相量分别为Δ₁＝π/2和Δ₂＝π/4，待求解的相位分布如图5所示；

首先，从第一幅干涉图和第二幅干涉图中分别选择第10，127，230行与列的数据，根据公式(5)～(8)得到每行数据对应的评价函数，并得到相应的评价函数曲线如图6所示，图中标出了从上述三行中计算出的全局移相量，图中还标出了曲线的最小值，也就是待求的全局移相量；同理，可以计算出三列数据的全局移相量；然后再根据式(4)进行平均，从而得到第一和第二幅干涉图之间全局移相量d₁＝1.57rad，再重复前述步骤，计算出第一和第三幅干涉图之间全局移相量d₂＝0.785rad，再利用公式(9)和(10)计算出压包相位分布，如图7所示；接着再利用相位解包方法进行相位展开操作，得到了如图8所示的解包相位分布；因为解包相位中还包含了线性载频部分，所以还需要使用平面拟合的方法计算出线性载频系数，将其从恢复相位中消去，最终的恢复相位结果如图9所示；

可见，所得的全局移相量的值一致、相位的恢复图与原设定图一致，因此，该模拟实验验证了本发明方法的准确性、可行性及合理性。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明，本文所定义一般原理可以在不脱离本发明的精神或范围的情况下，在其他实施例中实现。本发明的范围由权利要求及其等同物限定。

Claims (4)

1.一种从线性载频移相干涉图中恢复出真实相位的方法，其特征在于：包括如下步骤：

1)从干涉测量仪器的图像采集系统中获取至少三幅线性载频移相干涉图，且两幅图像的像素尺寸均为P行和Q列，对应像素点(x,y)的光强分布分别定义为M₁(x,y)、M₂(x,y)和M₃(x,y)，每两幅干涉图之间的全局移相量是未知的，根据线性载频干涉图的定义，上述干涉图光强分布用下式表示：

其中，a(x,y)表示背景不均匀分布，b(x,y)表示调制度，本方法中认为这两项分布在三幅移相干涉图中保持一致，

表示相位分布，k₁和k₂分别为x和y方向的角波数，代表了线性载频的大小，

为待求的相位分布，假设M₁(x,y)的全局移相量为0，则Δ₁表示M₂(x,y)相对于M₁(x,y)的全局移相量，Δ₂表示M₃(x,y)相对于M₁(x,y)的全局移相量；

2)选取其中的M₁(x,y)和M₂(x,y)，根据该两幅图像的像素尺寸合理地选择N行和T列的光强数据，每组光强数据分别记作

及

其中，上标n＝1,2…N、t＝1,2,…T分别表示行数据和列数据的分组序数，下标1和2表示分别取自第一和第二幅干涉图；

3)计算每一组行数据

和每一组列数据

对应的评价函数，利用评价函数求解每组行数据和列数据对应的全局移相量dⁿ和d^t；

4)根据下式计算出M₁(x,y)和M₂(x,y)之间全局移相量的平均值d₁，也就是公式(2)中Δ₁的计算值，

5)同理，采用步骤2)至步骤4)相同的流程，计算出M₁(x,y)和M₃(x,y)之间全局移相量的平均值d₂，第一幅移相干涉图的全局移相量设定为d₀＝0；

6)根据d₀、d₁、d₂、M₁、M₂及M₃，利用最小二乘原则，并结合公式计算得到包含线性载频的压包相位分布r(x,y)；

7)对步骤6)中得到的r(x,y)进行相位展开操作，得到包含线性载频的解包相位分布r’(x,y)；

8)对步骤7)中得到的包含线性载频的解包相位进行二维的一次项拟合，并消除线性载频，从而得到不包含线性载频的真实相位分布

2.根据权利要求1所述的从全局移相量未知的线性载频移相干涉图中恢复出真实相位的方法，其特征在于：所述步骤3)中求解每组行数据和列数据对应的全局移相量dⁿ和d^t过程为：

①选择一组行数据

进行计算，初始计算时n＝1；

②根据步骤①计算出的行数据分别利用下列公式计算出I_m(x)和I_p(x)：

③针对步骤②中得到的I_m(x)使用一维希尔伯特变换计算得到I_mc(x)：

其中，H[.]表示取一维希尔伯特变换；

④依据步骤②中得到的I_p(x)和步骤③中得到的I_mc(x)，利用下式计算评价函数：

F(δ′)＝∫_l|I_p(x)-I_mc(x)·cot(δ′/2)|dx， (8)

其中，l为该组数据的采样点数，δ’＝δ₀+iε，i表示步数，移相量的初值δ₀和步长ε依经验确定；

⑤依据步骤④中得到的评价函数绘制出评价函数的曲线，并找到最小值点，该点代表的全局移相量即为从这组数据中求得的全局移相量dⁿ＝δ’；

⑥判断是否已经按照上述流程处理完所有选择出的行数据组，如果没有完成，则n＝n+1，回到步骤①，代入下一组数据，重复步骤①至步骤⑤的流程，直至确定出所有的dⁿ；

⑦同理，采用步骤①至步骤⑥相同的流程处理完所有选择出的列数据组，确定出所有列数据对应的全局移相量d^t。

3.根据权利要求2所述的从全局移相量未知的线性载频移相干涉图中恢复出真实相位的方法，其特征在于：所述公式(8)中步数i的取值范围为200～400。

4.根据权利要求1所述的从全局移相量未知的线性载频移相干涉图中恢复出真实相位的方法，其特征在于：所述步骤6)中压包相位分布r(x,y)的计算过程为：

①利用最小二乘原则，将前述所得的d₀、d₁、d₂、M₁、M₂及M₃代入最小二乘方程(9)，计算出D₁、D₂、D₃，

其中，i＝1,2,3；

②将步骤①计算得到的D₂、D₃，代入下式即可计算出包含线性载频的压包相位分布r(x,y)：

Images