基于反余弦图和连续性判据从单幅干涉图恢复相位的方法

1.本发明属于光学精密测试技术领域，具体涉及一种基于反余弦图和连续性判据从单幅干涉图中恢复出相位的方法。


背景技术：

2.光学工业除了在我们日常生活中，占有极高的市场地位，提高人们的生活水平之外。高端光学工业还具有显著性的极端制造特征，即在极端条件下，制造极端尺度或极高功能器件或功能系统的大科学工程。其中，大型望远镜、光束控制、高精度投影曝光系统已经成为现代大系统工程的典型代表。这些系统中的光学元件对光学表面质量、粗糙度等指标要求远远超出一般光学元件，它们基本代表了光学元件测量、加工中的极限水平要求。检测是加工的基础。
3.不管我们随处可见的照相镜头，还是高精度实验室装配出来的投影曝光系统，光学元件的检测都是一个非常重要的环节。没有高精度的检测手段，就没有高精度的确定性加工。面形精度是光学元件的核心技术指标，能否达到该指标依赖于测量所能达到的精度。目前，一般光学系统的面形精度要求达到几十纳米甚至更高，一般的测量技术已不能满足当前的需要。而基于相干原理的干涉仪能够轻易达到波长量级的测量精度，满足光学元件的面形的常规测量要求。
4.干涉测量在光学检测领域中应用广泛，用这种方法可以实现高精度，无接触，和全场测量。这种方法将有用的物理信息编码在条纹相位中，操作者并不能直接得到所需的信息，往往只能得到由余弦函数调制的干涉图。可见，解相技术的研究对干涉测量起着至关重要的作用。为了求出相位，人们提出了许多可行的方法。
5.时间移相干涉(tps)是一种成熟的高精度解相技术，自其诞生以来，就广泛应用于各种应用场景之中。相移干涉法如今已成为各种高精度干涉仪的标准设置。然而移相法要求至少三幅图像进行求解，难以实现实时检测。而且移相干涉法对机器内部和外部条件的要求都很苛刻，环境中的机械振动，空气扰动等都会影响测试精度。为了实现动态检测，使用单幅图求解相位是一种诱人的尝试。通过加入空间载频，结合傅里叶分析的傅里叶(ft)方法，可以实现单幅图解相位，然而载波干涉需要在两路干涉波前之间加入倾斜或者离焦，当被测波前畸变较大时，会引入较大像差，从而影响检测精度。而且ft方法因为吉布斯现象的存在会在边缘处产生震荡引入较大误差。除此之外还有其它方法，诸如窗口傅里叶法(wft)，空间载波移相法(scps)等方法，但这些方法都有自己的缺点。而且求解出来的相位需要进一步进行解包裹操作。且当干涉图中出现闭合条纹时，则难以进行正确解相。
6.为了恢复单幅图求解闭合条纹的相位，正则化相位跟踪法(rpt)是一种较为有效的方法，其求解结果连续，边缘处无吉布斯现象。但其也有不少缺点，其处理时间长，精度不高，在极值点和鞍点容易出错，且检测结果尚依赖路径。


技术实现要素：

7.本发明要解决的技术问题为：克服现有技术的不足，提出一种单幅图求解相位的方法，这种方法不需要加入载波，可以求解闭合条纹。针对rpt求解慢和精度不高的问题，直接采用反余弦图进行作为基础，通过连续性判据求出真值。针对rpt路径依赖的问题，提出在连续性判据中加入二阶导数连续的约束。整个过程容易理解，易于实施。
8.为了实现所述目的，本发明所采取的技术方案为：一种基于反余弦图和连续性判据从单幅干涉图恢复相位的方法，该方法步骤如下：
9.(a)：对干涉图进行预处理，去除直流项，进行归一化，去除高频噪声，干涉图i 的表达式可以表示为其中，a(x,y)为直流项，可以通过希尔伯特变换或频域滤波等方法滤除，b(x,y)为调制项，为待求解的相位，x为待测相位点(x,y)平面直角坐标系的横轴坐标，y为待测相位点(x,y)平面直角坐标系的竖轴坐标，归一化指的是将b(x,y)置为1，经过预处理后的干涉图i’的表达式为：
[0010][0011]
(b)：计算i’每个点的反余弦值，反余弦值是指cos-1
(v)＝-iln(v+i(1-v2)
1/2
)， v为i’中的单个点的值，i为令反余弦值取值范围为[0,π]；所有反余弦值组成的二维矩阵称为反余弦图，表示为i”；
[0012]
(c)：通过判定算法，采用连续性判据，遍历整幅图像，恢复出相位真值，连续性判据u(z)用公式(1)可以表示为：
[0013][0013][0014]
反余弦图中的每一个值可设为w，计算w+2π(k+l)和-w+2π(k+l)，其中，k为整数，且|k|＝0,1,2......f，f为阶梯常数，由用户自己设置，一般取2可以满足要求，l为周期值，为上一个恢复相位点所在的周期，即上一个点的估计值对2π取整，将计算w+ 2π(k+l)和-w+2π(k+l)的计算值整合成一个向量，顺序任意，称为可能性向量，z为可能性向量的某个元素，计算可能性向量中每一个元素的u值，使u最小的z值即可认为是真实相位值；n(x,y)为待测相位点(x,y)的邻域；为邻域中的已估计的值，为邻域中各点xx方向的二阶导数，为待测相位点(x,y)的xx方向的二阶导数，为邻域中各点xy方向的二阶导数，为待测相位点(x,y)上xy 方向的二阶导数，为邻域中各点yy方向的二阶导数，为待测相位点(x,y) 上yy方向的二阶导数，m(x,y)是这样一个标志场：如果点m(x,y)的相位已经估计过，则 m(x,y)的值为1，否则为0；β，β2分别为控制一阶和二阶梯度平滑度的参数；m2(x,y)是这样一个标志场：对于能直接计算即xx方向二阶导数的点，其值为1，否则为零0；m3(x,y)是这样一个标志场：对于能直接计算即xy方向二阶导数的点，其值为1，否则为零0；m4(x,y)是这样一个标志场：对于能直接计算即 yy方向二阶导数的点，其值为1，否则为零0；m2(ε,η)为标志场m2(x,y)在待测相位点 (x,y)周围邻域的取值，m3(ε,η)为标志场m3(x,y)在待测相位点(x,y)周围
邻域的取值， m4(ε,η)为标志场m4(x,y)在待测相位点(x,y)周围邻域的取值，如果m2(ε,η)不为0，则在z值取定后可以计算，也就是说，对于不同的z值其二阶导数值不同；对于第一个点，取其值为反余弦值，之后可以按行遍历整幅图像，直到所有点都被估计过。
[0015]
进一步地，整个步骤(c)的流程可以描述如下：
[0016]
c(1)：输入反余弦图，由反余弦图用户自行计算并输入标志场m2(x,y),m3(x,y)和 m4(x,y)；
[0017]
c(2)：制定遍历的路径，路径用于遍历反余弦图，起始点指的是最先遍历，即遍历开始的那个点，下一点指的是按照遍历路径接下来应该处理的点，上一点指的是按照遍历路径之前处理的那个点；
[0018]
c(3)：起始点取反余弦值；
[0019]
c(4)：移至路径中的下一点；
[0020]
c(5)：获得可能性向量，之后计算可能性向量中各元素的u值；
[0021]
c(6)：将使u值最小的元素作为估计的相位真值；
[0022]
c(7)：计算此元素的l值；
[0023]
c(8)：此元素m(x,y)置为1；
[0024]
c(9)：判断标志场m(x,y)所有值是否全都为1，若是，进行c(10)；若否，回到c(4)；
[0025]
c(10)：输出所有空间位置上估计的相位真值。
[0026]
进一步地，c(8)具体包括；此元素空间位置(x,y)所对应的标志场m(x,y)上相同位置的标志场的值置为1，标志场在遍历开始前为跟反余弦图相同大小的全零矩阵，遍历起始点后，起始点位置对应的标志场的值置为1。
[0027]
本发明的原理在于：
[0028]
从单幅干涉图恢复出闭合条纹的相位，是基于以下两条假设：
[0029]
1.重构的相位所产生的干涉图要与原始干涉图保持一致。
[0030]
2.重构的相位需要足够光滑。
[0031]
换句话说，相位的恢复过程就是寻求一个能满足上述两个约束条件的相位的过程。对于第一个条件，使用反余弦图即可以令其满足，用反余弦图生成的干涉图和原始干涉图肯定是一致的。但反余弦图的相位被包裹在半个周期之中，毫无疑问是不连续的。本发明通过上面所述的连续性判据来保证相位的连续性。
[0032]
首先要生成可能性向量，可能性向量中元素的公式如步骤c所述，其是基于这样一种考虑，连续相位相邻两点之间的差距不会太大。度量两点之间差距的就是阶梯常数f。在前一个估计点的基础上，即前一个所在周期的基础上，本估计点所有的取值可能性为 w+2π(k+l)和-w+2π(k+l)，其中k为整数，且|k|＝0,1,2......f。l是前一个估计点的值对2π取整得到的，设其值为w-1
，为了稳定性，可以采用邻域中估计过的点的均值。在所有可能性中，有一个使得连续性判据最小的值，该值即可认为是真值。
[0033]
连续性判据的原理在于，对于连续相位，估计点的值和其周遭邻域的值一定是接近的，但若只采用一阶判据，可以看作是用平面进行拟合，这样的判据是路径相关的，于是加入二阶判据。二阶判据认为其二阶导数也是连续的，这样就是用曲面来进行拟合，可以稳定的恢复出相位，可以按行遍历整幅图像，省去路径规划的工作，减轻工作量。
[0034]
本发明与现有技术的相比优点在于：
[0035]
(1)用反向验证的方式，将优化问题变成了整数优化问题。通过列举所有可能性进行连续性验证，找到最优解，可以降低计算量，提高求解速度。
[0036]
(2)精度更高，因为反余弦值将细节部分的信息精确的计算了出来，只需要对其进行符号和周期的确定即可。
[0037]
(3)易于理解，容易操作。
[0038]
(4)路径无关，可以按行遍历，简化了流程。
[0039]
(5)无需加载波。
[0040]
(6)可以求解闭合条纹。
[0041]
(7)边缘光滑，无吉布斯效应。
附图说明
[0042]
图1为本发明基于反余弦图和连续性判据从单幅干涉图恢复相位的方法中步骤(c) 算法流程图；
[0043]
图2为原始干涉图；
[0044]
图3为反余弦图；
[0045]
图4为连续性处理后的干涉图。
具体实施方式
[0046]
为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。
[0047]
对干涉图进行归一化处理后，计算其反余弦图，按照如图1所示的流程进行优化重构，具体实施步骤为：
[0048]
(a)：对干涉图进行预处理，去除直流项，进行归一化，去除高频噪声，干涉图i 的表达式可以表示为其中，a(x,y)为直流项，可以通过希尔伯特变换或频域滤波等方法滤除，b(x,y)为调制项，为待求解的相位，x为待测相位点(x,y)平面直角坐标系的横轴坐标，y为待测相位点(x,y)平面直角坐标系的竖轴坐标，归一化指的是将b(x,y)置为1，经过预处理后的干涉图i’的表达式为：
[0049][0050]
可以得到如图2所示的原始干涉图。
[0051]
(b)：计算i’每个点的反余弦值，反余弦值是指cos-1
(v)＝-iln(v+i(1-v2)
1/2
)， v为i’中的单个点的值，i为令反余弦值取值范围为[0,π]；所有反余弦值组成的二维矩阵称为反余弦图，表示为i”；反余弦图如图3所示。
[0052]
(c)：通过判定算法，采用连续性判据，遍历整幅图像，恢复出相位真值。对于第一个点，取其值为反余弦值，之后可以按行遍历整幅图像，直到所有点都被估计过。
[0053]
具体的，通过判定算法，采用连续性判据，遍历整幅图像，恢复出相位真值，连续性判据u(z)用公式(1)可以表示为：
[0054][0055]
反余弦图中的每一个值可设为w，计算w+2π(k+l)和-w+2π(k+l)，其中，k为整数，且|k|＝0,1,2......f，f为阶梯常数，由用户自己设置，一般取2可以满足要求，l为周期值，为上一个恢复相位点所在的周期，即上一个点的估计值对2π取整，将计算w+ 2π(k+l)和-w+2π(k+l)的计算值整合成一个向量，顺序任意，称为可能性向量，z为可能性向量的某个元素，计算可能性向量中每一个元素的u值，使u最小的z值即可认为是真实相位值；n(x,y)为待测相位点(x,y)的邻域；为邻域中的已估计的值，为邻域中各点xx方向的二阶导数，为待测相位点(x,y)的xx方向的二阶导数，为邻域中各点xy方向的二阶导数，为待测相位点(x,y)上xy 方向的二阶导数，为邻域中各点yy方向的二阶导数，为待测相位点(x,y) 上yy方向的二阶导数，m(x,y)是这样一个标志场：如果点m(x,y)的相位已经估计过，则 m(x,y)的值为1，否则为0；β，β2分别为控制一阶和二阶梯度平滑度的参数；m2(x,y)是这样一个标志场：对于能直接计算即xx方向二阶导数的点，其值为1，否则为零0；m3(x,y)是这样一个标志场：对于能直接计算即xy方向二阶导数的点，其值为1，否则为零0；m4(x,y)是这样一个标志场：对于能直接计算即 yy方向二阶导数的点，其值为1，否则为零0；m2(ε,η)为标志场m2(x,y)在待测相位点 (x,y)周围邻域的取值，m3(ε,η)为标志场m3(x,y)在待测相位点(x,y)周围邻域的取值， m4(ε,η)为标志场m4(x,y)在待测相位点(x,y)周围邻域的取值，如果m2(ε,η)不为0，则在z值取定后可以计算，也就是说，对于不同的z值其二阶导数值不同；对于第一个点，取其值为反余弦值，之后可以按行遍历整幅图像，直到所有点都被估计过。
[0056]
其中，整个步骤(c)的流程可以描述如下：
[0057]
c(1)：输入反余弦图，由反余弦图用户自行计算并输入标志场m2(x,y),m3(x,y)和 m4(x,y)；
[0058]
c(2)：制定遍历的路径，路径用于遍历反余弦图，起始点指的是最先遍历，即遍历开始的那个点，下一点指的是按照遍历路径接下来应该处理的点，上一点指的是按照遍历路径之前处理的那个点；
[0059]
c(3)：起始点取反余弦值；
[0060]
c(4)：移至路径中的下一点；
[0061]
c(5)：获得可能性向量，之后计算可能性向量中各元素的u值；
[0062]
c(6)：将使u值最小的元素作为估计的相位真值；
[0063]
c(7)：计算此元素的l值；
[0064]
c(8)：此元素m(x,y)置为1；
[0065]
其中，c(8)具体包括；此元素空间位置(x,y)所对应的标志场m(x,y)上相同位置的标志场的值置为1，标志场在遍历开始前为跟反余弦图相同大小的全零矩阵，遍历起始点后，起始点位置对应的标志场的值置为1。
[0066]
c(9)：判断标志场m(x,y)所有值是否全都为1，若是，进行c(10)；若否，回到c(4)；
[0067]
c(10)：输出所有空间位置上估计的相位真值。
[0068]
恢复相位的干涉图正如图4所示。