基于波面面形的滤波计算

2016-02-11夏志林刘世杰张志刚

光学精密工程 2016年12期

关键词：球差波面面形

夏志林,曾翔,,刘世杰,周游,张志刚

(1.武汉理工大学材料科学与工程学院，湖北武汉 430070;2.中国科学院上海光学精密机械研究所强激光材料重点实验室，上海 201800)

基于波面面形的滤波计算

夏志林1,曾翔1,2,刘世杰2*,周游2,张志刚2

(1.武汉理工大学材料科学与工程学院，湖北武汉 430070;2.中国科学院上海光学精密机械研究所强激光材料重点实验室，上海 201800)

为提高光学干涉测量的精度，研究了用来消除干涉条纹图噪声的滤波算法。讨论了不同波面和不同干涉仪噪声下不同滤波方法的滤波效果。考虑干涉仪噪声对条纹图的影响与检测中的不同波面情况，基于泽尼克多项式原理用matlab仿真法建立了离焦、像散、彗差、球差4种不同的干涉条纹图，计算了滤波后不同波面面形的PV值和GRMS值差别的大小。结果表明：同种噪声与滤波器情况下，在离焦、像散、彗差、球差4种面型中，滤波后球差波面的PV差值最大，高达0.11λ，与在乘性噪声下其它波面的最高PV值0.08λ相比，滤波效果较差；滤波后像散波面的PV差值最小，除乘性噪声中值滤波情况外，PV差值均在0.05λ以下，最低达0.01λ，滤波效果较好。此外，彗差波面滤波情况要略优于离焦波面。滤波器中，低通滤波稳定性高，PV差值局限于0.01λ～0.04λ，滤波效果相对好；而空域滤波后波面的PV差值最高达0.09λ，但噪声针对性更好。

干涉测量；干涉条纹图；波面面形；滤波方式

1 引言

光学干涉测量是一种集光电探测、计算机技术和图像处理为一体的测量手段。光学干涉测量技术探测被测物表面特征的主要过程包括：利用光电探测器读取光强信息、利用干涉条纹图表征被测表面形貌信息、通过干涉条纹图的解包获取形貌信息。在测量过程中会不可避免地引入振动和噪声因素，这些因素将对干涉条纹图的解包等过程产生干扰，导致无法获取较为准确的结果。因此，如何有效去除条纹图的固有噪声减少测量误差已成为光学干涉测量技术的关键问题。目前已有相当多的滤波算法用于消除条纹图中的噪声，例如传统的空域、变换域滤波算法。这些算法虽然有着良好的噪声去除效率，但破坏了条纹图中某些高频有效信息。总的来说，滤波算法一般采取两个途径实现去噪，一个是直接对采集的干涉强度图进行滤波，再提取和反演出待测量信息；另一个是暂时对采集的干涉强度数据不进行滤波，而对由其反演来的包裹相位图进行滤波处理。两个方法的数据处理途径虽然不同，但对滤波的要求都是一致的，即都是在抑制噪声和保持信号保真度两方面达到一种完美的平衡。研究发现对不同面形的被测物表面，滤波算法对有效信息的破坏程度也有所不同。因此，本文研究了不同波面、不同噪声情况下，滤波算法的选取条件。

2 Zernike多项式表征不同面形

光学测试中，通常用多项式来表达波前数据。Zernike多项式就是其中一种，它由相同形式的许多项组成，与光学测试中观察到的各种像差具有相同的形式。Zernike多项式的指导思想是将一个任意波面看作由无穷多个基面的线性组合[1]，若用前n项Zernike多项式来代表波面，则对于这一波面相位W，有：

W=a0+a1z1+a2z2+…+aizi+…+anzn,

(1)

其中zi为定义在单位圆上的正交多项式,zi表示各基面，ai为相应基面的权因子。

在仿真实验中，Zernike多项式的以下特点让它容易用于不同面形波面的仿真[2-3]：

(1)被测光学器件或者光学系统一般都具有圆形光瞳或圆形通光孔，经过归一化后正好是单位圆，而Zernike多项式的单位圆具有正交性，这恰好满足圆形光瞳的特点，同时多项式的正交性使得拟合多项式的系数相互独立，避免了系数之间的相互影响。

(2)Zernike多项式所特有的旋转对称性,使得光学问题的求解过程具有良好的收敛性，易于计算。

(3)Zernike多项式与初级像差有着一定的对应关系。

z3=2x2+2y2-1

z4=y2-x2

z5=2xy，

z6=-2y+3y3+3x2y

z7=-2x+3xy2+3x3

z8=1-6y2-6x2+6y4+12x2y2+6x4

(2)

式(2)为笛卡尔坐标系中初级像差表达式，各项分别对应不同波面相位:

z3对应于离焦波面相位；z4对应于像散波面相位(0°&90°方向)；z5对应于像散波面相位(45°&-45°方向)；z6对应于彗差波面相位(x轴方向)；z7对应于彗差波面相位(y轴方向)；z8对应于球差波面相位。

采用Zernike多项式作为不同波面的相位表达式，在生成干涉图时，可以通过改变式(1)中各表达式的系数来生成不同的干涉图，或通过式(1)中各表达式系数的取舍来仿真不同面形的波面相位。

3 干涉仪噪声分析

3.1 噪声的分类

图像噪声是由于图像摄取和传输时受到的随机信号干扰产生的，它是妨碍人们接收信息的因素。光干涉图像中也存在噪声。

按产生的原因，大致可以将光干涉图噪声分为随机噪声和系统噪声两大类[4]。其中，随机性噪声通常与光和电的基本性质密切相关。例如，由于粒子出现的随机性运动而产生的散斑噪声，还有导体中的自由运动电子没有规则的随机运动形成的热噪声；系统噪声则是指系统以外的信号，如电磁波信号和过电源串进系统内部而产生的噪声，电气设备的电磁干扰、自然放电现象等引起的噪声。系统噪声经常体现为脉冲噪声，表现为突发性出现的幅度高而持续时间短的离散脉冲。

若按概率分布函数或概率密度分布函数来描述光干涉图噪声，其分为加性噪声、乘性噪声和脉冲噪声[5]。加性噪声往往与图像灰度值绝对大小无关，如运算放大器信号处理过程中产生的噪声，是图像信号转移过程中引入的信道噪声，在图像像素上表现为叠加关系。乘性噪声通常与图像信号强度有关，一般是由信道不理想引起的噪声，它们与信号是某种程度的相乘关系。

3.2 噪声对干涉条纹图的影响

在干涉仪系统中，不同噪声对干涉条纹图的影响形式不同，以上噪声与图像强度的关系如下：

I(x,y)=I0(x,y)+I1(x,y)Im(x,y)cosφ(x,y)+In(x,y).

(3)

上式描述了光干涉图像I(x,y)中不同形式的噪声因素。In(x,y)为单独项形式的加性噪声，由于光干涉图像I(x,y)中存在相位场余弦调制，条纹信号的频率信息与噪声信号一般混合在一起；而乘性噪声Im(x,y)是与信号幅值的绝对强度有关的噪声，一般与条纹光强分布呈乘性关系。

加性噪声在光干涉图像上一般表现为热噪声、散弹噪声。一般的热噪声从频域上看，具有均匀的功率谱密度，是高斯分布的白噪声的一种主要类型。在matlab程序计算中，通常利用高斯噪声进行仿真。乘性噪声在光干涉图像上则一般表现为相干性噪声和散斑噪声。相干性噪声是物体表面漫反射出的相干光在空间相遇发生干涉而在图像上形成的随机分布、且大小不一的亮点与暗斑。仿真乘性噪声可通过在matlab程序中添加乘性噪声实现。在光干涉图像中，条纹图像中条纹的幅值和背景光强度常常会发生变化，而这两者的变化程度总是不同的且没有明显规律。这和系统噪声一致。仿真中采用添加随机噪声、椒盐噪声进行模拟。

图1显示了像散波面条纹图受各种噪声影响的表现形式。

图1 无噪声及受不同噪声影响的干涉条纹图Fig.1 Interferogram affected by different types of noises

4 条纹图滤波方式

干涉图的滤波方法分为空域滤波和变换域滤波两种，其中空域滤波实际上就是在图像中借助模板对图像邻域进行操作。空域滤波器有很多种类，但在空域中都是利用模板的卷积实现的[6]。常用的是均值滤波和中值滤波。在变换域滤波中常用的是傅里叶变换，条纹图中的主要滤波方式为低通滤波。现简要介绍一下几种常用的滤波方式[7-9]。

4.1 均值滤波

均值滤波是以当前点为中心，选择一个圆形或方形的滤波窗口，用窗口内各点的平均值代替当前点的值。实际上就是对输出像素邻域进行平均化操作，再将平均值作为输出像素的灰度。滤波窗口一般选为3×3或5×5。本文仿真实验中采用5×5窗口的均值滤波方式。

作为一种线性滤波器，图像均值滤波算法简单，易于实现，且计算速度快。但窗口中灰度极高或极低的像素点对该方法的影响很大，易造成边缘模糊，不利于后期的相位解包裹处理。

4.2 中值滤波

中值滤波是基于排序统计理论的能有效抑制噪声的非线性处理技术。中值滤波的应用前提是图像具有一种特性：噪声经常是以孤立点的形式出现的，这些点对应的像素数很少，而图像则是由像素数较多、面积较大的小块构成的。

中值滤波的基本操作是：对滤波窗口内的像素值进行排序，如果有一个突出值，就会被排到最右侧或最左侧，最终输出的中值，一般不是噪声点，从而达到去噪的目的。滤波窗口内的中值输出不受少量冲激噪声的影响，可有效滤除脉冲噪声和椒盐噪声。仿真中采用5×5窗口，即进行25个像素点排序处理的中值滤波方式。

4.3 基于傅里叶变换的低通滤波

频域滤波是图像在频率域中进行的一种非常重要的处理手段[10]。在数字图像中，图像的边缘和噪声对应于傅里叶变换谱中的高频部分，因此通过使用低通滤波器在频域对这些高频成分的抑制作用，从而达到消除空域中图像噪声或对图像边缘进行平滑模糊处理的目的。低通滤波[11-12]的基本操作是：将空域上的图像变换到频域上滤除噪声点，再反变换回空域图像。

傅里叶滤波方法的原理：一幅大小为M×N的图像f(x,y)，经过二维快速傅里叶变换为：

(4)

频域空间滤波可以表示为：

G(u,v)=F(u,v)H(u,v),

(5)

式中：H(u,v) 为转移函数；G(u,v) 为对图像进行频域滤波后的傅里叶能量谱。

通过傅里叶反变换得到噪声去除后的图像f(x,y)：

(6)

值得注意的是，在数字图像处理中，图像的傅里叶变换一般取M=N。常见的频率低通滤波器有：理想低通滤波器、巴特沃斯低通滤波器、高斯低通滤波器。考虑到理想低通滤波的振铃效应严重，低通滤波中采用巴特沃斯低通滤波器和高斯低通滤波器。

5 仿真计算

5.1 仿真实验流程

仿真实验采用matlab软件实现，仿真过程中分别添加随机噪声、椒盐噪声、乘性噪声、高斯噪声，分别采取均值滤波、中值滤波、巴特沃斯低通滤波与高斯低通滤波进行滤波。

面形仿真：通过Zernike多项式仿真不同波面面形，如图2(彩图见期刊电子版)所示，利用不同面形的波面相位可仿真生成不同的干涉条纹图。

为简化计算，在不同方向选取同一种像差进行计算。在公式(1)中选取z3离焦、z4像散、z6彗差、z8球差作为计算波面。则可得相应的波面相位：

p=a·z3+b·z4+d·z6+f·z8 ，

(7)

式中a、b、d、f分别指波面系数。在仿真其中的一种波面时，为简化计算并保证与实际检测效果相接近，需保证无噪声时该波面PV值为3λ。例如：像散波面p=a·z3时，若保证波面PV值为3λ，则系数a=1.517 8。同理对于其它波面可得：像散系数b=1.523 7；彗差系数d=1.567 8；球差系数f=2.0。

图2 不同面形的相位图Fig.2 Phase maps of difference surface shapes

图3 滤波处理前与处理后结果图与原图相位差Fig.3 Phase difference between processed results and original diagram

仿真处理后，结果图与原图相位残差以及未经滤波处理的相位残差图如图3(彩图见期刊电子版)所示：从图中可以看出随机噪声的分布是无规律的，并可简单地观察滤波效果。

本文以不同面形的相位残差的PV值、RMS值与GRMS值作为该面形在各种噪声、滤波方式下的性能评定指标。

5.2 结果分析

5.2.1 不同滤波器的滤波结果

图4为离焦系数a=1.517 8情况下，离焦波面在仿真后的PV值、RMS、GRMS值。由图4可知，在不同噪声的影响下，PV值、RMS值、GRMS值都呈现相似的趋势：均值滤波与中值滤波区域，滤波情况变化较大，低通滤波情况较好且效果稳定；均值滤波与中值滤波的滤波效果好坏与噪声情况密切相关，在分别添加均值噪声和中值噪声时效果明显。因此，在波面检测中优先考虑低通滤波器、在已判断出噪声的情况下，可采用均值滤波与中值滤波，也可两者结合使用。

图4 添加不同噪声滤波后PV差值、RMS、GRMS值结果图Fig.4 PV differences, RMSs, GRMSs of different methods when image with different noises

5.2.2 不同波面下的滤波情况

仿真实验中，选取不同波面系数：离焦波面系数a=1.517 8；像散波面系数b=1.523 7；彗差波面系数d=1.567 8；球差波面系数f=2.0。截取512 pixel×512 pixel面形区域，原波面PV值均为3λ。仿真相位残差PV值结果如下：

从数据整体上来看，球差波面f滤波后PV差值最大，滤波效果最差；像散波面b滤波后PV差值最小，滤波效果较好；彗差波面d的滤波情况略优于离焦波面a。此外，除去特例，如中值滤波器过滤椒盐噪声的情况；各种噪声下，不同波面下几种滤波方法的滤波效果大致相同。均值滤波下各个面形的滤波效果基本相同，而球差波面不适用于中值滤波；低通滤波相对空域滤波稳定性较高。

图5 不同波面添加相同噪声经滤波后PV差值Fig.5 PV differences of filtered image with adding the same noise on different wave surfaces

6 结论

文章对不同波面和不同干涉仪噪声下，几种常用滤波算法的优劣进行了探讨。选取离焦、像散、彗差、球差4种波面面形进行matlab仿真研究。结果表明：(1)不同噪声下同种波面的滤波情况差别不大，基本遵循相同的规律，低通滤波的情况较好且效果稳定；均值滤波与中值滤波的性能与噪声情况密切相关，在分别添加均值噪声和中值噪声时效果明显。(2)相同噪声下不同波面的滤波情况存在差别，球差波面滤波后PV差值最大，最高可达0.11λ，与乘性噪声下其它波面最高PV值(0.08λ)相比，其滤波效果普遍低于其它波面的情况；像散波面滤波后PV差值最小，除乘性噪声中值滤波情况外，PV差值均在0.05λ以下，最低达0.01λ，滤波效果较好；此外，彗差波面的滤波情况要略优于离焦波面。

[1] 惠梅, 牛憨笨. 运用泽尼克多项式进行物面波前数据拟合[J]. 光子学报,1999,28(12):1113-1116. HUI M, NIU H B. A method of wavefront data fitting using Zernike polynomials[J].ChineseJournalofPolarScience,1999,28(12):1113-1116. (in Chinese)

[2] D马拉卡拉主编，白国强等译. 光学车间检测[M]. 北京：机械工业出版社, 1983. Editor-in-chief of D Malakala.OpticalShopInspection[M]. BAI G Q translated. Beijing: Machinery Industry Press, 1983.(in Chinese)

[3] 高芬, 蒋庄德, 李兵,等. 基于扩展平均的多步相移算法及误差抑制特性比较[J]. 光子学报, 2014, 43(4):0426001. GAO F, JIANG ZH D, LI B,etal.. Multi-step phase-shifting algorithm based on extended mean and error comparision[J].ChineseJournalofPolarScience, 2014, 43(4):0426001.(in Chinese)

[4] 章霄. 数字图像处理技术[M]. 北京：冶金工业出版社, 2005:122-124. ZHANG X.DigitalImageProcessingTechnology[M]. Beijing: Metallurgical Industry Press, 2005:122-124.(in Chinese)

[5] 张鲜妮. 基于噪声分布特征的干涉图像滤波新方法[J]. 北京测绘, 2014(2):1-4,16. ZHANG X N. A new method of interference image filter based on noise distribution characteristics[J].BeijingSurveyingandMapping, 2014(2):1-4,16.(in Chinese)

[6] 刘乾. 抗振动移相干涉测量算法与实验研究[D]. 绵阳：中国工程物理研究院, 2015. LIU Q.Experimentalstudyandanti-vibrationphaseshiftinginterferencemeasurementalgorithm[D]. Mianyang:China Academy of Engineering Physics, 2015.(in Chinese)

[7] 闪萨雷斯. 数字图像处理[M].第二版. 北京：电子工业出版社, 2007. SANCALES.DigitalImageProcessing[M].second edition. Beijing: Electronic Industry Press, 2007.(in Chinese)

[8] 李伟. 光干涉图像处理中的关键算法研究[D]. 西安:西安科技大学, 2013. LI W.Researchonkeyarithmeticinopticalinterferometryimageprocessing[D]. Xi'an:Xi'an Technology University, 2013. (in Chinese)

[9] 吕绪浩. 光学元件干涉图数据的采集和处理方法研究[D]. 长春:长春理工大学, 2014. LV X H.Studyondataacquisitionandprocessingmethodofopticalelementinterferogram[D]. Changchun:Changchun University of Science and Technology, 2014.(in Chinese)

[10] LIU Q, WANG Y, JI F,etal.. A three-step least-squares iterative method for tilt phase-shift interferometry[J].OpticsExpress, 2013, 21(24):29505-29515.

[11] BONE D J, BACHOR H A, SANDEMAN R J. Fringe-pattern analysis using a 2-D Fourier transform [J].AppliedOptics, 1986, 25(10): 1653-1660.

[12] RODDIER C, RODDIER F. Interferogram analysis using Fourier transform techniques [J].AppliedOptics, 1987, 26(9): 1668-1673.

Filter calculation based on wavefront figures

XIA Zhi-lin1,ZENG Xiang1,2,LIU Shi-jie2*,ZHOU You2,ZHANG Zhi-gang2

(1.SchoolofMaterialsScienceandEngineering,WuhanUniversityofTechnology,Wuhan430070,China; 2.ShanghaiInstituteofOpticsandFineMechanics,ChineseAcademyofSciences,KeyLaboratoryofMaterialsforHighPowerLaser,Shanghai201800,China) *Correspondingauthor,E-mail:shijieliu@siom.ac.cn

To improve the measuring accuracy of optical interferometry, the filtering methods to eliminate the noise in interference fringe patterns were researched. The results of different filtering methods were discussed in different wavefronts and different noises of the interferometer. By taking the effects of interferometric noise on fringe patterns and different wavefronts into account, four kinds of fringe patterns for defocus, astigmatism, coma, spherical aberration were established based on zernike polynomials by matlab simulation. The differences of PV and GRMS values before and after filtering were calculated for different wavefronts. The results on four kinds of wavefronts for defocus, astigmatism, coma, spherical aberration show that under the same kind of noise and filtering, the PV difference of the spherical aberration is largest, up to 0.11λ, which shows a poor filtering effect as compared with the 0.08λ from other wavefront differences under the multiplicative noise. Furthermore, the PV difference of the astigmatism wavefront is smallest, and it is down to 0.05λ excepting median filtering in the multiplicative noise, showing a better filter effect. Moreover, the result of coma wavefront filtering is slightly better than that of the defocus wavefront filtering. For the filter methods, low-pass filters offer high-stability and the differences of the PV is from 0.01λ to 0.04λ, showing a better filtering effect.

interferometry; interferometric fringe pattern; wavefront figure; filtering method

2016-10-13；

2016-11-17.

国家自然科学基金青年科学基金资助项目(No.11602280)

1004-924X(2016)12-3033-08

TH744.3；TN247

:Adoi:10.3788/OPE.20162412.3033