提取偏頭痛患者等的神經影像特征并進行識別模型的設計對相關疾病的輔助診斷具有重要意義。相較于常用的影像特征，本研究直接采用時間序列信號表征偏頭痛患者組和健康對照組的大腦功能狀態，可有效利用時間信息并減小分類模型訓練計算量。首先，本研究針對小樣本群體運用組水平獨立成分分析和字典學習劃分不同腦區后，提取區域平均時間序列信號；其次，將提取的時間序列平均劃分成多個子時間序列，以擴充模型輸入樣本；最后，使用雙向長短期記憶網絡對時間序列建模，學習每個時間序列內部的前后時序信息來刻畫周期性大腦狀態變化以提高偏頭痛的診斷準確率。研究結果顯示，偏頭痛患者組與健康對照組的分類準確率為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數據中捕獲與特征狀態相關的大腦動態變化，研究結果表明其在偏頭痛分類診斷任務中的可行性和優越性，并可改進后用于其他腦部疾病的分類診斷，具有潛在的臨床應用價值。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：孫昂主要負責實驗流程、數據記錄與分析、論文編寫以及算法程序設計；陳寧、何俐主要負責項目主持、協調溝通、理論指導和論文審閱；張俊然主要負責提供實驗指導，數據分析指導，論文審閱修訂。
倫理聲明：本研究通過了四川大學華西醫院醫學倫理委員會的審批。
