李芳 王明清 鄭明 盧葦 于慶南 賈燕 吳堅(jiān)

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

1 引 言

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

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

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

2 算法原理

2.1 理論模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 模擬分析和驗(yàn)證

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

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

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

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

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

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

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

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

圖6 實(shí)驗(yàn)結(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í)驗(yàn)測(cè)量和驗(yàn)證

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

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

4 結(jié) 論

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

參考文獻(xiàn)

[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)[王華英,劉飛飛,宋修法,廖薇,趙寶群,于夢(mèng)杰,劉佐強(qiáng)2013物理學(xué)報(bào)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é)報(bào) 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)[崔華坤,王大勇,王云新,劉長(zhǎng)庚,趙潔,李艷2011物理學(xué)報(bào)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)[王華英,劉飛飛,宋修法,廖微,于夢(mèng)杰,劉佐強(qiáng) 2013中國(guó)激光 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é)報(bào)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