夏志林,曾 翔,,劉世杰,周 游,張志剛

(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λ，濾波效果較好。此外，彗差波面濾波情況要略優(yōu)于離焦波面。濾波器中，低通濾波穩(wěn)定性高，PV差值局限于0.01λ～0.04λ，濾波效果相對好；而空域濾波后波面的PV差值最高達0.09λ，但噪聲針對性更好。

干涉測量；干涉條紋圖；波面面形；濾波方式

1 引 言

光學干涉測量是一種集光電探測、計算機技術和圖像處理為一體的測量手段。光學干涉測量技術探測被測物表面特征的主要過程包括：利用光電探測器讀取光強信息、利用干涉條紋圖表征被測表面形貌信息、通過干涉條紋圖的解包獲取形貌信息。在測量過程中會不可避免地引入振動和噪聲因素，這些因素將對干涉條紋圖的解包等過程產(chǎn)生干擾，導致無法獲取較為準確的結果。因此，如何有效去除條紋圖的固有噪聲減少測量誤差已成為光學干涉測量技術的關鍵問題。目前已有相當多的濾波算法用于消除條紋圖中的噪聲，例如傳統(tǒng)的空域、變換域濾波算法。這些算法雖然有著良好的噪聲去除效率，但破壞了條紋圖中某些高頻有效信息?？偟膩碚f，濾波算法一般采取兩個途徑實現(xiàn)去噪，一個是直接對采集的干涉強度圖進行濾波，再提取和反演出待測量信息；另一個是暫時對采集的干涉強度數(shù)據(jù)不進行濾波，而對由其反演來的包裹相位圖進行濾波處理。兩個方法的數(shù)據(jù)處理途徑雖然不同，但對濾波的要求都是一致的，即都是在抑制噪聲和保持信號保真度兩方面達到一種完美的平衡。研究發(fā)現(xiàn)對不同面形的被測物表面，濾波算法對有效信息的破壞程度也有所不同。因此，本文研究了不同波面、不同噪聲情況下，濾波算法的選取條件。

2 Zernike多項式表征不同面形

光學測試中，通常用多項式來表達波前數(shù)據(jù)。Zernike多項式就是其中一種，它由相同形式的許多項組成，與光學測試中觀察到的各種像差具有相同的形式。Zernike多項式的指導思想是將一個任意波面看作由無窮多個基面的線性組合[1]，若用前n項Zernike多項式來代表波面，則對于這一波面相位W，有：

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

(1)

其中zi為定義在單位圓上的正交多項式,zi表示各基面，ai為相應基面的權因子。

在仿真實驗中，Zernike多項式的以下特點讓它容易用于不同面形波面的仿真[2-3]：

(1)被測光學器件或者光學系統(tǒng)一般都具有圓形光瞳或圓形通光孔，經(jīng)過歸一化后正好是單位圓，而Zernike多項式的單位圓具有正交性，這恰好滿足圓形光瞳的特點，同時多項式的正交性使得擬合多項式的系數(shù)相互獨立，避免了系數(shù)之間的相互影響。

(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)中各表達式的系數(shù)來生成不同的干涉圖，或通過式(1)中各表達式系數(shù)的取舍來仿真不同面形的波面相位。

3 干涉儀噪聲分析

3.1 噪聲的分類

圖像噪聲是由于圖像攝取和傳輸時受到的隨機信號干擾產(chǎn)生的，它是妨礙人們接收信息的因素。光干涉圖像中也存在噪聲。

按產(chǎn)生的原因，大致可以將光干涉圖噪聲分為隨機噪聲和系統(tǒng)噪聲兩大類[4]。其中，隨機性噪聲通常與光和電的基本性質密切相關。例如，由于粒子出現(xiàn)的隨機性運動而產(chǎn)生的散斑噪聲，還有導體中的自由運動電子沒有規(guī)則的隨機運動形成的熱噪聲；系統(tǒng)噪聲則是指系統(tǒng)以外的信號，如電磁波信號和過電源串進系統(tǒng)內部而產(chǎn)生的噪聲，電氣設備的電磁干擾、自然放電現(xiàn)象等引起的噪聲。系統(tǒng)噪聲經(jīng)常體現(xiàn)為脈沖噪聲，表現(xiàn)為突發(fā)性出現(xiàn)的幅度高而持續(xù)時間短的離散脈沖。

若按概率分布函數(shù)或概率密度分布函數(shù)來描述光干涉圖噪聲，其分為加性噪聲、乘性噪聲和脈沖噪聲[5]。加性噪聲往往與圖像灰度值絕對大小無關，如運算放大器信號處理過程中產(chǎn)生的噪聲，是圖像信號轉移過程中引入的信道噪聲，在圖像像素上表現(xiàn)為疊加關系。乘性噪聲通常與圖像信號強度有關，一般是由信道不理想引起的噪聲，它們與信號是某種程度的相乘關系。

3.2 噪聲對干涉條紋圖的影響

在干涉儀系統(tǒng)中，不同噪聲對干涉條紋圖的影響形式不同，以上噪聲與圖像強度的關系如下：

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)是與信號幅值的絕對強度有關的噪聲，一般與條紋光強分布呈乘性關系。

加性噪聲在光干涉圖像上一般表現(xiàn)為熱噪聲、散彈噪聲。一般的熱噪聲從頻域上看，具有均勻的功率譜密度，是高斯分布的白噪聲的一種主要類型。在matlab程序計算中，通常利用高斯噪聲進行仿真。乘性噪聲在光干涉圖像上則一般表現(xiàn)為相干性噪聲和散斑噪聲。相干性噪聲是物體表面漫反射出的相干光在空間相遇發(fā)生干涉而在圖像上形成的隨機分布、且大小不一的亮點與暗斑。仿真乘性噪聲可通過在matlab程序中添加乘性噪聲實現(xiàn)。在光干涉圖像中，條紋圖像中條紋的幅值和背景光強度常常會發(fā)生變化，而這兩者的變化程度總是不同的且沒有明顯規(guī)律。這和系統(tǒng)噪聲一致。仿真中采用添加隨機噪聲、椒鹽噪聲進行模擬。

圖1顯示了像散波面條紋圖受各種噪聲影響的表現(xiàn)形式。

圖1 無噪聲及受不同噪聲影響的干涉條紋圖Fig.1 Interferogram affected by different types of noises

4 條紋圖濾波方式

干涉圖的濾波方法分為空域濾波和變換域濾波兩種，其中空域濾波實際上就是在圖像中借助模板對圖像鄰域進行操作?？沼驗V波器有很多種類，但在空域中都是利用模板的卷積實現(xiàn)的[6]。常用的是均值濾波和中值濾波。在變換域濾波中常用的是傅里葉變換，條紋圖中的主要濾波方式為低通濾波。現(xiàn)簡要介紹一下幾種常用的濾波方式[7-9]。

4.1 均值濾波

均值濾波是以當前點為中心，選擇一個圓形或方形的濾波窗口，用窗口內各點的平均值代替當前點的值。實際上就是對輸出像素鄰域進行平均化操作，再將平均值作為輸出像素的灰度。濾波窗口一般選為3×3或5×5。本文仿真實驗中采用5×5窗口的均值濾波方式。

作為一種線性濾波器，圖像均值濾波算法簡單，易于實現(xiàn)，且計算速度快。但窗口中灰度極高或極低的像素點對該方法的影響很大，易造成邊緣模糊，不利于后期的相位解包裹處理。

4.2 中值濾波

中值濾波是基于排序統(tǒng)計理論的能有效抑制噪聲的非線性處理技術。中值濾波的應用前提是圖像具有一種特性：噪聲經(jīng)常是以孤立點的形式出現(xiàn)的，這些點對應的像素數(shù)很少，而圖像則是由像素數(shù)較多、面積較大的小塊構成的。

中值濾波的基本操作是：對濾波窗口內的像素值進行排序，如果有一個突出值，就會被排到最右側或最左側，最終輸出的中值，一般不是噪聲點，從而達到去噪的目的。濾波窗口內的中值輸出不受少量沖激噪聲的影響，可有效濾除脈沖噪聲和椒鹽噪聲。仿真中采用5×5窗口，即進行25個像素點排序處理的中值濾波方式。

4.3 基于傅里葉變換的低通濾波

頻域濾波是圖像在頻率域中進行的一種非常重要的處理手段[10]。在數(shù)字圖像中，圖像的邊緣和噪聲對應于傅里葉變換譜中的高頻部分，因此通過使用低通濾波器在頻域對這些高頻成分的抑制作用，從而達到消除空域中圖像噪聲或對圖像邊緣進行平滑模糊處理的目的。低通濾波[11-12]的基本操作是：將空域上的圖像變換到頻域上濾除噪聲點，再反變換回空域圖像。

傅里葉濾波方法的原理：一幅大小為M×N的圖像f(x,y)，經(jīng)過二維快速傅里葉變換為：

(4)

頻域空間濾波可以表示為：

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

(5)

式中：H(u,v) 為轉移函數(shù)；G(u,v) 為對圖像進行頻域濾波后的傅里葉能量譜。

通過傅里葉反變換得到噪聲去除后的圖像f(x,y)：

(6)

值得注意的是，在數(shù)字圖像處理中，圖像的傅里葉變換一般取M=N。常見的頻率低通濾波器有：理想低通濾波器、巴特沃斯低通濾波器、高斯低通濾波器?？紤]到理想低通濾波的振鈴效應嚴重，低通濾波中采用巴特沃斯低通濾波器和高斯低通濾波器。

5 仿真計算

5.1 仿真實驗流程

仿真實驗采用matlab軟件實現(xiàn)，仿真過程中分別添加隨機噪聲、椒鹽噪聲、乘性噪聲、高斯噪聲，分別采取均值濾波、中值濾波、巴特沃斯低通濾波與高斯低通濾波進行濾波。

面形仿真：通過Zernike多項式仿真不同波面面形，如圖2(彩圖見期刊電子版)所示，利用不同面形的波面相位可仿真生成不同的干涉條紋圖。

為簡化計算，在不同方向選取同一種像差進行計算。在公式(1)中選取z3離焦、z4像散、z6彗差、z8球差作為計算波面。則可得相應的波面相位：

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

(7)

式中a、b、d、f分別指波面系數(shù)。在仿真其中的一種波面時，為簡化計算并保證與實際檢測效果相接近，需保證無噪聲時該波面PV值為3λ。例如：像散波面p=a·z3時，若保證波面PV值為3λ，則系數(shù)a=1.517 8。同理對于其它波面可得：像散系數(shù)b=1.523 7；彗差系數(shù)d=1.567 8；球差系數(shù)f=2.0。

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

圖3 濾波處理前與處理后結果圖與原圖相位差Fig.3 Phase difference between processed results and original diagram

仿真處理后，結果圖與原圖相位殘差以及未經(jīng)濾波處理的相位殘差圖如圖3(彩圖見期刊電子版)所示：從圖中可以看出隨機噪聲的分布是無規(guī)律的，并可簡單地觀察濾波效果。

本文以不同面形的相位殘差的PV值、RMS值與GRMS值作為該面形在各種噪聲、濾波方式下的性能評定指標。

5.2 結果分析

5.2.1 不同濾波器的濾波結果

圖4為離焦系數(shù)a=1.517 8情況下，離焦波面在仿真后的PV值、RMS、GRMS值。由圖4可知，在不同噪聲的影響下，PV值、RMS值、GRMS值都呈現(xiàn)相似的趨勢：均值濾波與中值濾波區(qū)域，濾波情況變化較大，低通濾波情況較好且效果穩(wěn)定；均值濾波與中值濾波的濾波效果好壞與噪聲情況密切相關，在分別添加均值噪聲和中值噪聲時效果明顯。因此，在波面檢測中優(yōu)先考慮低通濾波器、在已判斷出噪聲的情況下，可采用均值濾波與中值濾波，也可兩者結合使用。

圖4 添加不同噪聲濾波后PV差值、RMS、GRMS值結果圖Fig.4 PV differences, RMSs, GRMSs of different methods when image with different noises

5.2.2 不同波面下的濾波情況

仿真實驗中，選取不同波面系數(shù)：離焦波面系數(shù)a=1.517 8；像散波面系數(shù)b=1.523 7；彗差波面系數(shù)d=1.567 8；球差波面系數(shù)f=2.0。截取512 pixel×512 pixel面形區(qū)域，原波面PV值均為3λ。仿真相位殘差PV值結果如下：

從數(shù)據(jù)整體上來看，球差波面f濾波后PV差值最大，濾波效果最差；像散波面b濾波后PV差值最小，濾波效果較好；彗差波面d的濾波情況略優(yōu)于離焦波面a。此外，除去特例，如中值濾波器過濾椒鹽噪聲的情況；各種噪聲下，不同波面下幾種濾波方法的濾波效果大致相同。均值濾波下各個面形的濾波效果基本相同，而球差波面不適用于中值濾波；低通濾波相對空域濾波穩(wěn)定性較高。

圖5 不同波面添加相同噪聲經(jīng)濾波后PV差值Fig.5 PV differences of filtered image with adding the same noise on different wave surfaces

6 結 論

文章對不同波面和不同干涉儀噪聲下，幾種常用濾波算法的優(yōu)劣進行了探討。選取離焦、像散、彗差、球差4種波面面形進行matlab仿真研究。結果表明：(1)不同噪聲下同種波面的濾波情況差別不大，基本遵循相同的規(guī)律，低通濾波的情況較好且效果穩(wěn)定；均值濾波與中值濾波的性能與噪聲情況密切相關，在分別添加均值噪聲和中值噪聲時效果明顯。(2)相同噪聲下不同波面的濾波情況存在差別，球差波面濾波后PV差值最大，最高可達0.11λ，與乘性噪聲下其它波面最高PV值(0.08λ)相比，其濾波效果普遍低于其它波面的情況；像散波面濾波后PV差值最小，除乘性噪聲中值濾波情況外，PV差值均在0.05λ以下，最低達0.01λ，濾波效果較好；此外，彗差波面的濾波情況要略優(yōu)于離焦波面。

[1] 惠梅, 牛憨笨. 運用澤尼克多項式進行物面波前數(shù)據(jù)擬合[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]. 北京：機械工業(yè)出版社, 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] 章霄. 數(shù)字圖像處理技術[M]. 北京：冶金工業(yè)出版社, 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] 閃薩雷斯. 數(shù)字圖像處理[M].第二版. 北京：電子工業(yè)出版社, 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] 呂緒浩. 光學元件干涉圖數(shù)據(jù)的采集和處理方法研究[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

夏志林(1978-)，男，江蘇東臺人，研究員,2007年于中國科學院上海光學精密機械研究所獲得博士學位，現(xiàn)為武漢理工大學材料科學與工程學院教師，主要從事激光與物質相互作用方面的研究。E-Mail:xiazhilin@whut.edu.cn

劉世杰(1979-)，男，陜西人，副研究員，碩士生導師，2008年于中國科學院上海光學精密機械研究所獲得博士學位，主要從事光學檢測技術的研究。E-mail：shijieliu@siom.ac.cn