采用偏微分方程的磁共振弥散张量各向异性扩散滤波去噪

吴 锡周激流毕务忠谢明元罗代升

1(成都信息工程学院电子工程系，成都 610225)2(四川大学电子信息学院，成都 610065)

采用偏微分方程的磁共振弥散张量各向异性扩散滤波去噪

吴 锡1，2*周激流2毕务忠1谢明元1罗代升2

1(成都信息工程学院电子工程系，成都 610225)2(四川大学电子信息学院，成都 610065)

弥散张量磁共振成像(DTI)是无创研究大脑白质结构及其他人体纤维状组织结构的主要工具。由于合成DTI数据的弥散加权成像数据(DWI)易受噪声干扰，需要有效去噪以保证DTI数据精度和后续应用的实现。使用各向异性扩散滤波理论，综合考虑各方向通道DWI数据的几何形态和结构特点，重构其特征向量和特征值，获得统一平滑的结构张量，以期在有效去噪的基础上最大程度地保持DTI数据几何结构和特征。利用所提出方法在合成弥散张量数据上进行仿真，并在真实脑部DTI数据上进行实验。仿真和实验结果表明，该方法能有效减少噪声对DTI数据的影响，较之现有的时频分析去噪方法，可更准确地恢复DTI数据，减少主分量方向的偏差和部分各向异性值的误差。

弥散张量磁共振成像;各向异性扩散滤波;偏微分方程

Abstract:Diffusion tensor magnetic resonance imaging(DTI)is the main non-invasive utility to reveal the information of local diffusivity of white matter and other fibrous human tissues.Because the diffusion weighted image(DWI)which is used to acquire DTI is sensitive to noise，effective noise removal is required to improve the accuracy of the DTI data and its subsequent applications.Instead of denoising the DWI in different direction separately，an anistropic filtering method synthetically considering the structure information and characteristic of the DTI was proposed in this paper.The eigenvalue and eigenvector of the smoothing structure tensor were reconstructed to remove the noise and keep the structure characteristic simultaneously.Simulations using a synthetic DTI dataset and experiments using an in vivo brain DTI dataset were performed.The results demonstrated that the noise could be removed significantly and the direction bias of the main eigenvector and the error of fractional anisotropy were reduced effectively compared with the common methods.

Key words:diffusion tensor magnetic resonance imaging;anisotropic smoothing;partial differential equation

引言

弥散张量磁共振成像(diffusion tensor magnetic resonance imaging，DT-MRI或者 DTI)，是对大脑白质结构进行无创研究的主要工具［1］。利用脑白质中水分子弥散效应，对脑白质三维结构的每一体素使用一个3×3对称正定矩阵描述，该矩阵即为体素弥散张量。该矩阵通过对每一体素，采集6个非共面的扩散敏感梯度磁场方向的回波衰减信号测量值构成的弥散加权图像(diffusion weighted image，DWI)和一个不施加扩散敏感梯度磁场的MR信号参考测量值解得。在实际扫描过程中，一般采集多于6个的扩散敏感方向信号，在每个方向上采用更多的编码幅度，以此来抑制成像噪声。在得到多于6个的扩散编码方向的测量结果后，可采用最小二乘拟合方法求解弥散矩阵，获得DTI数据。DTI也同时大量应用在脑损伤、脑发育［2］及人体其他纤维状组织(如骨骼肌［3］、筋腱［4］等)研究。

由于DWI在快速成像时极易受系统噪声干扰，同时在扫描成像过程中由于呼吸、心跳等会产生生理性噪声，通过DWI获得的DTI数据受噪声干扰极大［5］，另外DTI的主要应用如纤维追踪成像等对于噪声极为敏感，少量随机噪声即会影响全部成像结果，因此DTI数据应用的首要工作就是去噪，以获得正确的弥散张量。

现有的DTI去噪方法一般分为3类:第一类是对DWI数据使用传统时域或者频域手段进行去噪［6］，但由于是分别对不同方向通道的DWI数据进行滤波，没有综合考虑DTI数据的形态和结构等信息，容易丢失信息;第二类是在DWI估计 DTI的过程中，对该过程进行正则化去噪［7］，但由于 DWI数据本身的复杂性，使该正则过程的适应受到很多约束;第三类是对 DTI数据进行去噪［8］，由于 DTI数据是二阶张量，将极大地增加去噪的计算量和计算复杂度。

本研究提出采用偏微分方程的DTI各向异性扩散去噪方法，综合考虑所有不同通道DWI数据所包含的DTI结构和形态信息，重构特征值和特征向量;对每个方向通道的DWI数据，使用统一平滑结构张量进行去噪，既能保持DTI的结构和边沿信息，同时又能有效去噪。下面，首先介绍采用偏微分方程的DTI各向异性扩散去噪的基本原理，然后阐述采用这种方法在人工合成数据集中的仿真试验和结果，最后介绍采用这种方法对真实脑部DTI数据进行去噪的实验和结果。

1 采用偏微分方程的DTI各向异性滤波去噪

采用偏微分方程的各向异性滤波，可用下述模型表示［9］为

式中，ρ为卷积核的尺度参数。

G的特征向量 v1，v2，v3及其对应的特征值λ1，λ2，λ3表征 ρ内灰度梯度的变化，最大特征值及其对应的特征向量代表灰度梯度变化最大的方向，反之亦然。当λ1≈λ2≈λ3≈0时，为各向同性结构(即平坦区域);当 λ1＞λ2≈λ3≈0时，则为各向异性结构(边缘区域)。结构张量T由灰度的梯度张量G按照一定比例重构特征值获得，有

式中，A为梯度张量G的各向异性系数，a为归一化参数，C为阈值参数。

在平坦区域，A≪C，λ3→a，使用各向同性滤波;在边缘区域，A≫C，λ3→1，使用各向异性滤波，在梯度小的方向进行大尺度滤波，即在边沿的切向方向增强滤波。

不同于传统DWI去噪仅考虑各独立方向通道的数据特点，对于不同方向通道的 DWI数据，综合考虑DTI结构信息，使用统一的平滑结构张量，使在一个DWI方向通道上失去的结构信息可以从其他DWI数据弥补，其梯度张量 G′可由各 DWI方向通道的梯度和与高斯核卷积获得，即

平滑结构张量T′由梯度张量G′重构获得。不同于传统各向异性滤波主要针对线状边沿，DTI由于其生理特点，其结构主要以带状分布为主，特别在脑白质纤维束内部，弥散信息变化较大，容易出现除传统各向同性结构(λ1≈λ2≈λ3≈0)和各向异性结构(λ1≈λ2＞λ3≈0)以外的特殊结构，其特征值满足 λ1≥λ2≥λ3＞0。针对脑白质 DTI数据的特殊性，在各向同性结构区域，希望在各个方向都有较强扩散去噪;在边缘区域，希望沿边沿切向方向有较强扩散去噪，垂直于边沿尽量保持不要扩散;最后对于脑白质带状纤维束内部进行相应小尺度各向异性扩散，以保持数据精度。

根据上述分析，对梯度张量G′进行重构，T′的特征向量为

T′的特征值重构为:

其中，λ1，λ2由下式获得，有式中 c(x，y，z)为权函数，且为的减函数，满足 0≤c(x，y，z)≤1，令其为c(x，y，z)

由上述特征值和特征向量，即可重构平滑结构张量T′。当位于各向同性结构时，各方向特征值均较大，保证各方向均进行较大的扩散滤波;当位于各向异性结构时，沿边缘方向特征值大，垂直边缘特征值小，保证边缘增强去噪;当位于脑白质纤维束内部区域时，特征值均较小，使用较小的扩散保持数据精度。

图1 合成DTI数据用本方法和Parker方法的去噪结果比较。(a)合成DTI三维数据;(b)合成DTI二维放大结构;(c)加入5%高斯噪声的结果;(d)加入10%高斯噪声的结果;(e)和(f)分别加入5%和10%高斯噪声的本方法去噪结果;(g)和(f)分别加入5%和10%高斯噪声Parker方法去噪结果Fig.1 Denoising results of the proposed method and Parker’s method using the synthetic phantom.(a)A 3D view of one slice of the synthetic dataset;(b)Enlarged 2D view of the boxed region in(a);(c)Result of adding 5%Gaussian noise;(d)Result of adding 10%Gaussian noise;(e)and(f)Result of the proposed method to 5%and 10%Gaussian noise;(g)and(h)Result of the Parker’s method to 5%and 10%Gaussian noise

2 仿真与实验结果

在合成数据集和脑部真实DTI数据集中，分别使用本方法和Parker的方法［6］进行去噪，并对结果进行比较。

首先，采用所提出的算法在合成张量数据上进行了仿真，实验平台为Matlab 7.0.4，Windows Vista操作系统，处理器为Intel Core2 Duo CPU(1.80 GHz)，系统内存4 GB。图1(a)为模拟数据集，扫描矩阵分辨率为64像素×64像素，层数为20层，合成纤维束张量的 FA为0.8，弥散度为2.1×10-5cm2/s。如图1(a)所示，合成弥散张量数据分为横向和纵向两组，其余部分为非弥散区域。图1(b)为图1(a)框中图像的放大结果，其中的直线段表征弥散张量，其长短和方向与该体素弥散张量主分量相同，非弥散区域用圆点表示。图1(c)和(d)分别为上述DTI数据集加入5%和10%高斯噪声的结果，由于噪声干扰DTI数据的弥散大小和方向都发生了变化，随着噪声增加，出现了明显的错误 DTI数据。图1(e)和(f)是使用本方法分别对受5%和10%高斯噪声干扰的DTI数据集去噪结果，图1(g)和(h)是使用Parker方法分别对受5%和10%高斯噪声干扰的DTI数据集去噪结果。

由图1(e)和(g)可知，对5%高斯噪声分别使用本方法和Parker方法，均可在一定程度上抑制噪声的影响。相比而言，使用本方法的去噪效果更明显，除少数体素外，DTI数据的弥散大小和方向基本与未受噪声干扰的DTI数据一致，特别是在非弥散区域，基本抑制了噪声影响，见图1(e)。由图1(f)和(g)可知，在噪声较大的情况下，两种方法均可有效去噪。相比而言，本方法的去噪结果其弥散大小和方向与未加噪数据基本吻合，见图1(f)。

DTI数据的主要应用之一为纤维追踪成像，主要使用DTI数据弥散方向信息进行追踪成像。因此，对DTI数据方向信息进行复原和保护是DTI去噪的主要目的。使用本方法，分别对加入5%和10%高斯噪声的DTI数据去噪，对方向信息的复原效果见图2，数据参数同图1。横轴为迭代次数，纵轴为所有体素DTI数据主分量方向与标准DTI数据主分量方向角度偏差的均值。

图2 使用本方法的DTI数据方向信息复原效果分析Fig.2 Variation of mean angular difference of the proposed method in different noise

图3 真实脑部DTI数据用本方法和Parker方法的去噪效果比较。(a)原始FA图;(b)加入5%高斯噪声的FA图;(c)加入10%高斯噪声的FA图;(d)和(e)用本方法加入5%和10%高斯噪声的去噪结果;(f)和(g)用Parker方法加入5%和10%高斯噪声的去噪结果Fig.3 Filtering results using the proposed method and Parker’s method in real human in vivo dataset.(a)One slice FA map of the in vivo dataset;(b)FA map of adding 5%Gaussian noise;(c)FA map of adding 10%Gaussian noise;(d)and(e)Result of the proposed method to 5%and 10%Gaussian noise;(f)and(g)Result of the Parker’s method to 5%and 10%Gaussian noise

由图2可知，较之标准的DTI数据集，被5%和10%高斯噪声干扰的DTI数据集其平均弥散方向分别有大约4.5°和9.5°的偏差。使用本方法进行去噪后，随着迭代次数增加，平均角度偏差逐渐减小，大约在第10个迭代周期后基本达到稳定，平均角度偏差分别下降至2°和5°。

此外，采用所提出的算法在人脑部DTI数据上进行了成像实验。人脑部 DTI数据使用 Philips Intera Achieva 3T MRI扫描仪获得，每次扫描矩阵分辨率为128像素×128像素，扫描层数为53层，层厚2 mm。弥散加权数据集由31个不同梯度方向数据(弥散加权值1 000 s/mm2)和1个未加权基本数据(弥散加权值 0 s/mm2)构成。使用 Philips Diffusion Registration Tool(Release 0.4)进行图像配准，消除偏移等影响。张量弥散模型使用传统的最小二乘模型［1］，并据此计算部分各向异性(fractional anisotropy，FA)。

FA是衡量DTI数据各向异性大小的主要指标，可客观描述体素弥散的大小程度。图3(a)为上述真实人体DTI数据中第25层的 FA图，图3(b)和(c)为分别在原始DWI数据中加入5%和10%高斯噪声后获得的DTI数据第25层FA图。由图可知，由于噪声影响，FA图的准确性下降。

使用本方法，分别对5%和10%的高斯噪声干扰数据去噪，获得的FA图见图3(d)和(e)。分别使用Parker方法，对5%和10%的高斯噪声干扰数据去噪，结果FA图。见图3(f)和(g)。由图可知，两种方法在不同噪声情况下，都可以在一定程度上去噪复原图像，相比而言，本方法可以在更大程度上去除噪声影响。图3(b)和(c)在胼胝体压部周围的非弥散区域，由于噪声影响产生一些错误的弥散信息，使用本方法根据平滑结构张量进行不同尺度的滤波，去噪效果明显，这些错误的弥散区域基本得到修正。同时，本方法在弥散区域根据白质结构自适应地进行各向异性滤波，在去除噪声的情况下尽可能地考虑结构对其的影响，由图3(d)和(e)可见，胼胝体压部的结构保持良好。

较之DTI数据的大小信息，去除噪声对DTI方向信息的影响更加重要。使用本方法，分别对加入5%和10%高斯噪声的DTI数据去噪，表1为对方向信息的复原效果(参数设置同上)，为各数据集体素弥散张量主分量方向与标准数据集弥散张量主分量方向偏离值的均值。

表1 本方法和Parker方法真实DTI数据集滤波去噪方向信息的复原效果比较Tab.1 Variation of mean angular differencein the proposed method and Parker’s

由表1可知，较之标准 DTI数据集，被5%和10%高斯噪声干扰的DTI数据集其主轴平均弥散方向分别有平均5.271°和9.198°的偏差。本方法和Parker方法均能降低噪声对弥散张量数据方向信息的影响，相比而言，本方法对方向信息复原效果更好，平均角度偏离值分别下降至2.087°和4.671°，较之Parker方法，复原结果方向信息更接近于标准数据集。

3 结论

使用偏微分方程的方法，对磁共振弥散张量成像的DWI数据进行各向异性扩散去噪。不同于现有方法对多方向通道的DWI数据独立去噪，本研究综合DTI的结构和数据特征，对各方向通道的DWI数据，按照各向异性扩散原理，重构统一的平滑结构特征张量，重新设计相应的特征值和特征向量。实现在非弥散的各向同性结构区域进行大尺度的各向同性扩散;在边缘区域，平行于边沿方向进行大尺度扩散去噪，垂直于边沿方向进行小尺度扩散;最后对于脑白质带状纤维束内部进行相应小尺度扩散，以保持数据精度。在合成数据和真实人脑部DTI数据集上进行去噪仿真和实验，较之现有方法效果更好。

附录

R方向张量扩散公式

张量矩阵D

部分各向异性FA

角度偏差均值Ψ

FA变化率均值ζ

［1］Basser P，Mattiello J，LeBihan D.MR diffusion tensor spectroscopy and imaging［J］.Biophysical Journal，1994，66(1):259-267.

［2］Susumu M，Zhang Jiangyang.Principles of diffusion tensor primer imaging and its applications to basic neuroscience research［J］.Neuron，2006，51(5):527-539.

［3］Heemskerk AM，Sinha TK，Wilson KJ，et al.Quantitative assessment of DTI-based muscle fiber tracking and optimal tracking parameters［J］.Magn Reson Med，2009，61(2):467-72.

［4］Fechete R，Demco DE，Eliav U，et al.Self-diffusion anisotropy of water in sheep Achilles tendon［J］.NMR in Biomedicine，2005，18(8):577-586.

［5］Anderson AW.Theoretical analysis of the effects of noise on diffusion tensor imaging［J］.Magnetic Resonance in Medicine，2001，46(6):1174-1188.

［6］Parker GJ，Schnabel JA，Symms MR，et al.Nonlinear smoothing for reduction of systematic and random errors in diffusion tensor imaging［J］.J.Magn Reson Imag，2000，11:702-710.

［7］白衡，高玉蕊，王世杰，等.DTI扩散张量的一种稳健估计方法［J］.计算机研究与发展，2008，45(7):1232-1238.

［8］Pennec X，Fillard P，Ayache N.A Riemannian framework for tensor computing［J］.Int J Comput Vis，2006，66:41-66.

［9］Weickert J.Coherence-enhancing diffusion filtering［J］.Int J Comput Vis，1999，31:111-127.

Diffusion Tensor Magnetic Resonance Imaging Denoising Using Anisotropic Partial Differential Function Filter

WU Xi1，2*ZHOU Ji-Liu2BI Wu-Zhong1XIE Ming-Yuan1LUO Dai-Sheng2
1(Department of Electronic Engineering，Chengdu University of Information Technology，Chengdu 610225，China)2(School of Electronics and Information Engineering，Sichuan University，Chengdu 610065，China)

TP391

A

0258-8021(2010)04-0481-05

10.3969/j.issn.0258-8021.2010.04.001

2010-01-12，

2010-05-12

四川省教育厅科技创新重大培育项目(09ZZ004)

*通讯作者。 E-mail:xiey@cuit.edu.cn