提取偏头痛患者等的神经影像特征并进行识别模型的设计对相关疾病的辅助诊断具有重要意义。相较于常用的影像特征，本研究直接采用时间序列信号表征偏头痛患者组和健康对照组的大脑功能状态，可有效利用时间信息并减小分类模型训练计算量。首先，本研究针对小样本群体运用组水平独立成分分析和字典学习划分不同脑区后，提取区域平均时间序列信号；其次，将提取的时间序列平均划分成多个子时间序列，以扩充模型输入样本；最后，使用双向长短期记忆网络对时间序列建模，学习每个时间序列内部的前后时序信息来刻画周期性大脑状态变化以提高偏头痛的诊断准确率。研究结果显示，偏头痛患者组与健康对照组的分类准确率为96.94%、曲线下面积为0.98，且计算时间相对较短。实验表明，本文方法具有较强的适用性，时序特征提取和双向长短期记忆网络模型结合能较好地用于偏头痛的分类诊断；这项工作为基于小样本的神经影像数据的轻量化诊断模型提供了新的思路，并有助于相关疾病神经鉴别机制的探索。
引用本文： 孙昂, 陈宁, 何俐, 张俊然. 基于小样本功能磁共振数据的偏头痛时序特征分类研究. 生物医学工程学杂志, 2023, 40(1): 110-117. doi: 10.7507/1001-5515.202206060 复制
0 引言
偏头痛是一种常见的慢性神经系统疾病，影响着全球约15%的人群，其典型症状有反复发作的中重度头痛、恶心、呕吐，以及对光和声音异常敏感等[1]。由于偏头痛的诊断是基于不同症状的组合，需排除可能的继发性头痛病因，目前尚缺乏可以明确诊断的生物学或影像学标志物。考虑到偏头痛疾病本身的特有症状，以及疼痛相关症状的重复发生有可能使得偏头痛患者发生特异性的脑功能改变，最近功能磁共振成像技术和相应数据分析方法的进步，使得从脑功能数据中提取偏头痛的生物标志物成为可能[2]。静息态功能磁共振成像（resting-state functional magnetic resonance imaging，rs-fMRI）是指受试者在不执行思考、运动等任务的情况下，通过检测血氧水平依赖（blood oxygen level dependent，BOLD）信号来表征大脑自发的神经活动的一种成像技术，是评估神经、精神疾病脑功能变化的有效技术手段之一[3]。
近年来，随着深度学习框架及其系列算法在神经影像数据分析处理领域的逐步深入，其在相关疾病分类诊断、疗效预测等方面已展现出巨大的潜力。目前，偏头痛的分类研究大多采用了基于图像的深度学习分析框架，即输入rs-fMRI图像给卷积神经网络（convolutional neural network，CNN）等深层卷积模型，并通过该模型进行特征提取和分类训练。如Yang等[4]提取了rs-fMRI数据的三种特征映射功能图（低频波动幅度、区域同质性和区域功能相关强度），并对特征映射功能图沿z轴和时间轴方向进行切片共生成了5 760张图像，采用了CNN模型其分类准确率最高接近99%。李翔等[5]提出了一种三维CNN无先兆偏头痛基础诊断（three dimensional CNN based diagnosis of migraine without aura，MwoA3D-Net）的智能辅助诊断算法，其最高诊断准确率也达到了98.40%。但是，基于图像输入的CNN通常对样本量和计算资源的要求较高，三维数据和大多数二维图像都比一维时间序列需要更多的计算资源和时耗。与深度学习应用中最常见的自然图像数据集相比，神经影像数据集通常样本量相对较少，大多局限在40～160例[6]；这就导致深层卷积网络等深度学习模型无法充分挖掘图像特征，容易出现过拟合等不足。具体到偏头痛研究领域，其影像数据相较于抑郁症等疾病样本量更少，因此需要新的特征提取方式和模型架构设计等方法层面的创新来解决偏头痛病患这种小样本影像数据的分类识别问题。
在特征提取方式设计上，本文拟直接采用时间序列及其变换作为模型训练输入数据，这样可以大幅减小计算资源开销和训练时耗，并且更适用于小样本数据的学习拟合。另外，直接针对时间序列进行处理也可以避免特征转换过程中可能出现的信息损失[7-9]。为了更好地提取小样本群体的脑区时间序列，本研究拟将组水平独立成分分析（group independent component analysis，Group ICA）方法和字典学习（dictionary learning，DictLearning）这两种数据驱动的方法结合在一起，以识别偏头痛的敏感脑区并研究它们的时序特征。其中，Group ICA可以将多变量信号分离为附加子分量，自动分离出rs-fMRI信号中相互独立的信号源，包括有意义的静息态功能网络和头动等噪声信号[5]；而DictLearning可以通过将整个大脑信号分解为多个各异的时间特征来估计静息态网络的字典，为观察到的样本提供有效的表示[10-12]。在模型架构设计上，为了直接从rs-fMRI时间序列中学习功能特征和大脑状态之间的映射关系，本研究拟采用长短期记忆（long short-term memory，LSTM）网络来捕捉大脑动态时序变化以用于偏头痛分类研究；该网络能够“记住”并加权时间序列或信号的过去值，通过从数据中学习来自适应地捕捉时间依赖性[13-14]。其中，双向长短期记忆（bi-directional long short-term memory，BiLSTM）网络能够挖掘隐藏在时间序列中的上下文信息[15]，与LSTM网络相比，更有机会捕捉到神经元活动的特异信息。此外，以往研究中提出rs-fMRI数据每分钟的BOLD信号对于每个受试者来说可能反映相似的神经元活动[16]，故本文对提取出的时间序列进行了裁剪，以扩充模型的输入序列。
相较于之前的偏头痛数据分类识别框架，本文在三个方面做出了改进：首先，本研究基于数据驱动方法来为小样本数据构建更合理的感兴趣区域（regions of interest，ROIs）划分，并基于ROIs提取时间序列。接着，本研究将所有受试者的输入时间序列裁剪为固定的子序列长度（68 s），扩充模型的输入序列以充分利用整个时间序列蕴含的动态信息。最后，本研究进一步使用BiLSTM网络在两个方向上访问远程上下文，从而捕捉神经元活动前后时间的依赖特性。目前，尚未见基于时间序列特征来进行偏头痛疾病识别的相关深度学习研究，本研究将有助于拓展这一领域的技术边界。
1 材料和方法
本文方法主要从特征提取方式和模型架构设计两个方面进行了改进，模型框架图如图1所示，主要包括rs-fMRI信号分解、提取区域平均时间序列、序列裁剪扩充、BiLSTM网络分析、归一化指数（softmax）函数分类。首先，利用Group ICA和DictLearning对rs-fMRI数据进行信号分解，提取出相应ROIs的区域平均时间序列。然后，将时间序列按时间点分段裁剪，扩充得到5倍样本量的子序列。接着，将子序列送入BiLSTM网络中进行模型训练，学习前向和后向信息。最后通过softmax函数完成二元分类，实现偏头痛患者组和健康对照组的识别与分类。
图1
偏头痛分类预测框架图
Figure1.
Framework of migraine classification and prediction
1.1 数据采集及预处理
1.1.1 受试者
本研究共招募了64名成年人参与rs-fMRI数据采集，如表1所示，分为两组：① 36名偏头痛患者，其中男性11名，女性25名，作为偏头痛患者组。偏头痛患者均是神经科医师根据《国际头痛疾病分类（第2版）》（international classification of headache disorders Ⅱ，ICHD-Ⅱ）[17]诊断标准进行评定入组。② 28名健康受试者，其中男性13名，女性15名，作为健康对照组。健康对照组年龄、性别和受教育年限与偏头痛患者组相匹配。两组受试者均排除了药物或酒精滥用、磁共振成像扫描禁忌、患有其他慢性疼痛疾病和神经发育障碍等情况。
表1
64名受试者的人口统计学信息（
）
Table1.
Demographic information of the 64 subjects（
）
所有的受试者年龄均在18～50岁之间，且均为右利手。研究中对所有受试者进行了24项汉密尔顿抑郁量表（24-Hamilton depression scale，24-HAMD）和14项汉密尔顿焦虑量表（14-Hamilton anxiety scale，14-HAMA）测试，相应结果如表1所示。四川大学华西医院磁共振研究中心对所有受试者进行了rs-fMRI脑部扫描，其中偏头痛患者在脑部扫描前至少72 h和扫描后48 h无发作。本研究经四川大学华西医院医学伦理委员会批准，并获得每位受试者的书面知情同意书。
1.1.2 rs-fMRI数据采集
本研究64名受试者的rs-fMRI数据通过梯度平面回波成像序列获取，其具体设置如下：重复时间/回波时间 = 2 000 ms/30 ms；翻转角 = 90°；切片厚度 = 5 mm；视野 = 240 mm × 240 mm；体素大小 = 3.75 mm × 3.75 mm × 5 mm；采集矩阵 = 64 × 64。所有受试者的扫描时长均为6 min（即360 s），对应180个时间点。
1.1.3 数据预处理
使用神经影像脑连接组分析软件图论网络分析工具箱（graph theoretical network analysis，GRETNA）对采集数据进行预处理。处理步骤包括：
（1）考虑到初始rs-fMRI信号不稳定和受试者前期的适应时间，去除每个功能数据的前10个时间点，剩余170个时间点（即340 s）用于后续处理；
（2）剩余的数据首先进行时间层校正，然后重新对齐以去除头动伪影的影响。将头动位移超过2.5 mm或者头动旋转角度超过2.5°的受试者剔除。本研究中所有受试者的头动位移和头动旋转角度均未超过上述标准。
（3）然后，使用蒙特利尔神经学研究所（Montreal Neurological Institute，MNI）的模板将得到的图像配准至3 mm体素的标准空间。
（4）最后对每个体素时间序列进行带通滤波（0.01～0.08 Hz），以减少低频漂移和高频生理噪声的影响。
1.2 rs-fMRI时间序列特征提取
预处理后的图像通常会进行脑区划分以降低后续分析的数据维度。脑网络划分的常用方法：一种是使用自动解剖标记（automated anatomical labeling，AAL）等固定图谱；一种是通过数据驱动的方式来确定大脑ROIs图谱。基于大脑固定图谱的方法在获取脑区平均时间序列时易受主观因素影响，可能存在一定的偏差，本文中选择两种数据驱动的方法来获取所需的ROIs时间序列。
考虑到本文的小样本群体，基于Group ICA的方法能够将模型应用于组数据[18]，在群体水平上对大脑进行分割，寻找在各个受试者之间较为一致的公共成分[19-20]。Group ICA方法中，先设Y是T×M的观测信号，即预处理后的rs-fMRI图像数据矩阵。其中T表示时间点的个数，M表示受试者大脑中体素的总数。假设每个体素的rs-fMRI信号是各个成分信号的线性混合，则Y可以分离成相互统计独立的成分，记P为N个分离的成分图，W是混合矩阵，其计算公式如式（1）所示：
 |
对于每个受试者的成分图空间分布未知的情况下，可利用盲源分离方法求取，记Ak为分离矩阵，它的列表示相应成分图的时间序列，k表示迭代次数，其计算公式如式（2）所示：
 |
Group ICA的成分估计过程需要进行跨个体的独立成分的融合，才能反映出一组rs-fMRI数据的整体特性。分离矩阵Ak需要在多个受试者中逐次迭代，直至算法收敛，从而得到各受试者一致的成分映射。
由于Group ICA需要遵守数据表示子空间上的正交性约束，这可能导致独立源数量受限于信号维度。针对这一不足，DictLearning可以通过构建一个原子或子空间的字典，为观察到的样本提供有效的稀疏表示，并有可能推导出稀疏信号的先验未知统计量[21-22]。同时，相关神经科学研究表明大脑神经活动中存在稀疏响应[11]，稀疏表示方法的内在性质可能与大脑神经活动的稀疏响应相契合。
在DictLearning方法中，给定rs-fMRI信号矩阵Y，维度为T×M，其中T是rs-fMRI时间序列点的总数，M是受试者大脑中体素的总数。Y中的rs-fMRI信号被建模为学习到的基础字典D，模型矩阵形式如式（3）所示：
 |
其中A保存稀疏表示的系数矩阵。具体来说，D中每个字典原子的信号代表特定脑网络的功能活动，而系数矩阵A中的向量代表对应脑网络的空间分布。字典D可用于划分ROIs，并提取相应的时间序列。
1.3 特征融合与数据扩增
Group ICA方法参照惯例在此提取了20个独立成分，每个成分涉及1～5个脑区，共计将大脑划分为58个ROIs。同样，DictLearning方法依据字典提取了20个原子，每个原子涉及2～8个脑区，共计将大脑划分为96个ROIs。两种方法在脑区划分中各有其优势和特点，本文结合了两种方法获得的154个ROIs并计算了各个ROIs的平均时间序列，生成维度为 (x, y, z)=(64, 170, 154) 的特征矩阵。其中x是样本数，y是时间点，z是ROIs数。接着，对每个ROIs时间序列应用了从-1～1的线性归一化。
对于每个受试者170个时间点（即340 s）的初始时间序列可被平均拆分成5段长度为34个时间点的子序列，每段代表68 s的成像时间。如图1所示，序列裁剪扩充后每一段68 s的子序列都可以作为模型的输入数据，这使模型数据集扩充到初始的5倍，总共达到320个输入序列。重构后的特征矩阵维度为 (320, 34, 154)，这将作为接下来BiLSTM分类网络的输入。
1.4 BiLSTM网络分类
本研究为了进一步学习每个时间序列内部的前后时序信息来刻画周期性大脑状态变化，在对rs-fMRI时间序列裁剪扩充的基础上采用了BiLSTM网络。如图1所示，输入序列 
经过LSTM层，学习前向和后向信息，为每个时间点返回相应的值 
。接着使用带有S型生长曲线（sigmoid）激活函数的密集层和带有softmax激活函数的第二个密集层作为二元分类器。模型使用二元交叉熵损失函数和具有默认参数值的自适应矩估计优化器进行训练。为了避免过拟合，本研究在所有隐藏层中引入随机失活（dropout）机制，随机失活50%神经元，在LSTM层中加入了正则化项L1和L2。实验设定了300次迭代训练，提前停止条件为验证损失连续在30次训练中小于0.01。此外，本研究通过将数据分为8份来执行交叉验证方法，以便受试者比例在所有测试中都大致相同。
2 结果
为了评估性能，实验中进行了8折交叉验证，计算了灵敏度、特异性和分类精度。如表2所示，Group ICA结合DictLearning的方法的敏感性、特异性和准确性均优于单一Group ICA或单一DictLearning。同时，本文绘制了受试者工作特性（receiver operating characteristic，ROC）曲线并计算了曲线下面积（area under curve，AUC），如图2所示。其中，曲线的横坐标为假阳性率（false positive rate，FPR），纵坐标为真阳性率（true positive rate，TPR）。ROC曲线所覆盖的区域面积AUC可以用来定量评价偏头痛患者组和健康对照组分类的准确性，ROC曲线越靠近平面的左上角，AUC的值越高，表明分类器分类效果越好。
表2
不同时间序列提取方法的分类结果
Table2.
Classification results of different time-series extraction methods
图2
受试者工作特征曲线
Figure2.
Operating characteristic curves of subjects
本文方法的平均AUC达到了0.98，并且显著优于传统的AAL固定图谱、单一Group ICA提取图谱和单一DictLearning提取图谱的结果。AAL的效果最差，这可能是由于样本特性与模板不适配，时间序列中存在噪声造成的。相比之下，本文方法结合了两种数据驱动的ROIs划分方法，可以根据数据特性学习内在的大脑状态转换。在8折交叉验证（ROC fold1～ROC fold8）下，每一折的AUC都在0.98～1.00范围内，并且在几种方法中数值波动最小，表现出对时间序列噪声良好的鲁棒性。
如表3所示，使用Group ICA结合DictLearning提取特征后，分类网络和输出单元数量的不同影响了分类精度。结果表明BiLSTM显著优于单向的LSTM，并且具有32个隐藏节点的BiLSTM准确率为96.94%，达到了最高的分类精度。
表3
LSTM与BiLSTM分类网络结果对比
Table3.
Comparison results of LSTM and BiLSTM classification networks
3 讨论
3.1 特征提取方式的对比分析
考虑到临床预测模型对部署难易度和计算效率的客观要求，本研究确立了以一维时间序列为核心的分析模式。传统rs-fMRI数据的时间序列分析，是通过皮尔逊相关等功能连接性度量来构建脑功能连接网络，以生成鉴别性特征用于后续分类。但这些特征构建方法以及相应的分类器没能充分考虑时间依赖性[7]，可能会影响模型解码性能。本研究前期在偏头痛患者的数据上做了预实验，通过皮尔逊相关构建了功能连接矩阵作为分类特征。与本文方法的最高分类精度准确率达96.94%相比，以功能连接矩阵作为分类特征准确率可达88%，则不算理想。Dvornek等[14]也在自闭症患者数据上使用时间序列特征与功能连接度量方法做了对比，时间序列特征同样展现了较优的分类性能。
在选择大脑ROIs图谱以提取时间序列过程中，本研究结合了两种数据驱动方法来划分ROIs，充分利用了体素间的空间联系，从而降低群体差异带来的影响。常用的固定脑图谱模板可能由于对大脑的划分精细度不同，或构建的群体基础不同，而对小样本数据的分析结果有较大的影响。如表2所示，Group ICA（58 ROIs）和DictLearning（96 ROIs）的指标均优于固定图谱AAL（116 ROIs），即使它们的脑区划分精细度并没有AAL高。这源于Group ICA具有良好的控制群体模型，可以分解出大脑信号中的背景噪声；而DictLearning则可以提取到自然稀疏的脑区划分，捕捉到数据的内在特征。Group ICA加DictLearning的方式结合了两种方法的优点，以结合方法来划分区域提取时序特征可以达到更好的分类精度。
3.2 不同模型框架的影响
LSTM是一种特殊类型的递归神经网络，内含的时间循环结构可以较好地刻画具有时间关联的序列数据，并学习到受试者大脑动态的模式信息。考虑到大脑可能会不断地使用上下文信息来指导更高层次的认知功能，而不是在时间序列结束时产生具有严格方向的输出，本文中使用了BiLSTM和LSTM进行了对比分析。
如表3所示，在相同隐藏节点设置下，BiLSTM准确率一直优于LSTM，表明了双向信息更有利于表征大脑信号。此外，隐藏节点的数量可能直接影响LSTM、BiLSTM网络的学习能力。本文比较了具有不同数量（8、16、32、64、128、154）隐藏节点输出单元模型的性能。在单层模型中，具有32个隐藏节点的BiLSTM达到了最高的分类精度。具有8个或16个隐藏节点的模型其性能有所降低，很可能是相应的隐藏节点输出单元太有限，无法存储足够的序列信息。具有64个以上隐藏节点的模型性能也略有下降，其原因可能是过拟合。如图3所示，在154个节点的参数设置下，训练准确度和验证准确度两条曲线整体波动较大，相距较远，说明模型的方差较大，模型有过拟合趋势。
图3
训练验证准确率变化曲线
Figure3.
Variation curves of training accuracy
3.3 局限与展望
本文提出的方法在对偏头痛患者的分类识别问题上取得了较好的实验结果，但仍需注意到研究中存在的一些局限性。首先，本研究在基于深度学习的框架中只包含了64名受试者，扩充后为320个样本，数据集依然很小。其次，尽管本研究成功提取了时序特征来提高分类器的性能，但还没有结合BOLD信号的特性来分析偏头痛患者的脑功能变化特征。
在针对偏头痛患者的神经影像学研究中，单个医院较难获得大量样本。研究界最近的趋势是提高神经影像数据的共享水平，这一趋势将会改善未来研究样本集的可及性，多中心样本的获取将有助于确定本文方法是否可以适用于更大的数据样本。此外，有研究指出偏头痛发作可能是对脑能量代谢和氧化应激失衡的反应[23]，进一步探求反映大脑能量代谢状态的时间序列特征与偏头痛疾病机制之间的关联模式或许可以成为更有价值的研究方向。在之后的工作中，基于本文提出的轻量化诊断模型和时序特征可以进一步探索相关疾病的神经鉴别机制，从而开展更具临床应用价值的研究。
4 总结
本研究直接使用rs-fMRI数据的时间序列信号作为特征输入，并进行模型架构设计来实现偏头痛患者组和健康对照组的高精度分类。除此之外，本研究还使用了一系列的优化策略来避免小数据集的脑影像噪声和潜在过度拟合，包括Group ICA和DictLearning相结合的时间序列脑区划分、基于裁剪的数据扩增、L1和L2正则化、dropout机制和交叉验证等。相较于传统模型，本文提出的模型充分考虑了时变特性和模型的预测效率，能快速有效地从rs-fMRI数据中捕获与特征状态相关的大脑动态变化，研究结果表明其在偏头痛分类诊断任务中的可行性和优越性，并可改进后用于其他脑部疾病的分类诊断，具有潜在的临床应用价值。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：孙昂主要负责实验流程、数据记录与分析、论文编写以及算法程序设计；陈宁、何俐主要负责项目主持、协调沟通、理论指导和论文审阅；张俊然主要负责提供实验指导，数据分析指导，论文审阅修订。
伦理声明：本研究通过了四川大学华西医院医学伦理委员会的审批。
0 引言
偏头痛是一种常见的慢性神经系统疾病，影响着全球约15%的人群，其典型症状有反复发作的中重度头痛、恶心、呕吐，以及对光和声音异常敏感等[1]。由于偏头痛的诊断是基于不同症状的组合，需排除可能的继发性头痛病因，目前尚缺乏可以明确诊断的生物学或影像学标志物。考虑到偏头痛疾病本身的特有症状，以及疼痛相关症状的重复发生有可能使得偏头痛患者发生特异性的脑功能改变，最近功能磁共振成像技术和相应数据分析方法的进步，使得从脑功能数据中提取偏头痛的生物标志物成为可能[2]。静息态功能磁共振成像（resting-state functional magnetic resonance imaging，rs-fMRI）是指受试者在不执行思考、运动等任务的情况下，通过检测血氧水平依赖（blood oxygen level dependent，BOLD）信号来表征大脑自发的神经活动的一种成像技术，是评估神经、精神疾病脑功能变化的有效技术手段之一[3]。
近年来，随着深度学习框架及其系列算法在神经影像数据分析处理领域的逐步深入，其在相关疾病分类诊断、疗效预测等方面已展现出巨大的潜力。目前，偏头痛的分类研究大多采用了基于图像的深度学习分析框架，即输入rs-fMRI图像给卷积神经网络（convolutional neural network，CNN）等深层卷积模型，并通过该模型进行特征提取和分类训练。如Yang等[4]提取了rs-fMRI数据的三种特征映射功能图（低频波动幅度、区域同质性和区域功能相关强度），并对特征映射功能图沿z轴和时间轴方向进行切片共生成了5 760张图像，采用了CNN模型其分类准确率最高接近99%。李翔等[5]提出了一种三维CNN无先兆偏头痛基础诊断（three dimensional CNN based diagnosis of migraine without aura，MwoA3D-Net）的智能辅助诊断算法，其最高诊断准确率也达到了98.40%。但是，基于图像输入的CNN通常对样本量和计算资源的要求较高，三维数据和大多数二维图像都比一维时间序列需要更多的计算资源和时耗。与深度学习应用中最常见的自然图像数据集相比，神经影像数据集通常样本量相对较少，大多局限在40～160例[6]；这就导致深层卷积网络等深度学习模型无法充分挖掘图像特征，容易出现过拟合等不足。具体到偏头痛研究领域，其影像数据相较于抑郁症等疾病样本量更少，因此需要新的特征提取方式和模型架构设计等方法层面的创新来解决偏头痛病患这种小样本影像数据的分类识别问题。
在特征提取方式设计上，本文拟直接采用时间序列及其变换作为模型训练输入数据，这样可以大幅减小计算资源开销和训练时耗，并且更适用于小样本数据的学习拟合。另外，直接针对时间序列进行处理也可以避免特征转换过程中可能出现的信息损失[7-9]。为了更好地提取小样本群体的脑区时间序列，本研究拟将组水平独立成分分析（group independent component analysis，Group ICA）方法和字典学习（dictionary learning，DictLearning）这两种数据驱动的方法结合在一起，以识别偏头痛的敏感脑区并研究它们的时序特征。其中，Group ICA可以将多变量信号分离为附加子分量，自动分离出rs-fMRI信号中相互独立的信号源，包括有意义的静息态功能网络和头动等噪声信号[5]；而DictLearning可以通过将整个大脑信号分解为多个各异的时间特征来估计静息态网络的字典，为观察到的样本提供有效的表示[10-12]。在模型架构设计上，为了直接从rs-fMRI时间序列中学习功能特征和大脑状态之间的映射关系，本研究拟采用长短期记忆（long short-term memory，LSTM）网络来捕捉大脑动态时序变化以用于偏头痛分类研究；该网络能够“记住”并加权时间序列或信号的过去值，通过从数据中学习来自适应地捕捉时间依赖性[13-14]。其中，双向长短期记忆（bi-directional long short-term memory，BiLSTM）网络能够挖掘隐藏在时间序列中的上下文信息[15]，与LSTM网络相比，更有机会捕捉到神经元活动的特异信息。此外，以往研究中提出rs-fMRI数据每分钟的BOLD信号对于每个受试者来说可能反映相似的神经元活动[16]，故本文对提取出的时间序列进行了裁剪，以扩充模型的输入序列。
相较于之前的偏头痛数据分类识别框架，本文在三个方面做出了改进：首先，本研究基于数据驱动方法来为小样本数据构建更合理的感兴趣区域（regions of interest，ROIs）划分，并基于ROIs提取时间序列。接着，本研究将所有受试者的输入时间序列裁剪为固定的子序列长度（68 s），扩充模型的输入序列以充分利用整个时间序列蕴含的动态信息。最后，本研究进一步使用BiLSTM网络在两个方向上访问远程上下文，从而捕捉神经元活动前后时间的依赖特性。目前，尚未见基于时间序列特征来进行偏头痛疾病识别的相关深度学习研究，本研究将有助于拓展这一领域的技术边界。
1 材料和方法
本文方法主要从特征提取方式和模型架构设计两个方面进行了改进，模型框架图如图1所示，主要包括rs-fMRI信号分解、提取区域平均时间序列、序列裁剪扩充、BiLSTM网络分析、归一化指数（softmax）函数分类。首先，利用Group ICA和DictLearning对rs-fMRI数据进行信号分解，提取出相应ROIs的区域平均时间序列。然后，将时间序列按时间点分段裁剪，扩充得到5倍样本量的子序列。接着，将子序列送入BiLSTM网络中进行模型训练，学习前向和后向信息。最后通过softmax函数完成二元分类，实现偏头痛患者组和健康对照组的识别与分类。
图1
偏头痛分类预测框架图
Figure1.
Framework of migraine classification and prediction
1.1 数据采集及预处理
1.1.1 受试者
本研究共招募了64名成年人参与rs-fMRI数据采集，如表1所示，分为两组：① 36名偏头痛患者，其中男性11名，女性25名，作为偏头痛患者组。偏头痛患者均是神经科医师根据《国际头痛疾病分类（第2版）》（international classification of headache disorders Ⅱ，ICHD-Ⅱ）[17]诊断标准进行评定入组。② 28名健康受试者，其中男性13名，女性15名，作为健康对照组。健康对照组年龄、性别和受教育年限与偏头痛患者组相匹配。两组受试者均排除了药物或酒精滥用、磁共振成像扫描禁忌、患有其他慢性疼痛疾病和神经发育障碍等情况。
表1
64名受试者的人口统计学信息（）
Table1.
Demographic information of the 64 subjects（）
所有的受试者年龄均在18～50岁之间，且均为右利手。研究中对所有受试者进行了24项汉密尔顿抑郁量表（24-Hamilton depression scale，24-HAMD）和14项汉密尔顿焦虑量表（14-Hamilton anxiety scale，14-HAMA）测试，相应结果如表1所示。四川大学华西医院磁共振研究中心对所有受试者进行了rs-fMRI脑部扫描，其中偏头痛患者在脑部扫描前至少72 h和扫描后48 h无发作。本研究经四川大学华西医院医学伦理委员会批准，并获得每位受试者的书面知情同意书。
1.1.2 rs-fMRI数据采集
本研究64名受试者的rs-fMRI数据通过梯度平面回波成像序列获取，其具体设置如下：重复时间/回波时间 = 2 000 ms/30 ms；翻转角 = 90°；切片厚度 = 5 mm；视野 = 240 mm × 240 mm；体素大小 = 3.75 mm × 3.75 mm × 5 mm；采集矩阵 = 64 × 64。所有受试者的扫描时长均为6 min（即360 s），对应180个时间点。
1.1.3 数据预处理
使用神经影像脑连接组分析软件图论网络分析工具箱（graph theoretical network analysis，GRETNA）对采集数据进行预处理。处理步骤包括：
（1）考虑到初始rs-fMRI信号不稳定和受试者前期的适应时间，去除每个功能数据的前10个时间点，剩余170个时间点（即340 s）用于后续处理；
（2）剩余的数据首先进行时间层校正，然后重新对齐以去除头动伪影的影响。将头动位移超过2.5 mm或者头动旋转角度超过2.5°的受试者剔除。本研究中所有受试者的头动位移和头动旋转角度均未超过上述标准。
（3）然后，使用蒙特利尔神经学研究所（Montreal Neurological Institute，MNI）的模板将得到的图像配准至3 mm体素的标准空间。
（4）最后对每个体素时间序列进行带通滤波（0.01～0.08 Hz），以减少低频漂移和高频生理噪声的影响。
1.2 rs-fMRI时间序列特征提取
预处理后的图像通常会进行脑区划分以降低后续分析的数据维度。脑网络划分的常用方法：一种是使用自动解剖标记（automated anatomical labeling，AAL）等固定图谱；一种是通过数据驱动的方式来确定大脑ROIs图谱。基于大脑固定图谱的方法在获取脑区平均时间序列时易受主观因素影响，可能存在一定的偏差，本文中选择两种数据驱动的方法来获取所需的ROIs时间序列。
考虑到本文的小样本群体，基于Group ICA的方法能够将模型应用于组数据[18]，在群体水平上对大脑进行分割，寻找在各个受试者之间较为一致的公共成分[19-20]。Group ICA方法中，先设Y是T×M的观测信号，即预处理后的rs-fMRI图像数据矩阵。其中T表示时间点的个数，M表示受试者大脑中体素的总数。假设每个体素的rs-fMRI信号是各个成分信号的线性混合，则Y可以分离成相互统计独立的成分，记P为N个分离的成分图，W是混合矩阵，其计算公式如式（1）所示：
|
对于每个受试者的成分图空间分布未知的情况下，可利用盲源分离方法求取，记Ak为分离矩阵，它的列表示相应成分图的时间序列，k表示迭代次数，其计算公式如式（2）所示：
|
Group ICA的成分估计过程需要进行跨个体的独立成分的融合，才能反映出一组rs-fMRI数据的整体特性。分离矩阵Ak需要在多个受试者中逐次迭代，直至算法收敛，从而得到各受试者一致的成分映射。
由于Group ICA需要遵守数据表示子空间上的正交性约束，这可能导致独立源数量受限于信号维度。针对这一不足，DictLearning可以通过构建一个原子或子空间的字典，为观察到的样本提供有效的稀疏表示，并有可能推导出稀疏信号的先验未知统计量[21-22]。同时，相关神经科学研究表明大脑神经活动中存在稀疏响应[11]，稀疏表示方法的内在性质可能与大脑神经活动的稀疏响应相契合。
在DictLearning方法中，给定rs-fMRI信号矩阵Y，维度为T×M，其中T是rs-fMRI时间序列点的总数，M是受试者大脑中体素的总数。Y中的rs-fMRI信号被建模为学习到的基础字典D，模型矩阵形式如式（3）所示：
|
其中A保存稀疏表示的系数矩阵。具体来说，D中每个字典原子的信号代表特定脑网络的功能活动，而系数矩阵A中的向量代表对应脑网络的空间分布。字典D可用于划分ROIs，并提取相应的时间序列。
1.3 特征融合与数据扩增
Group ICA方法参照惯例在此提取了20个独立成分，每个成分涉及1～5个脑区，共计将大脑划分为58个ROIs。同样，DictLearning方法依据字典提取了20个原子，每个原子涉及2～8个脑区，共计将大脑划分为96个ROIs。两种方法在脑区划分中各有其优势和特点，本文结合了两种方法获得的154个ROIs并计算了各个ROIs的平均时间序列，生成维度为 (x, y, z)=(64, 170, 154) 的特征矩阵。其中x是样本数，y是时间点，z是ROIs数。接着，对每个ROIs时间序列应用了从-1～1的线性归一化。
对于每个受试者170个时间点（即340 s）的初始时间序列可被平均拆分成5段长度为34个时间点的子序列，每段代表68 s的成像时间。如图1所示，序列裁剪扩充后每一段68 s的子序列都可以作为模型的输入数据，这使模型数据集扩充到初始的5倍，总共达到320个输入序列。重构后的特征矩阵维度为 (320, 34, 154)，这将作为接下来BiLSTM分类网络的输入。
1.4 BiLSTM网络分类
本研究为了进一步学习每个时间序列内部的前后时序信息来刻画周期性大脑状态变化，在对rs-fMRI时间序列裁剪扩充的基础上采用了BiLSTM网络。如图1所示，输入序列 经过LSTM层，学习前向和后向信息，为每个时间点返回相应的值 。接着使用带有S型生长曲线（sigmoid）激活函数的密集层和带有softmax激活函数的第二个密集层作为二元分类器。模型使用二元交叉熵损失函数和具有默认参数值的自适应矩估计优化器进行训练。为了避免过拟合，本研究在所有隐藏层中引入随机失活（dropout）机制，随机失活50%神经元，在LSTM层中加入了正则化项L1和L2。实验设定了300次迭代训练，提前停止条件为验证损失连续在30次训练中小于0.01。此外，本研究通过将数据分为8份来执行交叉验证方法，以便受试者比例在所有测试中都大致相同。
2 结果
为了评估性能，实验中进行了8折交叉验证，计算了灵敏度、特异性和分类精度。如表2所示，Group ICA结合DictLearning的方法的敏感性、特异性和准确性均优于单一Group ICA或单一DictLearning。同时，本文绘制了受试者工作特性（receiver operating characteristic，ROC）曲线并计算了曲线下面积（area under curve，AUC），如图2所示。其中，曲线的横坐标为假阳性率（false positive rate，FPR），纵坐标为真阳性率（true positive rate，TPR）。ROC曲线所覆盖的区域面积AUC可以用来定量评价偏头痛患者组和健康对照组分类的准确性，ROC曲线越靠近平面的左上角，AUC的值越高，表明分类器分类效果越好。
表2
不同时间序列提取方法的分类结果
Table2.
Classification results of different time-series extraction methods
图2
受试者工作特征曲线
Figure2.
Operating characteristic curves of subjects
本文方法的平均AUC达到了0.98，并且显著优于传统的AAL固定图谱、单一Group ICA提取图谱和单一DictLearning提取图谱的结果。AAL的效果最差，这可能是由于样本特性与模板不适配，时间序列中存在噪声造成的。相比之下，本文方法结合了两种数据驱动的ROIs划分方法，可以根据数据特性学习内在的大脑状态转换。在8折交叉验证（ROC fold1～ROC fold8）下，每一折的AUC都在0.98～1.00范围内，并且在几种方法中数值波动最小，表现出对时间序列噪声良好的鲁棒性。
如表3所示，使用Group ICA结合DictLearning提取特征后，分类网络和输出单元数量的不同影响了分类精度。结果表明BiLSTM显著优于单向的LSTM，并且具有32个隐藏节点的BiLSTM准确率为96.94%，达到了最高的分类精度。
表3
LSTM与BiLSTM分类网络结果对比
Table3.
Comparison results of LSTM and BiLSTM classification networks
3 讨论
3.1 特征提取方式的对比分析
考虑到临床预测模型对部署难易度和计算效率的客观要求，本研究确立了以一维时间序列为核心的分析模式。传统rs-fMRI数据的时间序列分析，是通过皮尔逊相关等功能连接性度量来构建脑功能连接网络，以生成鉴别性特征用于后续分类。但这些特征构建方法以及相应的分类器没能充分考虑时间依赖性[7]，可能会影响模型解码性能。本研究前期在偏头痛患者的数据上做了预实验，通过皮尔逊相关构建了功能连接矩阵作为分类特征。与本文方法的最高分类精度准确率达96.94%相比，以功能连接矩阵作为分类特征准确率可达88%，则不算理想。Dvornek等[14]也在自闭症患者数据上使用时间序列特征与功能连接度量方法做了对比，时间序列特征同样展现了较优的分类性能。
在选择大脑ROIs图谱以提取时间序列过程中，本研究结合了两种数据驱动方法来划分ROIs，充分利用了体素间的空间联系，从而降低群体差异带来的影响。常用的固定脑图谱模板可能由于对大脑的划分精细度不同，或构建的群体基础不同，而对小样本数据的分析结果有较大的影响。如表2所示，Group ICA（58 ROIs）和DictLearning（96 ROIs）的指标均优于固定图谱AAL（116 ROIs），即使它们的脑区划分精细度并没有AAL高。这源于Group ICA具有良好的控制群体模型，可以分解出大脑信号中的背景噪声；而DictLearning则可以提取到自然稀疏的脑区划分，捕捉到数据的内在特征。Group ICA加DictLearning的方式结合了两种方法的优点，以结合方法来划分区域提取时序特征可以达到更好的分类精度。
3.2 不同模型框架的影响
LSTM是一种特殊类型的递归神经网络，内含的时间循环结构可以较好地刻画具有时间关联的序列数据，并学习到受试者大脑动态的模式信息。考虑到大脑可能会不断地使用上下文信息来指导更高层次的认知功能，而不是在时间序列结束时产生具有严格方向的输出，本文中使用了BiLSTM和LSTM进行了对比分析。
如表3所示，在相同隐藏节点设置下，BiLSTM准确率一直优于LSTM，表明了双向信息更有利于表征大脑信号。此外，隐藏节点的数量可能直接影响LSTM、BiLSTM网络的学习能力。本文比较了具有不同数量（8、16、32、64、128、154）隐藏节点输出单元模型的性能。在单层模型中，具有32个隐藏节点的BiLSTM达到了最高的分类精度。具有8个或16个隐藏节点的模型其性能有所降低，很可能是相应的隐藏节点输出单元太有限，无法存储足够的序列信息。具有64个以上隐藏节点的模型性能也略有下降，其原因可能是过拟合。如图3所示，在154个节点的参数设置下，训练准确度和验证准确度两条曲线整体波动较大，相距较远，说明模型的方差较大，模型有过拟合趋势。
图3
训练验证准确率变化曲线
Figure3.
Variation curves of training accuracy
3.3 局限与展望
本文提出的方法在对偏头痛患者的分类识别问题上取得了较好的实验结果，但仍需注意到研究中存在的一些局限性。首先，本研究在基于深度学习的框架中只包含了64名受试者，扩充后为320个样本，数据集依然很小。其次，尽管本研究成功提取了时序特征来提高分类器的性能，但还没有结合BOLD信号的特性来分析偏头痛患者的脑功能变化特征。
在针对偏头痛患者的神经影像学研究中，单个医院较难获得大量样本。研究界最近的趋势是提高神经影像数据的共享水平，这一趋势将会改善未来研究样本集的可及性，多中心样本的获取将有助于确定本文方法是否可以适用于更大的数据样本。此外，有研究指出偏头痛发作可能是对脑能量代谢和氧化应激失衡的反应[23]，进一步探求反映大脑能量代谢状态的时间序列特征与偏头痛疾病机制之间的关联模式或许可以成为更有价值的研究方向。在之后的工作中，基于本文提出的轻量化诊断模型和时序特征可以进一步探索相关疾病的神经鉴别机制，从而开展更具临床应用价值的研究。
4 总结
本研究直接使用rs-fMRI数据的时间序列信号作为特征输入，并进行模型架构设计来实现偏头痛患者组和健康对照组的高精度分类。除此之外，本研究还使用了一系列的优化策略来避免小数据集的脑影像噪声和潜在过度拟合，包括Group ICA和DictLearning相结合的时间序列脑区划分、基于裁剪的数据扩增、L1和L2正则化、dropout机制和交叉验证等。相较于传统模型，本文提出的模型充分考虑了时变特性和模型的预测效率，能快速有效地从rs-fMRI数据中捕获与特征状态相关的大脑动态变化，研究结果表明其在偏头痛分类诊断任务中的可行性和优越性，并可改进后用于其他脑部疾病的分类诊断，具有潜在的临床应用价值。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：孙昂主要负责实验流程、数据记录与分析、论文编写以及算法程序设计；陈宁、何俐主要负责项目主持、协调沟通、理论指导和论文审阅；张俊然主要负责提供实验指导，数据分析指导，论文审阅修订。
伦理声明：本研究通过了四川大学华西医院医学伦理委员会的审批。
