CN112179269B - 一种基于Zernike多项式与WKF的相位解包裹方法 - Google Patents

Abstract

本发明公开了一种基于Zernike多项式与WKF的相位解包裹方法，通过Zernike多项式拟合的方法建立包裹相位的空间动态模型，并用WKF计算模型中待定的Zernike多项式系数。为了抑制大噪声点对WKF的影响，本发明通过计算包裹相位质量图的方法将高噪声点从后续的WKF计算中剔除，进一步提升稳定性，以应对WKF发散的极端情形。能够在高噪声环境下具有高的精度和稳定性。

Description

一种基于Zernike多项式与WKF的相位解包裹方法

技术领域

本发明属于干涉计量技术领域，具体涉及一种Zernike多项式拟合和包裹卡尔曼滤波算法(Wrapped Kalman Filter，WKF)的相位解包裹方法。

背景技术

干涉条纹图的相位恢复方法是保证干涉计量精度的关键技术，而相位解包裹又是大部分相位恢复方法的关键步骤。绝大部分的相位恢复方法都会引入反正切函数，使得求出的相位被包裹在(-π,π]之间，相位解包裹就是把求出的包裹相位进一步地还原成没有区间限制的真实相位。

目前，国内外对相位解包裹方法的研究主要可以分为两类：路径跟踪算法和路径无关算法。路径跟踪算法的特点是要满足积分结果与路径无关的条件，其又可分为枝切法、质量图导向法和掩膜割线法等。而路径无关算法的任务是求满足最小范数的解，是一种全局的相位解包裹算法，按照求解方式的不同又可分为最小二乘法、光强传输方程法和多项式拟合法。

然而，上述算法在大噪声下的鲁棒性都比较差。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种Zernike多项式拟合和包裹卡尔曼滤波算法(Wrapped Kalman Filter，WKF)的相位解包裹方法，把Zernike多项式拟合与WKF相结合，用WKF来估计包裹相位空间动态模型中Zernike多项式的系数，稳定适用于高噪声下高精度的相位解包裹。

本发明采用以下技术方案：

一种相位解包裹方法，包括以下步骤：

S1、建立基于Zernike多项式拟合的包裹相位空间动态模型，将Zernike多项式从极坐标转换为直角坐标系，计算Zernike多项式Z_i(x,y)，确定包裹相位图与Zernike多项式的关系，得到真实相位φ(x,y)包含未知数的表达式；

S2、计算包裹相位图的质量图，剔除高噪声点；

S3、根据步骤S1得到的真实相位φ(x,y)和步骤S2包裹相位图的质量图，利用包裹卡尔曼滤波算法确定Zernike多项式的系数

S4、用参数CH表示Zernike多项式高阶系数与低阶系数的比值，并将其与阈值cth比较，判断步骤S3中包裹卡尔曼滤波算法的计算结果是否发散；

若包裹卡尔曼滤波算法收敛，计算最终的真实相位φ(x,y)；

若包裹卡尔曼滤波算法发散，进行再包裹运算得到

重复步骤S2至S4直至包裹卡尔曼滤波算法收敛，计算出真实相位φ(x,y)，完成相位解包裹。

具体的，步骤S1中，真实相位φ(x,y)为：

其中，Z_i(x,y)为定义在单位圆内的正交Zernike多项式的第i阶，c_i为相对应的待定系数，

为Zernike多项式总共所用到的阶数。

进一步的，包裹相位与Zernike多项式的关系：

其中，

为包裹后的相位，

为包裹运算，

为Zernike多项式总共所用到的阶数。

具体的，步骤S2中，记录包裹相位图中非四周边界的任意一像素点(i,j)与其8邻域的一阶差分绝对值D_ij，确定点(i,j)的质量R(i,j)，然后剔除R<5的高噪声点，R为点的质量。

进一步的，点(i,j)的质量R(i,j)为：

R(i,j)＝T(D_ij(1),th)+T(D_ij(2),th)+...+T(D_ij(8),th)

其中，T(·,th)为阈值运算。

具体的，步骤S3具体为：

S301、初始化l＝0时对应的状态向量和其协方差矩阵P；

S302、对状态向量和其协方差矩阵P在第l个有效像素点进行先验估计；

S303、基于卡尔曼增益，更新状态向量和其协方差矩阵P的后验估计；

S304、重复步骤S302和S303直至所有的有效像素点都参加计算，得到最终的Zernike多项式系数

进一步的，步骤S303中，定义预先设定好的测量误差协方差矩阵Σ＝100。

具体的，步骤S4中，阈值cth＝0.1。

与现有技术相比，本发明至少具有以下有益效果：

本发明一种相位解包裹方法，建立基于Zernike多项式拟合的包裹相位空间动态模型，使得识别出的高噪声点可以不被不纳入WKF计算，抑制了大噪声下WKF的发散问题；并且通过再包裹的方法，进一步提升了本发明的噪声鲁棒性。本发明与基于差分Zernike多项式拟合的相位解包裹方法相比，不需要计算相位梯度和多项式的一阶导数，噪声鲁棒性更好；与其它基于非线性卡尔曼滤波的相位解包裹算法相比，稳定性更好且速度更快。

进一步的，通过Zernike多项式拟合的方式，把对真实相位的求解问题转化成对相对应Zernike多项式待待定系数的求解问题。

进一步的，通过包裹相位与Zernike多项式的关系，建立了包裹相位的空间动态模型，以便于接下来用WKF求解待定的Zernike多项式系数。

进一步的，通过质量图去除R<5的高噪声点，抑制了高噪声点对后续WKF算法稳定性的影响，整体上提升了本发明的噪声鲁棒性。

进一步的，质量评价指标R能够有效的检测出高噪声。

进一步的，通过噪声鲁棒性强的WKF算法计算出待定的Zernike多项式系数。

进一步的，根据多次试验总结的经验，合理的设置参数Σ＝100能够有效的避免本发明中WKF算法发散的问题。

进一步的，合理的设定阈值cth＝0.1，能够准确的判断WKF算法是否发散，进而决定是否需要采取再包裹策略。

综上所述，本发明把Zernike多项式与WKF算法相结合，并采用再包裹的策略，实现了对高噪声的光学包裹相位的稳定、快速、高精度地解调。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明的实施路线图；

图2为Matlab软件仿真下不同噪声水平包裹相位图的解包裹结果；

图3为本发明与差分Zernike多项式拟合解包裹方法在不同噪声水下解包裹精度的比较。

具体实施方式

请参阅图1，本发明一种Zernike多项式拟合和包裹卡尔曼滤波算法WKF的相位解包裹方法，包括以下步骤：

S1、建立基于Zernike多项式拟合的包裹相位空间动态模型，计算Zernike多项式Z_i(x,y)；

Zernike多项式拟合光学相位的方法，即用光学相位图中大量的像素点去拟合相对而言较少的Zernike多项式系数，是一种数据降维数的方法，其本质上是超定问题的求解。这就意味着，求解Zernike多项式系数的过程，并不一定要求所有的像素点都参加计算。在相位解包裹中引入Zernike多项式，其优点在于可以将包裹相位图中的高噪声点去除，使之不参与相位解包裹的计算，以提高算法的噪声鲁棒性。

在相位解包裹中，包裹相位和真实相位的关系表示为：

其中，φ(x,y)为真实相位，

为包裹后的相位，

为式(2)所示的包裹运算：

通常来说，由光学表面检测引入的光学相位分解成一连串正交Zernike多项式的和，如下式(3)所示：

其中，Z_i(x,y)为定义在单位圆内的正交Zernike多项式的第i阶，其在极坐标内计算公式如式(4)和(5)所示；c_i为相对应的待定系数，也是本发明最终要求的参数；

为Zernike多项式总共所用到的阶数，本发明中

的取值为36；

其中，k为，指数n和m分别为径向阶次和角频率，m≤n且n-|m|＝偶数；参数ρ(ρ＜1)和θ分别表示极坐标的半径和角度。

通过式(6)和(7)将Zernike多项式从极坐标转换为直角坐标系：

其中，M和N分别为相位图像素尺寸的长和宽。

由式(1)和(3)可得包裹相位与Zernike多项式的关系：

将式(8)看成包裹相位在空域内随坐标点变化的空间动态模型，只要求得该空间动态模型中的未知参数c_i，通过式(3)计算出最终的真实相位φ(x,y)。

本发明中，求解c_i用的是噪声鲁棒性强的WKF算法，其具体内容见步骤S3。

S2、计算包裹相位图的质量图，剔除高噪声点；

对于任何一个解包裹算法而言，过多的大噪声点都会降低该算法的精度和稳定性。为了应对这一问题，本发明把一些高噪声点从相位解包裹的过程中剔除，使之不参与Zernike多项式系数的计算；通过求包裹相位图的质量图筛选高噪声点。

质量图的计算方法如下：

对于包裹相位图中非四周边界的任意一像素点(i,j)(1＜i＜M，1＜j＜N，M和N分别为相位图像素尺寸的长和宽)，用向量D_ij来记录(i,j)与其8邻域的一阶差分绝对值，如式(9)所示：

进一步，则该点的质量R(i,j)表示为：

R(i,j)＝T(D_ij(1),th)+T(D_ij(2),th)+…+T(D_ij(8),th) (10)

其中，T(·,th)为阈值运算，其定义式为：

在本发明中，取th＝1.5。

R<5的点为高噪声点，不纳入相位解包裹的计算中，R表示该点的质量，总体上来讲，R越大的点，噪声越小。

S3、利用WKF确定Zernike多项式的系数；

卡尔曼滤波(KF)算法的特点是对噪声鲁棒性强，能够在高噪声(尤其是针对高斯噪声)的环境下精确地估计动态系统的未知参数。WKF是KF的一种变形拓展，它解决了动态模型中的包裹运算

的非线性问题，使得部分非线性模型(例如包裹相位空间动态模型)也能使用该方法来估计其系统参数。

在式(8)所表示的包裹相位空间动态模型中，用于表示模型未知参数的状态向量为：

C＝[c₁,c₂,…,c₃₆]^T (12)

进一步，状态向量的更新方程为：

C_l＝FC_l-1 (13)

其中，F为单位矩阵，C_l-1和C_l和分别表示第l-1和第l个有效像素点(去除高噪声点后，按从左到右从上到下的扫描顺序，包裹相位图有效孔径内第l个像素点)对应的状态向量。

进一步，由式(8)可得系统的观测方程为：

其中，

为第l个有效像素点的包裹相位；(x_l,y_l)为第l个有效像素点的坐标；η_l表示包裹相位

的噪声；Z_l为Zernike多项式在第l个有效像素点的值，其表达式如下式(15)所示：

Z_l＝[Z₁(x_l,y_l),Z₂(x_l,y_l),…,Z₃₆(x_l,y_l)] (15)

其中，Z_i(x_l,y_l)表示第i阶Zernike多项式在点(x_l,y_l)的值，可由式(4)(5)(6)(7)计算得出。

进一步，WKF算法的具体计算步骤如下所示：

S301、初始化l＝0时对应的状态向量和其协方差矩阵P：

其中，E为求数学期望的操作，T表示矩阵的转置，上标-和+分别表示相关参数为先验估计和后验估计，

和

的取值对WKF算法最终的结果影响非常小，在实际应用中可以根据经验在一个较大的范围内选取。

S302、状态向量和其协方差矩阵在第l个有效像素点的先验估计：

S303、基于卡尔曼增益，更新状态向量和其协方差矩阵的后验估计：

其中，K_l为卡尔曼增益；

代表包裹运算；I为大小为

的单位矩阵；Σ为预先设定好的测量误差协方差矩阵，其值选取得过大或者过小会使卡尔曼滤波算法发散。

经过反复实验，当Σ＝100时，本发明中的WKF算法比较稳定。

重复步骤S302和S303直至所有的有效像素点都参加计算，得到最终的Zernike多项式系数

S4、判断是否需要再包裹运算；

在某些包裹相位图噪声非常大的极端情况下，步骤S3中的WKF的算法会得到一个发散的结果，其表现为计算出的Zernike多项式高阶系数过大，与实际情况不符。用参数CH表示Zernike多项式高阶系数与低阶系数的比值，并将其与阈值cth(本发明中取cth＝0.1)比较，判断步骤S3中的WKF算法的计算结果是否发散，如式(23)(24)所示：

进一步，若WKF算法和大部分情况一样收敛，则可直接用式(3)计算出最终的真实相位φ(x,y)；

若WKF算法发散，则需要把原始的包裹相位

加上一个固定相移

(本发明中

)，进行再包裹运算得到

如下式(25)所示。

进一步，用

重复步骤S2、S3、S4，直至WKF算法收敛，计算出真实相位φ(x,y)。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

依托Matlab软件，模拟出原始的真实相位图和与之对应不同噪声水平下的包裹相位图；通过本发明中的相位解包裹算法将这些不同噪声水平下的包裹相位图一一还原成真实相位图。将最终的计算结果与原始的真实相位图相减，得到代表计算精度的残差图。

请参阅图2，为不同噪声水平下的包裹相位图、计算真实相位、残差图及其均方根误差emrs；从图中可以看出，当噪声在σ＝1.5的极端情况下，残差图和均方根误差依然很小，说明本发明的相位解包裹算法在高噪声下任然有比较好的鲁棒性。

图3为本发明与差分Zernike多项式拟合方法在不同噪声水平下解包裹精度的对比图。其中横坐标表示噪声水平(添加在包裹相位上高斯噪声的均方根值)，纵坐标表示解包裹相位的精度(残差图的均方根值)，实线代表本发明中的解包裹算法，虚线代表差分Zernike多项式解调算法。从图中可以看出，两种算法在噪声水平比较小时(σ<0.5)其计算精度都比较高；但随着噪声水平逐渐增大到一定程度时(σ>0.5),差分Zernike多项式解调算法的计算误差急剧上升，而相比之下本发明的计算误差始终保持在一个比较小的范围内。相对于差分Zernike多项式解调算法，本发明算法拥有更好的噪声鲁棒性和解包裹精度。

综上所述，本发明一种基于Zernike多项式与WKF的相位解包裹方法，将WKF算法与Zernike多项式相结合，能够针对高噪声的光学包裹相位进行可靠、快速、高精度的解调。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (6)

1.一种基于Zernike多项式与WKF的相位解包裹方法，其特征在于，包括以下步骤：

S1、建立基于Zernike多项式拟合的包裹相位空间动态模型，将Zernike多项式从极坐标转换为直角坐标系，计算Zernike多项式Z_i(x,y)，确定包裹相位图与Zernike多项式的关系，得到真实相位φ(x,y)包含未知数的表达式；

S2、记录包裹相位图中非四周边界的任意一像素点(i,j)与其8邻域的一阶差分绝对值D_ij，确定点(i,j)的质量R(i,j)，然后剔除R<5的高噪声点，R为点的质量，点(i,j)的质量R(i,j)为：

R(i,j)＝T(D_ij(1),th)+T(D_ij(2),th)+…+T(D_ij(8),th)

其中，T(·,th)为阈值运算；

S3、根据步骤S1得到的真实相位φ(x,y)和步骤S2得到的质量R(i,j)，利用包裹卡尔曼滤波算法确定Zernike多项式的系数

S4、用参数CH表示Zernike多项式高阶系数与低阶系数的比值，并将其与阈值cth比较，判断步骤S3中包裹卡尔曼滤波算法的计算结果是否发散；

若包裹卡尔曼滤波算法收敛，计算最终的真实相位φ(x,y)；

若包裹卡尔曼滤波算法发散，进行再包裹运算得到

重复步骤S2至S4直至包裹卡尔曼滤波算法收敛，计算出真实相位φ(x,y)，完成相位解包裹。

2.根据权利要求1所述的方法，其特征在于，步骤S1中，真实相位φ(x,y)为：

其中，Z_i(x,y)为定义在单位圆内的正交Zernike多项式的第i阶，c_i为相对应的待定系数，

为Zernike多项式总共所用到的阶数。

3.根据权利要求2所述的方法，其特征在于，包裹相位与Zernike多项式的关系：

其中，

为包裹后的相位，

为包裹运算，

为Zernike多项式总共所用到的阶数。

4.根据权利要求1所述的方法，其特征在于，步骤S3具体为：

S301、初始化l＝0时对应的状态向量和其协方差矩阵P；

S302、对状态向量和其协方差矩阵P在第l个有效像素点进行先验估计；

S303、基于卡尔曼增益，更新状态向量和其协方差矩阵P的后验估计；

S304、重复步骤S302和S303直至所有的有效像素点都参加计算，得到最终的Zernike多项式系数

5.根据权利要求4所述的方法，其特征在于，步骤S303中，定义预先设定好的测量误差协方差矩阵Σ＝100。

6.根据权利要求1所述的方法，其特征在于，步骤S4中，阈值cth＝0.1。

Images