扩散加权图像(DWI)能够通过计算得到扩散张量图像。扩散张量成像技术(DTI)是目前唯一非侵入活体获取脑白质结构的技术，从扩散张量图像中能够尽早发现急性脑梗塞症状，也能够揭示有关大脑神经的一些问题，且广泛的应用于人体大脑的神经纤维跟踪，该技术在医学上具有十分重要的意义。但在DWI图像的获取和传输的过程中不可避免的会受到一定程度的噪声干扰，从而影响到DTI数据的准确度以及后续的处理。因此，计算DTI数据之前对DWI图像进行降噪是非常重要的。

目前国内外关于DWI图像降噪的研究有许多。Zhang等人^[1]提出的使用高阶奇异值分解对DWI进行降噪，显著减少了在处理图像的过程中产生的伪条纹。Liu等人^[2]提出一种用于同时平滑和估计DTI图像的鲁棒变分法，该算法直接对DTI降噪，实现有效降噪的同时减少了计算时间，但是该算法不是直接处理DWI图像中的噪声，因而DTI数据存在可靠性问题。Haldar等人^[3]利用数据集中所观察到的高水平结构相似性，采用统计重构的方法，在提高信噪比的同时保持了DWI图像的分辨率。Bao等人^[4]通过定义了一个基于局部均值的相似性度量策略，提出将结构自适应稀疏的方法应用于DWI图像降噪，该算法能够较好的减少DWI图像中的噪声。Lam等人^[5-7]将降噪问题看成一个包含边缘和低秩模型的最大后验估计问题，通过秩和边缘的约束有效提高了DWI图像的信噪比，从而提高扩散参数估计精度且提高了执行效率。Becker等人^[8-9]提出在体素空间和扩散梯度空间中基于结构自适应平滑的DWI图像数据增强方法。McGraw等人提出通过加权的全变分(TV)范数最小化对扩散加权图像数据进行平滑，该算法较好地实现了对DWI图像的平滑^[10]。张相芬^[11]提出的基于各向异性扩散的DTI图像恢复，主要采用的是经典的各向异性(PM)算法，能够较好地保存图像的边缘。

以上算法均能够有效地改善DWI图像的质量，但是这些算法都没关注到DWI图像中的纹理细节信息。由于DWI图像的自相似性程度高，纹理和结构具有重复出现的特性并且纹理细节较多。而加权核范数降噪算法能够利用图像自相似性的特点，在整幅图像中寻找与局部块相似的图像块，将其堆叠在一起形成相似块矩阵并对相似块矩阵进行奇异值分解，根据奇异值的不同分配权值，再将相似块矩阵聚合得到估计原始图像从而达到降噪的目的。由于该算法是对整幅图像中每个局部块的相似块进行处理，因而可以更好地保留图像中的纹理细节。基于上述特点，本文将加权核范数降噪算法应用到DWI图像的降噪中，对DWI图像进行降噪的同时也能够更好地保存图像中的纹理细节。

扩散张量成像技术通过检测水分子在大脑中的扩散情况来获取其各向异性等信息，而人体大脑中水分子的扩散运动可以通过扩散张量来表示，其扩散状态可以间接的反映组织微观结构特点及其变化。DTI图像中的每一个体素都是由扩散张量$\mathit{\boldsymbol{D}}$组成的，扩散张量是一个对称的正定矩阵，一般用$\mathit{\boldsymbol{D}}$表示为

$ \mathit{\boldsymbol{D = }}\left[{\begin{array}{*{20}{c}} {{D_{xx}}}&{{D_{xy}}}&{{D_{xz}}} \\ {{D_{xy}}}&{{D_{yy}}}&{{D_{yz}}} \\ {{D_{xz}}}&{{D_{yz}}}&{{D_{zz}}} \end{array}} \right] $

(1)

扩散张量数据可以用于处理大脑中的许多问题如对大脑神经纤维进行跟踪等，在医学领域具有重要作用。DTI图像是根据Stejskal-Tanner公式对DWI图像进行计算得到的^[12-13]，而DWI图像是通过对普通核磁图像施加扩散敏感梯度脉冲得到的，它能够检测机体组织中水分子的扩散状态。由于扩散张量$\mathit{\boldsymbol{D}}$是对称矩阵，因而扩散张量$\mathit{\boldsymbol{D}}$实际只有6个独立分量，这6个独立分量表示的是水分子在不同方向上的扩散系数。所以令$\mathit{\boldsymbol{D = }}{\left[{\begin{array}{*{20}{l}} {{D_{xx}}}&{{D_{yy}}}&{{D_{zz}}}&{{D_{xy}}}&{{D_{xz}}}&{{D_{yz}}} \end{array}} \right]^{\text{T}}}$，具体计算公式为

$ \left\{ \begin{gathered} AD{C_i} =-\frac{1}{{{b_i}}}\ln \frac{{{S_i}}}{{{S_0}}} \hfill \\ \mathit{\boldsymbol{G}}*\mathit{\boldsymbol{D}} = ADC \hfill \\ \end{gathered} \right. $

(2)

式中，$\mathit{\boldsymbol{D}}$为扩散张量，$i$表示方向数且$i \geqslant 6$，$ADC$被称为表观扩散系数，${b_i}$为第$i$个方向的扩散敏感系数，用来表示扩散的加权程度。${S_0}$表示未施加敏感梯度脉冲的核磁图像，${S_i}$为施加了第$i$个方向的扩散敏感梯度脉冲的DWI图像。由于扩散张量$\mathit{\boldsymbol{D}}$有6个独立分量，所以要计算扩散张量就需要至少6个方程，也就是至少需要6个不同的方向。$\mathit{\boldsymbol{G}}$为扩散敏感梯度记为：$\mathit{\boldsymbol{G}} = \left[{\begin{array}{*{20}{l}} {x_i^2}&{y_i^2}&{z_i^2}&{2{x_i}{y_i}}&{2{x_i}{z_i}}&{2{y_i}{z_i}} \end{array}} \right]$，不同方向的扩散敏感梯度不同。

根据式(2)可知DTI图像的计算源于基准图像${S_0}$和加权的核磁图像${S_i}$。而这些图像在获取过程中均会受到噪声污染，且由于${S_i}$是施加了扩散敏感梯度的DWI图像，因此所含噪声更多图像更模糊，从而使得根据式(1)计算得到的DTI数据不可靠。并且由于式(2)是关于方程组的计算，这个方程组对噪声十分敏感^[14]，即使微小的偏差也会对扩散张量数据产生较大的影响。因此，对DWI进行有效降噪和对DTI的计算及其后续研究都具有非常重要的意义。

加权核范数降噪算法也可以称为加权核范数最小化算法(WNNM)是由Gu等人^[15]提出的，该算法利用了图像的自相似性的特点对图像进行降噪。而图像中普遍存在自相似性，即图像的纹理和结构具有重复出现的特性，且人的大脑图像在这一方面尤为突出。这种特性体现在图像中即为数据具有冗余性，这在图像处理上具有十分重要的意义。加权核范数降噪算法通过在整幅图像中寻找与图像局部块相似的块，将其堆叠形成相似块矩阵，对相似块矩阵进行处理来估计原始图像以达到降噪目的。由于图像块不仅能表示图像的灰度分布特性，还能表示图像的几何特性，因而用图像块的方式对图像进行处理能够较好的保存图像中的纹理细节信息，降噪效果较理想。

对于一幅含噪图像，其模型为

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{x}} + \mathit{\boldsymbol{n}} $

(3)

式中，$\mathit{\boldsymbol{y}}$为观察的含噪图像，$\mathit{\boldsymbol{x}}$为无噪的原始图像，$\mathit{\boldsymbol{n}}$是均值为0方差为${\sigma ^2}$的高斯噪声。将含噪图像$\mathit{\boldsymbol{y}}$分为若干个块，每个块称为局部块${y_j}$，对局部块使用块匹配的方法从整幅图像中找到与它相似的图像块，并将找到的所有相似块进行堆叠，形成一个相似块矩阵${\mathit{\boldsymbol{Y}}_j}$，则

$ {\mathit{\boldsymbol{Y}}_j} = {\mathit{\boldsymbol{X}}_j} + {\mathit{\boldsymbol{N}}_j} $

(4)

式中，${\mathit{\boldsymbol{X}}_j}$和${\mathit{\boldsymbol{N}}_j}$分别是原始图像的相似块矩阵与噪声图像块矩阵，$j$表示第$j$个局部块。由于图像本身的纹理和结构具有重复出现的特性，因此${\mathit{\boldsymbol{X}}_j}$可看做是一个低秩矩阵，但由于图像受到噪声${\mathit{\boldsymbol{N}}_j}$的影响，含噪图像的秩比原始图像的秩要高，因此通过低秩矩阵逼近的方法来得到原始图像。

文献[16]中证明了大多数的低秩矩阵都可以用核范数最小化(NNM)算法得到有效解决^[16]。NNM算法是通过对相似块矩阵的奇异值进行等量收缩来达到降噪的目的，但由于不同的奇异值具有不同的重要性，奇异值的大小表示所含图像信息的多少，而对奇异值等量收缩会导致许多图像的重要信息也被滤除，因此NNM算法在处理问题的时候不够灵活。于是使用加权的NNM算法即WNNM算法更为合理，通过对不同的奇异值分配不同的权值从而实现对图像的降噪。因此，式(4)可以使用WNNM算法从图像块${\mathit{\boldsymbol{Y}}_j}$中估计图像块${\mathit{\boldsymbol{X}}_j}$，再通过聚合所有降噪后的图像块，就可估计整个图像$ x$。通过使用噪声方差$\sigma _n^2$来标准化F范数的数据保真项$\left\| {{\mathit{\boldsymbol{Y}}_j}-{\mathit{\boldsymbol{X}}_j}} \right\|_{\text{F}}^2$，于是有能量函数

$ {\mathit{\boldsymbol{\hat X}}_j} = \arg \;\mathop {\min }\limits_{{\mathit{\boldsymbol{X}}_j}} \frac{1}{{\sigma _n^2}}\left\| {{\mathit{\boldsymbol{Y}}_j}-{\mathit{\boldsymbol{X}}_j}} \right\|_{\text{F}}^2 + {\left\| {{\mathit{\boldsymbol{X}}_j}} \right\|_{\omega, *}} $

(5)

式中，${\mathit{\boldsymbol{\hat X}}_j}$为估计得到的降噪图像块矩阵，${\left\| {{\mathit{\boldsymbol{X}}_j}} \right\|_{\omega, *}}$为加权的核范数矩阵，可表示为

$ {\left\| {{\mathit{\boldsymbol{X}}_j}} \right\|_{\omega, *}} = \sum\limits_i {{{\left| {{\omega _i}{\sigma _i}\left( \mathit{\boldsymbol{X}} \right)} \right|}_1}} $

(6)

式中，${\sigma _i}\left( {{\mathit{\boldsymbol{X}}_j}} \right)$表示为矩阵${\mathit{\boldsymbol{X}}_j}$的第$i$个奇异值。${\omega _i}$为权值函数，表示为第$i$个奇异值的权值，具体公式为

$ {\omega _i} = c\sqrt n /\left( {{\sigma _i}\left( {{\mathit{\boldsymbol{X}}_j}} \right) + \varepsilon } \right) $

(7)

式中, $n$为相似块的个数，$c$为正常数，$\varepsilon $是避免分母为0的常数。由于原始图像中相似块矩阵${\mathit{\boldsymbol{X}}_j}$的奇异值是未知的，因而要从含噪相似块矩阵${\mathit{\boldsymbol{Y}}_j}$的奇异值估计得到，计算公式为

$ {\hat \sigma _i}\left( {{\mathit{\boldsymbol{X}}_j}} \right) = \sqrt {\max \left( {\sigma _i^2\left( {{\mathit{\boldsymbol{Y}}_j}} \right)-n\sigma _n^2, 0} \right)} $

(8)

式中，$\sigma _n^2$为噪声方差，$n $为相似块个数，$\sigma _i^2\left( {{\mathit{\boldsymbol{Y}}_j}} \right)$表示图像块${\mathit{\boldsymbol{Y}}_j}$的第$ i$个奇异值。通过得到的权值就可以对原始图像块矩阵${\mathit{\boldsymbol{\hat X}}_j}$进行估计。

式(5)中的解${\mathit{\boldsymbol{\hat X}}_j}$可由${\mathit{\boldsymbol{\hat X}}_j} = \mathit{\boldsymbol{U}}{S_\omega }\left( \mathit{\boldsymbol{ \boldsymbol{\varSigma} }} \right){\mathit{\boldsymbol{V}}^{\text{T}}}$得到，$\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}$为奇异值矩阵，${S_\omega }\left( \mathit{\boldsymbol{ \boldsymbol{\varSigma} }} \right)$为软阈值函数，而通过对含噪图像的相似块矩阵${\mathit{\boldsymbol{Y}}_j}$进行奇异值分解即$\left[{\mathit{\boldsymbol{U}}, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}, \mathit{\boldsymbol{V}}} \right] = SVD\left( {\mathit{\boldsymbol{Y}}j} \right)$，可得到矩阵$\mathit{\boldsymbol{U}}$和矩阵$\mathit{\boldsymbol{V}}$，软阈值函数的作用是对奇异值进行收缩，其计算公式为

$ {S_\omega }{\left( \mathit{\boldsymbol{ \boldsymbol{\varSigma} }} \right)_{ii}} = \max \left( {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{ii}}-{\omega _i}, 0} \right) $

(9)

式中，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{ii}}$为奇异值矩阵$\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}$中的每个对角元素，$\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}$为第$i $个奇异值的权值。由式(7)和式(9)可知，因为相似块矩阵的奇异值越大所含图像的信息就越多，所以分配的权值就越小，收缩速度就减小，从而达到保留信息的作用；反之，奇异值越小所含信息越少，则分配的权值就越大，收缩的速度也就越大，滤除噪声的强度也更大。通过使用收缩后的奇异值矩阵来估计新的相似块矩阵${\mathit{\boldsymbol{\hat X}}_j}$，再将所有估计得到的相似块矩阵进行聚合处理，最终得到降噪图像。

WNNM算法利用了图像本身非局部块之间的相似性，通过对含噪图像中的每个块进行块匹配处理得到相似块矩阵，再根据奇异值的特点，对相似块矩阵的奇异值进行收缩处理得到降噪块估计，该算法能够更好的保持图像局部结构，也就是能够保留图像中较多的纹理细节信息同时也能够较大程度的减少噪声。而DWI图像中含有较多的纹理细节，且其自相似性程度高，噪声对DWI图像产生的破坏则会对DTI数据的准确性造成一定的影响。因此，本文根据DWI图像的特点提出将WNNM算法应用于DWI图像的降噪，该算法能够尽可能多的保存DWI图像中的纹理细节信息，且能够较好的减少噪声，具有较好的降噪效果，提高了数据的准确性。

本文通过模拟数据实验和真实数据实验，将WNNM算法与常用的DWI降噪算法如TV算法、各向同性算法、PM算法、Wiener算法^[17]以及纹理检测算法^[18]进行比较。

模拟数据实验通过建立模拟张量场对DWI图像进行加噪，再使用各种降噪算法对DWI图像进行降噪，观察其降噪效果。实验使用扩散加权方向为六个扩散方向，其方向编码如表 1所示。

所用的模拟数据是由7幅大小为128×128像素的DWI图像组成，其中包括一幅未加权的图像${S_0}$以及6幅根据表 1所示的扩散加权方向进行加权所得到的DWI图像${S_1}\sim{S_6}$。采用椭球模型来建立模拟张量场，但由于图像过大张量不便显示，于是选择了相同位置大小为15×15的张量视图。通过这7幅模拟数据计算得到的原始不含噪的张量场如图 1(a)所示。对这7幅图像${S_0}\sim{S_6}$加入均值为0标准差为$\sigma = 40$的高斯噪声，根据含噪的DWI图像计算得到含噪张量场如图 1(b)所示。然后再分别采用PM，TV，各向同性，纹理检测，wiener以及WNNM算法对DWI图像进行降噪，降噪后计算得到的张量场如图 1(c)-(h)所示。

从图 1(a)(b)可以看出，由于噪声的存在对张量的计算产生了较大的影响，使张量的大小形状及其主方向都发生了较大的改变。从图 1(c)-(h)中可以看出，使用WNNM算法降噪后计算得到的张量图与原始无噪张量图几乎趋于一致，其大小形状以及主方向基本未发生变化，降噪效果最好。通过对比其他算法可知，使用PM算法和各向同性扩散算法降噪后计算得到的张量图中，张量的大小形状与原始张量图差异较大。使用TV算法降噪后得到的部分张量主方向发生了改变。使用纹理检测算法降噪后得到的张量主方向基本没变，但其大小与原始张量仍存在少许差异。

为了客观的说明上述算法的降噪效果，本文采用了峰值信噪比(PSNR)，结构相似性(SSIM)和均方误差(MSE)作为评价标准。峰值信噪比表示的是最大可能信号和噪声之间的比值，该值越高，表明降噪图像所含的噪声越少且与原始图像差异越小；而结构相似性的值在0与1之间，当降噪图像与原始图像的结构细节相似度越高时，它的值就越接近1。在不同的噪声水平下，使用各类算法对模拟数据进行降噪后得到的降噪图像评价结果如表 2所示。

从表 2中可以看出，在同一噪声水平下使用WNNM算法降噪后得到的降噪图像PSNR值和SSIM值都最高，且高出其他降噪算法许多，MSE值也较小。并且，随着噪声水平的增加，噪声对图像的影响加大，而WNNM算法的各类评价值结果依然保持较高的水平，说明该算法的抗噪性能优于其他算法。

综合上述实验可知，通过使用WNNM算法对DWI图像进行降噪得到的客观评价结果比其他常用的降噪算法效果要好，说明WNNM算法较大程度地减少了DWI图像中的噪声，并且通过对降噪图像进行计算所得到的张量也与原始图像趋于一致，说明WNNM算法保存了尽更多的纹理细节特征，提高了张量计算的准确度。

真实数据的实验样本采用的是一套完整的人体大脑MRI图像，大小为128×128×58×7。其扫描层数为58层，每层有7幅图像，其中包含一幅未加权的MRI图像和施加了6个不同方向的梯度脉冲获得的6幅DWI图像，每幅图像的大小为128×128。扩散敏感系数$b$=1 000 s/mm²。同样使用各类降噪算法对DWI图像进行降噪，并对降噪后的图像进行神经纤维跟踪处理，通过跟踪结果来判断降噪效果。

实验对真实数据分别采用各降噪算法进行降噪，根据降噪图像计算得到DTI数据，然后采用经典的纤维连续跟踪法(FACT)算法^[19-20]对DTI数据进行神经纤维跟踪。FACT算法是神经纤维跟踪里最常用的算法。跟踪结果如图 2所示。

从图 2可以看出，使用WNNM算法进行降噪处理后得到的神经纤维图中的纤维密度明显较大，纤维普遍较长且平滑。从图 2的方框中可以看到，使用WNNM算法处理后得到的神经纤维密度最大，且长度普遍较长，而TV算法和PM算法跟踪到的该部分纤维较为稀疏，右边一簇神经纤维几乎没有跟踪到。通过对比其他算法可知，使用PM算法和各向同性算法进行降噪处理后跟踪得到的神经纤维中存在角度弯曲过大的纤维，这在人体大脑中是不可能存在的，而WNNM算法跟踪到的神经纤维中则没有这种弯曲角度过大的纤维，其纤维都较为平滑。

对于真实数据，要对神经纤维跟踪的效果进行好坏的判断，常用的方法是根据神经纤维数量和最长纤维长度来说明跟踪效果^[21]。但是由于图像中噪声的影响，会使较长的纤维被错分为几根较短的纤维，从而出现被跟踪出来的纤维又短又多的情况。因此以纤维数量和最长纤维长度来判断降噪效果并不可靠，本文提出用平均纤维长度来判断降噪效果，平均纤维长度越长则跟踪到的纤维更均匀平滑，能较好的说明降噪效果。平均纤维长度的计算为

$ av{e_{{\text{length}}}} = \frac{{\sum\limits_{i = 1}^n {length\left( i \right)} }}{n} $

(10)

式中，$n$表示为纤维数量，$length$表示每条纤维长度。

用平均纤维长度来对各算法降噪处理后跟踪出来的神经纤维进行统计，其结果如表 3所示。

从表 3的统计数据中可以看出，使用WNNM算法进行降噪后得到的降噪图像，其跟踪得到的神经纤维平均长度最长，说明得到的神经纤维长度普遍较长，则其有效性更高。而使用各向同性算法进行降噪后跟踪得到的神经纤维的最长纤维长度最长，但根据图 2(d)可以看到，该算法得到的神经纤维中有由于错误跟踪得到的神经纤维，因此进一步说明以最长纤维长度来判断其效果是不可靠的。

加权核范数降噪算法利用了图像本身的自相似性而形成的低秩的特性，通过对图像的相似块矩阵进行处理，得到最优化解，较好地保存了图像的局部结构及纹理细节同时也较大程度的减少了噪声。本文根据DWI图像本身的结构相似性程度高和纹理细节多的特点，提出了将加权核范数降噪算法应用于DWI图像的降噪，从而得到更加理想的降噪效果。本文采用模拟实验与真实数据的实验对该算法进行验证，并将该算法与目前DWI降噪中常用的算法进行比较。实验结果表明，本文算法取得的降噪效果更好，所得峰值信噪比较其他算法高出至少20 dB，结构相似性值要高0.2~0.5，且跟踪得到的神经纤维长度要长0.2~0.8, 。加权核范数降噪算法保留了DWI图像中较多的纹理细节信息，提高了DTI数据的准确度，使后续的处理更具有效性。本文对DWI图像的降噪是在建立2维数据的基础上，而DTI是3维的数据，所以我们将对图像的处理扩展到3维空间作为下一步的研究内容。