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

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

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

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

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分量。通過這些不同頻率的IMF進(jìn)行選取并重構(gòu)圖像，從而得到效果不同的圖像。比如需要對圖像進(jìn)行圖像退化從而降低噪聲的影響，則只需要選取低頻的IMF分量進(jìn)行重構(gòu)新的圖像；需要突出圖像細(xì)節(jié)從而增強(qiáng)圖像，則只需要選取高頻的IMF分量進(jìn)行重構(gòu)圖像。本文所提出的方法通過選取EMD分解出來的高頻IMF分量重構(gòu)出細(xì)節(jié)更加突出的圖像，隨后用Canny算子對重構(gòu)圖像進(jìn)行邊緣檢測，從而引導(dǎo)P-M方程實現(xiàn)更好的擴(kuò)散去噪性能。該散斑噪聲抑制方法的主要處理過程如下：

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

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

.（2）

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

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

.（3）

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

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

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

2.2EMD

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

（1）找出f(x,y)的極大值和極小值，通過插值法構(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ù)上述過程k次，直到滿足式（6）時方可停止循環(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á)式實際上是兩次相鄰迭代計算結(jié)果之間的標(biāo)準(zhǔn)差計算公式，主要反映了兩次迭代計算之間的偏差程度。當(dāng)SD≤α?xí)r，說明兩次相鄰迭代計算結(jié)果h1(k?1)(x,y)與h1k(x,y)之間相差很小，達(dá)到了可接受的計算精度，如果再繼續(xù)循環(huán)迭代，對于計算精度的提高作用不明顯。此時計算得到的h1k(x,y)是一個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）的計算過程之后，得到第2個IMF分量C2。

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

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

.（8）

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

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

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

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

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

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

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

，（9）

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

c(x)=exp(?xk2)

，（10）

c(x)=11+(xk)2

，（11）

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

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

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

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

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

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

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

div=14(cN×

?N

+cS×

?S

+cE×

?E

+cW×

?W)

.（12）

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

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

.（13）

3實驗

3.1實驗裝置及樣本

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

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

Fig.1SchematicdiagramandphysicaldiagramofdigitalrecordingopticalpathbasedonMach-Zehnderinterferometer.（a）Schematicdiagram；（b）Physicaldiagram.

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

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

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

Fig.2Samplesandholograms.（a）Sampleplate；（b）SampleA；（c）SampleB；（d）SampleC；（e）SampleD；（f）SampleE；（g）HologramAcorrespondingtosampleA；（h）HologramBcorrespondingtosampleB；（i）HologramCcorrespondingtosampleC；（j）HologramDcorrespondingtosampleD；（k）HologramEcorrespondingtosampleE.

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

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

首先，將全息圖進(jìn)行EMD，得到一系列從高頻到低頻的IMF分量，通過選取不同的IMF分量可以重構(gòu)出不同效果的圖像，圖3為任意選取全息圖B分解出的不同IMF分量進(jìn)行重構(gòu)的圖像。根據(jù)結(jié)果可知，C1重構(gòu)出來的圖像可以觀察到明顯的條紋，且條紋存在明顯的形變；C3重構(gòu)出來的圖像為圖案的輪廓信息，能夠大致識別圖案；C5只能重構(gòu)出非常粗略的部分輪廓信息；C7重構(gòu)出的圖像非常模糊，已無法分辨出圖像的形狀。因此，我們選取全息圖的最高頻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.3ReconstructionresultswithdifferentIMFcomponentsselected.（a）ReconstructedimagewithC1；（b）ReconstructedimagewithC3；（c）ReconstructedimagewithC5；（d）ReconstructedimagewithC7.

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

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

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

Fig.4Cannyoperatoredgedetectionresults.（a）EdgedetectionresultsofCannyoperatorforFig.3（a）；（b）EdgedetectionresultsofCannyoperatorforFig.2（h）.

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

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

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

Fig.5Hologramsprocessedbytheproposedmethod.（a）HologramA；（b）HologramB；（c）HologramC；（d）HologramD；（e）HologramE.

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

最后，對處理后的全息圖采用角譜衍射法［26］進(jìn)行數(shù)值重建，從再現(xiàn)像中恢復(fù)出包裹相位信息。采用基于離散余弦的最小二乘法［27］對包裹相位進(jìn)行展開，以此獲得真實相位。針對所搭建的基于Mach-Zehnder干涉儀的數(shù)字記錄光路，其相位誤差的主要來源是物光與參考光之間引入的傾斜誤差。采用最小二乘擬合法［28］進(jìn)行相位補(bǔ)償，最終得到的相位圖如圖6（a）~（e）所示。為了驗證本文所提出方法處理后恢復(fù)相位的準(zhǔn)確性，對未進(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個角上分布均勻，后者凸起部分雜亂分布；樣品C經(jīng)本文所提方法處理后的重建相位圖像和未經(jīng)任何處理的重建相位圖像分別為圖6（c）和圖6（h），前者呈規(guī)則的球型分布，后者圓球變形且在三維顯示中其頂部有不規(guī)則的凸起，由于高度變化較小，所以二維顯示對于不規(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)勢。

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

Fig.6Comparisonofthereconstructedphasemapsafterprocessingbytheproposedmethodandthereconstructedphasemapswithoutdenoisingholograms.（a）ReconstructedphasemapofhologramAafterprocessingbytheproposedmethod；（b）ReconstructedphasemapofhologramBafterprocessingbytheproposedmethod；（c）ReconstructedphasemapofthehologramCafterprocessingbytheproposedmethod；（d）ReconstructedphasemapofhologramDafterprocessingbytheproposedmethod；（e）ReconstructedphasemapofthehologramEafterprocessingbytheproposedmethod；（f）ReconstructedphasemapofhologramAwithoutdenoising；（g）ReconstructedphasemapofhologramBwithoutdenoising；（h）ReconstructedphasemapofhologramCwithoutdenoising；（i）ReconstructedphasemapofhologramDwithoutdenoising；（j）ReconstructedphasemapofhologramEwithoutdenoising.

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

3.3實驗對比

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

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

Fig.7ReconstructedphasemapsofhologramA，B，C，DandEafterprocessingbytheimprovedP-MequationmethodbasedontheCannyoperator，meanfiltermethodandmedianfiltermethod.（a）~（e）ReconstructedphasemapsafterprocessingbytheimprovedP-MequationmethodbasedonCannyoperator；（f）~（j）Reconstructedphasemapsafterprocessingbythemeanfilteringmethod；（k）~（o）Reconstructedphasemapsafterprocessingbythemedianfiltermethod.

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

4分析與討論

4.1降噪定量分析

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

SSIM的計算公式如式（14）所示：

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

，（14）

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

EPI的計算公式如式（15）所示：

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

，（15）

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

SSI的計算公式如式（16）所示：

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

，（16）

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

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

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

Tab.1CalculationresultsofSSIM，EPIandSSIfordifferentspecklesuppressionmethods

樣品方法SSIMEPISSI

樣品A本文所提方法0.96810.86390.7438

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

均值濾波法0.89880.82010.8087

中值濾波法0.91280.83230.8001

樣品B本文所提方法0.90130.94790.7537

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

均值濾波法0.81330.80660.8577

中值濾波法0.85080.83930.7810

樣品C本文所提方法0.96370.95060.6354

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

均值濾波法0.91410.91940.6920

中值濾波法0.94810.93690.6915

樣品D本文所提方法0.91260.84680.9222

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

均值濾波法0.77550.69940.9869

中值濾波法0.80310.77410.9637

樣品E本文所提方法0.94360.87310.8584

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

均值濾波法0.76800.70070.9109

中值濾波法0.87780.80600.8912

下載:導(dǎo)出CSV

從表1的定量分析計算結(jié)果可以得出，本文所提出方法對應(yīng)的SSIM、EPI的計算結(jié)果最大，SSI的計算結(jié)果最小。在5種樣品之中，以均值濾波對應(yīng)的SSIM為計算基準(zhǔn)，中值濾波法相對于均值濾波法的SSIM最大提高了14.2969%，最小提高了1.5576%，平均提高了5.5488%；基于Canny算子改進(jìn)P-M方程法相對于均值濾波法的SSIM最大提高了15.4297%，最小提高了2.4588%，平均提高了9.0262%；本文所提出方法相對于均值濾波法的SSIM最大提高了22.8646%，最小提高了5.4261%，平均提高了12.9000%。

以均值濾波對應(yīng)的EPI為計算基準(zhǔn)，中值濾波法相對于均值濾波法的EPI最大提高了15.0278%，最小提高了1.4876%，平均提高了6.6307%；基于Canny算子改進(jìn)P-M方程法相對于均值濾波法的EPI最大提高了23.9618%，最小提高了1.9997%，平均提高了12.4491%；本文所提出方法相對于均值濾波法的EPI最大提高了24.6064%，最小提高了3.3925%，平均提高了14.3861%。

以均值濾波的SSI為計算基準(zhǔn)，中值濾波法相對于均值濾波法的SSI最大降低了8.9425%，最小降低了0.0723%，平均降低了2.9183%；基于Canny算子改進(jìn)P-M方程法相對于均值濾波法的SSI最大降低了9.5488%，最小降低了2.5469%，平均降低了5.0893%；本文所提出方法相對于均值濾波法的SSI最大降低了12.1255%，最小降低了5.7635%，平均降低了8.1299%。

4.2相位截面曲線分析

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

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

Fig.8PhasemapsshowninFig.6（a），F(xiàn)ig.7（a），F(xiàn)ig.7（f）andFig.7（k）comparedwiththecross-sectionalcurvesoftheoriginalphasemapaty=711.

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

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

Fig.9PhasemapsshowninFig.6（b），F(xiàn)ig.7（b），F(xiàn)ig.7（g）andFig.7（l）arecomparedwiththecross-sectionalcurvesoftheoriginalphasemapatx=415.

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

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

Fig.10PhasemapsshowninFig.6（c），F(xiàn)ig.7（c），F(xiàn)ig.7（h）andFig.7（m）comparedwiththecross-sectionalcurvesoftheoriginalphasemapatx=589.

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

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

Fig.11PhasemapsshowninFig.6（d），F(xiàn)ig.7（d），F(xiàn)ig.7（i）andFig.7（n）comparedwiththecross-sectionalcurvesoftheoriginalphasemapaty=959.

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

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

Fig.12PhasemapsshowninFig.6（e），F(xiàn)ig.7（e），F(xiàn)ig.7（j）andFig.7（o）comparedwiththecross-sectionalcurvesoftheoriginalphasemapatx=1355.

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

從相位截面曲線變化來看，綠色虛線、藍(lán)色點、橙色