基于稀疏编码的弥散微循环模型参数估计神经网络

郑天舒，颜国辉，叶初阳，吴 丹

（1.浙江大学生物医学工程与仪器科学学院，杭州 310027；2.浙江大学医学院附属妇产科医院，杭州 310003；3.北京理工大学集成电路与电子学院，北京 100081）

引 言

弥散磁共振成像（Diffusion magnetic resonance imaging，dMRI）是一种重要的医学成像工具，它基于水分子在生物组织中受限弥散作用，可建立弥散信号的生物物理模型，无创地重建生物组织内的微观结构［1］。常用的表观弥散系数（Apparent diffusion coefficient，ADC）用一个单指数模型来计算，该单指数模型对脑卒中［2］、肿瘤等病理变化敏感［3-4］，但不能特异性地反映微观层面病理变化。近30年来，领域内提出了多种复杂的生物物理模型来表征特定的微结构特性，如体素内不相干运动模型（Intravoxel incoherent motion，IVIM）［5］、神经元定向弥散和密度成像（Neurite orientation dispersionand density imaging，NODDI）［6］、用于肿瘤细胞术的血管、细胞外和限制性扩散模型（Vascular，extracellular and restricted diffusion for cytometry in tumors，VERDICT）［7］，胞体和神经突密度成像（Soma and neurite density imaging，SANDI）［8］等。

准确估计微结构模型参数对诊断具有重要意义。然而，大多数前沿的dMRI 模型由多个房室组成，并且这些房室在数学上复杂且高度非线性。用传统的优化方法例如最小二乘法，对这些模型进行拟合，容易产生估计误差［9］。此外，从数据获取的角度来看，先进的dMRI 模型往往需要在q空间中获取多个b值和弥散梯度方向，这既耗时又容易受到运动伪影的影响。这对于运动的被试，例如腹部器官、胎儿和胎盘等，是一个不小的挑战。

为了减少这些估计误差和加速q空间采集，研究者们提出了许多方法。如Nedjati-Gilani等［10］首先提出了基于随机森林的方法来对Kärger模型进行拟合；同一时间，Alexander等［11］通过随机森林方法实现了在仅用单个b值的情况下，对NODDI 模型和球面平均技术（Spherical mean technique，SMT）模型的参数估计。

随着深度学习技术的发展，深度学习的技术为dMRI 模型拟合开辟了新的途径。Golkov 等［12］首先提出了基于q空间的深度学习方法，这是一种基于多层感知器（Multilayer perceptron，MLP）的方法，它利用q空间数据的一个子集来估计扩散峰度成像（Diffusion kurtosis imaging，DKI）参数。Barbieri 等［13］首次使用基于q空间的深度学习方法来对IVIM 模型进行参数估计。然而，传统的基于q空间的深度学习方法与生物物理模型无关，因此很难被解释。在深层神经网络中引入领域知识作为先验信息被认为是提高网络性能和可解释性的有效途径。Ye［14］首次将通过凸优化加速微结构成像（Accelerated microstructure imaging via convex optimization，AMICO）算法中字典构建的过程与深度学习结合，摒弃之前预先定义字典的方法，通过深度学习方法来得到字典，实现了对NODDI 模型参数的准确估计。

由于之前的基于q空间的深度学习方法［12］缺乏先验知识引导，受此基于稀疏先验的模型构建方法启发［14］，结合之前基于q空间的深度学习方法［12］可以学习到降采样的q空间信号与模型参数映射关系的思想。本文将两者融合，提出了一种基于稀疏编码深度学习网络（Sparse coding deep neural network，SCDNN）的IVIM 模型参数估计网络。该网络由稀疏表示编码器和参数映射解码器两部分构成，其中稀疏表示编码器和参数映射解码器，均通过数据驱动的方式学习得到。首先，假设IVIM 模型所需要的b值具有稀疏性，通过稀疏表示编码器，低维空间中稠密采集的dMRI 信号可以被转化为高维空间的稀疏离散信号；随后，通过参数映射解码器，稀疏离散信号可以与IVIM 模型参数进行体素对体素的映射。将这两部分结合，可以实现对IVIM 模型在q空间的降采样，在采集时间更少的情况下获得近似于相同质量的图像信息，具有较高的准确度和精确性。本文还将所拟的网络和其他已有的dMRI 参数估计网络在胎盘数据上进行了测试，发现基于稀疏编码深度学习网络的方法估计效果优于其他IVIM 模型估计方法，并且更具备模型的解释性。

1 理 论

本节将会简要介绍IVIM 模型、传统IVIM 模型参数的估计方法以及如何通过稀疏表示的方法来估计IVIM 模型的参数。

1.1 IVIM 模型

IVIM 模型由Le Bihan 等［5］提出，基于组织中和血液中水分子扩散率的不同建立双指数模型，可以同时获取组织（细胞、轴突等）中的水分子弥散信息和微循环（毛细血管、血小管）中的血流灌注信息。它已广泛应用于检测组织中的灌注情况，如大脑、肾脏、肝脏和胎盘等器官的灌注异常的检测［15-17］。然而IVIM 成像需要采集多b值（≥10 个b值）的扩散信号，采集较长，容易受到运动伪影的影响［18］。体部器官随呼吸运动伪影较大，特别地，胎盘的成像不仅受到呼吸的影响还受到胎儿运动的影响。

常见的IVIM 模型可以表示为

式中：b为施加的弥散梯度场强度，Sb为该b值下的信号值，S0为没有扩散加权时的信号值，f为微循环的体积分数，D为组织水分子的弥散系数，D*为微循环血液中水分子的伪弥散系数。

1.2 IVIM 模型的参数估计

常见的IVIM 模型参数估计方法有两类，分别为基于非线性最小二乘法的参数估计和基于贝叶斯的参数估计。对于胎盘而言，基于非线性最小二乘法的参数估计先通过分段方法，通过相对高b值（＞200 s/mm2）的数据，采用单指数模型估计D值以及组织中水分子组分对应的非扩散加权信号；通过相对低b值（＜200 s/mm2）的数据通过外插方法估计组织加血液对应非扩散加权信号，f值可通过进行估计；D*值可估计为D×10 得到f、D、D*的估算值，再以此为初始条件和上下限条件，通过建立非线性最小二乘拟合方法来对其中的f、D、D*来进行求解。贝叶斯估计则通过结合先验信息，通过贝叶斯公式，优化最大化后验概率的分布来求f、D、D*［19］。

1.3 IVIM 模型的稀疏表示

稀疏表示可以被认为是一种表示学习方法，它已被普遍用于磁共振（Magnetic resonance imaging，MRI）的加速问题中［20-21］。对于IVIM 模型，q空间的弥散信号可以用以下字典学习框架来表示。

式中：y=Sb/S0∈RK是一系列通过S0归一化后的信号Sb，其中K为该样本中的弥散梯度个数；φ∈RK×N是被定义的字典，其中N为稀疏系数向量的长度，在本文中，该字典通过神经网络学习得到；x为与字典φ相关的稀疏系数向量；η为与x对应的噪声向量。结合IVIM 模型，φ和x可被拓展为

式中：φD、φD*分别与离散化的D和D*相对应；x1-f对应于组织成分（D）的组成比例；xf为灌注成分（D*）的组成比例。其中为了满足各组分之和为1 的限制，稀疏系数向量x被归一化到［0，1］之间。

如果φ已经被预先定义好，则稀疏系数向量可以被通过求解如下的带有l1正则化的目标函数得到。

式中：β为一个常数，用来控制稀疏系数向量x的稀疏性，β越大时，稀疏系数向量x越稀疏。l1正则化的优化可以通过许多算法来解决。在MRI 领域，普遍使用迭代收缩阈值算法（Iterative shrinkage thresholding algorithm，ISTA）［22-23］。优化该目标函数可以通过如下的迭代过程展开。

式中：I为单位矩阵；HM为非线性算子，可定义为

通过如上步骤，可以获得稀疏表示的信号。随后，可以通过如下的线性组合得到模型估计参数。

2 方 法

本节将从如何构建基于IVIM模型的深度神经网络、训练所需样本以及训练方法3个方面来进行阐述。

2.1 网络构建

如前所述，基于IVIM 模型的深度神经网络由稀疏表示编码器和参数映射解码器构成，具体方法如图1 所示。

图1 基于IVIM 模型的SCDNN 结构示意图Fig.1 Structure of IVIM model based SCDNN

对于稀疏表示编码器，它将该接受q空间降采样的输入信号，该输入信号通过信号提取层后输入字典层A后得到其对应的输出信号，随后将这一输出信号输入连续的迭代单元，迭代单元依照可以通过对式（6）字典构建的方法进行设计，通过迭代过程实现字典的构建，其中红色虚线所示的即对应迭代过程。

每个迭代单元中输入信号依次经过阈值层和字典层B，阈值层为一非线性算子，对应于式（6）中的HM，字典层B 则对数据进行组合，之后将字典层B 输出的信号与该迭代单元的输入信号通过残差连接相叠加，作为通过该迭代单元后的输出信号并作为下一个迭代单元的输入信号；第n个迭代单元中通过对信号的处理实现阈值迭代算法在IVIM 模型中的通项表达式：xn+1=HM[φHy+(I-φHφ)xn]，其中φH和φ对应于输入的线性层和迭代单元中线性层所学习到的权重，φ为由离散化的组织水分子的弥散系数D和微循环血液中水分子的伪弥散系数D*所构成的字典向量，I是单位矩阵，HM为非线性算子，y为输入信号，xn和xn+1分别代表第n个迭代单元的输入信号和输出信号。如此进行反复迭代，得到通过所有迭代单元后的最终输出信号，将最终输出信号输入归一化层后，得到体积分数f，通过归一化后，使信号限制在0～1，满足实际生理状况；再将通过归一化层的输出分为等长的两部分，分别进行再次归一化后输入两个不同线性组合层进行信号组合，一个线性组合层输出由离散化的组织水分子的弥散系数D重新组合，另一个线性组合层输出由离散化的微循环血液中水分子的伪弥散系数D*重新组合。

2.2 数据及预处理

采用29 名孕期在13～37 周的被试数据中进行了数据集的构建。

2.2.1 序列参数

29 名磁共振扫描在通用电气SIGNA HDXT 1.5 T 磁共振扫描仪上进行，采用扩散加权的平面回波成像序列从母体矢状位采集胎盘图像：回波时间（Echo time，TE）/重复时间（Repetition time，TR）=76/3 000 ms，视野（Field of view，FOV）=320 mm×320 mm，平面内分辨率为1.25 mm×1.25 mm，层厚4 mm，共15 层，共采集10 个b值的数据，分别为0、10、20、50、80、100、150、200、300 和500 s/mm2。

2.2.2 配 准

对29 例被试中的每个胎盘数据进行自配准，以去除母体呼吸、胎儿运动导致的各个b值扩散加权图像间的位移。在配准过程中采用循环的配准方式操作6～10 次之间不等。对每一位被试，将其各个b值的扩散加权图像进行平均，得到平均模板。然后，通过刚体变换和仿射变换将各个b值的扩散加权图像与平均模板进行比对配准。将配准之后得到的数据再次进行平均，得到第二次的平均模板，如此循环6～10 次，选取其中位移最小的图像作为得到目标图像的配准结果。

2.2.3 金标准获取

将2.2.2 节中配准好的图像，通过贝叶斯方法进行参数估计［24］，得到对应的f、D、D*。由于通过贝叶斯方法得到的f、D、D*存在因噪声造成的误差，因此将对应的f、D、D*代入式（1）中，得到该组信号的无偏信号，并添加噪声。噪声的添加方式与文献［25］相同。

式中：Snoisy为数据集的输入信号；S为通过式（1）得到的信号；ξ1，ξ2～(0，σ2)，σ=S0/SNR，其中与文献［25］相同，S0=1，训练数据SNR=30。

2.3 训 练

将预处理好的29 例孕周在13～37 周的正常孕妇样本按照6∶4 的比例划分为训练集和测试集。其中训练集中有15 例样本（324 016 个体素，其中与文献［12］相同，90%用于训练，10%用于验证），测试中有9 例样本（248 059 个体素）。本文选择使用Adam 作为优化器，学习率设置为1×10-4，总计10 个历时，批次大小为128。在网络训练的前3 个历时中使用了余弦预热方法，学习率在3 个历时后衰减。

3 实 验

对于2.2 节中获取的数据集，采用基于IVIM 模型的深度神经网络进行了估计，并与其他用于dMRI 参数估计的方法进行了性能对比。同时，对基于IVIM 模型的深度神经网络参数设置进行了分析。

3.1 性能分析

对基于IVIM 模型的深度神经网络进行了性能分析。首先，对基于IVIM 模型的深度神经网络的在训练集和验证集上的损失函数衰减曲线分别进行了观察，相关结果如图2 所示。

图2 基于SCDNN 的IVIM 参数估计Fig.2 Parameter estimation of IVIM model based on SCDNN

3.1.1 收敛性分析

结合图2，可以观察到，第1 次迭代训练的损失为0.33 左右，验证集误差为0.27 左右，说明模型刚开始训练时，损失较大，准确率较低。随着迭代周期增加，训练集误差和验证集误差在5 个历时后趋于稳定，可以认为其收敛，此时得到的网络参数就是最终优化的网络参数指标。

3.1.2 结果分析

首先对于一例样本进行分析，结果如图3 所示。对于该例样本，计算各参数网络输出的结果与其对应的金标准之间的决定系数（Coefficient of determination）R2。R2可定义为

式中：yi为该参数的金标准；fi为该参数网络估计值；yˉ为该参数金标准的平均值。对于3 个参数f、D、D*，其对应的R2如图3（a）所示，分别为0.998 2、0.999 1、0.998 5。3 个参数的网络估计结果示意图及残差图（参数金标准-参数网络输出结果）如图3（b）所示。由图3（b）中的残差图可以看出，f、D、D*均可以达到较好的估计效果。

图3 神经网络输出与金标准的相关性图及神经网络输出的参数图和误差图Fig.3 Correlation between neural network output and gold standard and parameter maps and error maps of neural network output

在全部测试集上比较了5 种不同算法，包括两种传统优化方法（最小二乘法和贝叶斯估计）以及3种基于学习的方法（q空间MLP［12］、IVIM 网络（IVIM-NET）［13］和本文提出的基于SCDNN 的IVIM参数估计），结果如表1 所示。与传统优化方法最小二乘法和贝叶斯方法相比，基于学习的方法具有更高的估计精度。在基于学习的方法中，模型驱动方法的性能优于无先验信息的q空间的深度学习。

表1 比较5 种方法在使用减少的b 值（20、50、150、300、500 s/mm2）估计弥散微循环模型参数的均方误差Table 1 Mean squared error comparison of five algorithms using a subset of b-values（20，50，150，300，500 s/mm2）for IVIM parameter estimation

3.2 泛化性分析

本节对基于IVIM 模型的深度神经网络进行了泛化性分析。为了评估基于IVIM 模型的深度神经网络的泛化性，在另一台通用电气GE 750 W 的3.0 T 磁共振扫描仪上通过相同的扫描序列参数进行了数据采集。采用扩散加权的平面回波成像序列从母体矢状位采集胎盘图像：TE/TR=76/3 000 ms，FOV=320 mm×320 mm，平面内分辨率为1.25 mm×1.25 mm，层厚4 mm，共15 层，共采集10 个b值数据，分别为0、10、20、50、80、100、150、200、300 和500 s/mm2。与3.1.2 节相同，比较了所有的5 种算法，结果如表2 所示。从表2 中可以观察到，在同一品牌的3.0 T 数据上进行测试，发现其他基于学习的方法有不同程度上的误差升高，其中，无监督IVIM-NET 误差升高最多，本文方法具有最高的准确度，证明本文方法具有更好的泛化性。

表2 泛化性测试比较5 种方法在使用减少的b 值（20、50、150、300、500 s/mm2）在另一个中心采集到的数据估计弥散微循环模型参数的均方误差Table 2 Generalizability test on data acquired from another center，based on the mean squared error of five algorithms using a subset of b-values（20，50，150，300，500 s/mm2）for IVIM parameter estimation

3.3 网络参数分析

本节对SCDNN 的参数设置进行了分析。首先，对于不同b值的选取进行了实验，之后对于所建立字典大小进行了实验。

3.3.1b值组合分析

本节对于如何选取IVIM 模型训练时所需要的b值数量和b值组合进行了研究。在之前的实验中，选取了减少到5 的b值数量（20、50、150、300、500 s/mm2）来进行估计。在这里对5 个b值数量的不同组合进行了估计，同时也对不同数量的b值选取进行了估计。其中对数量为5 个，b值大小为20～500 s/mm2的5 种组合进行了估计，结果如表3 所示。

结合表3，可以发现，在b值数量为3 时，参数f、D、D*的估计结果与b值数量为5 或7 时相比较差。当b值数量为7 时，参数f、D、D*的估计结果与b值数量为5 时最优结果总体误差相接近，且并不优于b值数量为5 时的结果。因此，在最终的b值选取上选择5 个b值组合为20、50、150、300、500 s/mm2的降采样方式。

表3 不同数量及不同b 值组合得到的IVIM 参数估计值与金标准的均方误差Table 3 Mean squared errors of IVIM parameter estimation using different numbers of b-values and different combinations of b-values

3.3.2 字典大小分析

本节对字典大小为200～600 的结果进行了分析，结果如图4 所示。从图4 可以观察到，随着字典大小从200 到600 的增加，验证集均方误差逐渐降低，字典大小从600 增加800，验证集误差升高；并且随字典大小的增加，训练时间显著增长。因此选择字典大小为600。

图4 SCDNN 随字典大小变化验证集误差及训练时间示意图Fig.4 Changes of the validation error and training time of SCDNN with dictionary size

4 结束语

本文通过提出了一种基于稀疏编码深度学习网络的弥散微循环模型参数估计方法来估计IVIM 模型参数f、D、D*，该方法借助深度学习方法，结合弥散微循环模型的物理原理，可以通过少数b值的q空间信息拟合，有效缩短了采集数据所消耗的时间，并且能够比现有的参数估计方法有更小的误差，同时从多中心比较实验来看，该算法具有良好的泛化性能。