作者：翟天保Steven
版权声明：著作权归作者所有，商业转载请联系作者获得授权，非商业转载请注明出处

注：本文所讲内容为本人硕士毕业论文：《基于干涉图像质量分析的激光干涉仪抗振技术研究》，如有引用需要标注来源哈，如有疑问可以评论回复也可以邮箱1257147469@qq.com联系我。

前言

对干涉图组进行相位求解并完成解包裹后，可以得到一个倾斜的波面，从肉眼无法直观地看出波面的相关信息，因此需要通过Zernike多项式来拟合一个与所求波面最接近的理想波面，并根据拟合出的波面得到相关的系列像差信息和PV、RMS值，以此可以进一步完成光学系统成像质量的评价工作和面形的测量工作。

本文介绍了Zernike拟合的相关知识，包括Zernike多项式背景和发展、Zernike拟合具体实现过程、代码实现、参考文献、附录。

一、Zernike多项式背景和发展

通常光学成像系统在实际成像过程中，因系统本身的非理想特征，而使光线在经系统各面的传输过程中的路线出现偏差，进而形成各类像差，使结果显示出图像模糊、尺寸变化、形态异常等缺陷。像差是对光学系统成像质量的重要评价指标，对其校正是光学系统设计中非常关键的环节，一般来说，像差可分为球差、慧差、畸变、场曲、像散、位置色差和倍率色差七种，每种像差的形成原因和特征各不相同，因而可以根据其特征分析建立数学模型，并以此作为像差函数，该函数还可展开为幂级数或一组正交多项式。

泽尼克（Zernike）[1] [2]在利用刀口检验法和相衬法对一块圆形反射镜作测试时，推导出一套可表征不同像差函数的多项式，称作Zernike多项式，因其同像差多项式形式一致，常被用来拟合波面并对波前特性作分析。后来诸多学者在此基础上深入研究，提出了各类改进方案并进行了系列实验测试：中国科学院的惠梅[3]提出了一种基于Zernike多项式的波前数据拟合方法，通过阈值设定将误差点用掩膜过滤，有效提高测量精度；哈尔滨工业大学的刘剑锋等人[4]基于Zernike多项式推导出傅里叶变换公式，通过分析频域中信息以实现对面形误差的拟合；浙江大学的郑臻荣等人[5]在成像物镜设计中，应用Zernike多项式对自由曲面波像差进行分析，解决了自由曲面在成像系统设计中应用难度较大的问题；中国科学院的刘伟伟等人[6]结合三阶Zernike多项式，有效校正了全息阵列光镊中因算法导致的像差；苏州大学的王毅和余景池[7]根据Zernike36项系数，实现了对超精密车削系统的误差补偿，该方法效率高、稳定性好且准确；西班牙的I. Area等人[8]提出了一种通过差分算子和连接矩阵生成广义Zernike多项式的算法。

目前，Zernike多项式广泛应用于光学领域，特别是光学成像系统检测方面，是拟合波前的主要手段之一。在不同强度的振动环境下，拟合同一面形所得PV和RMS值的结果并不相同，根据这一特性可以判断抗振技术的能力强弱，当抗振效果好时，振动环境和无振动环境的结果相差不大。当然它也有一定的局限性：一方面，Zernike多项式无法直接计算出波面检测时的扰动信息，当外在振动过大时，其数据结果会随之出现不稳定的大区间波动情况，甚至可能出现无法拟合的情况；另一方面，它也无法有效描述某些车间加工中存在的机器制造误差。

综合来看，Zernike拟合是比较经典的波面拟合手段，不仅可以获得36项成像质量指标，也可以通过拟合后波面的PV、RMS值波动情况分析，来评价系统抗振性能等等。

二、Zernike拟合具体实现过程

在考虑光学系统像差时，像差函数因其光学结构不同而有所差异，与Zernike多项式相关的两个重要考虑方面分别是：1）系统的光瞳形状；2）是否为旋转对称结构。本文讲解的zernike拟合方法，是基于含圆形光瞳且不含旋转对称轴的结构，因此只说明对此类情况下的Zernike实现（其实只要把这个搞懂了，其他的类似拟合方法都清楚了，最初我研究Zernike拟合的时候，我们课题组纠结了许久它的实现逻辑，搞明白了发现还挺简单的。。。）。

在没有旋转轴对称的光学系统中，其像差函数不仅包含 $\cos (m \theta)$ 部分，还包含 $\sin (m \theta)$ 部分，将像差函数展开正交Zernike圆多项式 $\mathrm{Z}_{s}(\rho, \theta)$ 形式[9]：

$W(\rho, \theta)=\sum_{s} a_{s} \mathrm{Z}_{s}(\rho, \theta)$

(1)

其中， $W(\rho, \theta)$ 为拟合波面， $a_{s}$ 为多项式各项系数， $\rho$ 表示距离原点的长度， $\theta$ 表示相对原点的角度，多项式部分具体可分为三部分：

$\mathrm{Z}_{s}(\rho, \theta)=\left\{\begin{array}{c} \sqrt{2(n+1)} R_{n}^{m}(\rho) \cos (m \theta), m>0 \\ \sqrt{2(n+1)} R_{n}^{m}(\rho) \sin (m \theta), m<0 \\ \sqrt{n+1} R_{n}^{0}(\rho), m=0 \end{array}\right.$

(2)

一般n为7，即总项数为36。n值小的多项式靠前，n值一定时m值小的多项式靠前，按此顺序列出正交Zernike多项式 $\mathrm{Z}_{s}(\rho, \theta)$ 的36项表达式，如图1所示，该顺序同我分享的Matlab-zernike拟合程序里的顺序一致，多项式同《光学车间检测》这本书中第383-384页提到的一致，顺序可能有区别。

当然，在不同的应用和机构会采用不同的排列顺序和Zernike多项式，比如美国Zygo公司采取的就是如图2所示的多项式和顺序，许多公司也以Zygo作为标准（毕竟在这领域是老大哥。。40年发展太超前了，没办法，希望国内技术能早日超越它们，技术壁垒太难受了）。下图取自Zygo软件使用手册，目前我做相关研究用的也是这套多项式，和Zygo同一标准也是为了方便对比差异性。

继续讲，假设拟合波面 $W(\rho, \theta)$ 的尺寸为ROW×COL，ROW为行数，COL为列数，为了模拟像差函数的3D图形，需要先生成所需的网格数据（meshgrid操作）：设置一行x坐标，从-1到1，等距离分为COL份，将该行数据复制ROW行，建立X矩阵；设置一列y坐标，区间从1到-1，等距离分为ROW份，将该列数据复制COL列，建立Y矩阵。结合求得的X和Y矩阵，将其笛卡尔坐标系转换为极坐标系可得到两个矩阵数据信息R和Th，分别是 $\rho$ 和 $\theta$ 的网格数据，其对应公式为：

$\mathrm{R}(\mathrm{i}, \mathrm{j})=\sqrt{X(i, j)^{2}+Y(i, j)^{2}}$

$\operatorname{Th}(\mathrm{i}, \mathrm{j})=\arctan \left(\frac{Y(i, j)}{X(i, j)}\right)$

(3)

对每个像差函数，根据表1中对应的多项式表达式，再结合R和Th矩阵中各个像素点对应的 $\rho$ 和 $\theta$ 两项数据值，可得到此函数在对应点下的对应值，进而求得关于此函数尺寸大小为ROW×COL的矩阵信息，将矩阵数据以三维图展示，可直观看到其3D模拟样貌，36项像差的样貌图见文章最后的附录。

将上述的36项像差对应的矩阵进行变形，尺寸大小从ROW×COL变为（ROW*COL）×1，将原矩阵每行的数据按行序号顺序，依次以列数据形式放置在新矩阵内。之后，将36个变形后的矩阵按多项式顺序，合并在一起构成尺寸大小为（ROW*COL）×36的Zernike综合矩阵Z。将所需要拟合的解包裹后相位矩阵，同样变形为（ROW*COL）×1的新矩阵I。

此时，可通过矩阵求解计算Zernike36项系数，待解X矩阵大小为36×1：

ZX=I

(4)

式(4)的意义可以理解为：对整幅图上任何一个像素点而言，均满足该点对应的36项多项式值分别乘以对应的系数并累加，其结果约等同于该点在所需拟合的解包裹后相位图内的数值。

Zernike多项式的第一项为平移，第二三项分别是沿x和y方向的倾斜，将解包裹后相位矩阵减去Zernike前三项矩阵和对应系数的乘积和，即可得到不含倾斜量和平移量的真实面形图。对解包裹后相位信息经上述处理，得到的面形图如图3所示（对圆形干涉图外区域进行了屏蔽操作）。

对4-36项矩阵和对应系数的乘积进行求和，可得到不含倾斜量和平移量的理想面形图，如图4所示。

对图3所示的面形信息进行相关计算：先进行单位的转换，前文计算得到的结果均为相位信息（单位：rad），转为长度信息（单位：λ）；再计算PV值（面形最大值和最小值之差）和RMS值（面形的均方根值）。

综上所述，便是Zernike拟合的原理和其实现过程，讲的很细，大家不要嫌我啰嗦。。如果我当时研究这个的时候有人讲的细节一些，可能会给我节省许多时间。

三、代码实现

matlab代码见：

zernike拟合matlab程序_zernike多项式拟合-图像处理文档类资源-CSDN下载

C++代码考虑到商业保密性暂不公开分享，如有想法可评论或者邮箱联系我，一起交流学习~

四、参考文献

[1]Zernike F . Diffraction theory of the knife-edge test and its improved form, the phase-contrast method[J]. Monthly Notices of the Royal Astronomical Society, 1934, 94(2):377-384.

[2]F. Zernike. Beugungstheorie des schneidenver-fahrens und seiner verbesserten form, der phasenkontrastmethode[J]. Physica,1934(1):689-704

[3]惠梅,牛憨笨.运用泽尼克多项式进行物面波前数据拟合[J].光子学报,1999(12):3-5.

[4]刘剑峰,龙夫年,张伟,王治乐.基于泽尼克多项式进行面形误差拟合的频域分析[J].光学学报,2005(08):1062-1066.

[5]郑臻荣,孙旭涛,缪盈盈,顾培夫,刘旭.应用泽尼克多项式自由曲面的成像物镜设计[J].浙江大学学报(工学版),2008,42(12):2202-2206.

[6]刘伟伟,任煜轩,高红芳,孙晴,王自强,李银妹.泽尼克多项式校正全息阵列光镊像差的实验研究[J].物理学报,2012,61(18):524-531.

[7]王毅,余景池.基于泽尼克多项式系数的大相对孔径表面超精密车削误差的补偿方法[J].红外与激光工程,2012,41(03):724-728.

[8]A I A , B D K D , C E G . Recursive computation of generalised Zernike polynomials[J]. Journal of Computational and Applied Mathematics, 2017, 312(2):58-64.

[9]李萌阳,李大海,赵霁文,章辰,王琼华.基于方形域内标准正交矢量多项式的波前重建[J].光学学报,2014,34(07):144-150.

五、附录

总结

以上就是本文所讲的内容，简单介绍了Zernike拟合的相关知识。

如果文章帮助到你了，可以点个赞让我知道，我会很快乐~加油！

光学算法——Zernike拟合相关推荐

Matlab读Zygo干涉仪面形数据并进行37阶Zernike拟合（附Matlab代码）
本代码具有以下两个主要功能: 一.Matlab读Zygo干涉仪/轮廓仪.dat文件中的面形数据通过在Matlab中开发大量的反编译代码,我们成功地在Matlab中提取到Zygo .dat二进制文件中 ...
光学算法——经典枝切法（解包裹算法）
作者:翟天保Steven 版权声明:著作权归作者所有,商业转载请联系作者获得授权,非商业转载请注明出处注:本文所讲内容为本人硕士毕业论文:<基于干涉图像质量分析的激光干涉仪抗振技术研究> ...
【自动驾驶】车道线拟合算法---最小二乘法拟合直线
概览关于自动驾驶车道线拟合算法,常用的方法有B样条.三次样条插值.Ransac.最小二乘法等等. 但是针对于高精度地图的车道线拟合,由于车道线坐标点已知,所以不需要有控制点进行约束,那么B样条.贝塞 ...
数学建模常见算法：拟合算法
曲线拟合问题是指:已知平面上个点,,互不相同.寻求函数,使在某种准则下与所有数据点最为接近,即曲线拟合得很好. 从拟合的定义可以看出:拟合不需要使曲线一定经过给定的点,而插值要求多项式要经过所有样本 ...
建模算法笔记 | 拟合算法
在插值算法中,样本点太多,那么多项式的次数过高,会造成龙格现象.尽管可以用分段避免,但我们更想要获得一条确定的曲线,即使这条曲线不能经过每一个样本点,但只要保证误差足够小就可以. 1.最小二乘法 ar ...
zernike多项式波前拟合在屈光度测量的应用
之前有一个项目是做全自动验光仪的,用到了哈特曼原理,其实本质上就是zernike多项式波前拟合.关于zernike多项式,这里不做深入解析,只需要知道它是用来拟合曲面的函数,由于该函数在单位圆内每个多 ...
三维点集拟合：平面拟合、RANSAC、ICP算法
ACM算法分类:http://www.kuqin.com/algorithm/20080229/4071.html;CSDN容易吞图,不过编辑器里面图片还是显示的..... 一. 拟合一个平面 ...
Zernike函数拟合曲面--MATLAB实现
利用前36阶zernike函数拟合曲面: 脚本程序 clc;clear; load unwrap_ph.mat unwrap_ph=max(max(unwrap_ph))-unwrap_ph; unw ...
跟着清风学建模——拟合算法介绍及mathlab代码实现
目录拟合算法定义拟合好坏 mathlab代码部分: mathlab拟合工具使用拟合算法定义我们将随机的样本点为(xi,i),i=1,2,...,n,设置拟合曲线为y=kx+b.问题转为k和b值 ...

光学算法——Zernike拟合

前言