技术领域
[0001] 本
发明属于
信号处理技术改进领域，尤其涉及一种基于两步非线性移相的快速高精度相位提取方法。
背景技术
[0002] 精密零器件的数量级在微米甚至
纳米级别，在精密零器件的生产加工中，不可避免的存在各种误差，一些微小的形变将会导致零器件的功能受到影响，因此对零器件微观形貌的检测十分重要。形貌的检测可以分为
接触测量和非接触测量。接触测量使用机械探针，探针与物体接触会损坏物体，探针也需要经常更换，增加了成本。非接触测量具有效率高、
稳定性好、避免对零器件二次损伤的优点，受到人们的广泛关注。为了获得相位信息，好的相位提取
算法至关重要。经典的相位提取算法需要知道确定的移相值，然而在移相中，由于
光源不稳定、空气扰动以及
移相器的不稳定所引起的移相误差，会导致经典相位提取算法精度急剧下降。为了提高移相精度，克服平移误差(移相器未对准所导致)的影响，涌现了一大批相位提取算法，包括主成分分析法(PCA)，李萨如图形与椭圆拟合(LEF)，斯密特
正交与椭圆拟合(GS&LEF)，斯密特正交与最小平方
迭代(GS&LSI)，主成分分析与椭圆拟合(PCA&LEF)，椭圆拟合与最小平方迭代(LEF&LSI)等等。然而移相器的非线性以及非均匀性所导致的倾斜移相误差尚未克服，影响算法的精度。
[0003] 为了克服倾斜移相误差，研究者们提出了大量相位提取算法。2000年，陈等利用泰勒展开的方法来补偿随机倾斜移相误差，该算法只能处理小幅度的倾斜
相移。2008年，徐等提出利用改进的最小平方迭代算法提取倾斜移相值，再利用移相值推算出相位信息，此算法需要进行迭代，导致算法运行时间长。为了扩大使用范围，刘等提出利用三步迭代算法，可以提取大幅度倾斜移相，然而，这个算法耗时长，不利于动态测量。近期，WIELGUS用非线性误差最小化从两副图中提取倾斜移相值，可以获得较快的速度，然而这种算法需要提前滤除背景光强，这会给算法引入新的误差。而且这种算法只针对线性倾斜移相有效，对于非线性倾斜移相，此算法不能有效的提取出相位信息。
[0004] 虽然在人们提出的移相干涉相位提取算法中可以在倾斜移相的情况下提取出相位，但这些方法都各自存在不足，无法满足快速、高精度的从倾斜移相干涉图中精确的提取出相位。
发明内容
[0005] 本发明的目的在于提供一种基于两步非线性移相的快速高精度相位提取方法，旨在解决由移相器的非线性以及非均匀性引起的随机非线性移相误差、背景光强和调制幅度扰动等影响的问题。
[0006] 本发明是这样实现的，一种基于两步非线性移相的快速高精度相位提取方法，所述基于两步非线性移相的快速高精度相位提取方法包括以下步骤：
[0007] S1、利用两幅干涉图拟合出一个斜椭圆方程；
[0008] S2、利用最小二乘法求解斜椭圆方程的系数；
[0009] S3、根据求解的斜椭圆方程系数求解斜椭圆的参数；
[0010] S4、根据获取的斜椭圆参数将斜椭圆转换成标准的椭圆方程；
[0011] S5、通过标准的椭圆方程提取出相位信息。
[0012] 本发明的进一步技术方案是：所述步骤S1中两幅干涉图在存在背景光强、调制幅度扰动及随机非线性移相误差下表示为 其中，a1(x,y),a2(x,y)表示背景光强，b1(x,y),b2(x,y)表示调制幅度，δ(x,y)表示两幅干涉图之间的移相量。
[0013] 本发明的进一步技术方案是：所述步骤S1中在sin2φ+cos2φ＝1斜椭圆方程为I1表示第一幅移相干涉图，I2表示第二幅干移相涉图。
[0014] 本发明的进一步技术方案是：所述步骤S2中还包括以下步骤：
[0015] S21、对斜椭圆方程：
[0016] 进行化简得到a·I12+b·I1·I2+c·I22+d·I1+f·I2+g＝0；
[0017] S22、通过最小二乘法得到斜椭圆方程的系数
[0018] 其中，a表示I1的二次项系数，b表示I1,I2乘积项的系数，c表示I2的二次项系数，d表示I1的一次项系数，f表示I2的一次项系数，g表示常数项系数。
[0019] 本发明的进一步技术方案是：所述步骤S3中还包括如下步骤：
[0020] S31、根据椭圆和斜椭圆的关系得到方程
[0021]
[0022] S31、根据步骤S22中求得的斜椭圆系数通过斜椭圆由标准椭圆关系方程求得椭圆参数
[0023]
[0024]
[0025]
[0026] 其中，ax代表椭圆长半轴，ay代表椭圆短半轴，θ为旋转
角，x0，y0为平移量。
[0027] 本发明的进一步技术方案是：所述步骤S4中利用步骤S32中获取的椭圆参数将斜椭圆转换为一个标准椭圆方程
[0028] 本发明的进一步技术方案是：所述步骤S5中根据获取的标准椭圆方程，利用方程获取相位信息，其中，r＝ay/ax表示椭圆长半轴与短半轴的比值，X表示校正后的第一幅移相干涉图，Y表示校正后的第二幅移相干涉图。
[0029] 本发明的有益效果是：通过斜椭圆转换为标准椭圆的过程，校正了由光源不稳定、空气扰动以及移相器不稳定所带来的移相误差，提高算法的精度。不需要提前滤除
直流分量，利用非迭代的方法直接从干涉图中提取相位信息，提高算法的运行速度，可以实现快速的测量。
附图说明
[0030] 图1是本发明
实施例提供的基于两步非线性移相的快速高精度相位提取方法的流程示意图。
[0031] 图2是本发明实施例提供的计算机模拟存在扰动的相位提取结果图；
[0032] 图3是本发明实施例提供的计算机模拟干涉条纹数小于1的相位提取结果图。
[0033] 图4是本发明实施例提供的计算机模拟不同噪声下的相位提取结果图。
具体实施方式
[0034] 如图1所示，本发明提供的基于两步非线性移相的快速高精度相位提取方法，其详述如下：
[0035] 步骤S0、采集两幅干涉图，光源采用经过
准直后平行LED光源，利用单轴平晶将一束光分为有一定剪切量同方向的两束光，利用1/4波片以及检偏器进行随机移相，利用CCD工业相机采集两幅移相干涉图。
[0036] 步骤S1，利用两幅干涉图拟合出一个斜椭圆方程；拟合出一个斜椭圆。当存在背景光强和调制幅度扰动，并且存在随机非线性移相误差的情况下，两幅干涉图可以表示为：
[0037]
[0038] a1(x,y),a2(x,y)表示背景光强，b1(x,y),b2(x,y)表示调制幅度。δ(x,y)表示两幅干涉图之间的移相量。
[0039] 在接下来的推导中，我们省略(x,y)，对公式(1)进行化简运算，
[0040]
[0041] 得到以下结果：
[0042]
[0043] 因为sin2φ+cos2φ＝1，
[0044]
[0045] 步骤S2、利用最小二乘法求解斜椭圆方程的系数；观察公式(4)，用a来表示I1的二次项系数，用b来表示I1,I2乘积项的系数，用c来表示I2的二次项系数，用d来表示I1的一次项系数，用f来表示I2的一次项系数，用g来表示常数项系数我们可以得到：
[0046] a·I12+b·I1·I2+c·I22+d·I1+f·I2+g＝0 (5)
[0047] 其中，
[0048] 这些参数可以利用最小二乘法得到。对斜椭圆方程：进行化简得到a·
I12+b·I1·I2+c·I22+d·I1+f·I2+g＝0；通过最小二乘法得到斜椭圆方程的系数[0049] 其中，a表示I1的二次项系数，b表示I2的二次项系数，c表示I1,I2乘积项的系数，d表示I1的一次项系数，f表示I2的一次项系数，g表示常数项系数。
[0050] 令M＝[a,b,c,d,f,g]， 目标就是
[0051]
[0052] 其中，
[0053] 构造拉格朗日函数
[0054] L(M,λ)＝MTNNTM-λ(MTHM-1) (7)
[0055] 对其求导为零，可以得到
[0056]
[0057] 令S＝NNT,则SM＝λHM,通过求解广义
特征向量可以得到6个可能备选的M。利用MTHM＝1这个条件来筛选合格的M。
[0058] 接下来说明筛选条件。假设得到的特征值和特征向量对{λi,vi}，对于椭圆方程a·I12+b·I1·I2+c·I22+d·I1+f·I2+g＝0判别式为：
[0059]
[0060] 当Λ≠0，且ac-b2>0时为椭圆。
[0061] 步骤S3，根据求解的斜椭圆方程系数求解斜椭圆的参数；一个标准的椭圆参数方程如下所示
[0062]
[0063] 其中ax代表椭圆长半轴，ay代表椭圆短半轴。
[0064] 斜椭圆由标准椭圆旋转和平移得到，第一，只考虑旋转，旋转后的新椭圆可以由下式表示：
[0065]
[0066] 其中，x',y'表示经过旋转后的新椭圆参数，θ表示旋转角。
[0067] 第二，只考虑平移，平移后的新椭圆可以由下式表示：
[0068]
[0069] 其中，x”,y”表示经过平移后的新椭圆参数，x0，y0为平移量。
[0070] 因此，标准椭圆经过旋转和平移后的方程式可以表示为：
[0071]
[0072] 其中，x,y表示经过旋转和平移后的斜椭圆。
[0073] 用旋转平移后的x，y来表示标准椭圆X，Y，将斜椭圆转换成标准椭圆的推导过程如下：
[0074]
[0075] 将公式(14)带入公式 可以得到公式(15)
[0076]
[0077] 斜椭圆转换为标准椭圆的关键就是求解出椭圆参数ax，ay，θ，x0，y0；
[0078] 将公式(15)化简成公式(5)的形式，并将公式(15)的x用公式(5)中的I1代替，将公式(15)的y用公式(5)中的I2代替。对比公式(5)与公式(15)，则公式(5)的系数a，b，c，d，f，g可以表示为：
[0079]
[0080]
[0081]
[0082] 观察公式(16)，用a，b，c来表示d，f，g可以得到 利用这个方程可以求出x0,y0,ax,ay；从a·I12+b·I1·I2+c·I22+d·I1+f·I2+g＝0这个方程，可以写出倾角θ的表达式，因此我们可以得到椭圆的参数表达式如下：
[0083]
[0084]
[0085]
[0086] 步骤S4，根据获取的斜椭圆参数将斜椭圆转换成标准的椭圆方程；利用这些椭圆参数，斜椭圆是由标准椭圆旋转和平移得到，我们求出旋转角θ以及平移量x0,y0后，将公式(14)中的x用I1表示，y用I2表示，因此利用公式(14)我们可以把斜椭圆转换为一个标准椭圆，形式如下：
[0087]
[0088] 步骤S5，通过标准的椭圆方程提取出相位信息。
[0089] 标准椭圆方程参数的形式为： φ就是我们要求解的相位信息。
[0090]
[0091] 其中，
[0092] 因此利用标准椭圆方程，我们可以准确的提取出相位信息
[0093]
[0094] 为了验证该方法的有效性，用计算机对各种情况进行模拟仿真。
[0095] 案例1：在移相干涉测量中，存在光源不稳定、外界环境的扰动，这些会引起移相误差。因此将两幅干涉图的参数设置如下所示：
[0096] 背景光强：a1(x,y)＝1+0.1·exp(-0.02(x2+y2))，a2(x,y)＝0.9+0.21·exp(-0.02(x2+y2))
[0097] 调制幅度：b1(x,y)＝1+0.1·exp(-0.02(x2+y2))，b2(x,y)＝0.9+0.21·exp(-0.02(x2+y2))
[0098] 物体相位：φ(x,y)＝5·π·(x2+y2)
[0099] 两幅干涉图之间的移相差为：δ＝0.75+0.05·x2+0.15·y
[0100] 数值模拟生成512*512
像素的初始相位，如图2(a)所示，图2(b)是两幅随机干涉图中的一幅，图2(c)是两幅干涉图之间存在非线性误差的移相量，图2(d)是通过该发明提供的方法提取出相位，再经过相位解包裹运算计算,并进行倾斜调整后得到的物体表面相位信息。
[0101] 相位提取时间为：0.046847s。
[0102] 案例2：在移相干涉测量中，干涉条纹往往会小于1。然而在现有的很多相位提取算法中都要求干涉条纹数大于1才能有效提取相位。为了验证本发明在该情况的精确性，我们对干涉条纹数小于1的情况进行模拟仿真。
[0103] 背景光强：a1(x,y)＝1+0.1·exp(-0.02(x2+y2))，a2(x,y)＝0.9+0.21·exp(-0.02(x2+y2))
[0104] 调制幅度：b1(x,y)＝1+0.1·exp(-0.02(x2+y2))，b2(x,y)＝0.9+0.21·exp(-0.02(x2+y2))
[0105] 物体相位：φ(x,y)＝0.6·π·(x2+y2)
[0106] 两幅干涉图之间的移相差为：δ＝0.75+0.05·x2+0.15·y
[0107] 数值模拟生成512*512像素的初始相位，如图3(a)所示，图3(b)是两幅随机干涉图中的一幅，图3(c)是两幅干涉图之间存在非线性误差的移相量，图3(d)是通过该发明提供的方法提取出相位，再经过相位解包裹运算计算,并进行倾斜调整后得到的物体表面相位信息。
[0108] 相位提取时间为：0.044079s
[0109] 案例3：在移相干涉测量中，噪声不可避免。为了验证本发明的抗噪性，我们在不同
信噪比(SNR)进行仿真实验。
[0110] 背景光强：a1(x,y)＝1+0.1·exp(-0.02(x2+y2))，a2(x,y)＝0.9+0.21·exp(-0.02(x2+y2))
[0111] 调制幅度：b1(x,y)＝1+0.1·exp(-0.02(x2+y2))，b2(x,y)＝0.9+0.21·exp(-0.02(x2+y2))
[0112] 物体相位：φ(x,y)＝5·π·(x2+y2)
[0113] 两幅干涉图之间的移相差为：δ＝0.75+0.05·x2+0.15·y
[0114] 此外加入不同信噪比的高斯白噪声。SNR分别为：13dB，15dB,17dB,20dB,30dB。
[0115] 数值模拟生成512*512像素的初始相位，如图4(a)所示。图4(b)-4(f)是通过该发明提供的方法提取出相位，再经过相位解包裹运算计算出的物体相位信息，分别对应13dB，15dB,17dB,20dB,30dB信噪比下的相位提取结果。
[0116] 相位提取时间如表一所示：
[0117] 表一
[0118]SNR(dB) 13 15 17 20 30
Time(s) 0.053104 0.043537 0.045660 0.044094 0.043574
[0119] 案例1、案例2和案例3说明了该方法在不同条件下都能够快速、高精度的提取出干涉图的相位信息，具有较强的鲁棒性。
[0120] 通过斜椭圆转换为标准椭圆的这个过程，校正了由光源不稳定、空气扰动以及移相器不稳定所带来的移相误差，提高算法的精度。不需要提前滤除直流分量，利用非迭代的方法直接从干涉图中提取相位信息，提高算法的运行速度，可以实现快速的测量。
[0121] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何
修改、等同替换和改进等，均应包含在本发明的保护范围之内。
