基于高阶奇异值分解和Rician 噪声校正模型的扩散加权图像去噪算法

徐 朴，郭 莉，冯衍秋，张鑫媛

南方医科大学生物医学工程学院//广东省医学图像处理重点实验室//广东省医学成像与诊断技术工程实验室//粤港澳大湾区脑科学与类脑研究中心，广东 广州510515

扩散磁共振成像（dMRI）技术可通过检测活体组织内水分子的扩散运动状态来反映生物组织的微观结构变化［1］。扩散张量成像（DTI）技术可通过采集不同扩散编码方向的扩散加权（DW）图像进行张量成像，用于描述水分子扩散运动的三维形式［2］。临床上常用的DTI量化指标包括各向异性分数（FA）和平均扩散率（MD），这些参数可用于检测急性脑中风［3］、多发性硬化症［4］、癫痫［5］以及脑肿瘤［6］等脑疾病的细微结构改变。

与常规MR图像相比，DW图像信噪比较低，受噪声影响严重。为提高DW图像的质量，已有多种后处理图像去噪算法被相继提出［7-15］。然而磁共振噪声信号服从Rician［16,17］分布，常用的基于变换域的图像去噪方法主要是针对加性高斯噪声，如基于块匹配的三维滤波（BM3D）［18］等，因此不能直接应用于磁共振噪声图像上，目前常用的策略是先对图像进行方差稳定性变换（VST）［7,19］，得到噪声空间分布不变的图像，然后再利用基于变换域的去噪方法对该图像去噪，去噪性能依赖于VST算法的效果。针对扩散加权图像去噪，有研究提出了一种基于全局指导下的局部高阶奇异值分解（GLHOSVD）去噪方法［9］，其去噪性能处于国际领先水平，且能够降低条形伪影。然而，GL-HOSVD方法是对相似块组成的高维数组进行高阶奇异值分解（HOSVD）变换，因此在去噪同时不可避免的会引入条形伪影。此外，GL-HOSVD方法只适用于处理加性高斯噪声图像，因此也需在去噪之前对DW噪声图像进行VST变换。

本文中我们提出了一种联合HOSVD稀疏约束与Rician 噪声校正模型的去噪算法，该方法不需要借助VST变换，可直接处理DW噪声图像。此外，该算法不需要构造相似块组进行去噪，因此从根本上解决了条形伪影的问题。为验证本方法的有效性，我们分别进行了仿真数据和真实数据实验研究，将低秩+边缘约束（LR+Edge）［10］，GL-HOSVD［9］，BM3D［18］和非局部均值（NLM）［20,21］4种去噪算法作为实验对比方法，结果表明本方法无论是从定量还是定性方面均优于4种实验对比方法。

1 方法

1.1 高阶奇异值分解去噪原理

高阶奇异值分解是一种张量分解方法，可看成是对二维矩阵奇异值分解的一种推广，可处理任意阶数（N≥3）的张量［22］。已知三阶张量AI1×I2×I3的HOSVD分解可表示为：

已知U(1)，U(2)，U(3)，核心张量S可由以下公式得到：

假设一幅无噪图像具有结构稀疏性，即经过HOSVD分解后，可用较少的变换系数来表示。若A是受加性高斯噪声污染的张量数据，那么HOSVD核心张量S中较小的系数主要来自于噪声。因此可直接对S的系数进行阈值操作得到以达到去噪效果［23-25］。去噪后的张量数据可表示为

1.2 Rician噪声校正模型

磁共振成像中，采集到的原始数据是K空间数据，经傅里叶变换得到空间域的复数图像，再经模运算得到幅值图像，其噪声不再是简单的加性高斯噪声，噪声信号服从Rician分布，其概率密度函数表示为［26］：

其中x为无噪信号，y为噪声信号，σ为噪声信号在复数域中高斯噪声的标准差，I0为第1类0阶修正贝塞尔函数。在高信噪比区域，py(y|x,σ)近似均值为x，标准差为σ的高斯噪声，在低信噪比区域，py(y|x,σ)会产生偏移，噪声信号y的均值也会随之偏移理想信号x，且幅值图像的噪声波动大小与图像信号有关，是空间变化的。因此基于变换域的去噪方法不能直接用于磁共振图像上。目前现有方法是先将MR图像进行VST变换，使得MR幅值图像的噪声分布是空间不变的，然后再对其进行去噪。而本文利用了Rician噪声信号的期望，进行噪声偏差校正。Rician噪声信号y的期望可表示为［26］：

1.3 基于HOSVD稀疏约束和Rician噪声校正模型的去噪算法

本文采用Rician噪声信号期望进行噪声偏差校正，将其融合到HOSVD 去噪模型当中，从而不需要借助VST变换，可直接处理DW噪声图像。具体来说，本文所提去噪模型为：

其中Y为带有Rician噪声的DW图像，数据大小为M×N×L，其中M×N为二维DW图像大小，L为不同编码方向的DW图像的总数，X为待求解的无噪图像，νNLM（Y）为采用向量形式的非局部均值对Y进行滤波得到的预处理图像。

我们采用变量分离法求解以上去噪模型（公式（8）），即将去噪模型分解为以下两个子问题，采用迭代的形式进行交替求解，最终得到最优解，即最终的去噪图像

对于子问题1，假设S,U(1),U(2),U(3)均已知，通过l-BFGS优化算法得到X。对于子问题2，固定X，对X进行HOSVD变换和稀疏约束，得到U(1),U(2),U(3)。

其中公式（12）为硬阈值操作。本文采用经典的自适应系数收缩函数定义阈值

γ为噪声估计水平的预设参数，

本文所提去噪算法的具体实施步骤如下：

1.根据图像背景计算DW图像的高斯噪声σ，公式为Ybg为图像背景的灰度值；

2.初始化X(1)，公式为

3.针对X(k)，采用滑动窗形式，以步长为Nstep，选取大小为t×t×L的DW三维块Pi，对每个局部图像块Pi，进行HOSVD 分解（公式（11））以及硬阈值操作（公式（12）），得到

5.根据Z(k)计算vNLM的相似性权重，根据更新后的权重重新计算得到vNLM(Y)；

6.已知Z(k)或是S，U(1)，U(2)，U(3)，根据公式（9）得到X（k+1）；

7.若满足迭代停止条件，则退出，X（k+1）为最终的去噪图像，否则跳到步骤3，继续进行迭代。

1.4 实验数据

为了验证所提方法的有效性，我们分别进行了仿真数据和真实数据实验研究。

仿真数据实验中的DW 图像是由张量模型S0exp（-bgT Dg）产生的，S0为非扩散加权图像，g为单位扩散编码方向，D为扩散张量场，来自于生物医学信息学研究网络数据库中的一个成年大鼠数据，该数据是由一个b=0图像以及44个不同扩散编码方向b=2000 s/mm2的DW图像构成，图像大小为256×256。然后加入不同噪声水平的Rician噪声，公式如下：

其中X表示无噪的DW图像，N1，N2为标准差为σ的高斯噪声。经过平方-相加-开方运算，幅值噪声图像Y服从Rician分布。为了验证所提方法在不同噪声水平下是否有效，我们将仿真数据最大值归一化到1，然后仿真不同噪声水平的DW图像（σ=0.01，0.03，0.05，0.07和0.09）。

真实数据的采集设备是飞利浦3.0T TX MR扫描仪（Achieva,Philips Medical Systems,Best,The Netherlands），采用8通道脑部接收线圈，采集序列为4次激发的EPI序列。DW图像通过multiplexed-SENSE重建得到，没有进行加速采集。图像成像参数为：6个不同扩散编码方向b=1000 s/mm2的DW 图像，一个b=0 图像，TR/TE=3000/82 ms,层内分辨率=1 mm×1 mm，层厚=4 mm，矩阵大小=220×220，层数=14，为了得到高信噪比的参考图像，重复采集10次（NEX=10）。

1.5 定量指标

为了定量评价所提算法的有效性，本文采用了DW图像的峰值信噪比（PSNR），结构相似性测度（SSIM）和FA 图的均方根误差（FA-RMSE）分别对去噪后的DW图像和后续的量化参数FA 图进行定量分析，计算公式如下：

其中X为无噪DW图像为去噪DW图像，分别为图像的均值，分别为图像的标准差，c1，c2为常数，分别为无噪DW图像和去噪后的DW 图像得到的FA 参数图。M1和M2为参与计算PSNR和各向异性分数均方根误差（FA-RMSE）的所有像素值的个数。我们只采用前景区域的像素进行量化分析，从而排除背景干扰。PSNR值越大，SSIM值越接近于1说明去噪图像的质量越好，RMSE越小说明由去噪图像得到的量化参数FA越准确。

1.6 参数设置

与BM3D［18］，NLM等［20,21］去噪算法类似，本算法的去噪性能依赖于参数设置，主要包括以下5个参数：图像局部块大小t×t，步长Nstep，噪声标准差重估计尺度γ，基于Rician 噪声校正模型保真项的权重λ，以及参与vNLM计算的平滑尺度β。通过优化参数我们发现，当t固定不变时，步长Nstep越小，去噪效果越好，但计算时间越长，通过权衡计算时间和去噪效果，我们固定，t=60，Nstep=10。通过仿真实验，我们发现随着噪声水平变大，最优值λ也会随之变大。对于真实数据实验，我们通过目测去噪图像质量来调节β和λ。实验参数设置如表1所示。

表1 本文中仿真数据和真实数据实验中的具体参数取值Tab.1 Parameter values in simulation data and real data experiment

2 结果

2.1 仿真数据实验结果

图1为所提方法与4种去噪方法（NLM，BM3D，LR+Edge，GL-HOSVD）在不同噪声水平下的量化比较结果。从图1可知，无论是从去噪图像的PSNR，SSIM还是后续参数FA图的RMSE来看，所提方法和GL-HOSVD算法的去噪性能均显著优于NLM，BM3D和LR+Edge。此外，在较低的噪声水平下（1%和3%），所提方法和GL-HOSVD 的量化结果相似，进一步观察，我们会发现所提方法的PSNR 和SSIM 略高于GL-HOSVD，而FA-RMSE 要略低于GL-HOSVD。在噪声水平为7%和9%时，所提方法的PSNR 和SSIM 要明显高于GLHOSVD，在噪声水平为5%和7%时，所提方法的FARMSE要明显低于GL-HOSVD。图2显示了在所提出的算法中不使用vNLM与使用vNLM预处理的结果对比。结果显示，在所提去噪框架中引入vNLM滤波能够进一步提高去噪图像的PSNR以及降低后续量化参数FA图的RMSE，具有更高的去噪性能。图3给出了噪声水平为5%时，不同去噪算法得到的去噪图像、FA图以及带有彩色方向编码信息的FA图以及相应的误差图。从图3中可以看出，5种去噪方法都能有效降低图像噪声，恢复图像主要的结构信息，且5种方法都能够显著提高DW图像的PSNR以及降低后续量化参数FA图的误差值。此外，本文方法和GL-HOSVD的去噪效果要明显优于其它3种去噪方法，去噪同时能够更好的保留图像的细节信息，DW图像和FA图的边缘结构信息更加清晰。虽然从整幅图来看，本文方法和GL-HOSVD得到的去噪图像以及FA图的视觉效果相似，但从局部放大图（图4）来看，GL-HOSVD与BM3D得到的去噪图像中含有一些条形伪影（如箭头所示），并且这些条形伪影结构也会出现在后续的量化参数图中（如箭头所示），而本文所提方法不但能够得到高质量的去噪图像以及准确的FA图，并且不存在这种伪影结构。

图1 所提方法和4种去噪方法的量化对比（NLM，BM3D，LR+Edge和GL-HOSVD）Fig.1 Comparison of the proposed method with 4 denoising methods(NLM,BM3D,LR+Edge and GL-HOSVD).

图2 所提方法中不使用vNLM和使用vNLM预处理的结果对比Fig.2 Comparison of the results using the proposed method with and without vNLM preprocessing.

图3 不同去噪方法的仿真实验结果对比Fig.3 Comparison of different denoising algorithms on simulation data.

图4 仿真实验去噪结果的局部放大图Fig.4 Enlarged views of the denoised images for simulation data.

2.2 真实数据实验结果

图5为真实数据的去噪方法比较结果。4种去噪算法都能有效降低图像噪声，恢复图像结构信息。但其中NLM的去噪结果相比于NEX=10过平滑，模糊了部分细节，在局部放大图（图6）中，NLM存在一些条形伪影。BM3D保留了较多的边缘结构信息，但同样也在部分位置引入了条形伪影（图6箭头所示）。GL-HOSVD和本文所提方法的去噪图像与NEX=10的参考图像相似，甚至比NEX=10 的参考图像噪声还要小，FA 图也是如此。但是从图像的局部放大图来看（图6），GLHOSVD在去噪同时也引入了一些伪影（如箭头所示），但是本文方法得到的去噪图像以及FA图均没有伪影出现，该实验结果与仿真实验结果一致。

图5 不同去噪方法对真实数据进行处理的结果比较Fig.5 Comparison of different denoising algorithms on real data.

图6 真实数据去噪结果的局部放大图Fig.6 Enlarged view of the denoised images for real data.arrow:stripe artifact.

3 讨论

dMRI技术可通过采集不同编码方向的扩散加权信号，来探测组织内水分子扩散的运动状态，从而反应组织的微观结构变化，以及进行脑白质神经纤维束的构建［29］。如各向异性分数和表观扩散系数能够探测多种脑疾病的病理组织改变。通过dMRI技术得到的大脑神经纤维束能够指导临床中脑肿瘤手术方案的制订，以及进行术后的疗效评估［30］。但与常规T1和T2图像相比，DW图像受噪声影响严重。虽然多次采集平均方法可提高DW图像的信噪比，但是随着采集次数的增加，采集时间也会增加，从而增加了采集过程中运动伪影的可能性。

目前很多研究学者提出采用后处理去噪技术提高图像的信噪比。由于磁共振噪声服从Rician分布，当信噪比较高时，Rician噪声信号近似高斯噪声，期望近似等于理想信号，而当信噪比较低时，Rician噪声信号近似瑞利分布，期望有偏于理想信号值。NLM作为图像域的去噪算法通常需要联合Rician 噪声校正模型，NLM自适应差，去噪时容易出现过平滑的问题［8］。LR+Edge可以将Rician噪声融合到最大后验中，因此也可以解决Rician噪声校正的问题，然而LR+Edge需要足够鲁棒的噪声估计，较差的噪声估计将导致比较差的去噪效果［10］。MR图像在实部和虚部信号噪声服从高斯分布，在相位校正过程中针对复数域的实部和虚部去噪，最终可获得去噪后的幅值图像，然而在临床中得到的数据通常是重建后的幅值图像而不是实部和虚部图像，因此该方法在广泛应用于MR去噪方面有一定的局限性［11］。目前现有的基于变换域的图像去噪方法不能直接处理磁共振噪声图像，尤其是对于这种信噪比较低的DW图像。据本文所知，现有基于变换域的去噪方法均是先对MR 图像进行噪声方差处理，然后再进行图像去噪。GL-HOSVD 也是一种基于变换域的去噪方法［9］，GLHOSVD借鉴了BM3D［18］的思想，将结构相似的块放在一起，构成高维数据，然后进行HOSVD变换，因此去噪效果极佳，在扩散加权图像去噪方面处于国际领先水平。但是GL-HOSVD也需要结合VST变换使用，去噪效果高度依赖VST算法的性能。此外，通过寻找相似块构造结构更加稀疏的高维数组这类去噪方法，虽然能够得到较好的去噪效果，但是在噪声水平较大的时候，会引入伪影结构，尤其是对这种信噪比较低的DW 图像。虽然GL-HOSVD通过全局预滤波来指导后续的基于块匹配的HOSVD去噪，可以有效减少去噪过程中引入的伪影，但是并没有从理论上解决此类问题。

为了解决以上问题，我们提出一种新颖的HOSVD去噪方法，建立了基于HOSVD稀疏约束和Rician噪声校正模型的联合去噪方法，所提去噪模型可以直接处理Rician噪声，不需要借助VST算法。此外，本文假设，在噪声水平较大时，在图像均匀区域或是灰度变化比较缓慢的地方，噪声产生的灰度变化起主导作用，通过寻找相似块来构造高维数组，会将具有相似噪声模式的图像块聚在一起，在去噪过程中错把模式相似的噪声当作图像结构保留下来，形成条形伪影。此外，这些相似块存在高度重叠区域，重叠区域的噪声值相同，导致某些噪声值出现的概率变大，因此构成的高维数组中的噪声并不是严格的随机噪声，这也是导致伪影问题的一个可能因素。因此，本文所提方法对整幅DW图像或是每个局部DW图像单独进行HOSVD去噪，同一结构块中，不会重复出现模式相同的噪声，因此从根本上解决了伪影问题。此外，我们在Rician噪声期望保真项中引入了vNLM滤波，从而能够降低原始信号的噪声波动，提高模型拟合精度。本方法也可看成是一种双域的去噪模型，结合了空间域的非局部均值滤波和变换域的HOSVD去噪，这种双域的去噪模型能够优势互补，从而达到更好的去噪效果。

我们进行了仿真实验和真实实验验证，并且与4种典型的DW图像去噪算法进行了比较，其中GL-HOSVD算法的去噪性能处于国际领先水平。实验结果表明本方法能够有效降低图像噪声，保留图像细节结构，并且不会引入条形伪影，其去噪性能无论在定量还是定性方面都是最优的。与其他优化算法类似，本方法的高水平去噪性能依赖于参数的设置。由于PSNR/SSIM最高的DW图像产生的FA图的RMSE并不是最小的，因此我们同时考虑了DW 图像的PSNR，SSIM 和FA 图的RMSE（FA-RMSE）3个指标对所提算法进行参数优化。

本方法是以二维形式进行去噪的，可扩展到三维形式用以处理各向同性的DW数据。本方法同样适用于多通道采集得到的noncentral-χ分布的DW噪声图像，只需要将Rician噪声校正模型换成noncentral-χ分布下的噪声信号期望即可。目前，该方法只适用于空间不变的噪声水平，接下来我们可以研究如何对本方法进行改进，从而解决噪声水平随空间变化的DW图像。

综上所述，本文提出一种新颖的HOSVD 去噪方法，建立了基于HOSVD稀疏约束和Rician噪声校正的去噪模型。所提去噪方法可以直接处理带有Rician噪声的DW图像，不需要借助VST变换。此外，本方法还从根本上解决了去噪过程中引入伪影的问题。实验结果证明了该方法的有效性，其去噪性能比采用了VST变换以及结构相似块的GL-HOSVD方法更好。