余麗玲 陽 維* 馮衍秋 劉 閩 馮前進 陳武凡

1(南方醫(yī)科大學生物醫(yī)學工程學院，廣州 510515)

2(深圳出入境檢驗檢疫局工業(yè)品檢測技術(shù)中心，深圳 518067)

引言

MRI技術(shù)提供了豐富的、有價值的人體結(jié)構(gòu)信息。然而，MR圖像在采集的過程中經(jīng)常被隨機噪聲干擾，MR圖像中的噪聲一般認為服從Rician分布，是一種與信號相關(guān)的非加性噪聲。此外，隨著多線圈并行采集技術(shù)(敏感編碼、廣義自動校準部分并行采集)的采用［1］，MR圖像中的Rician噪聲水平在空間上可能不均勻。噪聲不僅會影響醫(yī)生的觀察，還會影響后續(xù)處理和分析的效果，如圖像分割、配準、可視化等，抑制MR圖像中的Rician噪聲是后續(xù)處理和分析的重要步驟。

現(xiàn)有的Rician噪聲抑制方法，通常假設(shè)噪聲水平在空間上為常數(shù)。經(jīng)過適當?shù)母脑?，非局部均值?］和 BM3D［3］算法，已被研究人員應(yīng)用于 Rician噪聲的抑制。其中，F(xiàn)oi等提出了利用方差穩(wěn)定變換(variance－stabilization transformations，VST)進行Rician噪聲抑制的方法［4］，利用估計的全局噪聲水平，對MR圖像幅值進行變換，使得變換后圖像中的噪聲與信號是獨立的，進而采用BM3D算法進行噪聲抑制，然而該方法僅適用于噪聲水平為常數(shù)的Rician噪聲圖像。Manjón等提出了依據(jù)局部Rician噪聲水平進行自適應(yīng)調(diào)整的非局部均值［5］(adaptive non－local means，ANLM)算法，用于抑制 MR圖像中空間變化的Rician噪聲，取得了較好的去噪效果，但該方法較為耗時，且需調(diào)整的參數(shù)較多。

針對空間變化的Rician噪聲，本研究的思路是:對空間變化的噪聲水平建模并進行估計，然后對圖像各處的幅值依據(jù)不同的局部噪聲水平進行方差穩(wěn)定變換，進而結(jié)合有效的去噪算法抑制這種空間變化的Rician噪聲。對于研究目的而言，其關(guān)鍵是如何有效估計MR圖像中不同位置的噪聲水平。

MR圖像中Rician噪聲水平的估計，較為常用的是基于背景區(qū)域的方法［6］。假設(shè)背景區(qū)域的信號值為零，采用最大似然或者幅值圖像的二階矩進行估計Rician噪聲的水平，但不適用于空間變化噪聲水平的估計。本研究提出了一種基于稀疏性約束的Rician噪聲水平場的估計方法和空間自適應(yīng)方差穩(wěn)定變換方法，結(jié)合BM3D算法，有效地實現(xiàn)了對空間變化Rician噪聲的抑制。

1 方差穩(wěn)定變換

1.1 Rician噪聲分布

MR成像過程中，原始K空間數(shù)據(jù)被復(fù)數(shù)高斯噪聲干擾。由于Fourier變換的正交性和線性，經(jīng)Fourier逆變換后圖像實部和虛部的噪聲仍呈高斯分布。這樣，MR幅值圖像為兩個獨立高斯隨機變量的平方根，其噪聲不再服從高斯分布，而是服從Rician分布。設(shè)M為MR幅值，則M服從分布［2］為

式中，S為無噪聲時的幅值，σ為實圖像和虛圖像中高斯噪聲的標準偏差，I0為零階第一類貝塞爾函數(shù)。Rician噪聲是信號相關(guān)的，既不是加性的也不是乘性的，針對加性高斯白噪聲的去噪方法不能直接應(yīng)用于MR圖像的去噪。

Rician噪聲抑制的一種可行方案，是利用Rician分布的最大似然在最大后驗(maximum a posteriori，MAP)框架下對 MR圖像進行恢復(fù)，但MAP框架一般僅可融入有限的先驗，而且其優(yōu)化過程相對耗時。另一種可行方案，是先將圖像進行變換，使得在變換域中噪聲與信號獨立，在變換域中對去噪后圖像進行逆變換，如小波變換、方差穩(wěn)定變換。本研究采用方差穩(wěn)定變換。方差穩(wěn)定變換可使圖像中的噪聲與信號近似不相關(guān)［4］，并且可通過相應(yīng)的逆變換得到信號的無偏估計。

1.2 方差穩(wěn)定變換

方差穩(wěn)定變換［4，7］是為了使異方差數(shù)據(jù)轉(zhuǎn)換為同方差的，更容易地運用標準方法來解決問題，去除Rician噪聲方差和無噪聲圖像幅值的依賴性。

圖1為所采用的方差穩(wěn)定變換及其對應(yīng)的逆變換，其中f為方差穩(wěn)定變換，D為去噪后逆變換之前的圖像，Vf為逆變換。對MR幅值圖像進行方差穩(wěn)定變換處理后，使得σf(θ)=std{f(M)|S，σ}=1，圖像各處噪聲的方差基本一致，這樣可將VST變換后的Rician噪聲作為加性噪聲進行抑制處理。對方差穩(wěn)定變換圖像進行噪聲抑制后，需經(jīng)過對應(yīng)的方差穩(wěn)定逆變換，才能得到無偏的去噪圖像。

圖1 方差穩(wěn)定變換及其逆變換。(a)方差穩(wěn)定變換;(b)方差穩(wěn)定逆變換Fig.1 VST and inverse VST.(a)VST;(b)Inverse VST

2 噪聲水平場估計

對于水平空間變化的Rician噪聲，依據(jù)局部的噪聲水平逐像素進行VST，即實現(xiàn)了空間自適應(yīng)VST?？臻g各處不同的噪聲水平在空間上構(gòu)成一個場，稱之為噪聲水平場(noise level field，NLF)。顯然，進行空間自適應(yīng)VST的關(guān)鍵在于NLF的準確估計。利用NLF進行空間自適應(yīng)VST和去噪的流程如圖2所示。主要包括4個步驟:估計NLF;利用NLF進行空間自適應(yīng) VST;對變換的圖像使用BM3D算法進行去噪;對去噪后圖像進行VST逆變換。

圖2 空間自適應(yīng)VST和去噪的流程圖Fig.2 The flow chart of spatially adaptive VST and denosing

2.1 噪聲水平場模型

采用多項式描述噪聲水平場［8］

式中，{at，l}為多項式系數(shù)，xi，yi為像素點 i的坐標，K為多項式階數(shù)。為了方便計算，{at，l}用列向量A表示，{}用列向量C表示，B表示估計的噪聲水平場。式(2)可以擴展到三維的NLF。

式中，λ為設(shè)定的正則化系數(shù)，B'為噪聲水平的局部估計值，其估計方法將在下一小節(jié)做詳細介紹。式(3)為l1正則化約束優(yōu)化問題，解決此類問題的優(yōu)化方法很多［9］，采用 L1_LS 算法［10］求解多項式系數(shù)A，進而由式(2)得到NLF。L1_LS算法結(jié)合了截斷牛頓內(nèi)點法、預(yù)條件共軛梯度算法，比普通的內(nèi)點法使用方向或共軛梯度法計算效率要高，求解時間要短。

2.2 局部Rician噪聲水平的估計方法

由式(3)可見，為了估計參數(shù)A，需要得到噪聲水平的局部估計。傳統(tǒng)的噪聲水平估計方法一般假設(shè)噪聲水平是空間不變的常數(shù)。對于MR圖像中的Rician噪聲，較為常用的是基于背景區(qū)域的估計方法，但這些方法對于空間變化的噪聲水平不再適用。假設(shè)圖像局部的噪聲水平為常數(shù)，采用修正的中值絕對偏差(median absolute deviation，MAD)估計方法估計圖像局部的噪聲水平［11］。

假設(shè)噪聲服從高斯分布，對局部噪聲水平進行MAD 估計［12－13］

式中，i，j表示圖像中的像素點，Ni，Nj表示點 i，點 j的鄰域，S(i)，S(j)是像素點的信號幅值，median(S)為中值濾波，為估計的噪聲標準偏差。與其它噪聲方差估計方法相比，MAD估計不但具有較強的魯棒性，而且運算量較小。當信噪比足夠大(＞5)時，Rician分布的標準偏差估計與高斯分布的趨于一致，無須進行修正;當信噪比較小時，式(4)的估計將會產(chǎn)生一定的偏差，須對估計的進行修正。這種修正方法是基于對Rician噪聲信噪比(signal noise ratio，SNR)的迭代估計［14］。用 MAD估計出的標準偏差初始化修正過程為

式中，θ為SNR值，ζ(θ)為修正因子，定義如下

式中，I0，I1分別是零階、一階修正的貝塞爾函數(shù)。修正因子ζ(θ)通過迭代估計直到收斂或者達到給定的迭代次數(shù)。使用|θt－θt－1|作為收斂條件，迭代過程可表示為

式中，〈S〉是所給像素信號幅值的平均值，t為迭代次數(shù)。在梯度大的像素點上估計的局部噪聲水平與實際噪聲水平相比，存在著較大偏差。因此，在式(3)估計噪聲水平場時，不使用梯度大和修正后信噪比仍然較小的局部噪聲水平估計。

2.3 空間自適應(yīng)方差穩(wěn)定變換去噪

利用估計的NLF對圖像中各像素進行方差穩(wěn)定變換

式中，β為可調(diào)節(jié)參數(shù)，B為估計的噪聲水平場，i，j表示圖像的位置。式(8)實際上相當于對MR圖像進行同態(tài)化和歸一化處理，使得整幅圖像的噪聲水平處處近似為1，噪聲與MR幅值和空間位置近似獨立。這樣，可采用針對加性高斯白噪聲的去噪方法進行處理，然后進行方差穩(wěn)定逆變換，最終得到無偏的去噪圖像。采用BM3D(BM4D)算法對方差穩(wěn)定變換后的圖像進行去噪。BM3D是當前公認對加性高斯白噪聲去噪性能良好的算法，它包含了非局部去噪的思想，同時又用到了變換域濾波的方法。通過空間自適應(yīng)VST，可有效利用BM3D的去噪能力，實現(xiàn)空間變化Rican噪聲的抑制。

3 實驗

為驗證所提出方法的有效性，進行仿真MR圖像實驗和真實MR圖像實驗，并與VST－BM3D、VSTBM4D和ANLM算法進行比較。為方便起見，本方法簡記為SAVST (spatially adaptivevariancestabilization transformations)。

3.1 仿真實驗

仿真實驗中，多項式的階數(shù)K設(shè)定為5，正則化系數(shù)λ為1，β為1.2。仿真的噪聲水平場由離散點三次方插值生成［5］，仿真的MR噪聲圖像的峰值信噪比(Peak signal to noise ratio，PSNR)為24 dB。

式中，RMSE為原始無噪聲圖像與去噪后圖像的均方差根誤差。

采用本方法估計噪聲水平場如圖3所示。由圖3可以看到，估計的噪聲水平場較好地逼近實際的噪聲水平場。噪聲水平場估計的精度采用均方根誤差(root mean square error，RMSE)和平均相對誤差(mean relative error，MRE)進行度量。表1列出對于3種不同的仿真噪聲水平場的估計精度，平均相對誤差小于0.2%。

圖3 仿真噪聲水平場的估計結(jié)果。(a)噪聲圖像;(b)仿真的噪聲水平場;(c)估計的噪聲水平場Fig.3 The estimated results of simulated NLF.(a)The noisy image;(b)The simulated NLF;(c)The estimated NLF

表1 噪聲水平場估計的精度Tab.1 The precision of NLF estimation

仿真實驗圖像數(shù)據(jù)來自BrainWeb。對T1加權(quán)圖像、T2加權(quán)圖像，添加仿真 Rician噪聲水平從1%到15%進行實驗，仿真的NLF如圖4(b)所示。采用PSNR評價去噪效果。

圖4 不同方法噪聲估計結(jié)果。(a)原始圖像;(b)仿真噪聲水平場;(c)噪聲圖像;(d)仿真的噪聲;(e)VSTBM3D的去噪效果;(f)VST-BM3D估計的噪聲;(g)本方法的去噪效果;(h)本方法估計的噪聲Fig.4 The estimated results of different methods.(a)Original image;(b)Simulated NLF;(c)Noisy image;(d)Simulated noise;(e)Denoised result of VST-BM3D;(f)Estimated noise of VST-BM3D;(g)Denoised result of the proposed method;(h)Estimated noise of the proposed method

圖4顯示了仿真噪聲水平為3%時，兩組圖像使用不同方法得到的去噪結(jié)果。其中，VST－BM3D使用噪聲水平為常數(shù)的方差穩(wěn)定變換去噪方法，圖4第一行為T1加權(quán)圖像的仿真實驗結(jié)果，第二行為T2加權(quán)圖像的仿真實驗結(jié)果。由圖4可見，本方法對于噪聲水平在空間上變化的Rician噪聲抑制效果較好，產(chǎn)生的模糊較小，并能保持原圖像的細節(jié)特征和邊緣信息。而VST－BM3D在信噪比較低的區(qū)域(如圖4(e)腦室區(qū)域)，噪聲未被有效抑制。對T1加權(quán)圖像、T2加權(quán)圖像采用本文方法(SAVST)和VST－BM3D方法進行去噪實驗，兩種去噪方法在不同噪聲水平下，對應(yīng)的去噪后PSNR值如圖5所示。其中，橫軸表示加噪水平，從1% ～15%，縱軸表示去噪后圖像的PSNR值。由圖可知，與VST－BM3D方法相比，本方法對應(yīng)的PSNR約提高了2 dB。

圖5 二維圖像去噪性能比較。(a)T1圖像;(b)T2圖像Fig.5 Comparison of denoising performance using different methods.(a)T1 image;(b)T2 image

本方法可應(yīng)用在三維MR圖像上，對同一噪聲水平場的三維腦部 MR圖像，采用 VST－BM4D、ANLM和本方法進行去噪實驗。圖6為在不同噪聲水平下，幾種方法去噪后PSNR的比較。從圖5和圖6可看出，本方法對噪聲水平空間變化的二維和三維MR圖像進行去噪，對應(yīng)的PSNR明顯高于VST－BM3D方法。ANLM結(jié)合了噪聲水平的局部估計，其去噪性能也略高于VST－BM3D方法，但相對耗時。對于圖像大小為181像素×217像素×181像素的腦部數(shù)據(jù)，本方法所需的運行時間約為340 s，而 VST－BM3D 方法為 378 s。

圖6 三維圖像去噪性能比較。(a)T1圖像;(b)T2圖像Fig.6 Comparison of denoising performance using different method.(a)T1 image;(b)T2 image

3.2 真實MR圖像實驗

所用真實的乳腺MR圖像，圖像大小為384像素×384像素，像素大小為1 mm×1 mm。真實乳腺MR圖像的噪聲水平場估計見圖7。

圖7 乳腺MR圖像噪聲水平場估計。(a)噪聲圖像;(b)估計的噪聲水平場Fig.7 The experimental results of breast MR/NLF estimation.(a)Noisy image;(b)Estimated NLF

由于真實噪聲圖像無對應(yīng)的參考圖像，為了評價去噪效果，采用Q度量［15］評價去噪圖像的圖像質(zhì)量。當參考圖像不能獲得的時候，Q度量可較好地評價圖像質(zhì)量。Q度量值越大，圖像質(zhì)量越好。為了使去噪效果達到最優(yōu)，利用Q度量調(diào)節(jié)式(8)中的參數(shù)β。圖8顯示了對兩幅乳腺MR圖像使用不同方法進行去噪的結(jié)果，每一行代表一幅含噪聲乳腺MR圖像的實驗結(jié)果，相應(yīng)的Q值如圖8(f)所示，可看出當參數(shù)β在1.2左右時，Q值達到最大值，去噪效果趨于最優(yōu)。本去噪方法對應(yīng)的Q值高于VST－BM3D，而且對于細節(jié)信息和邊沿結(jié)構(gòu)有更好的保存能力。

4 討論

本研究要解決的問題是抑制MR圖像中空間變化的Rician噪聲。提出的方法是通過估計Rician噪聲水平場和方差穩(wěn)定變換，將Rician噪聲轉(zhuǎn)化為常數(shù)水平的加性噪聲，然后使用一般的去噪算法實現(xiàn)對噪聲的抑制。本文選取的去噪算法為BM3D算法(因其性能良好)，其他針對加性高斯白噪聲的去噪算法，如非局部均值算法、高斯混合尺度模型(scale mixtures of Gaussians，GSM)算法［16］也可以替代BM3D算法。

圖8 兩幅真實MR圖像的處理結(jié)果(上行為一幅，下行為另一幅)。(a)噪聲圖像;(b)VST-BM3D的去噪效果;(c)(a)與(b)之間的殘差圖像;(d)本方法的去噪效果，(e)(a)與(d)之間的殘差圖像;(f)調(diào)節(jié)β和Q值的變化Fig.8 The experimental results of two real MR images(The top row concerns one image，the bottom row concerns the other one).(a)Noisy image;(b)Result of VST-BM3D;(c)Residual image between(a)and(b);(d)Result of our proposed method;(e)Residual image between(a)and(d);(f)The corresponding Q values with varying β

噪聲水平場的準確估計對于去噪效果有很大影響。文中，局部噪聲水平使用了修正的MAD方法進行估計。MAD估計對服從高斯分布的噪聲具有較高的精度，但當信噪比較低時，Rician分布不服從高斯分布，需要對估計的噪聲標準方差進行修正，才能得到無偏的噪聲標準方差，從而提高了噪聲水平場的估計精度。估計噪聲水平場時，相關(guān)參數(shù)的設(shè)置對估計精度有較大影響，當多項式階數(shù)K較小時，會得到偏低的估計，當K較大時，會得到偏高的估計;而正則化系數(shù)λ大小決定了噪聲水平場的平滑程度，需要依據(jù)經(jīng)驗或先驗知識進行調(diào)節(jié)，對于實際問題，可通過優(yōu)化Q度量確定最優(yōu)的λ。

5 結(jié)論

針對MR圖像中噪聲水平空間變化的Rician噪聲，本研究提出了一種噪聲水平場的估計方法，并用于抑制MR圖像中空間變化的Rician噪聲。通過Rician噪聲水平的局部估計和稀疏性約束模型估計噪聲水平場，然后對圖像各處的幅值依據(jù)不同的局部噪聲水平進行方差穩(wěn)定變換，使得噪聲與信號幅值和空間位置無關(guān)，然后利用BM3D算法的強大去噪能力，實現(xiàn)了對空間變化Rician噪聲的有效抑制。實驗結(jié)果表明，本方法能有效估計MR圖像中的Rician噪聲水平場，與使用全局噪聲水平估計的方差穩(wěn)定變換去噪方法、自適應(yīng)非局部均值去噪算法相比，本方法對應(yīng)的峰值信噪比和 Q度量均較高。

［1］陳武凡.并行磁共振成像的回顧、現(xiàn)狀與發(fā)展前景［J］.中國生物醫(yī)學工程學報，2005，24(6):649－654.

［2］Buades A，Coll B，Morel JM.A review of image denoising algorithm，with a new one［J］.SIAM Multiscale Modeling and Simulation，2005，4(2):490－530.

［3］Dabov K，F(xiàn)oi A，Katkovnik V，at al.Image denoising by sparse3D transform－domain collaborative filtering ［J］.IEEE Trans Image Process，2007，16(8):2080－2095.

［4］Foi A.Noise estimation and removal in MR imaging:the variance－stabilization approach ［C］// IEEE International Symposium on Biomedical Imaging.Chicago:IEEE，2011:1809－1814.

［5］Manjón JV，Coupé P，Martí－ Bonmatí L，et al.Adaptive nonlocal means denoising of MR images with spatially varying noise levels［J］.Magn Reson Imaging，2010，31(1):192 －203.

［6］Sijbers J，Poot D，Dekker AJ，et al.Automatic estimation of the noise variance from the histogram of a magnetic resonance image［J］.Phys Med Biol，2007，52(5):1335 －1348.

［7］FoiA. Optimization ofvariance－stabilizing transformations［OL］.http://www.cs.tut.fi/～foi/.

［8］Zheng Yuanjie.Estimation of image bias field with sparsity constraints［C］//Computer Vision and Pattern Recognition(CVPR).San Francisco:IEEE Computer Society，2010:255 －262.

［9］王浩.求解正則化問題的算法研究［D］.北京:北京航空航天大學，2010.

［10］Kim SJ，Koh K，Lustig M，at al.An interior－point method for large－scale l1－regularized leastsquares ［J］.IEEE Signal Processing，2007，1(4):606 －617.

［11］Rousseeuw PJ，Verboven S.Robust estimation in very small samples［J］.Comp Stat Data Anal，2002，40(4):741 －758.

［12］Maximov II，F(xiàn)arrher E，Grinberg F，at al.Spatially variable Rician noise in magnetic resonance imaging［J］.Med Image Anal，2012，16(2):536 －548.

［13］Coupé P，Manjón JV，Gedamu E，at al.Robust Rician noise estimation for MR images［J］.Med Image Anal，2010，14(4):483－493.

［14］Koay CG，Basser PJ.Analytically exact correction scheme for signal extraction from noisy magnitude MR signals［J］.J Magn Reson，2006，179(2):317－322.

［15］Zhu Xiang，Milanfar P. Automatic parameter selection fordenoising algorithms using a no－reference measure of image content［J］.IEEE Trans Image Process，2010，19(12):3116 －3132.

［16］Portilla J，Strela V，Wainwright M，et al.Image senoising using scale mixtures of gaussians in the wavelet domain ［J］.IEEE Trans on Image Process，2003，12(11):1338－1351.