Rician噪声水平场的估计及其在MR图像去噪中的应用

余丽玲 阳 维* 冯衍秋 刘 闽 冯前进 陈武凡

1(南方医科大学生物医学工程学院，广州 510515)

2(深圳出入境检验检疫局工业品检测技术中心，深圳 518067)

引言

MRI技术提供了丰富的、有价值的人体结构信息。然而，MR图像在采集的过程中经常被随机噪声干扰，MR图像中的噪声一般认为服从Rician分布，是一种与信号相关的非加性噪声。此外，随着多线圈并行采集技术(敏感编码、广义自动校准部分并行采集)的采用［1］，MR图像中的Rician噪声水平在空间上可能不均匀。噪声不仅会影响医生的观察，还会影响后续处理和分析的效果，如图像分割、配准、可视化等，抑制MR图像中的Rician噪声是后续处理和分析的重要步骤。

现有的Rician噪声抑制方法，通常假设噪声水平在空间上为常数。经过适当的改造，非局部均值［2］和 BM3D［3］算法，已被研究人员应用于 Rician噪声的抑制。其中，Foi等提出了利用方差稳定变换(variance－stabilization transformations，VST)进行Rician噪声抑制的方法［4］，利用估计的全局噪声水平，对MR图像幅值进行变换，使得变换后图像中的噪声与信号是独立的，进而采用BM3D算法进行噪声抑制，然而该方法仅适用于噪声水平为常数的Rician噪声图像。Manjón等提出了依据局部Rician噪声水平进行自适应调整的非局部均值［5］(adaptive non－local means，ANLM)算法，用于抑制 MR图像中空间变化的Rician噪声，取得了较好的去噪效果，但该方法较为耗时，且需调整的参数较多。

针对空间变化的Rician噪声，本研究的思路是:对空间变化的噪声水平建模并进行估计，然后对图像各处的幅值依据不同的局部噪声水平进行方差稳定变换，进而结合有效的去噪算法抑制这种空间变化的Rician噪声。对于研究目的而言，其关键是如何有效估计MR图像中不同位置的噪声水平。

MR图像中Rician噪声水平的估计，较为常用的是基于背景区域的方法［6］。假设背景区域的信号值为零，采用最大似然或者幅值图像的二阶矩进行估计Rician噪声的水平，但不适用于空间变化噪声水平的估计。本研究提出了一种基于稀疏性约束的Rician噪声水平场的估计方法和空间自适应方差稳定变换方法，结合BM3D算法，有效地实现了对空间变化Rician噪声的抑制。

1 方差稳定变换

1.1 Rician噪声分布

MR成像过程中，原始K空间数据被复数高斯噪声干扰。由于Fourier变换的正交性和线性，经Fourier逆变换后图像实部和虚部的噪声仍呈高斯分布。这样，MR幅值图像为两个独立高斯随机变量的平方根，其噪声不再服从高斯分布，而是服从Rician分布。设M为MR幅值，则M服从分布［2］为

式中，S为无噪声时的幅值，σ为实图像和虚图像中高斯噪声的标准偏差，I0为零阶第一类贝塞尔函数。Rician噪声是信号相关的，既不是加性的也不是乘性的，针对加性高斯白噪声的去噪方法不能直接应用于MR图像的去噪。

Rician噪声抑制的一种可行方案，是利用Rician分布的最大似然在最大后验(maximum a posteriori，MAP)框架下对 MR图像进行恢复，但MAP框架一般仅可融入有限的先验，而且其优化过程相对耗时。另一种可行方案，是先将图像进行变换，使得在变换域中噪声与信号独立，在变换域中对去噪后图像进行逆变换，如小波变换、方差稳定变换。本研究采用方差稳定变换。方差稳定变换可使图像中的噪声与信号近似不相关［4］，并且可通过相应的逆变换得到信号的无偏估计。

1.2 方差稳定变换

方差稳定变换［4，7］是为了使异方差数据转换为同方差的，更容易地运用标准方法来解决问题，去除Rician噪声方差和无噪声图像幅值的依赖性。

图1为所采用的方差稳定变换及其对应的逆变换，其中f为方差稳定变换，D为去噪后逆变换之前的图像，Vf为逆变换。对MR幅值图像进行方差稳定变换处理后，使得σf(θ)=std{f(M)|S，σ}=1，图像各处噪声的方差基本一致，这样可将VST变换后的Rician噪声作为加性噪声进行抑制处理。对方差稳定变换图像进行噪声抑制后，需经过对应的方差稳定逆变换，才能得到无偏的去噪图像。

图1 方差稳定变换及其逆变换。(a)方差稳定变换;(b)方差稳定逆变换Fig.1 VST and inverse VST.(a)VST;(b)Inverse VST

2 噪声水平场估计

对于水平空间变化的Rician噪声，依据局部的噪声水平逐像素进行VST，即实现了空间自适应VST。空间各处不同的噪声水平在空间上构成一个场，称之为噪声水平场(noise level field，NLF)。显然，进行空间自适应VST的关键在于NLF的准确估计。利用NLF进行空间自适应VST和去噪的流程如图2所示。主要包括4个步骤:估计NLF;利用NLF进行空间自适应 VST;对变换的图像使用BM3D算法进行去噪;对去噪后图像进行VST逆变换。

图2 空间自适应VST和去噪的流程图Fig.2 The flow chart of spatially adaptive VST and denosing

2.1 噪声水平场模型

采用多项式描述噪声水平场［8］

式中，{at，l}为多项式系数，xi，yi为像素点 i的坐标，K为多项式阶数。为了方便计算，{at，l}用列向量A表示，{}用列向量C表示，B表示估计的噪声水平场。式(2)可以扩展到三维的NLF。

式中，λ为设定的正则化系数，B'为噪声水平的局部估计值，其估计方法将在下一小节做详细介绍。式(3)为l1正则化约束优化问题，解决此类问题的优化方法很多［9］，采用 L1_LS 算法［10］求解多项式系数A，进而由式(2)得到NLF。L1_LS算法结合了截断牛顿内点法、预条件共轭梯度算法，比普通的内点法使用方向或共轭梯度法计算效率要高，求解时间要短。

2.2 局部Rician噪声水平的估计方法

由式(3)可见，为了估计参数A，需要得到噪声水平的局部估计。传统的噪声水平估计方法一般假设噪声水平是空间不变的常数。对于MR图像中的Rician噪声，较为常用的是基于背景区域的估计方法，但这些方法对于空间变化的噪声水平不再适用。假设图像局部的噪声水平为常数，采用修正的中值绝对偏差(median absolute deviation，MAD)估计方法估计图像局部的噪声水平［11］。

假设噪声服从高斯分布，对局部噪声水平进行MAD 估计［12－13］

式中，i，j表示图像中的像素点，Ni，Nj表示点 i，点 j的邻域，S(i)，S(j)是像素点的信号幅值，median(S)为中值滤波，为估计的噪声标准偏差。与其它噪声方差估计方法相比，MAD估计不但具有较强的鲁棒性，而且运算量较小。当信噪比足够大(＞5)时，Rician分布的标准偏差估计与高斯分布的趋于一致，无须进行修正;当信噪比较小时，式(4)的估计将会产生一定的偏差，须对估计的进行修正。这种修正方法是基于对Rician噪声信噪比(signal noise ratio，SNR)的迭代估计［14］。用 MAD估计出的标准偏差初始化修正过程为

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

式中，I0，I1分别是零阶、一阶修正的贝塞尔函数。修正因子ζ(θ)通过迭代估计直到收敛或者达到给定的迭代次数。使用|θt－θt－1|作为收敛条件，迭代过程可表示为

式中，〈S〉是所给像素信号幅值的平均值，t为迭代次数。在梯度大的像素点上估计的局部噪声水平与实际噪声水平相比，存在着较大偏差。因此，在式(3)估计噪声水平场时，不使用梯度大和修正后信噪比仍然较小的局部噪声水平估计。

2.3 空间自适应方差稳定变换去噪

利用估计的NLF对图像中各像素进行方差稳定变换

式中，β为可调节参数，B为估计的噪声水平场，i，j表示图像的位置。式(8)实际上相当于对MR图像进行同态化和归一化处理，使得整幅图像的噪声水平处处近似为1，噪声与MR幅值和空间位置近似独立。这样，可采用针对加性高斯白噪声的去噪方法进行处理，然后进行方差稳定逆变换，最终得到无偏的去噪图像。采用BM3D(BM4D)算法对方差稳定变换后的图像进行去噪。BM3D是当前公认对加性高斯白噪声去噪性能良好的算法，它包含了非局部去噪的思想，同时又用到了变换域滤波的方法。通过空间自适应VST，可有效利用BM3D的去噪能力，实现空间变化Rican噪声的抑制。

3 实验

为验证所提出方法的有效性，进行仿真MR图像实验和真实MR图像实验，并与VST－BM3D、VSTBM4D和ANLM算法进行比较。为方便起见，本方法简记为SAVST (spatially adaptivevariancestabilization transformations)。

3.1 仿真实验

仿真实验中，多项式的阶数K设定为5，正则化系数λ为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 仿真噪声水平场的估计结果。(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

仿真实验图像数据来自BrainWeb。对T1加权图像、T2加权图像，添加仿真 Rician噪声水平从1%到15%进行实验，仿真的NLF如图4(b)所示。采用PSNR评价去噪效果。

图4 不同方法噪声估计结果。(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%时，两组图像使用不同方法得到的去噪结果。其中，VST－BM3D使用噪声水平为常数的方差稳定变换去噪方法，图4第一行为T1加权图像的仿真实验结果，第二行为T2加权图像的仿真实验结果。由图4可见，本方法对于噪声水平在空间上变化的Rician噪声抑制效果较好，产生的模糊较小，并能保持原图像的细节特征和边缘信息。而VST－BM3D在信噪比较低的区域(如图4(e)脑室区域)，噪声未被有效抑制。对T1加权图像、T2加权图像采用本文方法(SAVST)和VST－BM3D方法进行去噪实验，两种去噪方法在不同噪声水平下，对应的去噪后PSNR值如图5所示。其中，横轴表示加噪水平，从1% ～15%，纵轴表示去噪后图像的PSNR值。由图可知，与VST－BM3D方法相比，本方法对应的PSNR约提高了2 dB。

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

本方法可应用在三维MR图像上，对同一噪声水平场的三维脑部 MR图像，采用 VST－BM4D、ANLM和本方法进行去噪实验。图6为在不同噪声水平下，几种方法去噪后PSNR的比较。从图5和图6可看出，本方法对噪声水平空间变化的二维和三维MR图像进行去噪，对应的PSNR明显高于VST－BM3D方法。ANLM结合了噪声水平的局部估计，其去噪性能也略高于VST－BM3D方法，但相对耗时。对于图像大小为181像素×217像素×181像素的脑部数据，本方法所需的运行时间约为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

由于真实噪声图像无对应的参考图像，为了评价去噪效果，采用Q度量［15］评价去噪图像的图像质量。当参考图像不能获得的时候，Q度量可较好地评价图像质量。Q度量值越大，图像质量越好。为了使去噪效果达到最优，利用Q度量调节式(8)中的参数β。图8显示了对两幅乳腺MR图像使用不同方法进行去噪的结果，每一行代表一幅含噪声乳腺MR图像的实验结果，相应的Q值如图8(f)所示，可看出当参数β在1.2左右时，Q值达到最大值，去噪效果趋于最优。本去噪方法对应的Q值高于VST－BM3D，而且对于细节信息和边沿结构有更好的保存能力。

4 讨论

本研究要解决的问题是抑制MR图像中空间变化的Rician噪声。提出的方法是通过估计Rician噪声水平场和方差稳定变换，将Rician噪声转化为常数水平的加性噪声，然后使用一般的去噪算法实现对噪声的抑制。本文选取的去噪算法为BM3D算法(因其性能良好)，其他针对加性高斯白噪声的去噪算法，如非局部均值算法、高斯混合尺度模型(scale mixtures of Gaussians，GSM)算法［16］也可以替代BM3D算法。

图8 两幅真实MR图像的处理结果(上行为一幅，下行为另一幅)。(a)噪声图像;(b)VST-BM3D的去噪效果;(c)(a)与(b)之间的残差图像;(d)本方法的去噪效果，(e)(a)与(d)之间的残差图像;(f)调节β和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分布不服从高斯分布，需要对估计的噪声标准方差进行修正，才能得到无偏的噪声标准方差，从而提高了噪声水平场的估计精度。估计噪声水平场时，相关参数的设置对估计精度有较大影响，当多项式阶数K较小时，会得到偏低的估计，当K较大时，会得到偏高的估计;而正则化系数λ大小决定了噪声水平场的平滑程度，需要依据经验或先验知识进行调节，对于实际问题，可通过优化Q度量确定最优的λ。

5 结论

针对MR图像中噪声水平空间变化的Rician噪声，本研究提出了一种噪声水平场的估计方法，并用于抑制MR图像中空间变化的Rician噪声。通过Rician噪声水平的局部估计和稀疏性约束模型估计噪声水平场，然后对图像各处的幅值依据不同的局部噪声水平进行方差稳定变换，使得噪声与信号幅值和空间位置无关，然后利用BM3D算法的强大去噪能力，实现了对空间变化Rician噪声的有效抑制。实验结果表明，本方法能有效估计MR图像中的Rician噪声水平场，与使用全局噪声水平估计的方差稳定变换去噪方法、自适应非局部均值去噪算法相比，本方法对应的峰值信噪比和 Q度量均较高。

［1］陈武凡.并行磁共振成像的回顾、现状与发展前景［J］.中国生物医学工程学报，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，Foi 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，Farrher 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.