結(jié)合線(xiàn)性回歸的離軸數(shù)字全息去載波相位恢復(fù)算法*

單明廣 劉翔宇 龐成 鐘志 于蕾 劉彬 劉磊?

1) (哈爾濱工程大學(xué)信息與通信工程學(xué)院，哈爾濱 150001)

2) (哈爾濱工程大學(xué)，先進(jìn)船舶通信與信息技術(shù)工信部重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150001)

為實(shí)現(xiàn)僅用一幅離軸數(shù)字全息圖便能直接恢復(fù)相位，提出一種利用空間載波相移技術(shù)(spatial carrier phase shift，SCPS)和線(xiàn)性回歸相結(jié)合的離軸數(shù)字全息去載波相位恢復(fù)算法.首先，利用SCPS 將一幅離軸數(shù)字全息圖分為四幅含有載波相移的全息圖，其中載波相移由沿行、列兩個(gè)方向的正交載波所引入;然后，將四幅載波相移全息圖作為輸入，將所求物體相位和兩個(gè)正交的載波作為未知量，結(jié)合最小二乘法和線(xiàn)性回歸同時(shí)求出載波和相位信息.相較于已有的去載波技術(shù)，本算法無(wú)需背景全息圖作為參考，便可準(zhǔn)確地去除載波，實(shí)現(xiàn)高質(zhì)量的相位重建.本文結(jié)合數(shù)值仿真和具體實(shí)驗(yàn)結(jié)果驗(yàn)證本算法的有效性和優(yōu)越性.

1 引言

數(shù)字全息技術(shù)[1-6]因其非接觸、分辨率高等優(yōu)點(diǎn)，在光學(xué)測(cè)量領(lǐng)域得到了廣泛的應(yīng)用.按照恢復(fù)算法，可將數(shù)字全息分為同軸數(shù)字全息[2，3]和離軸數(shù)字全息[4-6]兩類(lèi).其中，同軸數(shù)字全息可充分利用相機(jī)的空間帶寬利用率，具有較高的分辨率，但是通常需要犧牲視場(chǎng)利用率[3，7]或時(shí)間利用率[8，9]采集兩幅及以上的全息圖，并且需要復(fù)雜的相移設(shè)備引入相移.而離軸數(shù)字全息技術(shù)雖然對(duì)相機(jī)的空間帶寬利用率不高，但是僅需一幅全息圖便可恢復(fù)物體信息[10-12]，提高了系統(tǒng)的測(cè)量效率，也降低了環(huán)境對(duì)多幀測(cè)量的影響.

傅里葉變換法(Fourier transform algorithm，F(xiàn)TA)[10]是離軸數(shù)字全息最常用的相位恢復(fù)算法.該方法簡(jiǎn)單快速，但恢復(fù)過(guò)程中使用的帶通濾波會(huì)造成頻譜泄露，同時(shí)帶通濾波器截?cái)囝l譜的操作會(huì)在恢復(fù)結(jié)果引入吉布斯效應(yīng)，從而影響到恢復(fù)質(zhì)量;同時(shí)，F(xiàn)TA 需精確確定實(shí)像頻譜信息強(qiáng)度峰值的位置才能有效地去除載波，但由于全息圖的數(shù)字化特點(diǎn)，實(shí)像頻譜信息強(qiáng)度峰值很難恰好位于整數(shù)像素位置，從而造成載波殘留等影響恢復(fù)質(zhì)量的問(wèn)題.為去除載波的影響，目前最常用也是最有效的方法是雙曝光法(double exposure algorithm，DEA)[10，13]，即在實(shí)驗(yàn)之前采集一幅不含樣品的全息圖，以此記錄系統(tǒng)的詳細(xì)載波信息，通過(guò)簡(jiǎn)單的相減或相除便可準(zhǔn)確地去除載波，但是DEA 要求系統(tǒng)具有極高的穩(wěn)定性，保證系統(tǒng)參數(shù)在測(cè)量過(guò)程中保持不變;同時(shí)，某些反射式全息系統(tǒng)的無(wú)樣品全息圖很難采集，從而加大了實(shí)際使用的難度.為了實(shí)現(xiàn)亞像素級(jí)的載波估計(jì)，F(xiàn)an 等[14]提出了質(zhì)心法(spectrum centroid method，SCM)，通過(guò)計(jì)算峰值所在部分區(qū)域的質(zhì)心估算載波，但是SCM 的估算精度嚴(yán)重依賴(lài)所選區(qū)域，其估算精度在實(shí)際應(yīng)用中嚴(yán)重受限.Du 等[15]和Hincapié-Zuluaga 等[16]提出了補(bǔ)零法(zero padding，ZP)，通過(guò)對(duì)全息圖補(bǔ)零來(lái)提高載波估算精度，但是估算精度與計(jì)算量成正比，計(jì)算量巨大;同時(shí)，還需反復(fù)嘗試才能找到合適的補(bǔ)零倍數(shù)，使得該算法在實(shí)際應(yīng)用中難以取得良好效果.近期，基于全息圖翻轉(zhuǎn)變換[17-19]的去載波算法得到了極大的發(fā)展，但是該類(lèi)算法對(duì)物體的形狀和在全息圖中的分布位置有著嚴(yán)格的要求，并不具有普適性，制約了該類(lèi)算法的適用范圍.

為了避免帶通濾波操作對(duì)相位恢復(fù)造成影響，空間載波相移技術(shù)(spatial carrier phase shift，SCPS)被應(yīng)用到離軸全息的相位恢復(fù)領(lǐng)域.SCPS 是指將一幅離軸全息圖在空域沿一個(gè)方向進(jìn)行像素移位，獲取兩幅以上全息圖，并將提取出的全息圖按照相移全息圖對(duì)待，使用同軸數(shù)字全息的恢復(fù)算法計(jì)算出含有載波和物體相位的相位信息，再去除其中載波信息，最終得到準(zhǔn)確的物體相位信息.Guo 等[20]在SCPS 的基礎(chǔ)上，提出了采用最小二乘的相位恢復(fù)算法，該算法首先選取局部區(qū)域并估算其所含有的載頻信息，再通過(guò)最小二乘法計(jì)算物體相位，但是該算法的恢復(fù)質(zhì)量與局部區(qū)域的位置選擇和尺寸密切相關(guān)，而如何準(zhǔn)確地選擇合適的區(qū)域限制了該算法使用.Stykyu 和Patorski[21]指出該方法并沒(méi)有解決背景和調(diào)制度造成的誤差，提取的相位是不穩(wěn)定的.Xu 等[22]針對(duì)Guo 等提出算法的缺點(diǎn)，提出了基于最小二乘的空間載波相移算法，提取了高精度的相位，并且提高了算法的穩(wěn)定性.Huang等[23]提出在頻譜中初步估算峰值信息以減少迭代次數(shù).Liu 等[24]基于SCPS，利用利薩如圖形和橢圓擬合實(shí)現(xiàn)了相位恢復(fù)，提升了算法的抗噪能力和恢復(fù)質(zhì)量.由此可知，相較于相移法，SCPS 僅需一幅全息圖便可完成相位恢復(fù)，并可避免吉布斯效應(yīng)和頻譜泄露等問(wèn)題.但是上述基于SCPS 的相位恢復(fù)算法仍需要準(zhǔn)確的載波信息才能恢復(fù)出物體相位信息.綜上，F(xiàn)TA 和SCPS 在相位恢復(fù)時(shí)，均面臨著如何準(zhǔn)確、有效地去除載波的問(wèn)題.

因此，基于SCPS，本文提出一種結(jié)合線(xiàn)性回歸的離軸數(shù)字全息去載波相位恢復(fù)算法，實(shí)現(xiàn)從一幅離軸全息圖同時(shí)直接恢復(fù)出載波與物體相位.本算法將物體相位、載波信息和全息圖的直流量與調(diào)制量均作為未知量，結(jié)合線(xiàn)性回歸等技術(shù)，實(shí)現(xiàn)載波和物體相位的同時(shí)直接恢復(fù).本文將給出本算法的原理及計(jì)算過(guò)程，通過(guò)仿真比較幾種相位恢復(fù)算法的性能，并使用實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證本算法的有效性和優(yōu)越性.

2 理論方法

一般情況下，離軸數(shù)字全息圖可以用下式表示:

其中A(x，y)，B(x，y)分別表示全息圖的直流量和調(diào)制量;φ(x，y)代表待測(cè)物體相位信息;kx，ky分別代表全息圖中沿x方向、y方向的載波頻率.

在實(shí)際全息圖中，相較于載波，直流量A、調(diào)制量B和相位信息φ的變化非常緩慢.將I(x，y)在空域上進(jìn)行移位，即下移、右移以及右下移各一個(gè)像素，從而得到四幅相移全息圖:

其中由于直流量、調(diào)制量和相位變化緩慢，A2，A3，A4和A近似相等，B2，B3，B4和B近似相等，φ2，φ3，φ4和φ近似相等.因此，經(jīng)過(guò)空間移位，從一幅離軸全息圖提取出四幅含有載波相移的全息圖，其中相移值分別為0，kx，ky，kx+ky.將(2)式整理為一般表達(dá)式，如下:

式中，下標(biāo)m=1，2，3，4 代表圖像索引，上標(biāo)t代表理論值，δxm和δym代表沿x，y兩個(gè)方向的相移量.設(shè)沿行、列方向的載波相移分別為Δym(y)，Δxm(x)，可以得到:

設(shè)坐標(biāo)(x，y)的總載波相移為Δm(x，y)，為了進(jìn)一步闡述載波相移和沿x，y兩個(gè)方向的載波和相移之間的關(guān)系，Δm(x，y)可由下式表達(dá):

在完成上述基本分析之后，接下來(lái)將詳細(xì)地描述算法流程，該算法的每一次迭代都包含三個(gè)步驟:

步驟1利用行、列的載波相移計(jì)算相位φ

定義a(x，y)A(x，y)，b(x，y)B(x，y)cos[φ(x，y)]，c(x，y)-B(x，y)sin[φ(x，y)]，則(3)式可改寫(xiě)為

步驟1 中，假設(shè)kxm，kym和相移δxm，δym是已知的，因此Δm(x，y)是已知的，而a(x，y)，b(x，y)，c(x，y)是所要求的量，根據(jù)最小二乘法[25]，可得

其中M代表干涉圖的數(shù)量;Δm代表第m張干涉圖的載波相移;Im代表第m張全息圖的實(shí)際強(qiáng)度.

則物體相位信息可通過(guò)下式計(jì)算:

步驟2用物體相位和行載波相移確定列載波相移，進(jìn)而使用線(xiàn)性回歸算法計(jì)算列載波kx

由(3)式和(4)式可知，列載波相移Δxm(x)、行載波相移Δym(y)分別代表著x，y方向的載波相移量，二者相互獨(dú)立，且具有各自對(duì)應(yīng)的線(xiàn)性函數(shù).步驟2 將利用相位φ(x，y)和行載波相移Δym(y)計(jì)算列載波相移Δxm(x).

式中，Y代表總圖像矩陣第x列的長(zhǎng)度，即圖像矩陣的行高.

則第x列載波相移Δxm(x)可以通過(guò)下式計(jì)算:

在離軸全息圖中，載波頻率通常都較大，上式通過(guò)反正切求取Δxm(x)必然出現(xiàn)包裹.又因?yàn)棣m(x)是關(guān)于x的線(xiàn)性函數(shù)，因此，可以通過(guò)線(xiàn)性回歸[26]得到載波頻率kxm，即

最后，對(duì)所求得的kxm求取均值，得到列載波頻率kx.

步驟3用物體相位和列載波相移確定行載波相移，進(jìn)而使用線(xiàn)性回歸算法計(jì)算行載波ky

X代表總圖像矩陣第y行的長(zhǎng)度，即圖像矩陣的列寬.第y行載波相移Δym(y)可以通過(guò)下式?jīng)Q定:

與步驟2 同理，對(duì)求得的Δym(y)通過(guò)線(xiàn)性回歸得到載波頻率kym，最后對(duì)kym求取均值，得到列載波頻率ky.

上述三個(gè)步驟構(gòu)成了一次迭代過(guò)程.在迭代次數(shù)達(dá)到預(yù)設(shè)值或是恢復(fù)結(jié)果達(dá)到預(yù)設(shè)精度，則計(jì)算結(jié)束.本文中，將迭代過(guò)程中前后兩次求得的kx和ky的誤差值作為迭代是否停止的標(biāo)準(zhǔn)，計(jì)算過(guò)程如下:

其中閾值δ是預(yù)設(shè)值;i代表迭代次數(shù).在本文仿真及實(shí)驗(yàn)中δ被設(shè)置為10—6.

為了更清楚地描述恢復(fù)過(guò)程，每次迭代的流程如圖1 所示.

圖1 算法的基本流程Fig.1.Flowchart of the proposed algorithm.

最后，對(duì)本算法進(jìn)行總結(jié)，算法流程如下:

1) 使用SCPS 技術(shù)，從一幅離軸全息圖中提取四幅載波相移全息圖;

2) 將兩個(gè)方向的載波相移作為已知量，通過(guò)最小二乘法求出物體相位;

3) 將物體相位、行載波相移作為已知量，通過(guò)最小二乘法和線(xiàn)性回歸求出列載波相移和列載波;

4) 將物體相位、列載波相移作為已知量，通過(guò)最小二乘法和線(xiàn)性回歸求出行載波相移和行載波;

5) 重復(fù)步驟(2)—(4)直至滿(mǎn)足迭代停止條件.

當(dāng)滿(mǎn)足迭代條件計(jì)算結(jié)束時(shí)，便可同時(shí)得到載波信息和物體相位信息.由上述算法描述可知，不同于目前已有的最小二乘相位恢復(fù)算法，本算法利用線(xiàn)性回歸技術(shù)建立了載波相移與載波之間的聯(lián)系，實(shí)現(xiàn)了待測(cè)物體相位信息和載波信息的同時(shí)直接獲取.

3 仿真實(shí)驗(yàn)

為了驗(yàn)證本文所提算法的有效性和優(yōu)越性，首先利用仿真對(duì)本算法進(jìn)行驗(yàn)證，同時(shí)與FTA，SCM，ZP 和DEA 進(jìn)行比較，仿真過(guò)程中使用波長(zhǎng)為632.8 nm 的光源.首先，仿真生成一個(gè)光程差為100 nm 的相位型半球，像素大小為256×256，像素尺寸為4.4 μm;其次，引入沿x方向載波kx=0.76，沿y方向載波ky=0.64，從而生成一幅離軸全息圖，如圖2(a)所示.分別使用FTA，SCM，10 倍ZP，20 倍ZP，DEA 和本算法對(duì)該全息圖進(jìn)行相位恢復(fù)，得到6 個(gè)恢復(fù)結(jié)果如圖2(b)—(f)所示.由仿真結(jié)果可知，6 個(gè)恢復(fù)結(jié)果中，僅有本算法和DEA 能去除載波，得到準(zhǔn)確的相位恢復(fù)結(jié)果.

圖2 (a) 光程差100 nm 相位型半球生成的離軸全息圖;(b)—(f) 使用FTA，SCM，10 倍ZP，20 倍ZP，DEA 和本算法的恢復(fù)結(jié)果Fig.2.(a) Off-axis hologram of a phase hemisphere with 100 nm optical path difference;(b)—(f) retrieved phase maps by FTA，SCM，ZP with 10 times zero-padding，ZP with 20 times zero-padding，DEA and the proposed algorithm.

經(jīng)多次仿真實(shí)驗(yàn)可知，本算法對(duì)于載波數(shù)值的計(jì)算精度可達(dá)到10—7級(jí);而ZP 法欲達(dá)到10—7的精度，需要將圖像進(jìn)行107倍的補(bǔ)零，數(shù)據(jù)量已經(jīng)遠(yuǎn)遠(yuǎn)超過(guò)普通計(jì)算機(jī)的內(nèi)存容量，致使該算法無(wú)法運(yùn)行;同時(shí)，SCM 法的定位精度嚴(yán)重依賴(lài)于所選區(qū)域，很難在實(shí)際應(yīng)用中取得良好的效果.因此，雖然本算法的迭代過(guò)程耗時(shí)較長(zhǎng)，但是不需要提前采集背景全息圖就能準(zhǔn)確地去除載波，具有極高地實(shí)用價(jià)值.

為了驗(yàn)證本算法的抗噪能力，向生成全息圖中加入均值為0、方差為0.005 的高斯白噪聲，進(jìn)行仿真驗(yàn)證.鑒于DEA 算法的優(yōu)越性，僅使用其和本算法對(duì)含噪聲的全息圖進(jìn)行恢復(fù)，得到恢復(fù)結(jié)果分別如圖3(a)和圖3(b)所示.為了便于分析比較，求出恢復(fù)結(jié)果相較于初始值的殘差，如圖3(c)和圖3(d)所示，并將殘差圖的峰谷值(peak-valley，PV)和標(biāo)準(zhǔn)差(standard deviation，SD)標(biāo)注在殘差圖的右上角.從圖中可以看出，相較于DEA，本算法的PV 和SD 更小，具有更好的恢復(fù)質(zhì)量.

圖3 確定高斯白噪聲影響下的(a)本算法和(b)DEA 的恢復(fù)結(jié)果，以及(c)本算法和(d)DEA 的殘差圖Fig.3.Retrieved phase maps from the hologram with Gaussian white noise by (a) the proposed algorithm and (b)DEA，and the corresponding residue maps by (c) the proposed algorithm and (d) DEA.

為了更深入地對(duì)比本算法和DEA 的抗噪能力，使用本算法和DEA 對(duì)含有0 均值和不同標(biāo)準(zhǔn)差高斯白噪聲的全息圖進(jìn)行相位恢復(fù)，分別計(jì)算兩種算法恢復(fù)結(jié)果與初始值間殘差的PV 和SD，獲得結(jié)果如圖4 所示.從圖4 中可以看出，隨著噪聲的放大，兩種算法恢復(fù)結(jié)果的PV 和SD 都在增大，標(biāo)志著恢復(fù)質(zhì)量也在逐漸降低.但是，即使噪聲方差增至0.025 的情況下，本算法仍能較高質(zhì)量的恢復(fù)出物體相位信息，從而證明了本算法具有較強(qiáng)的抗噪干擾能力.

圖4 不同噪聲情況下本算法和DEA 算法對(duì)應(yīng)殘差值的PV 和SDFig.4.PV and SD of the residue maps by the proposed algorithm and DEA with different variance.

由此得知，本文提出的算法在不需要任何背景全息圖的條件下，僅僅利用單幅離軸數(shù)字全息圖即可在去除載波的同時(shí)，實(shí)現(xiàn)待測(cè)相位準(zhǔn)確直接恢復(fù).

4 實(shí)驗(yàn)結(jié)果及分析

為了進(jìn)一步驗(yàn)證所提算法的有效性，利用FTA，10 倍ZP，SCM，DEA 和本算法分別對(duì)實(shí)際獲取的一幅全息圖進(jìn)行恢復(fù).在實(shí)驗(yàn)中，使用分光瞳離軸數(shù)字全息系統(tǒng)[27]測(cè)量一塊反射式的硅片樣品，所測(cè)得的全息圖如圖5(a)所示，相位恢復(fù)結(jié)果分別如圖5(b)—(f)所示.從圖5(b)—(f)中的恢復(fù)結(jié)果可知，如同仿真結(jié)果，F(xiàn)TA，10 倍ZP 和SCM 均無(wú)法準(zhǔn)確實(shí)現(xiàn)相位恢復(fù)，而只有DEA 和本算法得到了較為準(zhǔn)確的恢復(fù)結(jié)果.為了進(jìn)一步驗(yàn)證本算法的恢復(fù)質(zhì)量，對(duì)圖5(e)和圖5(f)實(shí)線(xiàn)處做一維剖面，其剖面數(shù)據(jù)如圖6 所示.本算法高、低位置的標(biāo)準(zhǔn)差分別為0.0613 rad 和0.1218 rad，而DEA 高、低位置的標(biāo)準(zhǔn)差分別為0.1051 rad 和0.1303 rad.由此可見(jiàn)，本算法在不需要額外記錄背景全息圖的情況下，就能達(dá)到和DEA 相近的恢復(fù)質(zhì)量.

圖5 硅片實(shí)驗(yàn)結(jié)果 (a)全息圖;利用(b) FTA，(c)10 倍ZP，(d) SCM，(e) DEA 和(f)本算法的恢復(fù)結(jié)果Fig.5.Experimental results for silicon wafer:(a) Hologram;retrieved phase maps by (b)FTA，(c)ZP with 10 times zero-padding，(d) SCM，(e) DEA and (f) the proposed algorithm.

圖6 圖5(e)和圖5(f)中白線(xiàn)所標(biāo)剖面數(shù)據(jù)Fig.6.1D phase profile along the white lines in Fig.5(e)and Fig.5(f).

上述實(shí)驗(yàn)已經(jīng)驗(yàn)證了本算法在邊緣突變物體上的有效性，接下來(lái)將對(duì)酒精蒸發(fā)過(guò)程中的一幀數(shù)據(jù)進(jìn)行處理，以驗(yàn)證本算法對(duì)于邊緣連續(xù)的物體仍具有有效性.圖7(a)為酒精蒸發(fā)過(guò)程中采集到的一幀數(shù)據(jù)，圖7(b)—(f)是分別使用FTA，10 倍ZP，SCM，DEA 和本算法的恢復(fù)結(jié)果.可見(jiàn)，相比于其他算法，本算法和DEA 仍能更準(zhǔn)確的恢復(fù)出物體相位信息.為了更直觀地展現(xiàn)本算法的恢復(fù)質(zhì)量，將圖7(e)和圖7(f)中白線(xiàn)標(biāo)記的剖面數(shù)據(jù)提取出來(lái)，如圖8 所示.從圖7(b)—(f)所示的恢復(fù)結(jié)果可以看出，本算法可以實(shí)現(xiàn)連續(xù)形貌高質(zhì)量恢復(fù).

圖7 酒精實(shí)驗(yàn)結(jié)果 (a)全息圖;利用(b) FTA，(c)10 倍ZP，(d) SCM，(e) DEA 和(f)本算法的恢復(fù)結(jié)果Fig.7.Experimental results for alcohol:(a) Hologram;retrieved phase maps by (b)FTA，(c)ZP with 10 times zero-padding，(d) SCM，(e) DEA and (f) the proposed algorithm.

圖8 圖7(e)和(f)中白線(xiàn)標(biāo)注一維剖面圖Fig.8.1D phase profiles along the white lines in Fig.7(e)and Fig.7(f).

5 結(jié)論

針對(duì)離軸數(shù)字全息相位恢復(fù)存在的問(wèn)題，基于SCPS，本文基于SCPS 提出了一種結(jié)合線(xiàn)性回歸的直接提取載波與物體相位的相位恢復(fù)方法.該算法利用SCPS，將一幅離軸全息圖分解成4 幅含有載波相移的全息圖，再利用線(xiàn)性回歸建立起載波和相移之間的關(guān)系，最后通過(guò)迭代過(guò)程恢復(fù)出物體相位信息和沿x，y方向的載波信息.仿真與實(shí)驗(yàn)表明，提出的方法不需要背景全息圖作參考，即可精確地提取載波信息和物體相位信息，且具有與最常用的DEA 相近的恢復(fù)效果.該方法不僅不需要先驗(yàn)信息，而且相較于以往的去載波方法具有更高的精度，對(duì)于離軸數(shù)字全息的實(shí)際應(yīng)用具有重要的意義.