一種有效解決離軸數(shù)字全息相圖傾斜畸變的數(shù)字參考平面方法?

李芳 王明清 鄭明 盧葦 于慶南 賈燕 吳堅

(北京航空航天大學(xué)物理科學(xué)與核能工程學(xué)院,北京 100191)

1 引 言

離軸數(shù)字全息由于參考光的傾斜入射可以使零級衍射、重建像及共軛像在波場的重建過程中彼此分離,從而方便地在成像中消除零級衍射和共軛像的影響而獲得廣泛的應(yīng)用[1,2].其中,平面光波由于計算簡單,是離軸數(shù)字全息中廣泛使用的一種參考光波形[3?5].然而,這種離軸成像方式由于引入了參考波的附加載波頻率會導(dǎo)致重建像產(chǎn)生傾斜相位畸變[6].針對此問題最直接的解決方案就是設(shè)法測得平面參考波的附加載波頻率或傾角,然后將其代入重建算法而消除傾斜畸變.這種方法的核心問題在于如何準確地測得平面參考光波的附加載波頻率或傾角值.因為到目前為止,通過實驗方法很難準確地測出參考光波的附加載波頻率值或傾角.

為解決這個問題,國內(nèi)外的研究者提出了多種解決方案,大致可以分為雙曝光法[7]、數(shù)字參考波法[8?13]、頻譜中心法[14?16]及自動修正算法[17?21].其中,雙曝光法通過分別記錄包含物體和不含物體的兩幅雙光束干涉全息圖并進行相減以對重建相圖的傾斜畸變做出修正.由于該方法需要記錄兩幅全息圖和要求系統(tǒng)具有高穩(wěn)定性,因此難以用于動態(tài)成像.而數(shù)字參考波法是通過建立和開展數(shù)字參考波的計算以確定其載波頻率,從而達到對重建相圖的傾斜畸變修正的目的.該方法需要一些特殊的條件,例如在全息圖平面的物波場需為緩變場、針對微小物體或平坦物體成像.除此以外,像素大小對條紋分辨率的影響也限制了獲取參考波載波頻率的精確性,因此該方法具有一定的局限性.頻譜中心法是將全息圖變換到頻域,通過尋找和定位物體頻譜的最大值點并移動到坐標中心而消除重建相圖傾斜畸變的方法.這種方法易受像素分辨率及噪聲的影響而難以準確定位物體頻譜的最大值,因此重建像仍會存在一定的殘留畸變.自動修正算法通過在波場的重建計算中引入相位掩模來修正相位畸變.該方法更多地用于透鏡導(dǎo)致的二次相位畸變的修正,是一種理論上的理想模型,如果存在噪聲和環(huán)境擾動的影響,則會出現(xiàn)較大誤差.

為了找到能更好地解決離軸數(shù)字全息成像中相圖傾斜畸變修正問題的方法,本文提出了一種基于無透鏡離軸數(shù)字全息成像的傾斜相位畸變修正的數(shù)字參考平面算法.對該方法的論述包括:1)相圖的解包裹重建和抑噪處理;2)從相圖選點以構(gòu)建一個能準確反映相圖傾斜的數(shù)字參考平面;3)建立數(shù)字參考平面參量與平面參考波載波頻率的數(shù)學(xué)關(guān)系;4)以構(gòu)建的數(shù)字參考平面為判據(jù)進行傾斜相位畸變修正的迭代計算;5)實驗驗證.

2 算法原理

2.1 理論模型

離軸數(shù)字全息成像中,物平面、全息圖平面以及像平面的坐標關(guān)系如圖1所示.

圖1 離軸數(shù)字全息成像中的物平面、全息圖平面以及像平面的坐標關(guān)系Fig.1.Coordinate systems of object,hologram and reconstructed image planes for digital o ff-axis holography.

物光波和參考光波在全息圖平面干涉產(chǎn)生的全息圖強度可表達為

其中,和分別表示物光波和參考光波的復(fù)振幅;||2+||2表示零級衍射強度;分別表示再現(xiàn)像及其共軛像強度.對于無透鏡離軸數(shù)字全息系統(tǒng),在全息圖平面(z=0)處的平面參考光波復(fù)振幅(x,y)可以表示為

其中,R(x,y)是參考光波的實振幅,(fx,fy)是平面參考光波分別在x和y方向上的實際載波頻率.

通過對參考光波傾角的適當設(shè)計可以方便地使(1)式中的三項在成像平面彼此完全分離,如此零級衍射和共軛像可以通過簡單的頻域濾波消除.因此,將(2)式代入(1)式并去除零級斑和共軛項,則濾波后的全息圖強度可表示為

依據(jù)傅里葉光學(xué)和基爾霍夫衍射理論,在像平面重建的物波場可以寫為

式中,i/λ是系數(shù);λ是光波長;r是全息圖平面和像平面上任意兩點間的距離;是用于圖像重建計算的參考波復(fù)振幅,其表達式為

其中R′(x,y)是實振幅,是參考波分別在x′和y′方向上的名義載波頻率值,它等于平面參考光波載波頻率的設(shè)計值.由于實際記錄全息圖的平面參考光波的載波頻率難以準確測知,因此用于圖像重建的參考波的名義載波頻率并不嚴格地等于參考光波的實際載波頻率(fx,fy),即設(shè)計值不等于實際值.因此,兩個參考波的載波頻率與(fx,fy)會存在一定的頻差,這個差值將導(dǎo)致重建相圖的傾斜畸變.將(3)和(5)式代入(4)式得到重建像的光場為

其中,根據(jù)傅里葉光學(xué)理論,由于積分中第一項的相位因子2π(?fxx+?fyy)是線性相因子,因此重建的圖像會發(fā)生一定的偏轉(zhuǎn)而偏離圖像平面中心,即導(dǎo)致圖像的傾斜相位畸變.根據(jù)(6)式,顯然如果?fx和?fy等于零,則重建的圖像會準確地呈現(xiàn)在平面中心而沒有傾斜畸變.由于在數(shù)字全息中用于圖像重建的參考波是以數(shù)字波函數(shù)的形式體現(xiàn),因此可以通過修正該參考波的名義載波頻率值使(?fx,?fy)=0,即達到與全息圖記錄使用的參考光波實際載波頻率一致.

為實現(xiàn)這個目標,我們構(gòu)建了一種能夠準確反映圖像重建中傾斜相位畸變程度的數(shù)字參考平面.具體方法是首先獲得三維重建相圖并做抑噪處理以減小修正誤差.然后在相圖的平坦區(qū)域內(nèi)隨機選取三個點的空間坐標值以建立數(shù)字參考平面方程:

其中,(a,b,c)是反映數(shù)字參考平面傾斜程度的系數(shù),它對應(yīng)于重建圖像的傾斜.數(shù)字參考平面的法向量指示了傾斜的方向和大小,可以表示為

其中,(i,j,k)對應(yīng)于(x′,y′,z′)方向的單位向量;系數(shù)(a,b,c)的值可以通過計算n的行列式來獲得.由于(6)式中的線性相因子引起的圖像偏轉(zhuǎn)傾斜和(7)式中差分系數(shù)(a,b)表示的圖像平面傾斜指示的是同樣的結(jié)果,因此以下的對應(yīng)關(guān)系能夠被建立:

以上關(guān)系顯示,如果系數(shù)(a,b)的值等于零,則顯然(?fx,?fy)也等于零,即圖像重建的數(shù)字平面參考波的名義載波頻率將等于全息圖記錄的實際平面參考光波的載波頻率值(fx,fy).在這種情況下,實際平面參考光波的載波頻率值(fx,fy)也得以準確地確定.基于這種考慮,在隨后對名義載波頻率值的修正計算中,將數(shù)字參考平面參量(a,b)值的變化設(shè)為判定?fx和?fy是否為零的判據(jù).修正名義載波頻率值和傾斜相位畸變的迭代計算過程如圖2所示.

圖2 修正數(shù)字參考波名義載波頻率和傾斜相位畸變的迭代計算流程Fig.2.Flow-chart for correcting nominal carrier frequencies,of the numerical reference wave and tilt phase distortion.

2.2 模擬分析和驗證

對該方法的計算機模擬實驗結(jié)果如圖3所示,圖3(a)是原始物體.模擬實驗中,引入高斯噪聲作為對激光散斑噪聲的模擬[22].假設(shè)用于全息圖記錄的斜入射平面參考光波的載波頻率設(shè)計值為實際值為fx=890 cm?1,fy=890 cm?1.因此,用于全息圖重建的數(shù)字參考波的名義載波頻率值為與參考光波的實際載波頻率不等.利用菲涅耳重建算法[23]和PUMA相位解包裹技術(shù)[24]獲得的重建相圖如圖3(b)所示,其結(jié)果不僅顯示了明顯的相位傾斜畸變,而且還有噪聲引起的相位畸變.因此,為了有效地修正相圖的傾斜畸變,在使用平穩(wěn)的PUMA解包裹算法的基礎(chǔ)上,做了進一步的抑噪處理.首先通過雙邊濾波[25]以有效保留去噪后的圖像高對比度邊緣特征信息,然后進一步通過包含小波收縮的短時傅里葉變換以更好地保留圖像紋理特征等低對比度信息[26].抑噪處理的基本原理可以在數(shù)學(xué)上表示為

其中,Sp,f是噪聲傅里葉系數(shù);Fp是所有小波收縮信號的平均值;Kp,f是與噪聲標準差σn有關(guān)的收縮系數(shù);p表示噪聲被抑制后的結(jié)果.然而,由于殘余噪聲的不確定性,很難確定已解包裹相位的噪聲標準差σn.對此,我們使用噪聲標準差的評估公式,該值表示為

其中,是高頻子頻帶的正交小波轉(zhuǎn)換系數(shù),Median表示在計算中的中值,0.6745是基于經(jīng)驗公式獲得的參數(shù)[27].經(jīng)上述抑噪處理的結(jié)果如圖3(c)所示.解包裹后的重建相圖中的噪聲得到有效抑制,為隨后的數(shù)字參考平面構(gòu)建和傾斜畸變的修正提供了保障.

圖3 仿真計算結(jié)果 (a)原始物體;(b)包含噪聲的PUMA解包裹相圖;(c)抑噪處理后的相圖;(d)構(gòu)建的數(shù)字參考平面;(e)傾斜相位畸變修正后的相圖;(f)傾斜相位畸變的誤差分析Fig.3.Simulation results:(a)Original object;(b)reconstructed phase map with noise using PUMA phase unwrapping technique;(c)phase map with noise removed;(d)constructed numerical reference plane;(e)phase map with tilt phase distortion corrected;(f)error analysis of the tilt phase distortion correction.

根據(jù)上述的數(shù)字參考平面構(gòu)建原理,在圖3(c)的平坦區(qū)域選擇三個點的坐標,利用(7)式建立數(shù)字平面方程,獲得反映圖像傾斜程度的數(shù)字參考平面,如圖3(d)所示.顯然,在這種情況下(a,b)值不等于零.然后,利用圖2的流程,對參考波的名義載波頻率進行修正的迭代計算,計算的誤差值控制在|a|<10?5和|b|<10?5.當頻率值為時,傾斜相位畸變修正的結(jié)果如圖3(e)所示.修正的誤差如圖3(f)所示.此時,(a,b)值為 (?5.4976×10?6,?5.8756×10?6).重建的相圖(如圖3(b)或圖3(c))的傾斜畸變得到了很好的修正.同時,由(9)式和(a,b)值可知,(?fx,?fy)足夠小,因此獲得平面參考光波的實際載波頻率為

作為比較,圖4和表1給出了采用雙曝光法[7]、數(shù)字參考波法[13]、頻譜中心法[15]和線性擬合法[18]獲得的相圖傾斜畸變修正的仿真計算結(jié)果.圖4(a)的雙曝光法結(jié)果顯示,相圖的傾斜畸變能夠得到很好的修正,但這僅僅是一種理想情況.由于該方法需要記錄兩幅全息圖,因此不能有任何擾動的影響,也就很難在實際環(huán)境中獲得理想結(jié)果.圖4(b)是使用數(shù)字參考波法得到的計算結(jié)果,由于該算法僅在振幅域?qū)D像進行處理,因此在相位域上仍會存在一定的誤差和殘余畸變,即相圖的傾斜畸變未能完全修正.而圖4(c)顯示的頻譜中心法計算結(jié)果是通過選取頻譜最大值點并將其移至頻譜中心實現(xiàn)相圖畸變修正.該方法在像素分辨率及噪聲的影響下,僅能定位物頻譜的最大值點,參考波傾角或載波頻率很小時,重建的相圖仍會存在一些殘留畸變而不能完全修平.圖4(d)顯示了采用線性擬合法獲得的相圖傾斜畸變修正結(jié)果.顯然該方法的結(jié)果優(yōu)于數(shù)字參考波法和頻譜中心法的結(jié)果.由于線性擬合法是從相圖中做直線擬合,即使取(x,y)雙方向做直線擬合,仍然存在所取的直線區(qū)域過于狹小而難以準確表征整個二維平面,因此重建相圖在使用該方法后仍然存在少量的殘余畸變.

圖4 不同方法的仿真結(jié)果對比 (a)雙曝光法;(b)數(shù)字參考波法;(c)頻譜中心法;(d)線性擬合法Fig.4.Comparison of the simulation results from different methods:(a)Double exposure technique;(b)numerical reference wave technique;(c)spectrum-centering algorithm;(d)linearly fitting method.

表1列出了圖3和圖4結(jié)果的量化對比分析.由于a和b是表征圖像傾斜畸變修正程度的關(guān)鍵特征參量,因此對圖4中不同方法得到的圖像傾斜畸變修正結(jié)果,通過(7)式建立對應(yīng)的平面方程可求得其(a,b)值.結(jié)果顯示,本文提出的方法在修正相圖傾斜畸變上更加精確.同時,基于圖4的幾種方法獲得的實際平面參考光波的載波頻率值與已知的890 cm?1比較也存在一定的誤差.其中,雙曝光法無法通過圖像傾斜畸變的修正結(jié)果推算出實際平面參考光波的載波頻率.

表1 不同方法的圖像傾斜畸變修正結(jié)果Table 1.Results of image tilt distortion correction with different methods.

圖5 離軸數(shù)字全息的實驗光路圖Fig.5.Experimental setup for digital o ff-axis holography.

圖6 實驗結(jié)果 (a)重建的分辨率板局部振幅像;(b)解包裹和去噪后的局部分辨率板相圖;(c)構(gòu)建的數(shù)字參考平面;(d)傾斜畸變修正后的分辨率板相圖Fig.6.Experimental results:(a)Reconstructed amplitude image of a partial resolution plate;(b)unwrapped and de-noised phase map of the partial resolution plate;(c)constructed numerical reference plane;(d)phase map of the resolution plate with the tilt phase distortion corrected.

3 實驗測量和驗證

離軸數(shù)字全息實驗光路如圖5所示.其中,波長λ=632.8 nm的激光束經(jīng)擴束系統(tǒng)和分束鏡形成物光束和參考光束;被測物為三維分辨率板;用于記錄全息圖的電荷耦合器(CCD)傳感器像素數(shù)1280×1024,像素尺寸5.2μm×5.2μm.傾斜的平面參考光名義載波頻率設(shè)置為

實驗結(jié)果如圖6所示.其中圖6(a)是利用菲涅耳算法重建的分辨率板局部振幅像,其結(jié)果包含了散斑等噪聲信息;圖6(b)是通過PUMA解包裹算法結(jié)合雙邊濾波及小波收縮短時傅里葉變換處理后得到的物體相圖.顯然,在重建物體三維圖像的同時,實驗中的噪聲得到了有效的消除.但由于參考光的傾斜入射引入了附加的載波頻率,重建相圖出現(xiàn)一定的傾斜畸變.在相圖的平坦部分取三點構(gòu)建表征圖像傾斜程度的數(shù)字參考平面,如圖6(c)所示.此時,(a,b)值明顯不等于零.利用圖2的方法對數(shù)字參考波的名義載波頻率進行修正迭代計算,實驗單次迭代計算時間為7.8 s,迭代計算步長d的初始值為10,為提高計算精度,設(shè)置步長d=0.1×d.當計算的載波頻率值為時,獲得(a,b)=(2.43946×10?6,?7.43256×10?6),計算完成.其中,迭代計算33次,用時258 s.此時,圖6(b)的傾斜相圖得到了非常好的修正,如圖6(d)所示.因此,獲得的實際參考光波的載波頻率值為fx=?291.47 cm?1,fy=?137.43 cm?1.誤差在10?6量級.

4 結(jié) 論

在離軸數(shù)字全息中,使用傾斜的平面參考波以消除成像中的零級衍射和共軛像是一種簡捷和廣泛使用的方法,但由于傾斜平面參考波引入了附加的載波頻率并很難通過實驗測量準確地獲得傾斜參考波的附加載波頻率值,因此會導(dǎo)致重建的物體相圖存在一定的傾斜畸變而無法完全修正的問題.對此本文提出了一種數(shù)字參考平面算法用以解決這一問題.該算法利用重建相圖的平坦區(qū)域選點構(gòu)建一個能準確表征相圖傾斜的數(shù)字參考平面,并建立該平面參量與參考波附加載波頻率的數(shù)學(xué)關(guān)系和作為隨后相圖畸變修正計算的判據(jù).同時,所提出的方法不僅使用了平穩(wěn)的PUMA相位解包裹重建算法,而且結(jié)合了雙邊濾波及小波收縮短時傅里葉變換的抑噪處理,最大限度地在消除噪聲對修正傾斜相位畸變影響的同時獲得較好的圖像分辨率.因此該方法在環(huán)境和系統(tǒng)噪聲的影響下仍然有效,不僅能很好地實現(xiàn)對傾斜相位畸變的準確修正,而且能準確地獲得傾斜平面參考波的附加載波頻率.模擬分析和實驗結(jié)果均證明了該方法的有效性.對離軸數(shù)字全息獲得準確的重建圖像具有重要意義.

參考文獻

[1]Wang H Y,Liu F F,Song X F,Liao W,Zhao B Q,Yu M J,Liu Z Q 2013Acta Phys.Sin.62 024207(in Chinese)[王華英,劉飛飛,宋修法,廖薇,趙寶群,于夢杰,劉佐強2013物理學(xué)報62 024207]

[2]Gu T T,Huang S J,Yan C,Miao Z,Chang Z,Wang T Y 2015Acta Phys.Sin.64 064204(in Chinese)[谷婷婷,黃素娟,閆成,繆莊,常征,王廷云 2015物理學(xué)報 64 064204]

[3]Wu Y C,Wu X C,Yao L C,Xue Z L,Wu C Y,Zhou H,Cen K 2017Fuel195 12

[4]Yuan C J,Zhong L Y,Wang Y P,Xu L X,Qian X F 2004Laser Technol.28 482

[5]Rong L,Xiao W,Pan F,Liu S,Li R 2010Chin.Opt.Lett.8 653

[6]Cuche E,Marquet P,Dahlgren P,Depeursinge C 2000Interferometry in Speckle Light(Berlin: Springer-Verlag)pp213–218

[7]Ferraro P,de Nicola S,Finizio A,Coppola G,Grilli S,Magro C,Pierattini G 2003Appl.Opt.42 1938

[8]Liebling M,Blu T,Unser M 2004J.Opt.Soc.Am.A21 367

[9]Colomb T,Cuche E,Charrière F,Kühn J,Aspert N,Frédéric M,Pierre M,Christian D 2006Appl.Opt.45 851

[10]Miccio L,Al fieri D,Grilli S,Ferraro P 2007Appl.Phys.Lett.90 041104

[11]Cui H K,Wang D Y,Wang Y X,Liu C G,Zhao J,Li Y 2011Acta Phys.Sin.60 044201(in Chinese)[崔華坤,王大勇,王云新,劉長庚,趙潔,李艷2011物理學(xué)報60 044201]

[12]Belashov A V,Petrov N V,Semenova I V,Vasyutinskii O S 2014J.Phys.C536 012003

[13]Pang T 2015J.Mod.Opt.62 816

[14]Cuche E,Marquet P,Depeursinge C 2000Appl.Opt.39 4070

[15]Cui H,Wang D,Wang Y,Zhao J,Zhang Y 2011Opt.Commun.284 4152

[16]Liu Y,Wang Z,Li J,Gao J,Huang J 2016Opt.Laser Eng.86 115

[17]Nguyen T,Bui V,Lam V,Raub C B,Chang L C,Nehmetallah G 2017Opt.Express25 15043

[18]Kim D C,Cho H J,Shin S,Jung W,Yu Y H 2009J.Opt.Soc.Korea13 451

[19]Wang H Y,Liu F F,Song X F,Liao W,Yu M J,Liu Z Q 2013Chin.J.Lasers40 196(in Chinese)[王華英,劉飛飛,宋修法,廖微,于夢杰,劉佐強 2013中國激光 40 196]

[20]Liu S,Xiao W,Pan F 2014Opt.Laser Technol.57 169

[21]Zhang D,Fan J,Zhao H,Lu X,Liu S,Zhong L 2014Optik125 5148

[22]Wang D Y,Wang Y X,Guo S,Rong L,Zhang Y Z 2014Acta Phys.Sin.63 154205(in Chinese)[王大勇,王云新,郭莎,戎路,張亦卓2014物理學(xué)報63 154205]

[23]Schnars U,Falldorf C,Watson J,Jüptner W 2015Digital Holography and Wavefront Sensing(Berlin:Heidelberg)pp39–68

[24]Bioucas-Dias J M,Valadao G 2007IEEE Trans.Image Process.16 698

[25]Durand F,Dorsey J 2002Acm.T.Graphic21 257

[26]Allen J B 1977IEEE.T.Acoust.Speech.25 235

[27]Donoho D L,Johnstone J M 1994Biometrika81 425