信号自适应衰减的多壳扩散磁共振成像方法

罗伶俐，王远军

(上海理工大学 医学影像工程研究所，上海 200093)

1 引 言

扩散磁共振成像(Diffusion Magnetic Resonance Imaging，dMRI)技术被广泛应用于诸如阿尔茨海默症、精神分裂症、脑损伤等诸多大脑疾病的研究中.dMRI的基础原理是由于扩散加权MR信号对器官组织中的内生水分子的随机运动十分敏感，这其中应用最为广泛的扩散加权磁共振成像(Diffusion-weighted MRI，DW-MRI)是使用了两种在180°附近的射频脉冲自旋回波，其磁场梯度的方向和幅度在MR序列上均相同[1]，若在两种梯度中氢核的位置不一致时，其磁矩会发生净位移.若一群随机运动的氢核在扩散过程中呈现的相位不一致，会导致整体信号的衰减，进而生成扩散MR信号.

此处用E(q)表示归一化后的扩散MR信号，即是在波矢量q处所测的dMRI信号与不施加扩散灵敏梯度场时的原始MR信号之间的比值.从数学上讲，E(q)与一个重要指标：总体平均扩散传播算子(Ensemble Average diffusion Propagator，EAP)有关，其重建对于dMRI意义重大：不仅能由EAP推导出诸如用于描述纤维束方向特征的方向分布函数(Orientation Distribution Function，ODF)，其他相关细胞大小、方向，及跨膜交换的信息特征同样可由此推导.

EAP的一种重建方法需要一个将粒子扩散过程与表示扩散信号的函数关联起来的理想模型，但前提是信号的函数表示够精确.这其中最经典的成像方法为扩散张量成像(Diffusion Tensor Imaging，DTI)[2]，DTI方法假设E(q)是一个以原点为中心的高斯函数，但这种过于简化的假设对于含有复杂纤维结构(如交叉、相接的纤维束)的体素，建模局限性很大.相关研究表明，在当前dMRI成像分辨率的条件下，超过90%的体素存在复杂的多纤维群[3].为突破DTI的这一局限性，许多基于高角度分辨率扩散成像(High Angular Resolution Diffusion Imaging，HARDI)的方法被提出，这其中包括扩散光谱成像(Diffusion Spectrum Imaging，DSI)技术[4].DSI需在q空间中密集的笛卡尔网格点上进行数据采样，在得出表示归一化信号的连续函数E(q)后，对其进行傅里叶逆变换可得到相应的EAP[5,6]：

(1)

其中P(r)表示运动粒子净位移为r的概率.q为q空间中的波矢量.由于扩散信号是扫描采样所得数据，因此EAP的求解问题简化成了如何根据所采样的散点信号数据，估计出表示信号的连续函数E(q).

但过长的采样时间使得DSI在临床上并不实用.为解决这一问题，近年来的改进方法主要通过建立合适的信号模型以减少信号采样，或用适当的插值基函数表示q空间中的扩散信号.例如，球极傅里叶基函数[7]，贝塞尔傅里叶基函数[8]，以量子力学简单谐波振荡器哈密顿量(亦称为埃尔米特函数)，其本征函数的线性组合表示三维q空间中的扩散信号的MAP[9]方法.而MAPL[10]则是在MAP-MRI的基础上，利用信号的拉普拉斯范数以对所求解的系数进行优化.但由于更多纤维组织的潜在结构信息在高b值条件下更易显现，需将成像条件由单壳推广到多壳条件(即多个b值条件下)时，以上方法均忽略了dMRI信号随b值增大而衰减的相关性质，该性质为研究白质微结构提供了重要信息.Rathi等人[11]以球型脊波作为dMRI信号的插值基函数，用径向函数项作为约束条件规范了整体信号的径向衰减性，并结合信号函数的全变分(Total Variation，TV) 范数优化，该方法在稀疏的多壳采样条件下实现了不错的EAP重建效果.Ning等人[12]则提出用径向基函数(Radial Basis Function，RBF)作为信号的插值基函数，该方法将dMRI信号表示为以q空间中不同位置为中心的各向异性高斯基函数的线性组合，并构建相应的约束条件以保证信号的衰减性，该方法也取得了理想的EAP重建效果并呈现了较强的鲁棒性.然而上述方法中约束信号衰减的条件均为构建信号插值基函数的外部条件，并没有将信号的衰减性考虑进基函数的构建中.为进一步提高EAP的重建精度，同时提升重建算法的计算效率，本文在以RBF作为频域信号插值基函数方法的基础上，添加了信号自适应衰减项以对信号衰减性进行建模，这不仅保留了原方法能够明确估算EAP的二阶、四阶矩张量的优势，且就体模数据的各项重建指标的结果来看，该方法能在提升计算效率的同时，进一步提高了EAP的重建精度.

2 径向基函数在dMRI中的应用

对于径向基函数φ(x)，x∈d，可被写作φ(‖x‖)，‖x‖表示由原点到x点的欧式距离，即该函数在点x处的值只取决于‖x‖.有时也可用马氏距离替代欧氏距离[13].RBF通常被写作式(2)形式以进行函数估计：

(2)

其中s(x)作为待估计函数，可被表示成N个中心点在cn处，对应权重为ωn的RBF的线性组合.相关RBF的基础理论[14]表明：当RBF的数量N足够大，且中心点cn的选取适当时，则RBF几乎可以逼近任意连续函数.

(3)

(4)

由于dMRI信号具有偶对称性质，即E(q)=E(-q)，因此对偶项的权重保持一致：

(5)

式(5)中的第一项v0φ0(q)，可使该模型更精确地估计各方向信号衰减基本呈各向同性的信号(如大脑灰质和脑脊液区域)，而剩余项则可被视作对信号衰减呈各向异性部分的信号估计.从而更全面地对整体dMRI信号特征进行建模.此外，用于插值dMRI 频域信号所用的RBF是一种常见的高斯核函数φ(x)=exp(-σ‖x‖2)，而使用高斯核函数作为插值基函数的优势在于：由于高斯核函数经傅里叶变换之后依旧是高斯核函数，因此可更简便解析地计算后续EAP、ODF以及其他相关成像统计标量.

3 EAP及基于EAP计算的相关指标

3.1 EAP的计算

由于将RBF作为q空间信号的插值基函数，则所估计的信号函数表达式E(q)是一系列高斯核函数的线性组合，因此对其进行傅里叶逆变换后所得的EAP也可对应地被写作相关基函数的线性组合：

(6)

(7)

3.2 方向分布函数(ODF)

ODF常被用于可视化纤维组织的方向特征，可由EAP通过式(8)中的积分公式进行计算：

(8)

(9)

3.3 均方位移(Mean-Squared-Displacement，MSD)及平均四阶位移(Mean Fourth order Displacement，MFD)

由于潜在纤维组织的重要特征可通过对EAP的高阶矩张量计算所得，用高斯核RBF作为信号插值基函数的一大优势在于：可解析地计算EAP的高阶矩统计量，这对于发掘白质结构微小异常的研究非常重要.因此衍生出了基于二阶矩张量及四阶矩张量计算的成像统计标量：均方位移(MSD)，及平均四阶位移(MFD).

其中MSD正比于测量时间内的水分子平均扩散量，MFD则对EAP的各向异性部分更敏感，因此可捕捉到更多潜在纤维组织的相关信息.二者计算方法为：

(10)

(11)

3.4 回归原点概率(Return-to-the-Origin Probability，RTOP)、回归轴向概率(Return-to-the-Axis Probability，RTAP)

用于衡量在两个扩散灵敏梯度场间水分子净位移的概率是dMRI中的重要成像指标，其中回归原点概率(RTOP)的取值由P(0)决定，由式(6)、式(7)可知，RTOP的计算公式为：

(12)

(13)

4 基于RBF的dMRI信号自适应衰减模型

4.1 信号自适应衰减项

在Ning等人提出用RBF作为dMRI信号插值基函数的方法中，约束dMRI信号衰减性质的方法是：构建约束矩阵B={φKq-φK(q+1)|q=1，2，…，Q}，对于所求系数v，通过保证Bv>0以确保信号随b值增大而衰减的性质.然而该约束矩阵的引入占据了整个求解框架90%以上的计算时间.因此为进一步提高重建算法的计算效率，本文在原有径向插值基函数上添加了信号自适应衰减项以对信号幅度的衰减率进行建模：

(14)

(15)

信号自适应衰减项H(q)=exp(-αq2)，α>0，H(q)∈[0，1]对多壳dMRI信号进行单调递减的指数衰减建模，这一项的引入可确保在高b值的低信噪比条件下整体信号依旧呈衰减趋势.当α=0时，式(14)则是针对单壳条件成像的表达式.此外整个求解框架除了计算系数v，还有对H(q)中衰减系数α的估算.而少量的待估参数确保了整个算法的稳定性，以及更好的鲁棒性.

4.2 参数估算

4.2.1 扩散张量计算

在计算插值基函数之前，须先用DTI方法求解二阶扩散张量：

(16)

a)对式(16)两边取对数后，用加权线性回归方法估算等式右边的权重D，由此减少线性求解过程中噪声带来的影响[19].

b)为保证张量的正定性，对于上一步所求的权重进行平方根分解(又称Cholesky分解)：D=U0TU0，其中U0为上三角阵.

c)将U0作为估算初值，基于Levenberg-Marquadt非线性拟合算法拟合式(17)以计算上三角阵U，最终扩散张量由D0=UTU计算所得.

(17)

4.2.2 构建基函数矩阵及约束条件

但当所采集的信号数量较少时(例如小于基函数数量N的情况)，那么方程Av=E就成了一个欠定问题，因此需对所求解的系数向量进行相关正则化以减小误差，结合常用的正则化方法：l2范数正则(又称蒂霍诺夫正则化)，目标函数构建为：

min‖Av-E‖2+λ‖v‖2

(18)

不过基于压缩感知理论[19]：对于欠定问题Av=E，假设向量v的元素分布足够稀疏，则由Av还原出的信号精度越高.因此l1范数正则化也常被用于求解该类线性逆问题.此处设η为待估计信号与实际信号之间的噪声，则目标函数为：

(19)

4.2.3 参数求解

由于在本实验中除了求解系数向量v，还需估算信号的衰减系数α，因此本文采用共轭梯度算法来交替地估算v和α，算法流程为：

算法1.求解系数向量v与衰减系数α

输入：基函数矩阵A，信号向量E

输出：系数向量v，衰减系数α

a) 初始化α

b) whileα>0

c) 根据目标函数式(18)或式(19)求解系数向量v

d) 根据‖Av-E‖2对α求导，使用梯度下降法更新α

e) ifα<0 或‖v‖1不再下降

f) break;

g) end

而对于式(19)中结合系数向量l1正则化的目标函数，本实验中则运用交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)[20]，可有效地将该目标函数分解成一系列更简单的优化问题，此处将辅助变量Z引入式(19)，将其写成另一种形式：

(20)

该算法求解过程中的每一次迭代均通过交替地更新参数v与Z来最小化以下拉格朗日增广函数：

(21)

其中ρ1，ρ2，ρ3分别为对应约束项的惩罚因子，以保证所求解在可行域中.λ1，λ2，λ3则分别是约束条件AZ-E=0，cTZ-1=0，v-Z=0的乘子项.设在第i+1次迭代中，各参数的更新公式为：

(22)

(23)

乘子项的更新迭代公式为：

(24)

(25)

(26)

5 实验方法与结果

本实验中所用的扫描数据来源于图1中的球型物理体模[21]，该体模与人脑内部组织具有相似的扩散特性.该体模上有两道大小为1×0.7cm2，交叉角度为45°的沟槽.

(a) Spherical phantom (b) Baseline image of the dataset

体模数据由西门子3T扫描仪以2mm×2mm×7mm的空间分辨率对体模进行扫描所得.对于五个不同的b值，即b={1000，2000，3000，4000，5000}，每个b值条件下在81个梯度方向下扫描十次，并取十次扫描数据的均值，并用该数据计算所得的成像标量(如MSD等)，作为后续实验数据所对比的金标准.而用于实验的测试数据则是在同扫描条件下，对于单个b值分别取K个梯度方向，K={15，20，25，30，40，45，50，55，60}，按照扫描参数b值的高低分类，实验数据被分为b={1000，3000}、b={3000，5000}这两组.

5.1 实验参数设置

关于算法1：α的初始值设为10-6，而梯度下降法中用于更新α的单次下降步长d=2×10-10.对于结合系数l2正则化的实验方法，设正则化参数λ=0.00025，单次迭代次数为1000次.对于结合ADMM算法的系数l1正则化参数求解方法，设各约束项的惩罚因子为ρ1=ρ2=ρ3=1×10-4，用于判断是否终止迭代的极小值设为ε1=ε2=ε3=1×10-8，单次迭代次数为300次.

5.2 相关评估指标

5.2.1 归一化均方差(Normalized Mean-Squared Error，NMSE)

(27)

(28)

同样，对于RTAP，MSD，MFD等标量的NMSE计算方式与式(28)一致，评估方法也一致：即更低的NMSE意味着更高的重建精度.

5.2.2 估算角度(Estimated angle，EA)

每个体素的纤维主方向主要是由该处的ODF的峰值所对应的方向所决定的.纤维束追踪技术往往根据ODF的峰值来还原大脑中神经纤维束的连接结构，因此ODF的精确重建对于大脑白质中神经连接性的相关研究至关重要.本实验中所用体模的纤维交叉角度为45°，根据双纤维交叉处的每个体素对应的ODF，计算平均EA的公式为：

(29)

其中ΩD表有双峰值ODF的体素集合，且|ΩD|表示该体素集的体素总个数.dx，1，dx，2则为位于x处的体素，其ODF双峰值所对应的两个单位矢量.显然当EA越接近实际值45°时，ODF的重建精度越高.

5.3 实验结果分析

5.3.1 ODF与估算角度

三种方法所估算的ODF可视化结果见图2，图2左侧一列是根据金标准数据计算所得ODF，图2右侧一列是根据在K=60，b={1000，3000}条件下(常用的临床扫描方式)采样数据所估计的ODF.

可见三种方法均能根据金标准数据成功地检测出体素中的交叉纤维.然而三种方法对ODF的重建效果均随着采样数据的减少而逐渐变差：不难看出在交叉区域，三种方法均有未检测到的交叉纤维.图2(d)，图2(f)中未显示左上角和右下角区域是因为在这些各向同性的区域处所估算的ODF出现了负值.

图2 ODF可视化Fig. 2 Visualization of ODF

估算角度结果见图3，可见无论在低b值还是高b值条件下，三种方法均在K=15处取得了最大值，这是由于采样数据过稀疏所导致的.当b={1000，3000}时，由原始方法，l1，l2正则化方法所求EA与实际值的平均角度误差分别是：1.29°，0.46°与1.05°.在高b值处，三种方法所得EA与实际值的平均角度误差分别是：0.43°，0.39°与1.14°.因此结合两种结果来看，l1正则化方法在估算纤维交叉角度问题上更胜一筹.而在高b值处由于信号整体信噪比降低，l1正则化方法使得整体EA的平均值小于45°，但是整体误差相较低b值时有所降低.

图3 (a)b={1000，3000}(b)b={3000，5000}时的估算角度Fig.3 Estimated Angle on b-value shells with (a)b={1000，3000}(b)b={3000，5000}

5.3.2 各标量的NMSE

重建信号对应的NMSE见图4，在b={1000，3000}时，l1，l2正则化方法对应的NMSE相较原方法分别下降了76.4%和64.9%，而在b={3000，5000}时，二者相较原方法NMSE分别下降了65.8%和68.1%，可见改进方法大大提升了信号的重建精度.

图4 重建信号的NMSEFig.4 NMSE of Signal on b-value shells

b={1000，3000}条件下各标量评估结果见图5，可见两种改进方法重建的MFD对应的NMSE远低于原方法，l2正则化方法对应的NMSE相较原方法降低了32%，效果略优于l1正则化方法(二者相差近10.2%).且l2正则化方法取得了最佳的MSD的重建效果，而l1正则化方法对MSD的重建效果却不尽人意，随着K的增大，NMSE显著提高，在K=60处取得了最大值，相较原方法提升了44%.但MSD的整体NMSE保持在0.0012以下，依旧保留了较高的重建精度.

图5 b={1000, 3000}时(a)MFD (b)MSD (c)RTAP (d)RTOP的NMSEFig.5 NMSE of (a)MFD (b)MSD (c)RTAP (d)RTOPon b-value shells with b={1000,3000}

就重建RTAP、RTOP而言，改进方法在K=15处NMSE取得了峰值，甚至RTAP对应的NMSE最大值超过了原方法，这是由于梯度方向过少，采样数据过稀疏所导致的高误差.当K>20，改进方法所重建RTAP对应的NMSE明显降低，l1，l2正则化方法下的NMSE相较原方法分别下降了69.6%及60%.不过对于RTOP的重建，即使在K=15处，基于l1，l2正则化方法的NMSE相较原方法也分别降低了38.5%与44.3%，整体NMSE远远低于原方法(二者分别降低了88.2%与86.4%).

b={3000，5000}条件下各标量评估结果见图6，由图6(c)可知改进方法在K≤30时RTAP的重建效果不如原始方法，可见改进方法在稀疏采样的条件下对RTAP的重建并不理想.但当K>30时，基于l1，l2正则化的改进方法对应的NMSE相较原方法分别下降了29.7%与42.7%.此外改进方法下其他指标对应NMSE的数据趋势与低b值条件下类似，均达到了相对理想的重建效果，且l1正则化方法对MSD的重建效果相较低b值时显著提升，l1，l2约束方法对应的NMSE相较原方法平均分别下降了84.5%和79%，重建效果显著提升.但在高b值条件下标量的整体NMSE高于低b值同条件下的NMSE，这是由于随着b值增大，信噪比随之降低，噪声对信号的影响程度变大所导致的.

图6 b={3000, 5000}时(a)MFD (b)MSD (c)RTAP (d)RTOP的NMSEFig.6 NMSE of (a)MFD (b)MSD (c)RTAP (d)RTOP on b-value shells with b={3000,5000}

本文所有实验数据的运行条件如下：中央处理器为Intel(R)Core(TM) i7-7700 CPU，运行内存与主频分别为8GB和3.60GHz，测试平台是Windows 8环境下的MATLAB 2018b.三种方法的计算时间对比见图7.在b={1000，3000}，同采样数据条件下，l1，l2正则化方法所需平均计算时间分别是原方法的6%与33.8%，大大提升了计算效率.而在b={3000，5000}的同数据计算条件下，l1，l2正则化方法所需计算时间更短，平均计算时间分别为2.72秒及8.32秒.

图7 (a) b={1000,3000} (b) b={3000,5000}时的计算时间Fig.7 Caculating time on b-value shells with (a)b={1000,3000}(b)b={3000, 5000}

结合两种方法对应的NMSE表现来看，在低b值条件下l2正则化方法的重建效果相较l1正则化方法更稳定，能够有效提升各标量重建精度，因此更适用于低b值条件下的成像.而在高b值条件下，虽然两种改进方法对各标量重建精度的提升相差无几，但l1正则化方法对稀疏采样数据的重建相较l2正则化方法略胜一筹，而在采样数据充分(K≥30)的条件下，选取l1正则化方法可兼顾计算效率及各标量重建精度.

6 讨论与总结

本文在用高斯核RBF作为频域dMRI信号插值基函数的方法基础上，引入了信号自适应衰减项，以实现在多壳条件下对每个体素进行信号衰减建模，同时改进了扩散张量的求解方式，运用基于l1，l2正则化的方法分别求解参数以进行对比.该改进算法不仅保留了原方法中以RBF做基函数的优势，即以稳定解析的方式计算对受限扩散信号更敏感的高阶统计量(如MSD、MFD等)，并保证了ODF的计算方式在任何纤维束追踪算法下的普遍适用性，还能在提升整体重建精度的同时大大缩短所需计算时间，减少了原方法中相关约束矩阵带来的计算负担.

但是本文中的实验数据只局限于体模仿真数据，没有将实际扫描过程中的噪声因素引入到重建过程中，且因缺少实际脑成像数据，无法分析该算法对脑成像的实际影响.同时在稀疏采样条件下该改进方法就个别成像指标的重建效果仍有待提升.此外，由于dMRI技术的普遍缺陷：采集信号所需的长时间扫描仍是dMRI领域的一大难题，因此为减少扫描时间，提升dMRI的临床实用性，如何在多壳条件下根据少量梯度方向下采集的稀疏数据恢复频域信号，并进一步重建出相关成像指标是接下来的工作.