林盼明， 羅 都， 李 軍， 黃佐華

(華南師范大學(xué)物理與電信工程學(xué)院， 廣州 510006)

相襯法是常見的定量相位分布成像方法,在透明生物醫(yī)學(xué)細(xì)胞顯微成像、透明材料、微透鏡檢測、干涉計(jì)量、流場密度分布測量等領(lǐng)域有廣泛應(yīng)用[1]. 在全場的真實(shí)相位分布重構(gòu)中，常常涉及相位的解包裹過程. 實(shí)際上，相位解包裹計(jì)算已被廣泛應(yīng)用于干涉合成孔徑雷達(dá)、基于結(jié)構(gòu)光的三維測量、光學(xué)干涉測量與成像等測量方法中[2-4]. 在從干涉圖獲取相位的過程中，一般要用反正切范數(shù)提取相位，致使連續(xù)的相位被截?cái)啵涤虮幌薅ㄔ?-π,π]，形成包裹相位. 故需要對包裹相位進(jìn)行解包裹，將截?cái)嘞辔恢貥?gòu)，以求得真實(shí)的連續(xù)相位[5]. 但實(shí)際情況下，包裹相位由于欠采樣、含有噪聲、局部陰影等影響因素，常會出現(xiàn)殘差點(diǎn)區(qū)域，影響解包裹的精度，給解包裹增加了難度[6]，也間接影響真實(shí)相位重構(gòu)的精確度.

現(xiàn)有的相位解包裹算法眾多[7-15]，主要有2種：路徑跟蹤算法和路徑無關(guān)算法. 其中，路徑跟蹤算法主要有枝切法、質(zhì)量圖導(dǎo)向法和掩膜割線法等[16-18]；路徑無關(guān)算法主要是尋求滿足最小范數(shù)的解，如基于離散余弦變換的最小二乘法、基于快速傅里葉變換的最小二乘法等[19]. 嚴(yán)重欠采樣和較強(qiáng)噪聲對解包裹來說是兩個(gè)難題，錢曉凡等[4]提出了橫向剪切最小二乘法(LSBLS)，郭媛等[20]提出了四向橫向剪切最小二乘法(FLSBLS). 這些算法可以恢復(fù)欠采樣相位圖，但仍然無法有效求解嚴(yán)重欠采樣的包裹相位圖.

近年來，自適應(yīng)閾值在圖像分割、去噪、去霧、相位轉(zhuǎn)換等方面應(yīng)用廣泛. SERT等[21]基于結(jié)構(gòu)光系統(tǒng)提出采用自適應(yīng)閾值相位轉(zhuǎn)換三級校準(zhǔn)算法進(jìn)行三維物體的測量. 在圖像處理方面，有研究者提出基于自適應(yīng)閾值移除背景的邊緣投影輪廓測定法[22]，也有將自適應(yīng)閾值用于細(xì)胞圖像處理[23-24]及圖像去噪聲[25]等. 目前，自適應(yīng)閾值在圖像解包裹方面尚未被廣泛應(yīng)用.

針對常見解包裹方法存在的不足，本文提出一種基于自適應(yīng)閾值翻轉(zhuǎn)的相位解包裹算法(SATR). 對于包裹相位被截?cái)嗟亩喾迥Ｐ?，采用程序自?dòng)尋找最佳閾值，確定相位凹陷處的邊界，得到一個(gè)只有凹陷處值的矩陣，再對該矩陣進(jìn)行翻轉(zhuǎn)操作，替換原矩陣中凹陷處的值，從而實(shí)現(xiàn)真實(shí)相位的還原. 通過程序自動(dòng)尋找自適應(yīng)閾值，能有效避免手動(dòng)設(shè)置閾值，提供正確閾值選擇，提高了相位還原的效率及準(zhǔn)確性，有較強(qiáng)的抗干擾能力.

1 研究方法

1.1 相襯成像的基本原理

相襯成像的基本思想：在含有被測物體相位信息的物光波中，通過改變零頻分量或高頻分量的相位，將物光波的相位分布調(diào)制成可以觀察的強(qiáng)度分布. 設(shè)樣品為純相位型物體，振幅為A的平行光正入射到物平面上，則物光場的復(fù)振幅可表示為：

f(x,y)=Aeiφ(x,y).

(1)

若忽略相干光學(xué)系統(tǒng)有限孔徑的影響，其頻譜為：

(2)

如果在頻譜面上加澤尼克相板，其空間濾波函數(shù)為：

(3)

其中,t與φd分別為相板的振幅透過率與相移量. 加相板后，像面的光強(qiáng)分布數(shù)學(xué)表達(dá)式可寫為[26]：

I=A2{t2+2[1-cosφ(x,y)-tcosφd+

tcos(φd-φ(x,y))]}.

(4)

要求解相位分布φ(x,y)，還需要知道不放置樣品和相板時(shí)，像面本底光強(qiáng)度I0的大小. 考慮到對實(shí)際樣品成像時(shí)，多次移動(dòng)相板會增加許多無謂的工作量，因此在不移動(dòng)相板的情況下，直接移開樣品或讓樣品中不含物體信息的區(qū)域光場通過. 此時(shí)，像面本底光強(qiáng)度可表示為：

I0=A2t2.

(5)

由I(x,y)和I0可計(jì)算出被測樣品的相位主值分布：

(6)

其中,B=(1+t2-2tcosφd)1/2，β=arctan(tsinφd/(tcosφd-1)). 根據(jù)式(6)可以確定0～π范圍內(nèi)的相位值，因此，當(dāng)被測物體相位變化較小時(shí)，即其相位分布變化范圍在0～π時(shí)，I(x,y)能直接反映物體的真實(shí)相位分布；如果物體相位超過π，拍攝2幅相襯圖像(即本底光強(qiáng)圖I0及相襯包裹干涉圖I(x,y))，由于相位被包裹，不能直接求解，這就需要利用上述反余弦關(guān)系，通過解包裹才能重構(gòu)被測相位物體的真實(shí)相位分布.

1.2 自適應(yīng)閾值翻轉(zhuǎn)相位解包裹算法

在自適應(yīng)閾值翻轉(zhuǎn)相位解包裹算法中，采用最佳閾值的方法確定包裹相位凹陷邊界、找到凹陷范圍并進(jìn)行翻轉(zhuǎn)操作. 基于翻轉(zhuǎn)法的解包裹方法要求處理包裹相位分布中的包裹層疊部分，將其移到真實(shí)區(qū)間，因此，如何尋找凹陷、確定凹陷處的邊界以及如何實(shí)現(xiàn)自動(dòng)翻轉(zhuǎn)是翻轉(zhuǎn)算法解包裹的核心.

自適應(yīng)閾值翻轉(zhuǎn)算法的基本原理是通過尋找閾值(閾值設(shè)為W，0.9≤W≤1)作為獲取凹陷邊界的依據(jù)(此處共有11個(gè)閾值). 在提取包裹相位矩陣中，相位大于Wφmax像素并將該像素區(qū)域相位設(shè)置為1，其中φmax為輸入相位矩陣的相位最大值，故標(biāo)記為1的部分為凹陷部分，而其余像素設(shè)置為0；找到相位凹陷后，對包裹相位凹陷部分進(jìn)行自動(dòng)翻轉(zhuǎn)替換；將翻轉(zhuǎn)替換后的矩陣記為φ00，進(jìn)行濾波平滑處理后將矩陣記為φ00T，保存并比較每一個(gè)閾值對應(yīng)的φ00和φ00T的相對誤差，找到最小誤差對應(yīng)的相位矩陣和閾值并輸出，此時(shí)的閾值為最佳閾值，輸出濾波平滑后的相位矩陣φ00T則為最佳解包裹相位分布. 自適應(yīng)閾值翻轉(zhuǎn)算法解包裹的程序框圖如圖1所示，其中，實(shí)現(xiàn)自動(dòng)翻轉(zhuǎn)替換算法的實(shí)現(xiàn)過程可用圖2所示的程序框圖來體現(xiàn)，其基本原理：對程序找到的只有凹陷處相位分布的矩陣進(jìn)行判斷. (1)若該矩陣為非零矩陣，則對該矩陣進(jìn)行翻轉(zhuǎn)，翻轉(zhuǎn)后用得到的矩陣替代原相位矩陣φ00中凹陷處的位置，最后進(jìn)行判斷;(2)若該矩陣為零矩陣，則終止循環(huán)，輸出所得到的經(jīng)過翻轉(zhuǎn)替代的矩陣φ00，否則將繼續(xù)循環(huán)直到找不到凹陷為止. 另外，翻轉(zhuǎn)的對稱軸在第一次為π，第二次為2π，第三次為3π，以此類推.

圖1 自適應(yīng)閾值翻轉(zhuǎn)算法(SATR)解包裹的程序框圖Figure 1 The block diagram of self-adaptive threshold reversal algorithm for phase unwrapping

圖2 自動(dòng)翻轉(zhuǎn)替換算法的程序框圖Figure 2 The block diagram of automatic flip-to-replace algorithm

2 數(shù)值仿真結(jié)果

2.1 4種相位解包裹算法精度的比較

構(gòu)建一個(gè)512 px×512 px且2倍峰值的原始相位分布，其最大相位為16.2 rad. 在不考慮噪聲的條件下，分別采用LS、LSBLS以及FLSBLS與本文SATR算法對無噪聲的二維包裹相位分布進(jìn)行解包裹.

計(jì)算LS、LSBLS、FLSBLS和SATR算法的真實(shí)相位與解包裹得到相位的均方根誤差(RMSE)分別為0.749 8、0.749 8、0.752 1和0.119 8 rad，本文SATR算法對相位恢復(fù)精確度明顯高于其他3種算法，相對相位反演精度為99.26%. 算法使用凹陷相位翻轉(zhuǎn)進(jìn)行相位恢復(fù)，每次翻轉(zhuǎn)替換得到的相位都由前一次翻轉(zhuǎn)決定，而對下一次翻轉(zhuǎn)無影響. 本算法通過設(shè)置程序?qū)ふ议撝礧，將查找到的相位大于Wφmax像素作為凹陷部分，隨后對相位凹陷部分進(jìn)行自動(dòng)翻轉(zhuǎn)替換. 將翻轉(zhuǎn)替換后的矩陣與經(jīng)過濾波平滑處理后的矩陣作相對誤差分析，找到最小誤差對應(yīng)的相位矩陣和閾值并輸出，此時(shí)的閾值為最佳閾值，通過這種方法得到恢復(fù)相位的誤差也較小.

2.2 抗噪性能比較

在實(shí)際應(yīng)用中，圖像都會存在一定的噪聲. 在原始相位分布圖像中加入不同程度的高斯噪聲后，對512 px×512 px且2倍峰值的相位分布取絕對值，然后進(jìn)行解包裹. 加入噪聲后，先對原始相位進(jìn)行濾波預(yù)處理. 解包裹后，要對恢復(fù)相位進(jìn)行均值濾波處理，得到最終的重構(gòu)相位分布. 表1為基于離散余弦變換的(DCT)最小二乘法與SATR算法解包裹的相位誤差的性能比較. 圖3是加入信噪比為3%的高斯噪聲后得到上述2種方法解包裹的實(shí)驗(yàn)結(jié)果及比較.

表1 2種相位解包裹算法的抗噪性能比較Table 1 The comparison of the anti-noise performance of two phase unwrapping algorithms

由表1和圖3可以看出：對相位加噪聲后，DCT解包裹的平均誤差和均方根誤差比SATR算法的都要大，SATR算法的解包裹精度優(yōu)于DCT，相對相位反演精度為96.40%～98.10%；SATR算法的解包裹誤差會隨著噪聲的增大而增大，但SATR算法的相位誤差始終小于DCT的解包裹誤差. 但總體上，大的噪聲越大一定程度上會影響解包裹的結(jié)果. 通過以上對比可知，SATR算法可以在一定程度上減小噪聲對相位解包裹的影響，提高了解包裹精度.

圖3 含噪聲包裹圖的2種相位解包裹方法的結(jié)果比較Figure 3 The comparison of the results of two phase unwrapping methods for enveloping diagrams with noise

此外，實(shí)驗(yàn)還表明：經(jīng)典的相位解包裹算法[27]由于迭代速度較慢而不適合網(wǎng)格大于64 px×64 px的相位解包裹運(yùn)算，而提出的閾值翻轉(zhuǎn)算法能夠適合計(jì)算量較大、網(wǎng)格大小為1 501 px×1 501 px的相位解包裹運(yùn)算，且解包裹精度可達(dá)96%.

3 實(shí)驗(yàn)驗(yàn)證

3.1 系統(tǒng)設(shè)置

為了驗(yàn)證SATR算法的實(shí)際應(yīng)用效果，搭建基于相板調(diào)制的定量相襯成像實(shí)驗(yàn)光路(圖4). 從半導(dǎo)體激光器(波長635 nm)發(fā)出的光束經(jīng)可調(diào)衰減器NF、擴(kuò)束鏡、準(zhǔn)直透鏡形成平行光波照明樣品，透鏡L2和L3(焦距均為180 mm)組成4f系統(tǒng)對樣品進(jìn)行成像，在透鏡L3的后焦面上產(chǎn)生相襯干涉圖樣. 在透鏡L3的成像面上放置CCD數(shù)字?jǐn)z像機(jī)(型號Guppy PRO F-125B，像元尺寸為3.75 μm)，對相襯干涉圖像圖樣進(jìn)行采集. 實(shí)驗(yàn)采用振幅透過率t為0.3、相移量為1.1π的自制相板.

圖4 相襯成像實(shí)驗(yàn)光路Figure 4 The experimental optical path of phase contrast imaging

3.2 實(shí)驗(yàn)分析

選用微透鏡陣列作為待測相位物體，其中單元微透鏡球面的厚度約3.52 μm，用波長為635 nm的激光平行照射該微透鏡陣列時(shí)，其球面頂部到底部的相位差約15.88 rad(5.055π).

實(shí)驗(yàn)系統(tǒng)像面上加相板后的光強(qiáng)分布如圖5所示，圖像噪聲比較嚴(yán)重. 由圖5A、B可分別得到微透鏡的背景光強(qiáng)矩陣I0和相襯光強(qiáng)矩陣I，從而計(jì)算得到一個(gè)微透鏡的包裹相位分布.

圖5 微透鏡陣列光強(qiáng)分布圖Figure 5 The profile of micro lens array intensity

將微透鏡的背景光強(qiáng)矩陣I0和相襯光強(qiáng)矩陣I代入反演公式(6)，得到二維包裹相位圖(圖6A)和三維包裹圖(圖6B)，在解包裹之前先將反余弦包裹轉(zhuǎn)化為反正切包裹分布，再采用自適應(yīng)閾值翻轉(zhuǎn)解包裹算法對其進(jìn)行處理，經(jīng)過濾波后得到還原相位(圖6C). 另采用橫向剪切最小二乘法對其進(jìn)行解包裹，經(jīng)濾波處理后得到圖6D. 通過對比圖6C和D可知，這2種算法均可恢復(fù)被測物的基本相位分布，但SATR算法解包裹得到的重構(gòu)相位分布更完整，LSBLS算法雖可以再現(xiàn)微透鏡相位的大體輪廓，但由于誤差傳遞的影響比較大，在峰頂出現(xiàn)了缺陷.

圖6 SATR算法與LSBLS算法還原微透鏡單元的結(jié)果對比Figure 6 The comparison of the results of reducting elements of micro lens between SATR and LSBLS

SATR反演得到單元微透鏡的相位，其最大值約16.35 rad，厚度約3.62 μm，與微透鏡陣列參數(shù)對比，存在約0.10 μm的絕對誤差，相對誤差為2.84%；LSBLS反演得到的單個(gè)微透鏡的相位最大值約15.22 rad，厚度約3.37 μm，與微透鏡陣列參數(shù)對比，存在約0.15 μm的絕對誤差，相對誤差為4.16%. 可見，SATR算法應(yīng)用于微透鏡單元解包裹得到的結(jié)果誤差較小，能更好地還原包裹相位，其解包裹相位與原始相位基本保持一致，是一種可行且精確度高的實(shí)用算法.

4 結(jié)論

在相位解包裹中，自適應(yīng)閾值相位翻轉(zhuǎn)恢復(fù)相位方法最大的優(yōu)點(diǎn)是利用自適應(yīng)閾值確定相位翻轉(zhuǎn)邊界，能夠以最小的相對誤差求解出完整的真實(shí)相位，從而提高了相位解包裹的精確度. 通過與幾種常見解包裹算法進(jìn)行模擬計(jì)算和實(shí)驗(yàn)驗(yàn)證效果的對比，結(jié)果表明：該算法不僅可靠、抗干擾性強(qiáng)，而且適合于大計(jì)算量及任意相位周期的解包裹，具有應(yīng)用價(jià)值.