CN112066909B - 一种基于倾斜平面高精度提取的抗振动干涉测量方法 - Google Patents

Abstract

本发明公开了一种基于倾斜平面高精度提取的抗振动干涉测量方法，包括以下步骤：对振动环境下采集的干涉图进行Fourier变换，并在频域中提取+1级峰值坐标；以步骤1得到的峰值坐标为中心，细化其周围频谱数据，获得新的峰值坐标，重复该过程直至获得亚像素精度峰值坐标；利用所述亚像素精度峰值坐标求解倾斜相位平面；利用最小二乘法求解最终相位分布信息。本发明无需对干涉仪进行硬件改动，具有较强的实用性和适用性，有望为振动环境下的移相干涉测量提供一种低成本、高精度的解决方案。

Description

一种基于倾斜平面高精度提取的抗振动干涉测量方法

技术领域

本发明属于光学检测领域，特别涉及一种基于倾斜平面高精度提取的抗振动干涉测量方法。

背景技术

移相干涉测量技术PSI(Phase-shifting interferometry)由于具有高精度的优点，因此被广泛应用于光学元件面形测量领域。为保证PSI的测量精度，在移相过程中必须确保移相量为预设的常量。然而在实际测量中，外部环境和机电结构的振动会导致移相量产生随机误差，造成同一帧干涉图内各像素点处的移相值不一样，使得测量结果精度下降。

为了解决振动带来的问题，国内外学者提出了一些处理倾斜误差的相位提取算法，主要分为迭代法、空间载波法和时域信号分析法。迭代法效率较低，且一般需要对初始移相量进行预估，对于倾斜移相的非线性解耦合，现有算法存在误差较大或优化速度慢等不足；空间载波法基于空间载波干涉图检测移相量，然而引入较大的空间载波将引起一定的系统误差，测量前需对其系统误差进行校正，在复杂波面情况下可能会导致算法失败；时域信号分析法需要大量干涉图，其抗振性能取决于移相速度，对硬件要求较高，且测量时间较长可能会引入新的误差源。

发明内容

本发明的目的在于提供一种基于倾斜平面高精度提取的抗振动干涉测量方法。

实现本发明目的的技术解决方案为：一种基于倾斜平面高精度提取的抗振动干涉测量方法，所述方法包括以下步骤：

步骤1，对振动环境下采集的干涉图进行Fourier变换，并在频域中提取+1级峰值坐标；

步骤2，以步骤1得到的峰值坐标为中心，细化其周围频谱数据，获得新的峰值坐标，重复该过程直至获得亚像素精度峰值坐标；

步骤3，利用所述亚像素精度峰值坐标求解倾斜相位平面；

步骤4，利用最小二乘法求解最终相位分布信息。

进一步地，步骤1所述对振动环境下采集的干涉图进行Fourier变换，并在频域中提取+1级峰值坐标，具体包括：

步骤1-1，对于振动环境下采集的干涉图，其光强表达式为：

I_n(x,y)＝a(x,y)+b(x,y)cos(φ(x,y)+p_n(x,y)),n＝1,2,...N

式中，a(x,y)为背景强度，b(x,y)为调制幅度，φ(x,y)为待测相位，p_n(x,y)为振动倾斜相位平面；

对干涉图I_n(x,y)进行Fourier变换获得频谱F(u,v)：

式中，A(u,v)为干涉图背景光强频谱，

为-1级旁瓣，

为+1级旁瓣，C为c(x,y)的Fourier变换，

φ_n(x,y)表示第n帧干涉图的待测相位，C^*为C的共轭，f_x和f_y为频谱正负一级峰值坐标，k_n为第n幅干涉图载频系数；

步骤1-2，从F(u,v)中滤除干涉图背景光强频谱A(u,v)，并取正频率部分获得频谱F′(u,v)，提取F′(u,v)的峰值坐标(f_x,f_y)：

进一步地，步骤2所述以步骤1得到的峰值坐标为中心，细化其周围频谱数据，获得新的峰值坐标，重复该过程直至获得亚像素精度峰值坐标，具体过程包括：

步骤2-1，令r＝1，进行第一次迭代；

步骤2-2，对干涉图I(X,Y)进行补零延拓，以使整个矩阵的横向和纵向长度都变为原来的K倍，获得新的矩阵I(X′,Y′)：

式中，X′＝[0,1,...,KM-1]^T，Y′＝[0,1,...KN-1]^T，M、N分别为干涉图在x、y方向的像素数；

步骤2-3，对I(X′,Y′)中的([f_x-0.5,f_x+0.5]，[f_y-0.5,f_y+0.5])范围作离散余弦变换DCT，获得频谱F(U_r,V_r)；

步骤2-4，从频谱F(U_r,V_r)中提取+1级旁瓣F′(U_r,V_r)，并提取新的峰值坐标

对

除以补零倍数K，获得具有K倍精度的峰值坐标：

步骤2-5，判断

的精度是否达到亚像素，若没达到，则令r＝r+1，返回执行步骤2-2，直至满足条件，否则输出

进一步地，步骤3所述利用所述亚像素精度峰值坐标求解倾斜相位平面，具体包括：

步骤3-1，假设

和k_n分别为第n幅干涉图的x方向载频系数、y方向载频系数和相对载频系数，利用所述

求解

和k_n，计算公式为：

将

代入步骤1-1的F(u,v)，求得k_n：

式中，Arg(·)为计算相位角函数；

步骤3-2，定义振动倾斜相位平面p_n(x,y)为：

p_n(x,y)＝α_nx+β_ny+γ_n

式中，α_n、β_n分别为倾斜相位平面在x方向和y方向的倾斜系数，γ_n为平动移相量；

根据所述干涉图载频系数计算α_n、β_n和γ_n，计算公式为：

式中，

和k₁为第1幅干涉图载频系数；

由此获得振动倾斜相位平面p_n(x,y)。

进一步地，步骤4所述利用最小二乘法求解最终相位分布信息，具体包括：

步骤4-1，令c₁＝a，c₂＝bcosφ，c₃＝-bsinφ，同时省略(x,y)，将步骤1-1光强公式改写为：

I_n＝c₁+c₂ cos p_n+c₃ sin p_n

步骤4-2，令实际测量得到的干涉图为I′_n，引入残差平方和函数：

式中，N表示参与计算的干涉图帧数；

步骤4-3，求解c₁、c₂、c₃：

令：

求解得C＝A^-1B，其中：

C＝[c₁ c₂ c₃]^T

由上述最小二乘法求解出系数c₁、c₂、c₃，可得待测相位结果为φ＝tan^-1(-c₃/c₂)。

本发明与现有技术相比，其显著优点为：1)可用于求解具有倾斜移相且条纹数较少的干涉图相位分布；干涉图条纹数较少(一般为三到四根)，基本可以忽略由于引入载频导致的系统回程误差；2)无需对干涉仪做硬件改动，测量成本较低且方法简单易行，具有很强的实用性和适用性。

下面结合附图对本发明作进一步详细描述。

附图说明

图1为本发明基于倾斜平面高精度提取的抗振动干涉测量方法流程图。

图2为局部频谱数据细化过程示意图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

在一个实施例中，结合图1，提供了一种基于倾斜平面高精度提取的抗振动干涉测量方法，所述方法包括以下步骤：

步骤1，对振动环境下采集的干涉图进行Fourier变换，并在频域中提取+1级峰值坐标；

步骤2，以步骤1得到的峰值坐标为中心，细化其周围频谱数据，获得新的峰值坐标，重复该过程直至获得亚像素精度峰值坐标；

步骤3，利用所述亚像素精度峰值坐标求解倾斜相位平面；

步骤4，利用最小二乘法求解最终相位分布信息。

进一步地，在其中一个实施例中，步骤1所述对振动环境下采集的干涉图进行Fourier变换，并在频域中提取+1级峰值坐标，具体包括：

步骤1-1，对于振动环境下采集的干涉图，其光强表达式为：

I_n(x,y)＝a(x,y)+b(x,y)cos(φ(x,y)+p_n(x,y)),n＝1,2,...N

式中，a(x,y)为背景强度，b(x,y)为调制幅度，φ(x,y)为待测相位，p_n(x,y)为振动倾斜相位平面；

将上式改写为复数形式：

式中，

和k_n为干涉图的载频系数；

对干涉图I_n(x,y)进行Fourier变换获得频谱F(u,v)：

式中，A(u,v)为干涉图背景光强频谱，

为-1级旁瓣，

为+1级旁瓣，C为c(x,y)的Fourier变换，

φ_n(x,y)表示第n帧干涉图的待测相位，C^*为C的共轭，f_x和f_y为频谱正负一级峰值坐标，k_n为第n幅干涉图载频系数；

步骤1-2，从F(u,v)中滤除干涉图背景光强频谱A(u,v)，并取正频率部分获得频谱F′(u,v)，提取F′(u,v)的峰值坐标(f_x,f_y)：

进一步地，在其中一个实施例中，步骤2所述以步骤1得到的峰值坐标为中心，细化其周围频谱数据，获得新的峰值坐标，重复该过程直至获得亚像素精度峰值坐标，具体过程包括：

步骤2-1，令r＝1，进行第一次迭代；

步骤2-2，由于离散傅里叶变换(DCT)的性质，步骤1所得峰值坐标(f_x,f_y)为整数精度。为了提高(f_x,f_y)的定位精度，对干涉图I(X,Y)进行补零延拓，以使整个矩阵的横向和纵向长度都变为原来的K倍，获得新的矩阵I(X′,Y′)：

式中，X′＝[0,1,...,KM-1]^T，Y′＝[0,1,...KN-1]^T，M、N分别为干涉图在x、y方向的像素数；

补零延拓后的频谱为：

其中，U′＝[0,1,...,KM-1]^T，V′＝[0,1,...,KN-1]^T。经过补零后的频谱中峰值位置为(f′_x,f′_y)，根据离散Fourier变换的性质，由于整个矩阵的横向和纵向长度都变为原来的K倍，此时(f′_x,f′_y)达到了1/K的像素精度，有：

(f′_x,f′_y)＝(f_x,f_y)/K

步骤2-3，为保证计算精度，延拓倍数K要足够大。然而当K取值较大时，会使得计算量极大增加。由于补零后的I(X′,Y′)中除了M×N个点外其余部分均为0，为减少计算量，在Fourier计算时可忽略(M-1)N个零值点。根据矩阵分块计算理论，上述补零延拓后的频谱F(U′,V′)可以简化为：

步骤2-4，对I(X′,Y′)中的([f_x-0.5,f_x+0.5]，[f_y-0.5,f_y+0.5])范围作离散余弦变换DCT，获得频谱F(U_r,V_r)(因为+1级旁瓣峰值位置(f_x,f_y)满足

和

为进一步减少计算量，只需要计算这部分频谱数据F(U_r,V_r)即可)；

令：

将这两个公式代入F(U′,V′)中，获得：

通过矩阵Fourier变换的性质省去了大量冗余数据的计算，极大提高了计算速度。快速频谱细化的具体过程如图2所示，可以看出最终计算三个大小为(K+1)×M，M×N和N×(K+1)的矩阵连乘，计算量大大降低；

步骤2-5，从频谱F(U_r,V_r)中提取+1级旁瓣F′(U_r,V_r)，并提取新的峰值坐标

对

除以补零倍数K，获得具有K倍精度的峰值坐标：

步骤2-6，判断

的精度是否达到亚像素，若没达到，则令r＝r+1，返回执行步骤2-2，直至满足条件，否则输出

这里，判断条件也可以设置为迭代次数是否达到预设次数阈值。

进一步地，在其中一个实施例中，步骤3所述利用所述亚像素精度峰值坐标求解倾斜相位平面，具体包括：

步骤3-1，假设

和k_n分别为第n幅干涉图的x方向载频系数、y方向载频系数和相对载频系数，利用所述

求解

和k_n，计算公式为：

将

代入步骤1-1的F(u,v)，求得k_n：

式中，Arg(·)为计算相位角函数；

步骤3-2，定义振动倾斜相位平面p_n(x,y)为：

p_n(x,y)＝α_nx+β_ny+γ_n

式中，α_n、β_n分别为倾斜相位平面在x方向和y方向的倾斜系数，γ_n为平动移相量；

根据所述干涉图载频系数计算α_n、β_n和γ_n，计算公式为：

式中，

和k₁为第1幅干涉图载频系数；

由此获得振动倾斜相位平面p_n(x,y)。

进一步地，在其中一个实施例中，步骤4所述利用最小二乘法求解最终相位分布信息，具体包括：

步骤4-1，令c₁＝a，c₂＝bcosφ，c₃＝-bsinφ，同时省略(x,y)，将步骤1-1光强公式改写为：

I_n＝c₁+c₂ cos p_n+c₃ sin p_n

步骤4-2，令实际测量得到的干涉图为I′_n，引入残差平方和函数：

式中，N表示参与计算的干涉图帧数；

步骤4-3，求解c₁、c₂、c₃：

令：

求解得C＝A^-1B，其中：

C＝[c₁ c₂ c₃]^T

由上述最小二乘法求解出系数c₁、c₂、c₃，可得待测相位结果为φ＝tan^-1(-c₃/c₂)。

本发明利用Fourier变换将振动环境下采集的干涉图变换到频域，再采用频域细分操作对峰值坐标进行亚像素精度定位从而求解出振动倾斜平面，最后利用最小二乘求出待测相位分布。不需要对初始移相量进行预估，也无需考虑背景光强和调制度是否一致，操作方便简单，具有很强的实用性和适用性；此外，通过迭代矩阵Fourier变换实现对峰值局部数据的计算，快速实现0.001像素的定位精度，极大保证了振动倾斜平面的计算精度。总的来说，本发明无需硬件改动，且计算精度较高，有望为动态环境下的移相干涉测量提供一种低成本、高精度的解决方案。

以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (2)

1.一种基于倾斜平面高精度提取的抗振动干涉测量方法，其特征在于，所述方法包括以下步骤：

步骤1，对振动环境下采集的干涉图进行Fourier变换，并在频域中提取+1级峰值坐标；具体包括：

步骤1-1，对于振动环境下采集的干涉图，其光强表达式为：

I_n(x,y)＝a(x,y)+b(x,y)cos(φ(x,y)+p_n(x,y)),n＝1,2,...N

式中，a(x,y)为背景强度，b(x,y)为调制幅度，φ(x,y)为待测相位，p_n(x,y)为振动倾斜相位平面；

对干涉图的光强表达式I_n(x,y)进行Fourier变换获得频谱F(u,v)：

式中，A(u,v)为干涉图背景光强频谱，

为-1级旁瓣，

为+1级旁瓣，C为c(x,y)的Fourier变换，

φ_n(x,y)表示第n帧干涉图的待测相位，C^*为C的共轭，f_x和f_y为频谱正负一级峰值坐标，k_n为第n幅干涉图载频系数；

步骤1-2，从F(u,v)中滤除干涉图背景光强频谱A(u,v)，并取正频率部分获得频谱F′(u,v)，提取F′(u,v)的峰值坐标(f_x,f_y)：

步骤2，以步骤1得到的峰值坐标为中心，细化其周围频谱数据，获得新的峰值坐标，重复该过程直至获得亚像素精度峰值坐标；具体过程包括：

步骤2-1，令r＝1，进行第一次迭代；

步骤2-2，对干涉图I(X,Y)进行补零延拓，以使整个矩阵的横向和纵向长度都变为原来的K倍，获得新的矩阵I(X′,Y′)：

式中，X′＝[0,1,...,KM-1]^T，Y′＝[0,1,...KN-1]^T，M、N分别为干涉图在x、y方向的像素数；K＝10^r，r为正整数；

步骤2-3，对I(X′,Y′)中的([f_x-0.5,f_x+0.5]，[f_y-0.5,f_y+0.5])范围作离散余弦变换DCT，获得频谱F(U_r,V_r)；

步骤2-4，从频谱F(U_r,V_r)中提取+1级旁瓣F′(U_r,V_r)，并提取新的峰值坐标

对

除以补零倍数K，获得具有K倍精度的峰值坐标：

步骤2-5，判断

的精度是否达到亚像素，若没达到，则令r＝r+1，返回执行步骤2-2，直至满足条件，否则输出

步骤3，利用所述亚像素精度峰值坐标求解倾斜相位平面；具体包括：

步骤3-1，假设

和k_n分别为第n幅干涉图的x方向载频系数、y方向载频系数和相对载频系数，利用所述

求解

和k_n，计算公式为：

将

代入步骤1-1的F(u,v)，求得k_n：

式中，Arg(·)为计算相位角函数；

步骤3-2，定义振动倾斜相位平面p_n(x,y)为：

p_n(x,y)＝α_nx+β_ny+γ_n

式中，α_n、β_n分别为倾斜相位平面在x方向和y方向的倾斜系数，γ_n为平动移相量；

根据所述干涉图载频系数计算α_n、β_n和γ_n，计算公式为：

式中，

和k₁为第1幅干涉图载频系数；

由此获得振动倾斜相位平面p_n(x,y)；

步骤4，利用最小二乘法求解最终相位分布信息。

2.根据权利要求1所述的基于倾斜平面高精度提取的抗振动干涉测量方法，其特征在于，步骤4所述利用最小二乘法求解最终相位分布信息，具体包括：

步骤4-1，令c₁＝a，c₂＝bcosφ，c₃＝-bsinφ，同时省略(x,y)，将步骤1-1光强公式改写为：

I_n＝c₁+c₂cosp_n+c₃sinp_n

步骤4-2，令实际测量得到的干涉图为I_n′，引入残差平方和函数：

式中，N表示参与计算的干涉图帧数；

步骤4-3，求解c₁、c₂、c₃：

令：

求解得C＝A^-1B，其中：

C＝[c₁ c₂ c₃]^T

由上述最小二乘法求解出系数c₁、c₂、c₃，可得待测相位结果为φ＝tan^-1(-c₃/c₂)。

Images