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

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

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

2) (哈爾濱工程大學,先進船舶通信與信息技術工信部重點實驗室,哈爾濱 150001)

數字全息技術是目前應用最廣泛的定量相位成像技術之一,但是當測量相位較大的物體時,需要解包裹算法才能計算出正確的相位信息.目前,已有的解包裹算法均面臨計算量巨大、耗時的問題.為了解決上述問題,本文基于傅里葉變換相位恢復算法,利用復振幅相位信息的完整性,提出針對薄相位和連續大相位的基于數字差分-積分的快速直接解包裹算法.該算法首先通過基本的傅里葉變換相位恢復算法操作后,得到含有物體完整相位的復振幅信息;隨后,從中提取兩幅子復振幅信息,并將二者相除,再相位提取出其中信息,便可得到一個物體真實相位差分信息;最后,沿差分方向對提取的相位差分信息進行積分,便可得到解包裹后的相位信息.同時通過仿真與具體實驗對該算法進行了驗證.結果表明,本文算法可以實現快速準確的解包裹相位直接恢復.

1 引言

數字全息技術是目前應用最廣泛的定量相位成像技術之一,它具有無需標記、大視場、非接觸和無損等優點.目前為止,國內外眾多課題組針對數字全息技術開展了大量的研究工作,并取得了豐富的成果[1-11].其中,相位重建算法[12,13]是數字全息領域的一個極重要的研究方向.

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

為了解決數值解包裹算法計算量大、耗時長的問題,本文提出了一種針對薄相位和連續大相位的基于數字差分-積分的快速直接解包裹算法.該算法利用基于傅里葉變換相位恢復算法在逆傅里葉變換所得二維復數矩陣中含有完整物體相位信息的特點,通過數字位移提取和復數相除等操作,得到一幅含有物體相位差分信息的復振幅矩陣,求取其相位信息便可得到無包裹的相位差分信息,最后,沿差分方向對提取的相位差分信息進行積分,便可得到解包裹后的相位信息.該算法避免了局部或全局優化過程,因此在保證解包裹質量的同時,大幅降低了計算量,提高了重建速度.

2 理論方法

假設采集到一幅離軸數字全息圖,表達如下式:

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

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

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

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

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

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

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

從第2 列像素開始提取尺寸與II′相同的復振幅信息

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

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

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

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

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

步驟1二維傅里葉變換.對全息圖進行二維傅里葉變換,得其頻譜圖.

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

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

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

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

3 仿真結果及分析

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

圖2 (a) Peak 值相位型物體及其(b)離軸全息圖;(c)傅里葉變換法的恢復結果;使用(d)枝切法和 (e)最小二乘法的解包裹結果;(f)本文所提算法的恢復結果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 個恢復出的相位轉換成對應的光程差,并提取黑色虛線標注位置的數據展示在圖3 中,三者數據完全重合.經計算3 種恢復算法恢復結果標準差均為0.5371 nm.

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

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

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

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

4 實驗結果及分析

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

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

將兩個恢復結果中白色虛線標注的相位信息提取出來,計算出對應的高度信息,一維高度剖面如圖5 所示,二者恢復結果完全重合.高度剖面的高度差為583.24 nm,與原子力顯微鏡測量結果基本吻合.此外,一維高度剖面的標準差(standard deviation,SD)分別為20.997 nm 和20.997 nm.本實驗驗證了該技術可以實現精確測量.

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

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

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

通過上述兩個實驗驗證了所提算法能夠對具體實驗中所采集到富含噪聲的全息圖進行有效正確的相位解包裹,充分驗證了本文算法的有效性和抗噪能力.

5 結論

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