數(shù)字差分-積分快速相位解包裹算法研究*

王子碩 劉磊2)? 劉晨博 劉珂 鐘志2) 單明廣2)?

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

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

數(shù)字全息技術(shù)是目前應(yīng)用最廣泛的定量相位成像技術(shù)之一,但是當(dāng)測(cè)量相位較大的物體時(shí),需要解包裹算法才能計(jì)算出正確的相位信息.目前,已有的解包裹算法均面臨計(jì)算量巨大、耗時(shí)的問(wèn)題.為了解決上述問(wèn)題,本文基于傅里葉變換相位恢復(fù)算法,利用復(fù)振幅相位信息的完整性,提出針對(duì)薄相位和連續(xù)大相位的基于數(shù)字差分-積分的快速直接解包裹算法.該算法首先通過(guò)基本的傅里葉變換相位恢復(fù)算法操作后,得到含有物體完整相位的復(fù)振幅信息;隨后,從中提取兩幅子復(fù)振幅信息,并將二者相除,再相位提取出其中信息,便可得到一個(gè)物體真實(shí)相位差分信息;最后,沿差分方向?qū)μ崛〉南辔徊罘中畔⑦M(jìn)行積分,便可得到解包裹后的相位信息.同時(shí)通過(guò)仿真與具體實(shí)驗(yàn)對(duì)該算法進(jìn)行了驗(yàn)證.結(jié)果表明,本文算法可以實(shí)現(xiàn)快速準(zhǔn)確的解包裹相位直接恢復(fù).

1 引言

數(shù)字全息技術(shù)是目前應(yīng)用最廣泛的定量相位成像技術(shù)之一,它具有無(wú)需標(biāo)記、大視場(chǎng)、非接觸和無(wú)損等優(yōu)點(diǎn).目前為止,國(guó)內(nèi)外眾多課題組針對(duì)數(shù)字全息技術(shù)開(kāi)展了大量的研究工作,并取得了豐富的成果[1-11].其中,相位重建算法[12,13]是數(shù)字全息領(lǐng)域的一個(gè)極重要的研究方向.

數(shù)字全息技術(shù)按照重建算法的區(qū)別可分為同軸數(shù)字全息[14-16]和離軸數(shù)字全息[9-11]兩類(lèi)全息技術(shù).同軸數(shù)字全息一般需要采集多幅同軸全息圖消除共軛像,因此耗時(shí)較長(zhǎng),而且對(duì)實(shí)驗(yàn)裝置要求很高,實(shí)驗(yàn)過(guò)程中實(shí)現(xiàn)難度較大.而離軸數(shù)字全息技術(shù)僅需一張離軸數(shù)字全息圖便可完成相位恢復(fù),使其成為最具效率的數(shù)字全息技術(shù).學(xué)者曾提出一系列基于傅里葉變換的相位恢復(fù)算法[12,13]對(duì)離軸數(shù)字全息進(jìn)行快速且準(zhǔn)確相位恢復(fù).基于傅里葉變換的相位恢復(fù)算法通常都包含二維傅里葉變換、帶通濾波、二維逆傅里葉變換和求取相位等步驟.其中帶通濾波會(huì)損失很多高頻信息,降低了成像的分辨率.因此,也有學(xué)者研究同-離軸混合數(shù)字全息恢復(fù)方法,以兼顧離軸數(shù)字全息和同軸數(shù)字全息的優(yōu)勢(shì)[17,18].但是,無(wú)論采用哪種數(shù)字全息技術(shù)進(jìn)行測(cè)量,由于相位的周期性,僅能提取包裹的相位信息.因此,那些算法恢復(fù)出的相位存在著包裹問(wèn)題,需要解包裹算法對(duì)包裹相位進(jìn)行解包裹才能得到準(zhǔn)確的相位信息.目前,解包裹算法主要分為路徑跟蹤法[19-21]、最小范數(shù)法[22-24]和基于深度學(xué)習(xí)方法[25-27]三大類(lèi).其中,枝切法是最經(jīng)典的路徑跟蹤解包裹算法,雖然已有若干改進(jìn)的基于枝切法的相位解包裹算法被提出,但是這些解包裹算法易受局部噪聲影響,會(huì)出現(xiàn)拖影的解包裹失敗的現(xiàn)象.而最小二乘算法利用全局最小化來(lái)進(jìn)行相位解包裹,該方法可以得到整幅解包裹結(jié)果,但是會(huì)平滑相位信息,降低重建結(jié)果的精度.最小二乘算法可分為無(wú)加權(quán)值[22]和有加權(quán)值[23]兩種,其中權(quán)值代表噪聲對(duì)像素的干擾程度,權(quán)值越大,干擾程度越低;反之亦然.為了提升解包裹的精度,Wei 等[24]將枝切法和最小二乘法相結(jié)合,提出了留數(shù)點(diǎn)校準(zhǔn)最小二乘解包裹算法,該方法利用留數(shù)點(diǎn)作為最小二乘算法的加權(quán)值,對(duì)含噪聲相位圖表現(xiàn)出卓越的抗噪性,但需多次迭代運(yùn)算才可求取最優(yōu)結(jié)果,運(yùn)算時(shí)間較長(zhǎng).需要注意的是,不論是路徑跟蹤法還是最小范數(shù)法,重建過(guò)程都很復(fù)雜,計(jì)算量很大,無(wú)法滿足快速相位重建的實(shí)際需求.此外,基于深度學(xué)習(xí)的相位解包裹算法發(fā)展迅速,但是這類(lèi)算法不僅需要事先采集處理大量數(shù)據(jù)作為數(shù)據(jù)集,還需要搭建網(wǎng)絡(luò)進(jìn)行長(zhǎng)時(shí)間的訓(xùn)練才能生效,增加了算法實(shí)現(xiàn)的難度,也很難對(duì)不同全息系統(tǒng)采集到的數(shù)據(jù)實(shí)現(xiàn)高質(zhì)量的解包裹.

為了解決數(shù)值解包裹算法計(jì)算量大、耗時(shí)長(zhǎng)的問(wèn)題,本文提出了一種針對(duì)薄相位和連續(xù)大相位的基于數(shù)字差分-積分的快速直接解包裹算法.該算法利用基于傅里葉變換相位恢復(fù)算法在逆傅里葉變換所得二維復(fù)數(shù)矩陣中含有完整物體相位信息的特點(diǎn),通過(guò)數(shù)字位移提取和復(fù)數(shù)相除等操作,得到一幅含有物體相位差分信息的復(fù)振幅矩陣,求取其相位信息便可得到無(wú)包裹的相位差分信息,最后,沿差分方向?qū)μ崛〉南辔徊罘中畔⑦M(jìn)行積分,便可得到解包裹后的相位信息.該算法避免了局部或全局優(yōu)化過(guò)程,因此在保證解包裹質(zhì)量的同時(shí),大幅降低了計(jì)算量,提高了重建速度.

2 理論方法

假設(shè)采集到一幅離軸數(shù)字全息圖,表達(dá)如下式:

其中,A(x,y)表示全息圖所含直流量,B(x,y) 表示全息圖調(diào)制量,φ(x,y)表示物體相位信息,k1和k2表示沿x和y方向的載波.

對(duì)全息圖做二維傅里葉變換、帶通濾波和逆傅里葉變換,表達(dá)式如下:

其中,FT[·]表示二維傅里葉變換操作,IFT[·]表示二維逆傅里葉變換操作,BPF 表示帶通濾波器.已有的傅里葉變換的相位恢復(fù)算法是直接提取復(fù)數(shù)矩陣的相位信息,得到包裹的相位信息,計(jì)算過(guò)程如下:

其 中,exp[i(k1x+k2y)] 是對(duì)無(wú)樣品全息圖進(jìn)行(2)式操作得到的系統(tǒng)所含載波信息和背景相位;ANGLE(·)表示取相位操作.

根據(jù)傅里葉變換的基本特性可知,在(2)式求取的復(fù)振幅數(shù)據(jù)所含相位信息應(yīng)為完整的物體相位信息,但是由于相位的周期性和相位提取算法的局限性,(3)式僅能恢復(fù)出(-π,π]的包裹相位φ′(x,y),必須使用解包裹算法才能完成相位重建.但是,目前已有的解包裹算法大多計(jì)算量大、耗時(shí)長(zhǎng).

為了簡(jiǎn)化解包裹過(guò)程,提高相位重建速度,在經(jīng)(2)式求出復(fù)振幅信息exp{i[φ(x,y)+k1x+k2y]}

后,并不直接求取其相位信息,而是從中提取兩幅具有一個(gè)像素位移量的尺寸略小的子復(fù)振幅信息: 從第1 列像素開(kāi)始提取得到復(fù)振幅信息

從第2 列像素開(kāi)始提取尺寸與II′相同的復(fù)振幅信息

并將I′(x,y) 與II′′(x,y) 相除,可得含待測(cè)相位沿y方向差分相位信息 ?yφ(x,y),計(jì)算如下:

對(duì)(4)式所得復(fù)振幅信息求取其相位信息:

其中exp(-ik2)是對(duì)一幅無(wú)樣品全息進(jìn)行相同操作得到的一幅參考復(fù)振幅矩陣.由于相位的周期性和相位提取算法的局限性,所求?yφ(x,y)也處于(-π,π],因此只要所求相位φ(x,y)相位中相鄰相位間差值處于(-π,π],?yφ(x,y)便能表達(dá)真實(shí)的相位變化,因此通過(guò)將(5)式計(jì)算所得?yφ(x,y)沿y方向積分,便可得真實(shí)的相位信息.

根據(jù)上述的分析,本文算法的計(jì)算流程(圖1)總結(jié)如下:

圖1 快速解包裹算法流程圖Fig.1.Algorithm flow diagram for fast phase unwrapping.

步驟1二維傅里葉變換.對(duì)全息圖進(jìn)行二維傅里葉變換,得其頻譜圖.

步驟2帶通濾波、逆傅里葉變換.對(duì)步驟1中得到的傅里葉頻譜進(jìn)行帶通濾波、逆傅里葉變化,得到相位中僅含物體相位信息的復(fù)振幅矩陣.

步驟3提取生成兩幅新的具有數(shù)字位移的子復(fù)振幅信息.

步驟4復(fù)振幅相除及相位提取.將步驟2 和步驟3 中得到的復(fù)振幅矩陣相除,得到一幅新的復(fù)振幅信息,并求出其相位信息.

步驟5積分操作.對(duì)步驟4 中得到的相位信息沿步驟3 中像素平移方向反向積分,即可得到解包裹的相位信息.

3 仿真結(jié)果及分析

為了驗(yàn)證所提方法的可行性,基于CPU i7-87 50H@2.2 GHz,32 GB 內(nèi)存的硬件平臺(tái)對(duì)該方法進(jìn)行計(jì)算機(jī)模擬驗(yàn)證.仿真生成尺寸為500 pixel×501 pixel 的Peak 的相位信息,如圖2(a)所示.使用波長(zhǎng)為532 nm 的光源,生成離軸全息圖如圖2(b)所示.使用傳統(tǒng)的復(fù)振幅相除的傅里葉變換相位恢復(fù)算法,得到恢復(fù)相位如圖2(c)所示.可見(jiàn),使用傳統(tǒng)算法恢復(fù)出的相位包裹,必須使用解包裹算法才能得到準(zhǔn)確的解包裹信息.分別使用枝切法和最小二乘法對(duì)圖2(c)進(jìn)行解包裹,得到解包裹相位為圖2(d)和圖2(e)所示.然后,使用本文算法得到恢復(fù)相位如圖2(f)所示.可見(jiàn),本文算法恢復(fù)出來(lái)的相位已經(jīng)是解包裹相位.為了顯示統(tǒng)一,僅在文中顯示尺寸為500 pixel×500 pixel 的區(qū)域.

圖2 (a) Peak 值相位型物體及其(b)離軸全息圖;(c)傅里葉變換法的恢復(fù)結(jié)果;使用(d)枝切法和 (e)最小二乘法的解包裹結(jié)果;(f)本文所提算法的恢復(fù)結(jié)果Fig.2.(a) Phase object with peak values and (b) its off-axis hologram;(c) the retrieved result of the Fourier transform;the unwrapped results by (d) branch-cut and (e) least square;(f) the retrieved result of this proposed algorithm.

將3 個(gè)恢復(fù)出的相位轉(zhuǎn)換成對(duì)應(yīng)的光程差,并提取黑色虛線標(biāo)注位置的數(shù)據(jù)展示在圖3 中,三者數(shù)據(jù)完全重合.經(jīng)計(jì)算3 種恢復(fù)算法恢復(fù)結(jié)果標(biāo)準(zhǔn)差均為0.5371 nm.

圖3 圖2(d)–(f)中黑色虛線所標(biāo)剖面數(shù)據(jù)Fig.3.One-dimensional (1D) phase profile along the black dashed lines in Figs.2 (d)–(f).

此外,為了更全面評(píng)價(jià)所提算法在相位恢復(fù)方面的速度優(yōu)勢(shì),使用枝切法、最小二乘法和該算法對(duì)多種尺寸(256 pixel×256 pixel,512 pixel×512 pixel 和1024 pixel×1024 pixel)的全息圖進(jìn)行恢復(fù),所能達(dá)到的恢復(fù)速度如表1 所列.在處理尺寸為256 pixel×256 pixel 的全息圖時(shí),相較于枝切法和最小二乘法,本文算法分別提升了7 倍和2 倍,并且隨著全息圖尺寸的增大,速度提升越大.可見(jiàn)本文算法在保證恢復(fù)質(zhì)量的同時(shí),大幅縮短計(jì)算時(shí)間.

表1 三種算法比較Table 1. Comparison of three algorithms.

上述數(shù)值仿真實(shí)驗(yàn)驗(yàn)證了本文所提算法的有效性,不僅可以直接得到解包裹的真實(shí)相位信息,而且恢復(fù)質(zhì)量和精度也與傳統(tǒng)算法一致.此外,由于僅需一幅數(shù)字位移的復(fù)振幅矩陣和一維積分等操作,整個(gè)算法所需時(shí)間遠(yuǎn)低于目前已有的單波長(zhǎng)解包裹算法.

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

使用課題組最近提出的離軸數(shù)字全息系統(tǒng)[28]采集全息圖,其中使用波長(zhǎng)為632.8 nm 的 He-Ne激光器作為光源,并使用 1280 pixel×1024 pixel(像素尺寸為 4.8 μm×4.8 μm)的 CCD 相機(jī)記錄全息圖.為了驗(yàn)證該技術(shù)的有效性,首先測(cè)量高度約為 580.22 nm(由BRUKER 原子力顯微鏡測(cè)量所得)的相位臺(tái)階物體,該相位臺(tái)階由折射率為1.5168 的BK7 玻璃制成.在本次實(shí)驗(yàn)中,采集的離軸全息圖如圖4(a)所示,使用基于傅里葉變換的相位恢復(fù)算法進(jìn)行相位重建,恢復(fù)出的相位圖如圖4(b)所示.由于引起光程差小于使用波長(zhǎng),所以并未引起包裹.使用本文所提的算法,得到恢復(fù)結(jié)果如圖4(c)所示.

圖4 (a) 全息圖;(b)傅里葉變換法恢復(fù)結(jié)果;(c)本文算法恢復(fù)結(jié)果Fig.4.(a) A hologram;(b) the retrieved result by Fourier transform;(c) the retrieved result by the proposed algorithm.

將兩個(gè)恢復(fù)結(jié)果中白色虛線標(biāo)注的相位信息提取出來(lái),計(jì)算出對(duì)應(yīng)的高度信息,一維高度剖面如圖5 所示,二者恢復(fù)結(jié)果完全重合.高度剖面的高度差為583.24 nm,與原子力顯微鏡測(cè)量結(jié)果基本吻合.此外,一維高度剖面的標(biāo)準(zhǔn)差(standard deviation,SD)分別為20.997 nm 和20.997 nm.本實(shí)驗(yàn)驗(yàn)證了該技術(shù)可以實(shí)現(xiàn)精確測(cè)量.

圖5 圖4(b)和圖4(c)中白色虛線所標(biāo)剖面數(shù)據(jù)Fig.5.1D phase profile along the black dashed lines in Fig.4(b) and (c).

為了進(jìn)一步驗(yàn)證本文算法對(duì)連續(xù)大相位物體的相位重建能力,使用本文算法對(duì)酒精蒸發(fā)過(guò)程中某一時(shí)刻的狀態(tài)進(jìn)行處理.圖6(a)顯示了記錄的酒精蒸發(fā)過(guò)程中的全息圖之一,圖6(b),(c)顯示了經(jīng)過(guò)枝切法和最小二乘解包裹后的兩個(gè)相位信息,圖6(d)顯示了本文算法直接恢復(fù)的相位信息.為了更好地驗(yàn)證本文算法的有效性,將圖6(b)–(d)中白色虛線標(biāo)注的數(shù)據(jù)提取出來(lái),展示在圖7中,二者高度一致.

圖7 圖6(b)–(d)中白色虛線所標(biāo)剖面數(shù)據(jù)Fig.7.1D phase profile along the black dashed lines in Fig.6 (b)-(d).

通過(guò)上述兩個(gè)實(shí)驗(yàn)驗(yàn)證了所提算法能夠?qū)唧w實(shí)驗(yàn)中所采集到富含噪聲的全息圖進(jìn)行有效正確的相位解包裹,充分驗(yàn)證了本文算法的有效性和抗噪能力.

5 結(jié)論

為了解決單波長(zhǎng)數(shù)字全息解包裹算法計(jì)算量大、速度慢等問(wèn)題,本文利用復(fù)振幅中相位信息完整性的特點(diǎn),通過(guò)數(shù)字位移結(jié)合復(fù)振幅相除,提出數(shù)字差分-積分的快速解包裹算法.本文從數(shù)值仿真和具體實(shí)驗(yàn)兩方面驗(yàn)證了所提算法,不僅恢復(fù)質(zhì)量和精度與已有算法高度一致,而且恢復(fù)速度也得到了大幅提升.并且,本文算法簡(jiǎn)單便于實(shí)現(xiàn),為數(shù)字全息的快速、高質(zhì)量重建提供有利的條件.