變分模態(tài)分解與稀疏SURE的電子圖像噪聲抑制

李 慶,Steven Y.Liang,2

(1.東華大學 機械工程學院,上海 201620； 2.佐治亞理工學院 喬治-伍德拉夫機械工程學院，佐治亞州 亞特蘭大 30332-0405)

在微觀機械制造領(lǐng)域如納米切削、微磨削、激光光整加工等，電子背散射衍射(electron backscatter diffraction，EBSD)作為一種新的微觀圖像采集與分析技術(shù)，其集成了顯微組織分析與晶體學圖像分析等模塊，在晶體材料 (如單晶鋁、單晶銅) 的亞微米級尺度對晶粒織構(gòu)與微觀特性分析中 (如晶粒尺寸、取向、滑移位錯、不同相分布、失效機理研究和應變評估) 發(fā)揮了重要作用[1-3]. 然而，電子圖像在采集、傳輸以及儲存的過程中不可避免的受到外界噪聲 (如入射電子束加速電壓、不穩(wěn)定電流、取向成像步長過大等) 的干擾，導致圖像保真度差，影響了材料微觀特性分析.在工程應用中，由于實時在線采集的需要，且為了防止材料晶粒發(fā)生馳豫現(xiàn)象，一般實驗條件不允許圖像重新采集，因此，對已污染微結(jié)構(gòu)圖像進行復原是微觀機械制造領(lǐng)域中一個亟待解決的難題.

一般地，電子微結(jié)構(gòu)圖像復原表現(xiàn)在兩個方面：有效去除干擾噪聲以及保留或增強微結(jié)構(gòu)圖像的邊緣與紋理等固有結(jié)構(gòu)特征.傳統(tǒng)的圖像去噪方法可分為空間域去噪和變換域去噪.1)空間域去噪方法.如偏微分方程方法、維納濾波方法、空間自適應濾波器等[4-6]，其本質(zhì)是考慮待估像素點與其鄰域相鄰像素點的加權(quán)平均來代替待估像素點的真實值，但該方法易造成邊緣信息模糊以及紋理、邊界丟失現(xiàn)象.2)變換域去噪方法.如正交基變換方法(如小波基等)、緊框架變換方法(如Ridgelet變換、Contourlet變換)等[7-8]，正交基變換與多尺度緊框架變換去噪方法雖具有良好的時頻域局部分析能力，但在圖像邊緣與細節(jié)處易產(chǎn)生平滑現(xiàn)象，且這兩種方法均沒有考慮圖像結(jié)構(gòu)紋理特征與冗余信息，無法對含有線、面等奇異特性的圖像進行表征. 因此，電子微結(jié)構(gòu)圖像去噪的難點在于：如何在抑制噪聲的同時避免圖像邊緣被平滑且保留圖像固有的結(jié)構(gòu)與紋理特征.

稀疏表示方法利用過完備冗余字典取代傳統(tǒng)的正交基與多尺度緊框架在稀疏域表示圖像，過完備冗余字典充分考慮了微結(jié)構(gòu)圖像的冗余信息，可以更好地描述圖像的線及邊界等奇異特性，特別適合處理含有邊界與紋理區(qū)域的電子圖像.然而，冗余字典的設計是圖像稀疏表示的關(guān)鍵環(huán)節(jié)，字典的選擇直接影響稀疏向量的迭代計算與算法收斂性.因此，如何選擇合適的冗余字典是本文研究的問題之一.

文獻[9-10]提出了小波收縮去噪算法，其利用收縮函數(shù)對圖像小波域系數(shù)進行取舍，該方法的核心表現(xiàn)為：1)收縮曲線的選取;2)最優(yōu)閾值的選擇.目前，對于閾值的選取有Minimaxi閾值方法、Neighshrink閾值方法[11]以及SURE閾值方法[12].然而Neighshrink閾值與Minimaxi閾值多趨向于過扼殺小波變換系數(shù)，造成高頻信息流失；SURE閾值多趨向于過保留小波變換系數(shù)，造成去噪效果不明顯.

考慮到電子圖像結(jié)構(gòu)紋理信息與噪聲信息的頻率成分具有不確定性，且噪聲干擾成分很可能掩蓋圖像特征信息，導致在去噪時，圖像結(jié)構(gòu)信息成分與噪聲信息成分無法分離或相混疊，影響了圖像特征信息的提取與分析.

對于冗余字典的設計問題，Haar小波是具有對稱性和緊支撐特性的正交小波，且具有最優(yōu)的時空域分辨率，故選擇Harr小波冗余字典輔助圖像稀疏表示.對于最優(yōu)化的閾值計算問題，提出了一種基于Stein無偏風險估計準則的自適應閾值選擇方法，該方法可計算得到最優(yōu)化的閾值，有效平衡閾值選擇的過大與過小問題，提高了圖像信噪比.變分模態(tài)分解方法作為一種新的自適應信號與圖像分解技術(shù)，該方法實質(zhì)為多個自適應維納濾波器，利用迭代搜尋求解約束變分模型，實現(xiàn)圖像結(jié)構(gòu)紋理部分與噪聲成分的有效分離，最終得到若干個帶通分量.因此，在圖像去噪前，將二維變分模態(tài)分解方法引入到微結(jié)構(gòu)圖像分解中，實現(xiàn)圖像固有特征信息與噪聲的初步分離.

基于上述分析，提出了一種基于變分模態(tài)分解與Stein無偏風險估計方法相結(jié)合的電子圖像噪聲抑制方法，以鋁合金、雙相鋼與鈦合金Ti6Al4V 3種材料的晶粒取向EBSD圖像為研究對象，利用提出方法對模擬含噪圖像進行去噪，同時與一些主流的去噪方法如小波閾值方法、Neigh-Shrink方法、稀疏SURE方法以及KSVD方法進行比較分析，實驗結(jié)果證實了本文提出算法的優(yōu)越性.

1 變分模態(tài)分解

不同于傳統(tǒng)經(jīng)驗模態(tài)分解、局部均值分解與局部特征尺度分解等方法所使用的循環(huán)篩選與剝離的方式獲取內(nèi)稟模態(tài)分量或乘積模態(tài)分量，變分模態(tài)分解方法利用交替方向乘子算法(alternate direction method of multipliers,ADMM)不斷搜尋約束變分模型的最優(yōu)解，實現(xiàn)圖像的自適應分解，其整體框架是一求解變分問題，每個模態(tài)的估計帶寬與頻率中心在迭代求解變分模型的過程中不斷更新，最終將原始圖像自適應剖分為若干個模態(tài)函數(shù)之和[13].一維信號變分模態(tài)分解方法的約束變分問題可表達為

式中：f為原始信號；{uk}、{ωk}分別為分解得到的第k個模態(tài)的時域信號和中心頻率；δ(t)為Dirac函數(shù)；k為模態(tài)總數(shù)；符號*為卷積.

推廣至二維圖像分解領(lǐng)域，對應的約束變分問題可變換為

(1)

為求解上述變分約束模型的最優(yōu)解，引入二次罰函數(shù)項參數(shù)α和Lagrange乘子λ(t),變分約束模型(1)可簡化為

利用ADMM算法求解上述變分約束問題的最優(yōu)解[15-16]，可將圖像自適應分解為k個窄帶IMF分量.具體算法步驟如下[13-14].

步驟2執(zhí)行循環(huán),n=n+1.

其中， (1+sgn(ω·ωk))uk(ω)=

步驟5更新λ(ω),即

步驟6重復步驟 2～步驟5，直至滿足迭代停止條件，即

結(jié)束循環(huán).

為了有效地分離各個幾何圖案，且避免算法運算時間較長，將懲罰函數(shù)項參數(shù)α設為5 000，分量個數(shù)k設為5，利用二維變分模態(tài)分解算法(Bi-variational mode decomposition，BVMD)對合成重疊圖像進行分解.圖2(a)～圖2(e)為分解得到的5個二維內(nèi)稟模態(tài)分量，圖2(f)為重構(gòu)圖像，從圖中可以看出BVMD方法分解得到的模態(tài)分量可顯示清晰的橢圓圖案與矩形圖案.為了比較分解效果，利用二維經(jīng)驗模態(tài)分解方法(Bi-empirical mode decomposition,BEMD)對原始合成重疊圖像進行分解[17-18].圖3為利用BEMD方法分解得到的BIMF分量，可發(fā)現(xiàn)BEMD方法僅得到1個BIMF分量，即各個分量分解失敗.比較圖2、3可知， BEMD分解算法對于該重疊圖像完全沒有分解，即分解產(chǎn)生了嚴重的模態(tài)混疊現(xiàn)象；而BVMD方法不僅能有效去除偽分量，且每個IMF分量均為某一尺度范圍內(nèi)的的模態(tài)，彼此之間沒有模態(tài)混疊，實現(xiàn)了對圖像的多尺度表征.因此，二維變分模態(tài)分解方法分解效果優(yōu)于二維經(jīng)驗模態(tài)分解方法.

圖1 原始合成重疊圖像[13-14]

2 稀疏Stein無偏風險估計建模與推導

設y為含噪污染圖像，x為理想無污染圖像，A為冗余字典，圖像稀疏表示可通過求解x的歐幾里得范數(shù)表達為

(2)

即

將該解帶入到約束條件Ax=y中，可得：

圖2 二維變分模態(tài)分解方法分解得到的BIMF圖像[13-14]

圖3 二維經(jīng)驗模態(tài)分解方法分解得到的BIMF圖像

(3)

采用自適應閾值濾波算法，利用Haar小波構(gòu)造字典A，去噪過程可由下式表示：

式中:A+為字典A的Moore-Penrose逆;Sλ(z)為閾值標量算子;λ為閾值，即：當|z|<λ，Sλ(z)=0，否則，Sλ(z)=z. Haar小波冗余字典的構(gòu)造過程如下：首先，對Haar小波進行2層分解，即利用核函數(shù)[+1,+1]/2與[+1,-1]/2的二維濾波器進行水平濾波，然后，利用同樣核函數(shù)對2個輸出進行垂直濾波，這樣得到4副同樣大小的圖像，之后再對低通-低通(LL)圖像進行相同的處理，最終得到7副同樣大小的圖像，字典冗余度即為7∶1[19].

(4)

式中W為對角矩陣，對角線原子的范數(shù)為‖ai‖2.

展開上述表達式，得到:

(5)

可知式(5)第2項為未知量y0的函數(shù)，將y0=y-v帶入式(5)第2項，可得:

E[h(y,θ)T·y0]=

E[h(y,θ)T·y]-E[h(y,θ)T·v].

設h(y,θ)的第i元素為hi(y,θ)，噪聲v的第i元素為vi，計算:

E[h(y,θ)T·y0]=

(6)

E[hi(y,θ)·vi]=

(7)

E[hi(y,θ)·vi]=

綜上所述，最終最小均方誤差為

2σ2·h(y,θ)]+Const.

(8)

(9)

(10)

(11)

求導可得：

(12)

由于硬閾值收縮曲線不具有連續(xù)可導性質(zhì)，為了保證收縮曲線連續(xù)且具有無窮階導數(shù)，利用偶數(shù)K對硬閾值收縮曲線進行平滑逼近.本文選取K=20所對應的閾值曲線，該收縮曲線連續(xù)且具有無窮階導數(shù)，便于數(shù)值計算[19].

將式(10)～式(12)帶入式(9)中，可得：

(13)

利用A1=AW-1,A2=AW,式(13)第3項可簡化為

(14)

3 結(jié)果與分析

實驗采用牛津儀器公司(http://www.ebsd.cn)公開提供的鋁合金、雙相鋼與鈦合金Ti6Al4V 3種金屬材料的電子微結(jié)構(gòu)圖像進行去噪分析.首先對材料進行去油污、打磨，然后利用丙酮溶液浸泡、機械-化學拋光等技術(shù)對表面進行預處理，最后利用電子背散射衍射儀 (AZtecHKL系統(tǒng)) 采集得到鋁合金，雙相鋼與鈦合金Ti6Al4V 3種樣品晶粒取向成像圖.由于電子背散射衍射儀在圖像采集系統(tǒng)內(nèi)部可能由電子束電壓或電流不穩(wěn)定等原因會引起噪聲，而這種噪聲可能是高斯噪聲或斑紋噪聲等多個噪聲疊加而成，故本實驗在理想晶粒取向圖像中加入高斯噪聲和Speckle斑紋噪聲，模擬含噪晶粒取向成像圖.圖4(a)為理想鋁合金晶粒取向圖像，圖4(b)為疊加了噪聲標準差為30的高斯噪聲與Speckle斑紋噪聲的灰度圖像，圖像大小為391×391，可知被污染的圖像質(zhì)量嚴重下降、邊緣特征較模糊.

圖4 鋁合金原始晶粒取向圖像與加噪圖像

Fig.4 Original crystal orientation image of aluminium alloy and its noisy image

利用BVMD方法對含噪圖像進行分解，模態(tài)個數(shù)k設為5，懲罰函數(shù)項參數(shù)α設為5 000.圖5為分解得到的5個BIMF模態(tài)分量結(jié)果，由圖5分解效果可看出BVMD方法能很好地克服模態(tài)混疊，實現(xiàn)了對加噪圖像的多尺度表征，其中，前4個模態(tài)分量為高頻分量或者噪聲信息，而BIMF5模態(tài)分量包含了原始圖像的固有特征信息，即BVMD方法實現(xiàn)了圖像固有特征信息與干擾噪聲的初步分離.由于BIMF1～BIMF4模態(tài)分量干擾斑紋過多，特征頻率成分與原始固有特征無關(guān)，因此選取BIMF5模態(tài)分量，利用稀疏SURE方法進行去噪分析.

為驗證本文提出算法的有效性，分別利用小波閾值去噪方法、Neigh鄰域閾值收縮法(Neigh shrink)[22],單一稀疏SURE方法以及文獻[23-24]提出的KSVD方法進行對比分析.字典A由冗余 Haar小波字典構(gòu)成，字典冗余度為7，收縮函數(shù)為圖5所示，偶數(shù)K取20，采用歸一化均方誤差(normalized mean squared error,NMSE)、峰值信噪比(PSNR)與結(jié)構(gòu)相似度(structural similarity,SSIM)參數(shù)作為判定標準，衡量以上方法的去噪復原效果.各參數(shù)定義如下：

圖5 二維變分模態(tài)分解得到的BIMF分量 (算法執(zhí)行時間639.200 995 s)

選擇BIMF5模態(tài)分量作為稀疏SURE輸入圖像，圖6為PSNR參數(shù)隨閾值λ變化而變化的趨勢，圖中虛線為含噪圖像對無污染圖像求得的PSNR值，可看出本文提出的方法(即BVMD+稀疏SURE)得到的PSNR曲線與真實PSNR曲線十分吻合.因此，可知BVMD+稀疏SURE曲線的峰值點(80,23.593 5)所對應的去噪效果最佳，最佳閾值為80.

圖6 稀疏收縮PSNR與閾值λ的變化關(guān)系(峰值點：80,23.593 5)

Fig.6 Relationship between sparse contraction PSNR and thresholdλ(maximum point: 80, 23.593 5)

圖7(a)～圖7(d)分別為利用小波閾值去噪方法、Neigh鄰域閾值收縮法、單一稀疏SURE方法、KSVD方法與本文方法得到的去噪效果圖，對于小波閾值去噪算法：本文利用小波軟閾值去噪算法進行去噪，軟閾值函數(shù)為

圖7 各類方法得到的鋁合金晶粒取向圖像去噪效果

表1 各類方法的去噪性能指標對比

觀察5種方法降噪圖像的視覺質(zhì)量，小波閾值去噪算法生成的圖像比較模糊，依然殘存大量的噪聲斑點，去噪效果差；Neigh-Shrink方法提高了圖像的去噪效果，但在圖像邊緣和細節(jié)處產(chǎn)生一定程度的平滑現(xiàn)象，使得輪廓不夠清晰，原因在于該方法進行閾值處理是考慮其周圍系數(shù)的分布特點；單一稀疏SURE方法與本文方法生成圖像較好的保存了紋理細節(jié)，沒有偽像，各類評價指標很接近，但單一稀疏SURE方法得到的圖像整體比較平滑，晶界部分較模糊，而本文方法較好的克服了這一點，晶界線條更為清晰，視覺效果有明顯改善，這驗證了BVMD方法在微結(jié)構(gòu)圖像預分解起到了重要作用.在客觀衡量指標上，由表1可看出，本文提出的方法得到的峰值信噪比最高，達到23.593 5 dB，也突破了稀疏SURE收縮曲線得到最大值(23.445 1 dB)，比Neigh-Shrink方法高0.39 dB，比KSVD方法高2.895 dB，比小波閾值去噪算法高3.07 dB.

為了考察本文方法相對于其他幾種方法對噪聲的魯棒性，選取標準差分別為10,20,30,…,90,100高斯噪聲與Speckle斑紋噪聲的10副微結(jié)構(gòu)晶粒取向成像圖作為去噪對象，采用峰值信噪比PSNR作為判定標準.圖8顯示了各類方法的去噪效果PSNR隨噪聲標準差變化的趨勢，可以看出，本文方法的PSNR相對于其他3種算法具有明顯的優(yōu)勢，尤其是當噪聲強度增大時優(yōu)越性更加明顯；小波閾值方法的去噪效果最差，這是由于閾值選擇缺乏自適應性，獲取的小波系數(shù)與真實小波系數(shù)有偏差，造成重構(gòu)圖像精度降低，產(chǎn)生振鈴與偽Gibbs效應，導致圖像模糊；當噪聲標準差大于70時，Neigh-Shrink方法的去噪效果迅速降低，這是由于Neigh-Shrink方法是依賴小波變換所得細節(jié)子帶系數(shù)與其鄰域圖像結(jié)構(gòu)信息之間相關(guān)分布關(guān)系，噪聲級數(shù)變大，鄰域圖像結(jié)構(gòu)信息附帶的噪聲信息增多，導致去噪效果降低；另外，從整體來看，KSVD算法隨噪聲標準差變化幅度下降較快，在噪聲方差較低時，KSVD冗余字典可有效刻畫微結(jié)構(gòu)圖像的細節(jié)信息，重構(gòu)精度高，但隨著噪聲方差進一步增大，KSVD字典過于冗余，字典含有較多無用信息，圖像邊緣細節(jié)信息無法表達，去噪效果大大降低.相比單純稀疏SURE方法，本文提出的方法由于預先對含噪圖像利用BVMD分解，去除了高頻無用噪聲，圖像噪聲抑制效果優(yōu)于單純稀疏SURE方法.

圖8 各類方法對晶粒取向圖像去噪PSNR值隨標準差-Sigma變化曲線

Fig.8 Variation of crystal orientation image denoising PSNR value with noise standard variance-Sigma under different methods

為了進一步驗證提出方法的普適性，分別引入牛津儀器(http:www.ebsd.cn)公開的雙相鋼樣品的EBSD相分布微結(jié)構(gòu)圖像 (奧氏體和鐵素體) 與鈦合金Ti6Al4V材料微結(jié)構(gòu)EBSD圖像再次進行去噪分析，在理想微結(jié)構(gòu)圖像中分別加入噪聲標準差為30的高斯噪聲和Speckle斑紋噪聲，圖9(a)、圖9(b)分別為雙相鋼樣品的EBSD相分布微結(jié)構(gòu)圖像與噪聲污染圖像，圖10(a)、圖10(b)分別為鈦合金Ti6Al4V材料微結(jié)構(gòu)圖像與噪聲污染圖像.同樣利用BVMD分解方法、構(gòu)造Harr小波冗余字典以及Stein無偏風險估計方法進行去噪，最終兩種材料的去噪圖像分別如圖9(c)、圖10(c)所示，由于篇幅所限，此處省略BVMD分解結(jié)果.對比去噪前與去噪后的圖像可知，雙相鋼材料圖像的峰值信噪比由加噪前的18.598 1 dB提升至25.784 8 dB，鈦合金Ti6Al4V材料圖像的峰值信噪比由加噪前的18.569 9 dB提升至24.587 5 dB，噪聲圖像邊緣特征模糊現(xiàn)象得到改善，去噪后的圖像邊界紋理信息更加清晰，視覺效果大大提升.圖11、圖12分別為利用小波軟閾值方法、Neigh-Shrink方法、稀疏SURE方法與KSVD方法對雙相鋼材料與鈦合金材料進行去噪得到的去噪圖像對比圖，結(jié)果進一步驗證了本文所提出算法的優(yōu)越性.

圖9 雙相鋼材料的EBSD微結(jié)構(gòu)圖像、噪聲污染圖像與本文方法去噪圖像

圖10 鈦合金Ti6Al4V材料微結(jié)構(gòu)EBSD圖像、噪聲污染圖像與本文方法去噪圖像

圖11 各類對比方法得到的雙相鋼材料晶粒取向圖像去噪結(jié)果

圖12 各類對比方法得到的鈦合金 Ti6Al4V材料晶粒取向圖像去噪結(jié)果

4 結(jié) 論

1)針對電子微結(jié)構(gòu)圖像的噪聲干擾問題，提出了一種基于變分模態(tài)分解與稀疏Stein無偏風險估計方法相結(jié)合的去噪方法，具體以鋁合金，含奧氏體和鐵素體的雙相鋼與鈦合金Ti6Al4V的背散射衍射含噪模擬圖像為例進行降噪分析.

2)利用二維變分模態(tài)分解方法可有效地將含噪電子微結(jié)構(gòu)圖像按照不同尺度、不同頻率分離為固有特征信息成分與高頻噪聲信息成分.

3)推導了基于稀疏Stein無偏風險估計的目標優(yōu)化收縮曲線，并用黃金分割法搜索得到全局最佳自適應閾值，該閾值的選取接近真實峰值信噪比曲線的最大點.圖像稀疏表示過程中利用了Haar小波冗余字典，其能充分考慮微結(jié)構(gòu)圖像的結(jié)構(gòu)冗余模式，可以更好地刻畫圖像的邊界、紋理等固有特性.

4)結(jié)果表明,無論在客觀去噪指標還是視覺效果方面，相比其他去噪方法，本文提出的方法均顯示較優(yōu)越的去噪性能. 以鋁合金為例，證明了當噪聲標準差為30時，該方法去噪結(jié)果的峰值信噪比PSNR值高于Neigh-Shrink去噪方法0.39 dB，比KSVD方法高2.895 dB，高于小波閾值去噪算法3.07 dB，較單純稀疏SURE方法更大程度提高了圖像視覺清晰度；并隨著噪聲標準差級別的增大，圖像去噪魯棒性較高，為海量的微結(jié)構(gòu)圖像去噪提供了新的思路.

5)在BVMD預分解圖像的過程中，模態(tài)個數(shù)k與二次罰函數(shù)項α的選擇具有一定的經(jīng)驗性，且交替方向乘子算法求解變分約束問題最優(yōu)解的循環(huán)計算，Haar小波冗余字典的構(gòu)造，導致整體抑噪運行速率較慢，仍具有較大的改進空間.