基于經(jīng)驗(yàn)?zāi)B(tài)分解的激光散斑噪聲抑制方法

定量相位成像（QPI）［1-2］是一種新興的無(wú)標(biāo)記光學(xué)成像方法，可以實(shí)現(xiàn)反射樣品三維形貌的可視化，測(cè)量透明和半透明樣品的內(nèi)部結(jié)構(gòu)和折射率分布，其在透明生物樣品上的應(yīng)用是當(dāng)前生物光學(xué)成像的研究熱點(diǎn)之一。定量相位成像技術(shù)主要依賴全息［3］和干涉技術(shù)記錄透明物體的相位信息，或是記錄光場(chǎng)通過(guò)透明物體的衍射強(qiáng)度圖，從記錄的強(qiáng)度圖中提取物體的相位。它可以與壓縮感知［4］、深度學(xué)習(xí)［5-7］、自動(dòng)對(duì)焦［8］、拓?fù)湔{(diào)制［9］或補(bǔ)償技術(shù)［10］相結(jié)合，實(shí)現(xiàn)更高精度的測(cè)量結(jié)果。當(dāng)用激光照射待測(cè)樣品時(shí)，物體表面相對(duì)于激光波長(zhǎng)量級(jí)而言相當(dāng)粗糙，由于激光的高相干性，各個(gè)物點(diǎn)的散射光之間可能發(fā)生干涉，從而產(chǎn)生呈顆粒狀無(wú)規(guī)則分布的散斑噪聲。其存在會(huì)顯著降低成像結(jié)果的分辨率、減少圖像的信噪比，無(wú)法滿足高精度相位重建的要求，因此有必要對(duì)激光散斑噪聲進(jìn)行抑制。

目前，常用的散斑噪聲抑制的方法主要分成兩類［11］，一類是添加或替換硬件降低光源的相干性，如使用LED等低相干性的光源［12］等。Farrokhi［13］等人采用一動(dòng)一靜雙散射片系統(tǒng)來(lái)獲得高速散斑場(chǎng)照明，以抑制定量相位成像中的相干散斑。李煊［14］等人提出一種基于旋轉(zhuǎn)雙散射片的全息散斑噪聲抑制方法，其光學(xué)信號(hào)的非相關(guān)周期大，為獲得大量非相關(guān)全息圖求均值去噪提供了依據(jù)?？傮w來(lái)說(shuō)，這類硬件方法對(duì)光程差進(jìn)行精確控制的要求更高，增加了搭建光路的復(fù)雜度和難度。另外一類則是通過(guò)圖像處理手段實(shí)現(xiàn)散斑噪聲抑制。Uzan等人［15］對(duì)位于圖像塊中心的有噪聲像素采用非局部均值（NLM）濾波，對(duì)散斑噪聲起到了較好的抑制效果。但該方法在去噪的同時(shí)平滑了部分邊緣處的信息，會(huì)導(dǎo)致相位出現(xiàn)誤差。劉吉等人［16］提出一種基于正余弦分解結(jié)合自適應(yīng)全變分的去噪方法，該方法提高了圖像的峰值信噪比，減少了去噪后相位圖的波動(dòng)性。但該方法中自適應(yīng)參數(shù)的選取是圖像去噪和保留邊緣信息的關(guān)鍵，需要通過(guò)縮小自適應(yīng)參數(shù)的取值范圍來(lái)達(dá)到最優(yōu)去噪的重建效果。Montresor等人［17］應(yīng)用殘差網(wǎng)絡(luò)（DNCNN）對(duì)數(shù)字全息干涉測(cè)量的相位圖進(jìn)行去噪，具有非高斯統(tǒng)計(jì)和非平穩(wěn)特性，并表現(xiàn)出空間相關(guān)長(zhǎng)度和出色的去噪性能。但該方法需要使用大量的含散斑噪聲的全息圖數(shù)據(jù)集對(duì)多等級(jí)神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，實(shí)驗(yàn)條件復(fù)雜。牛瑞等人［18］提出利用二維高斯窗口對(duì)包裹相位進(jìn)行散斑噪聲抑制，不僅保留了包裹相位的跳變邊緣，還提高了相位重建的精度。但該方法中二維高斯窗口大小以及閾值會(huì)影響相位重建精度，需要研究最佳窗口大小及閾值。吳育民等人［19］提出一種基于Canny算子改進(jìn)各向異性擴(kuò)散（P-M）方程的全息散斑抑制方法，在去噪的同時(shí)有利于保留圖像的細(xì)節(jié)信息。但該方法的再現(xiàn)圖像仍然存在輕微的二次干涉條紋，且圖像相較于原始圖像而言存在一定的模糊。這類圖像處理手段雖然在定量相位成像中達(dá)到了一定的散斑抑制效果，但對(duì)于相位重建精度還有待提高。

針對(duì)以上問(wèn)題，本文提出了一種結(jié)合經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）的Canny算子改進(jìn)P-M方程的定量相位成像散斑噪聲抑制方法。首先對(duì)所記錄的全息圖進(jìn)行EMD得到一系列頻率由高到低的固有模態(tài)函數(shù)（IMF）分量，僅用高頻的IMF分量重構(gòu)圖像以增強(qiáng)圖像細(xì)節(jié)信息。隨后引入Canny算子對(duì)細(xì)節(jié)信息突出的重構(gòu)圖像進(jìn)行邊緣檢測(cè)，從而更好地控制擴(kuò)散去噪過(guò)程，提高了P-M方程抑制散斑噪聲的性能，并通過(guò)實(shí)驗(yàn)驗(yàn)證了該散斑噪聲抑制方法的可行性。

2 基本原理

2.1　結(jié)合EMD的Canny算子改進(jìn)P-M方程法在全息散斑噪聲抑制中的應(yīng)用

在全息圖中的散斑噪聲是一種乘性噪聲［20-21］，其數(shù)學(xué)模型如式（1）所示：

H'(x,y)=H(x,y)N(x,y)

 ， （1）

其中：H'(x,y)為電荷耦合器件（CCD）所記錄的圖像，其包含著散斑噪聲；H(x,y)為不含噪聲的理想全息圖；N(x,y)為散斑噪聲。

結(jié)合EMD的Canny算子改進(jìn)P-M方程法首先將圖像H'(x,y)經(jīng)EMD處理，得到一系列頻率從高到低的IMF分量。通過(guò)這些不同頻率的IMF進(jìn)行選取并重構(gòu)圖像，從而得到效果不同的圖像。比如需要對(duì)圖像進(jìn)行圖像退化從而降低噪聲的影響，則只需要選取低頻的IMF分量進(jìn)行重構(gòu)新的圖像；需要突出圖像細(xì)節(jié)從而增強(qiáng)圖像，則只需要選取高頻的IMF分量進(jìn)行重構(gòu)圖像。本文所提出的方法通過(guò)選取EMD分解出來(lái)的高頻IMF分量重構(gòu)出細(xì)節(jié)更加突出的圖像，隨后用Canny算子對(duì)重構(gòu)圖像進(jìn)行邊緣檢測(cè)，從而引導(dǎo)P-M方程實(shí)現(xiàn)更好的擴(kuò)散去噪性能。該散斑噪聲抑制方法的主要處理過(guò)程如下：

（1）對(duì)圖像H'(x,y)進(jìn)行EMD分解，式（1）可以改寫(xiě)為：

H'(x,y)=(∑i=1nCi+r)N(x,y)

 . （2）

由于全息圖的細(xì)節(jié)信息主要包含在高頻分量C1中，因此只使用高頻分量C1來(lái)重構(gòu)新的全息圖以突出細(xì)節(jié)信息，式（2）可以改寫(xiě)為：

H''(x,y)=C1N(x,y)

 . （3）

（2）對(duì)重構(gòu)的圖像H''(x,y)進(jìn)行Canny邊緣檢測(cè)，得到邊緣檢測(cè)的結(jié)果C(x,y)，并記錄上閾值k1。

（3）根據(jù)邊緣檢測(cè)結(jié)果引導(dǎo)各向異性擴(kuò)散方程的擴(kuò)散程度，進(jìn)行去噪，得到散斑抑制后的全息圖H'''(x,y)。

（4）隨后對(duì)散斑抑制后的全息圖H'''(x,y)進(jìn)行數(shù)值重建，從再現(xiàn)像中提取出相位信息進(jìn)行定量相位成像。由于此時(shí)的相位信息包裹在（-π，π）之中，需要使用相應(yīng)的算法進(jìn)行相位展開(kāi)以獲得真實(shí)的相位信息。最后根據(jù)相位誤差的來(lái)源，采用相對(duì)應(yīng)的相位補(bǔ)償算法對(duì)誤差進(jìn)行補(bǔ)償。

2.2　EMD

二維經(jīng)驗(yàn)?zāi)B(tài)分解的目的是將一個(gè)二維信號(hào)f(x,y)分解成一系列頻率從高到低的IMF分量，高頻部分包含物體細(xì)節(jié)信息，低頻部分包含物體輪廓信息，具體過(guò)程如下：

（1）找出f(x,y)的極大值和極小值，通過(guò)插值法構(gòu)造出極大值曲面u(x,y)和極小值曲面v(x,y)，極大值曲面和極小值曲面的均值為：

e1(x,y)=u(x,y)+v(x,y)2

 . （4）

均值e1(x,y)與f(x,y)的差記為h1(x,y)，有：

h1(x,y)=f(x,y)?e1(x,y)

 . （5）

（2）將h1(x,y)作為新的輸入重復(fù)上述過(guò)程k次，直到滿足式（6）時(shí)方可停止循環(huán)。其中，SD的值一般在0.1~0.5之間，α一般為0.2。

SD=∑x=1X∑y=1Y∣∣h1(k?1)(x,y)?h1k(x,y)∣∣2h1k2(x,y)≤α

 . （6）

值得注意的是，SD的表達(dá)式實(shí)際上是兩次相鄰迭代計(jì)算結(jié)果之間的標(biāo)準(zhǔn)差計(jì)算公式，主要反映了兩次迭代計(jì)算之間的偏差程度。當(dāng)SD≤α?xí)r，說(shuō)明兩次相鄰迭代計(jì)算結(jié)果h1(k?1)(x,y)與h1k(x,y)之間相差很小，達(dá)到了可接受的計(jì)算精度，如果再繼續(xù)循環(huán)迭代，對(duì)于計(jì)算精度的提高作用不明顯。此時(shí)計(jì)算得到的h1k(x,y)是一個(gè)IMF分量（C1），包含著圖像f(x,y)的最高頻信息。

（3）令r1(x,y)表示圖像f(x,y)去除最高頻信息C1之后剩余的部分，則有：

r1(x,y)=f(x,y)?C1

 . （7）

（4）令r1(x,y)作為新的待分析圖像，并重復(fù)步驟（1）和步驟（2）的計(jì)算過(guò)程之后，得到第2個(gè)IMF分量C2。

（5）重復(fù)上述過(guò)程n次，當(dāng)rn(x,y)和Cn小于預(yù)定誤差或者rn(x,y)為單調(diào)函數(shù)時(shí)，則不能再?gòu)膱D像f(x,y)中提取出IMF分量。此時(shí)，圖像f(x,y)如式（8）表示：

f(x,y)=∑i=1nCi+r

 . （8）

2.3　Canny算子改進(jìn)P-M方程法

基于偏微分方程的圖像去噪方法得到了廣泛的應(yīng)用，從線性均勻擴(kuò)散最終發(fā)展到各向異性擴(kuò)散。在各向異性擴(kuò)散中，應(yīng)用最為廣泛的是Perona和Malik共同提出的Perona-Malik（P-M）方程［22］，許多方法都是基于P-M方程發(fā)展起來(lái)的［23-24］。

然而傳統(tǒng)的P-M方程法存在以下問(wèn)題［25］：

（1）由于梯度算子?無(wú)法去除大的孤立噪聲點(diǎn)，在大的噪聲點(diǎn)處的梯度值會(huì)很大，同時(shí)擴(kuò)散系數(shù)c與梯度值成反比，這會(huì)導(dǎo)致此處的擴(kuò)散系數(shù)c的值變小，達(dá)不到擴(kuò)散去噪的效果。

（2）梯度算子?對(duì)噪聲的敏感度高，抗噪性能不強(qiáng)，不能識(shí)別偽邊緣。

為了克服上述問(wèn)題，通過(guò)引入Canny邊緣檢測(cè)算子改進(jìn)P-M方程從而克服梯度算子?抗噪能力不強(qiáng)的問(wèn)題［19］。其數(shù)學(xué)模型如式（9）所示：

?I(x,y,t)?t=div[(1?K) ?c(|?I(x,y,t)|)??I(x,y,t)]

 ， （9）

其中：t為時(shí)間，I(x,y,0)為初始狀態(tài)，?為梯度算子，div為散度算子，K為Canny算子。值得注意的是，當(dāng)該像素點(diǎn)是圖像邊緣時(shí)，K趨近于1，這會(huì)使得P-M方程的擴(kuò)散速度減慢，有利于保留圖像細(xì)節(jié)信息；否則，當(dāng)K趨近于0時(shí)，會(huì)使得P-M方程的擴(kuò)散速度加快，從而更好地去除噪聲。c為擴(kuò)散系數(shù)，用于控制P-M方程的擴(kuò)散程度，常用的計(jì)算公式有兩種，分別如式（10）和式（11）所示：

c(x)=exp(?xk2)

 ， （10）

c(x)=11+(xk)2

 ， （11）

在式（10）和式（11）中，k為控制擴(kuò)散過(guò)程的一個(gè)系數(shù)。

使用Canny算子K改進(jìn)P-M方程法的優(yōu)勢(shì)在于：一是Canny算子K相比于梯度算子?對(duì)噪聲的敏感度更低，抗噪能力更好。二是Canny邊緣檢測(cè)中的雙閾值篩選一般涵蓋圖像70%的非邊緣像素點(diǎn)，令擴(kuò)散系數(shù)c中的參數(shù)k等于上閾值k1，顯然這是一個(gè)合理的選擇。該方法具體過(guò)程如下所示：

（1）對(duì)圖像g(x,y)進(jìn)行Canny邊緣檢測(cè)，得到邊緣檢測(cè)的結(jié)果C(x,y)，以及記錄上閾值k1。

（2）按上、下、左、右4個(gè)方向求解圖像g(x,y) 4個(gè)方向上的梯度，即?N、?S、?E和?W。

（3）令擴(kuò)散系數(shù)c中的參數(shù)k等于上閾值k1，同時(shí)結(jié)合4個(gè)方向的梯度求解得到cN、cS、cE和cW。

（4）由于Canny邊緣檢測(cè)結(jié)果C(x,y)中只有0和1兩種值，當(dāng)C(x,y)=1時(shí)，令Canny算子K=0.01；當(dāng)C(x,y)=0時(shí)，令Canny算子K=0.99。

（5）散度算子div可以按照式（12）計(jì)算：

div=14(cN×

?N

+cS×

?S

+cE×

?E

+cW×

?W)

 . （12）

（6）圖像g(x,y)經(jīng)過(guò)Canny算子改進(jìn)P-M方程法處理后，可以表示為：

g'(x,y)=g(x,y)+(1?K)div

 . （13）

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

3.1　實(shí)驗(yàn)裝置及樣本

本文搭建了一套基于Mach-Zehnder干涉方法的數(shù)字記錄光路，如圖1所示。將一個(gè)波長(zhǎng)為632.8 nm的He-Ne激光器作為光源，其發(fā)出的光首先經(jīng)過(guò)由一塊焦距為15 mm的平凹透鏡和一塊焦距為250 mm的平凸透鏡組成的擴(kuò)束準(zhǔn)直器，隨后被分束鏡BS1分成兩束等光強(qiáng)的光。一束經(jīng)過(guò)樣品之后被反射鏡M1反射到分束鏡BS2，稱之為物光；另外一束光則直接被反射鏡M2反射到BS2上，稱之為參考光。最后，物光和參考光在BS2處合成一束光，在CCD（BFLY-U3-23S6M-C，相機(jī)分辨率為1 920×1 200，像素大小為5.86 μm）的記錄面上干涉形成全息圖，并由CCD記錄下來(lái)。

圖1  基于Mach-Zehnder干涉儀的數(shù)字記錄光路原理圖及實(shí)物圖。（a）原理圖；（b）實(shí)物圖。

Fig.1  Schematic diagram and physical diagram of digital recording optical path based on Mach-Zehnder interferometer. （a） Schematic diagram； （b） Physical diagram.

下載: 原圖 | 高精圖 | 低精圖

在如圖2（a）所示的相位型樣品板中，其以玻璃材料為基底，表面圖案具有一定的厚度，當(dāng)光透射樣品板時(shí)，圖案的厚度會(huì)引起干涉條紋發(fā)生形變。我們選取圖2（b）~（f）所示的5種不同的樣本置于光路之后，得到5張對(duì)應(yīng)的全息圖，分別如圖2（g）~（k）所示。

圖2  樣品及全息圖。（a）樣品板；（b）樣品A；（c）樣品B；（d）樣品C；（e）樣品D；（f）樣品E；（g）樣品A對(duì)應(yīng)的全息圖A；（h）樣品B對(duì)應(yīng)的全息圖B；（i）樣品C對(duì)應(yīng)的全息圖C；（j）樣品D對(duì)應(yīng)的全息圖D；（k）樣品E對(duì)應(yīng)的全息圖E。

Fig.2  Samples and holograms. （a） Sample plate； （b） Sample A； （c） Sample B； （d） Sample C； （e） Sample D； （f） Sample E； （g） Hologram A corresponding to sample A； （h） Hologram B corresponding to sample B； （i） Hologram C corresponding to sample C； （j） Hologram D corresponding to sample D； （k） Hologram E corresponding to sample E.

下載: 原圖 | 高精圖 | 低精圖

3.2　結(jié)合EMD的Canny算子改進(jìn)P-M方程法的相位重建

首先，將全息圖進(jìn)行EMD，得到一系列從高頻到低頻的IMF分量，通過(guò)選取不同的IMF分量可以重構(gòu)出不同效果的圖像，圖3為任意選取全息圖B分解出的不同IMF分量進(jìn)行重構(gòu)的圖像。根據(jù)結(jié)果可知，C1重構(gòu)出來(lái)的圖像可以觀察到明顯的條紋，且條紋存在明顯的形變；C3重構(gòu)出來(lái)的圖像為圖案的輪廓信息，能夠大致識(shí)別圖案；C5只能重構(gòu)出非常粗略的部分輪廓信息；C7重構(gòu)出的圖像非常模糊，已無(wú)法分辨出圖像的形狀。因此，我們選取全息圖的最高頻IMF分量，即C1進(jìn)行圖像重構(gòu)。

圖3  選取不同IMF分量的重構(gòu)結(jié)果。（a）選取C1的重構(gòu)圖像；（b）選取C3的重構(gòu)圖像；（c）選取C5的重構(gòu)圖像；（d）選取C7的重構(gòu)圖像。

Fig.3  Reconstruction results with different IMF components selected. （a） Reconstructed image with C1； （b） Reconstructed image with C3； （c） Reconstructed image with C5； （d） Reconstructed image with C7.

下載: 原圖 | 高精圖 | 低精圖

其次，用Canny算子對(duì)全息圖進(jìn)行邊緣檢測(cè)。圖4所示為對(duì)圖3（a）和圖2（h）進(jìn)行邊緣檢測(cè)的結(jié)果。通過(guò)對(duì)比可知，前者的邊緣檢測(cè)結(jié)果更好。

圖4  Canny算子邊緣檢測(cè)結(jié)果。（a） Canny算子對(duì)圖3（a）的邊緣檢測(cè)結(jié)果；（b） Canny算子對(duì)圖2（h）的邊緣檢測(cè)結(jié)果。

Fig.4  Canny operator edge detection results. （a） Edge detection results of Canny operator for Fig.3（a）； （b） Edge detection results of Canny operator for Fig.2（h）.

下載: 原圖 | 高精圖 | 低精圖

我們用Canny算子改進(jìn)P-M方程法對(duì)5張由最高頻IMF分量重構(gòu)的圖像進(jìn)行去噪，結(jié)果如圖5所示。

圖5  經(jīng)過(guò)所提出方法處理后的全息圖。（a）全息圖A；（b）全息圖B；（c）全息圖C；（d）全息圖D；（e）全息圖E。

Fig.5  Holograms processed by the proposed method. （a） Hologram A； （b） Hologram B； （c） Hologram C； （d） Hologram D； （e） Hologram E.

下載: 原圖 | 高精圖 | 低精圖

最后，對(duì)處理后的全息圖采用角譜衍射法［26］進(jìn)行數(shù)值重建，從再現(xiàn)像中恢復(fù)出包裹相位信息。采用基于離散余弦的最小二乘法［27］對(duì)包裹相位進(jìn)行展開(kāi)，以此獲得真實(shí)相位。針對(duì)所搭建的基于Mach-Zehnder干涉儀的數(shù)字記錄光路，其相位誤差的主要來(lái)源是物光與參考光之間引入的傾斜誤差。采用最小二乘擬合法［28］進(jìn)行相位補(bǔ)償，最終得到的相位圖如圖6（a）~（e）所示。為了驗(yàn)證本文所提出方法處理后恢復(fù)相位的準(zhǔn)確性，對(duì)未進(jìn)行任何去噪處理的原始全息圖進(jìn)行相位成像，結(jié)果如圖6（f）~（j）所示。樣品A經(jīng)本文所提方法處理后的重建相位圖像和未經(jīng)任何處理的重建相位圖像分別為圖6（a）和圖6（f），前者相位分布更加規(guī)則，后者存在明顯的相位形變；樣品B經(jīng)本文所提方法處理后的重建相位圖像和未經(jīng)任何處理的重建相位圖像分別為圖6（b）和圖6（g），前者凸起部分在4個(gè)角上分布均勻，后者凸起部分雜亂分布；樣品C經(jīng)本文所提方法處理后的重建相位圖像和未經(jīng)任何處理的重建相位圖像分別為圖6（c）和圖6（h），前者呈規(guī)則的球型分布，后者圓球變形且在三維顯示中其頂部有不規(guī)則的凸起，由于高度變化較小，所以二維顯示對(duì)于不規(guī)則凸起不是很明顯，但能看出其存在明顯形變；樣品D經(jīng)本文所提方法處理后的重建相位圖像和未經(jīng)任何處理的重建相位圖像分別為圖6（d）和圖6（i），兩者較為接近，但可以看出前者中心部分的相位細(xì)節(jié)更加突出；樣品E經(jīng)本文所提方法處理后的重建相位圖像和未經(jīng)任何處理的重建相位圖像分別為圖6（e）和圖6（j），前者凸起的條形分布更加規(guī)則，后者條形呈鋸齒形狀。因此，根據(jù)成像結(jié)果可知，經(jīng)本文所提方法處理后恢復(fù)的相位分布輪廓更加光滑，其在定量相位成像中具有一定的優(yōu)勢(shì)。

圖6  經(jīng)過(guò)所提出方法處理后的重建相位圖和未去噪全息圖的重建相位圖對(duì)比。（a）全息圖A經(jīng)過(guò)所提出方法處理后重建的相位圖；（b）全息圖B經(jīng)過(guò)所提出方法處理后重建的相位圖；（c）全息圖C經(jīng)過(guò)所提出方法處理后重建的相位圖；（d）全息圖D經(jīng)過(guò)所提出方法處理后重建的相位圖；（e）全息圖E經(jīng)過(guò)所提出方法處理后重建的相位圖；（f）未去噪全息圖A的重建相位圖；（g）未去噪全息圖B的重建相位圖；（h）未去噪全息圖C的重建相位圖；（i）未去噪全息圖D的重建相位圖；（j）未去噪全息圖E的重建相位圖。

Fig.6  Comparison of the reconstructed phase maps after processing by the proposed method and the reconstructed phase maps without denoising holograms. （a） Reconstructed phase map of hologram A after processing by the proposed method； （b） Reconstructed phase map of hologram B after processing by the proposed method； （c） Reconstructed phase map of the hologram C after processing by the proposed method； （d） Reconstructed phase map of hologram D after processing by the proposed method； （e） Reconstructed phase map of the hologram E after processing by the proposed method； （f） Reconstructed phase map of hologram A without denoising； （g） Reconstructed phase map of hologram B without denoising； （h） Reconstructed phase map of hologram C without denoising； （i） Reconstructed phase map of hologram D without denoising； （j） Reconstructed phase map of hologram E without denoising.

下載: 原圖 | 高精圖 | 低精圖

3.3　實(shí)驗(yàn)對(duì)比

為了驗(yàn)證所提出方法的有效性，將所提出的方法與基于Canny算子改進(jìn)P-M方程法、均值濾波法以及中值濾波法進(jìn)行比較。所記錄的5種全息圖經(jīng)過(guò)基于Canny算子改進(jìn)P-M方程法、均值濾波法和中值濾波法處理后，經(jīng)過(guò)角譜衍射法數(shù)值重建［26］、基于離散余弦變換的最小二乘法相位展開(kāi)［27］以及最小二乘擬合法［28］相位補(bǔ)償之后得到的相位圖如圖7所示。相對(duì)于未去噪聲的相位圖，其相位分布均勻且輪廓光滑，在成像質(zhì)量方面得到了較大的提升。通過(guò)將其與圖6（a）~（e）相比可知，基于Canny算子改進(jìn)P-M方程法處理后重建的相位圖像與本文所提出方法處理后重建的相位圖像更為相似，但后者相位輪廓分布更加規(guī)則且細(xì)節(jié)形狀突出。均值濾波和中值濾波在一些樣品中存在丟失細(xì)節(jié)等問(wèn)題，不利于高精度的恢復(fù)物體的相位。例如，在圖7（k）中，橢圓形凸起部分存在明顯的左低右高；圖7（g）和圖7（l）中凸起小圖案的左下角存在明顯缺失；圖7（i）雖然整體上光滑，但其中間部分已經(jīng)基本看不出具體輪廓；圖7（j）和圖7（o）的凸起條形存在明顯的鋸齒形狀。

圖7  全息圖A、B、C、D、E經(jīng)過(guò)基于Canny算子改進(jìn)P-M方程法、均值濾波法與中值濾波法處理后重建的相位圖。（a）~（e）經(jīng)過(guò)基于Canny算子改進(jìn)P-M方程法處理后重建的相位圖；（f）~（j）經(jīng)過(guò)均值濾波法處理后重建的相位圖；（k）~（o）經(jīng)過(guò)中值濾波法處理后重建的相位圖。

Fig.7  Reconstructed phase maps of hologram A， B， C， D and E after processing by the improved P-M equation method based on the Canny operator， mean filter method and median filter method. （a）~（e） Reconstructed phase maps after processing by the improved P-M equation method based on Canny operator； （f）~（j） Reconstructed phase maps after processing by the mean filtering method； （k）~（o） Reconstructed phase maps after processing by the median filter method.

下載: 原圖 | 高精圖 | 低精圖

4 分析與討論

4.1　降噪定量分析

為了對(duì)降噪后的重建相位圖做定量分析，本文采用結(jié)構(gòu)相似性（SSIM）［29］、邊緣保持指數(shù)（EPI）［30］和散斑抑制指數(shù)（SSI）［31］作為定量分析的指標(biāo)。

SSIM的計(jì)算公式如式（14）所示：

SSIM=l(x,y)?c(x,y)?s(x,y)

 ， （14）

其中：l(x,y)表示圖像的亮度，c(x,y)表示圖像的對(duì)比度，s(x,y)表示圖像的結(jié)構(gòu)。SSIM的值越大，說(shuō)明圖像失真的程度越小。當(dāng)兩張圖像相同時(shí)，SSIM的值等于1。其中，參考圖像是通過(guò)多幅全息圖疊加平均法獲得。

EPI的計(jì)算公式如式（15）所示：

EPI=∑i=1m∑j=1n|GR11?GR21|after∑i=1m∑j=1n|GR12?GR22|before

 ， （15）

其中：GR11和GR21表示圖像抑制后兩個(gè)相鄰的像素點(diǎn)，GR12和GR22表示圖像抑制前兩個(gè)相鄰的像素點(diǎn)。邊緣保持指數(shù)的最大值為1，最小值為0。邊緣保持指數(shù)越高，則說(shuō)明經(jīng)過(guò)處理之后對(duì)原圖像的細(xì)節(jié)信息保持得越好。

SSI的計(jì)算公式如式（16）所示：

SSI=var(Iafter)????????√mean(Iafter)?mean(Ibefore)var(Ibefore)????????√

 ， （16）

其中：Iafter表示經(jīng)過(guò)散斑抑制后的圖像；Ibefore表示散斑抑制前的圖像，即含噪的圖像；var(I)?????√表示圖像I(x,y)的標(biāo)準(zhǔn)差；mean(I)表示圖像I(x,y)的均值。通常，SSI小于1，散斑被抑制。SSI越小，則說(shuō)明散斑抑制能力越強(qiáng)。

5種全息圖經(jīng)過(guò)所提出方法、基于Canny算子改進(jìn)P-M方程法、均值濾波法以及中值濾波法處理之后重建相位圖對(duì)應(yīng)的SSIM、EPI和SSI計(jì)算結(jié)果如表1所示。

表1  不同散斑抑制方法的SSIM、EPI、SSI計(jì)算結(jié)果

Tab.1  Calculation results of SSIM， EPI and SSI for different speckle suppression methods

樣品 方法 SSIM EPI SSI

樣品A 本文所提方法 0.968 1 0.863 9 0.743 8

基于Canny算子改進(jìn)P-M方程法 0.920 9 0.836 5 0.786 1

均值濾波法 0.898 8 0.820 1 0.808 7

中值濾波法 0.912 8 0.832 3 0.800 1

樣品B 本文所提方法 0.901 3 0.947 9 0.753 7

基于Canny算子改進(jìn)P-M方程法 0.878 4 0.925 1 0.775 8

均值濾波法 0.813 3 0.806 6 0.857 7

中值濾波法 0.850 8 0.839 3 0.781 0

樣品C 本文所提方法 0.963 7 0.950 6 0.635 4

基于Canny算子改進(jìn)P-M方程法 0.951 1 0.942 3 0.658 3

均值濾波法 0.914 1 0.919 4 0.692 0

中值濾波法 0.948 1 0.936 9 0.691 5

樣品D 本文所提方法 0.912 6 0.846 8 0.922 2

基于Canny算子改進(jìn)P-M方程法 0.893 3 0.833 0 0.930 8

均值濾波法 0.775 5 0.699 4 0.986 9

中值濾波法 0.803 1 0.774 1 0.963 7

樣品E 本文所提方法 0.943 6 0.873 1 0.858 4

基于Canny算子改進(jìn)P-M方程法 0.886 5 0.868 6 0.887 7

均值濾波法 0.768 0 0.700 7 0.910 9

中值濾波法 0.877 8 0.806 0 0.891 2

下載: 導(dǎo)出CSV

從表1的定量分析計(jì)算結(jié)果可以得出，本文所提出方法對(duì)應(yīng)的SSIM、EPI的計(jì)算結(jié)果最大，SSI的計(jì)算結(jié)果最小。在5種樣品之中，以均值濾波對(duì)應(yīng)的SSIM為計(jì)算基準(zhǔn)，中值濾波法相對(duì)于均值濾波法的SSIM最大提高了14.296 9%，最小提高了1.557 6%，平均提高了5.548 8%；基于Canny算子改進(jìn)P-M方程法相對(duì)于均值濾波法的SSIM最大提高了15.429 7%，最小提高了2.458 8%，平均提高了9.026 2%；本文所提出方法相對(duì)于均值濾波法的SSIM最大提高了22.864 6%，最小提高了5.426 1%，平均提高了12.900 0%。

以均值濾波對(duì)應(yīng)的EPI為計(jì)算基準(zhǔn)，中值濾波法相對(duì)于均值濾波法的EPI最大提高了15.027 8%，最小提高了1.487 6%，平均提高了6.630 7%；基于Canny算子改進(jìn)P-M方程法相對(duì)于均值濾波法的EPI最大提高了23.961 8%，最小提高了1.999 7%，平均提高了12.449 1%；本文所提出方法相對(duì)于均值濾波法的EPI最大提高了24.606 4%，最小提高了3.392 5%，平均提高了14.386 1%。

以均值濾波的SSI為計(jì)算基準(zhǔn)，中值濾波法相對(duì)于均值濾波法的SSI最大降低了8.942 5%，最小降低了0.072 3%，平均降低了2.918 3%；基于Canny算子改進(jìn)P-M方程法相對(duì)于均值濾波法的SSI最大降低了9.548 8%，最小降低了2.546 9%，平均降低了5.089 3%；本文所提出方法相對(duì)于均值濾波法的SSI最大降低了12.125 5%，最小降低了5.763 5%，平均降低了8.129 9%。

4.2　相位截面曲線分析

為了驗(yàn)證降噪后定量相位重建的精確性，在如圖6（a）、圖7（a）、圖7（f）和圖7（k）所示相位圖的y=711處（圖6（a）中紅色虛線處）做截面曲線對(duì)比分析，同時(shí)與原始相位圖比較，結(jié)果如圖8所示；在如圖6（b）、圖7（b）、圖7（g）和圖7（l）所示相位圖的x=415處（圖6（b）中紅色虛線處）做截面曲線對(duì)比分析，同時(shí)與原始相位圖比較，結(jié)果如圖9所示；在如圖6（c）、圖7（c）、圖7（h）和圖7（m）所示相位圖的x=589處（圖6（c）中紅色虛線處）做截面曲線對(duì)比分析，同時(shí)與原始相位圖比較，結(jié)果如圖10所示；在如圖6（d）、圖7（d）、圖7（i）和圖7（n）所示相位圖的y=959處（圖6（d）中紅色虛線處）做截面曲線對(duì)比分析，同時(shí)與原始相位圖比較，結(jié)果如圖11所示；在如圖6（e）、圖7（e）、圖7（j）和圖7（o）所示相位圖的x=1 355處（圖6（e）中紅色虛線處）做截面曲線對(duì)比分析，同時(shí)與原始相位圖比較，結(jié)果如圖12所示。其中，原始相位圖的獲取采用了多幅全息圖疊加平均法，通過(guò)記錄6幅同一場(chǎng)景下的全息圖，對(duì)6張?jiān)佻F(xiàn)像疊加取平均后就獲得相對(duì)真實(shí)的相位分布。該相位分布輪廓較為光滑，用其來(lái)評(píng)價(jià)不同散斑抑制方法性能的高低。

圖8  圖6（a）、圖7（a）、圖7（f）和圖7（k）所示相位圖與原始相位圖在y=711處的截面曲線對(duì)比。

Fig.8  Phase maps shown in Fig.6（a）， Fig.7（a）， Fig.7（f） and Fig.7（k） compared with the cross-sectional curves of the original phase map at y=711.

下載: 原圖 | 高精圖 | 低精圖

圖9  圖6（b）、圖7（b）、圖7（g）和圖7（l）所示相位圖與原始相位圖在x=415處的截面曲線對(duì)比。

Fig.9  Phase maps shown in Fig.6（b）， Fig.7（b）， Fig.7（g） and Fig.7（l） are compared with the cross-sectional curves of the original phase map at x=415.

下載: 原圖 | 高精圖 | 低精圖

圖10  圖6（c）、圖7（c）、圖7（h）和圖7（m）所示相位圖與原始相位圖在x=589處的截面曲線對(duì)比。

Fig.10  Phase maps shown in Fig.6（c），F(xiàn)ig. 7（c）， Fig.7（h） and Fig.7（m） compared with the cross-sectional curves of the original phase map at x=589.

下載: 原圖 | 高精圖 | 低精圖

圖11  圖6（d）、圖7（d）、圖7（i）和圖7（n）所示相位圖與原始相位圖在y=959處的截面曲線對(duì)比。

Fig.11  Phase maps shown in Fig.6（d）， Fig.7（d）， Fig.7（i） and Fig.7（n） compared with the cross-sectional curves of the original phase map at y=959.

下載: 原圖 | 高精圖 | 低精圖

圖12  圖6（e）、圖7（e）、圖7（j）和圖7（o）所示相位圖與原始相位圖在x=1 355處的截面曲線對(duì)比。

Fig.12  Phase maps shown in Fig.6（e）， Fig.7（e）， Fig.7（j） and Fig.7（o） compared with the cross-sectional curves of the original phase map at x=1 355.

下載: 原圖 | 高精圖 | 低精圖

從相位截面曲線變化來(lái)看，綠色虛線、藍(lán)色點(diǎn)、橙色點(diǎn)線、紫色雙點(diǎn)線與紅色實(shí)線的變化趨勢(shì)基本一致，但是綠色虛線與紅色實(shí)線的偏差要比藍(lán)色點(diǎn)、橙色點(diǎn)線、紫色雙點(diǎn)線與紅色實(shí)線的偏差要小，具體的偏差如表2所示。其中，綠色虛線、藍(lán)色點(diǎn)、橙色點(diǎn)線、紫色雙點(diǎn)線與紅色實(shí)線的最大偏差分別為偏差1、偏差2、偏差3、偏差4；綠色虛線、藍(lán)色點(diǎn)、橙色點(diǎn)線、紫色雙點(diǎn)線與紅色實(shí)線的最小偏差分別為偏差5、偏差6、偏差7、偏差8；綠色虛線、藍(lán)色點(diǎn)、橙色點(diǎn)線、紫色雙點(diǎn)線與紅色實(shí)線的平均偏差分別為偏差9、偏差10、偏差11、偏差12。

表2  綠色虛線、藍(lán)色點(diǎn)、橙色點(diǎn)線和紫色雙點(diǎn)線與紅色實(shí)線的最大、最小、平均偏差

Tab.2  Maximum， minimum and average deviation of the green dashed line， the blue point， the orange point line and the purple double point line from the red solid line

偏差 樣品A 樣品B 樣品C 樣品D 樣品E

偏差1 1.867 1 0.440 3 0.084 6 0.735 5 0.387 5

偏差2 2.770 2 0.465 3 0.228 8 0.825 1 0.420 9

偏差3 4.620 1 0.708 4 0.148 6 0.927 5 1.208 3

偏差4 4.321 8 0.467 9 0.144 5 0.840 5 1.106 5

偏差5 0.000 1 0.0001 0.000 1 0.000 1 0.000 1

偏差6 0.0001 0.000 1 0.000 1 0.000 2 0.000 1

偏差7 0.000 2 0.000 2 0.000 1 0.001 3 0.000 2

偏差8 0.000 2 0.000 1 0.000 1 0.000 2 0.000 1

偏差9 0.278 5 0.053 2 0.0240 0.113 9 0.067 6

偏差10 0.342 0 0.061 9 0.029 6 0.147 3 0.086 0

偏差11 0.802 9 0.071 6 0.045 3 0.266 5 0.140 1

偏差12 0.616 8 0.068 1 0.034 1 0.177 9 0.095 9

下載: 導(dǎo)出CSV

綜合上述結(jié)果可以得出，在定量相位成像中，本文所提出的方法——結(jié)合EMD的Canny算子改進(jìn)P-M方程法相對(duì)于基于Canny算子改進(jìn)P-M方程法、均值濾波法和中值濾波法的SSIM、EPI與SSI計(jì)算結(jié)果以及截面曲線對(duì)比結(jié)果均更好。本文所提出方法具備更好的全息散斑噪聲抑制效果，誤差更小，相位重建精度更高。

5 結(jié)論

本文提出一種結(jié)合EMD的Canny算子改進(jìn)P-M方程的散斑噪聲抑制方法，該方法用來(lái)抑制在定量相位成像中高相干性光源產(chǎn)生的散斑噪聲對(duì)相位重建精度的影響。通過(guò)EMD獲得高頻分量的IMF重構(gòu)出細(xì)節(jié)突出的圖像，使得Canny算子檢測(cè)更加準(zhǔn)確，從而更好地控制P-M方程的擴(kuò)散去噪過(guò)程，在去噪的同時(shí)又保留高精度的相位信息。文中對(duì)本文所提出方法、基于Canny算子改進(jìn)P-M方程法、均值濾波法和中值濾波法降噪后的重建相位圖用SSIM、EPI和SSI進(jìn)行了定量分析。其中，以均值濾波為參考，本文所提方法去噪后的重建相位圖像的SSIM提高了12.900 0%，EPI提高了14.386 1%，SSI降低了8.129 9%。通過(guò)對(duì)不同方法降噪處理后重建的相位圖像進(jìn)行截面曲線分析，經(jīng)對(duì)比發(fā)現(xiàn)，本文所提方法與原始相位偏差更小。實(shí)驗(yàn)結(jié)果表明，本文所提方法在性能上優(yōu)于其他3種降噪方法，有較好的散斑抑制效果，相位重建精度更高。