李 翔,魏本征,吳宏赟,李徐周,洪雁飛,叢金玉
1.山東中醫藥大學 智能與信息工程學院,濟南250355
2.山東中醫藥大學 醫學人工智能研究中心,山東 青島266112
3.山東中醫藥大學 青島中醫藥科學院,山東 青島266112
4.山東中醫藥大學附屬醫院 腦病科,濟南250014
5.山東青年政治學院 信息工程學院,濟南250103
偏頭痛是一種反復發作的、可致殘的原發性頭痛,已成為嚴重危害人類健康的腦部疾病之一。其臨床表現為頭部一側或兩側搏動性的劇烈頭痛,對外界的運動、視覺等刺激敏感,可伴有畏光、嘔吐等癥狀[1]。2015年世界衛生組織的一項研究結果表明,偏頭痛已經成為世界第三大流行疾病,也是全球第六大致殘疾病,影響了約15%的世界人口[2]。在中國,偏頭痛的發病率約為9.3%[3]。除高發病率外,偏頭痛患者在頭痛發作時還會出現暴躁、易怒等現象,嚴重干擾患者的正常生活[4]。臨床中常見的偏頭痛類型以無先兆偏頭痛為主,其診斷主要由臨床醫生通過患者的臨床表現,根據《國際頭痛疾病分類(第3 版)》(International Classification of Headache Disorders,ICHD-III)[5]進行。有效和精確的診斷是無先兆偏頭痛治療至關重要的一步,由于缺乏明顯的前驅癥狀,傳統的診斷方法在鑒別無先兆偏頭痛與其他原發性或繼發性頭痛時充滿了挑戰[6]。因此,為提高無先兆偏頭痛的診斷準確率,降低誤診率,設計自動化的臨床輔助診斷系統具有重要的臨床應用價值。
近年來,隨著MRI設備的快速更新以及神經影像分析方法的迅速發展,研究人員發現實現某一功能的不同腦區存在一致波動的低頻信號,這些具有一致波動低頻信號的不同腦區構成了靜息態腦網絡(resting-state brain network,RSN)[7]。而依靠血氧水平依賴信號來反映神經元活動的靜息態功能磁共振(resting-state functional MRI,Rs-fMRI),尤其適合用于觀察無先兆偏頭痛RSN的異常情況。RSN 的異常可反映出患者大腦功能區的功能連接異常。當前,功能連接的分析主要有基于模型驅動的方法和基于數據驅動的方法[8],這些方法可從不同角度探究無先兆偏頭痛患者不同大腦區域或組織的功能連接異常。雖然無先兆偏頭痛的發病原因和病理機制尚不明確,但已有研究結果初步表明無先兆偏頭痛患者某些與疼痛有關的RSN 存在異常,如默認模式網絡(default mode network,DMN)[9]、執行控制網絡(executive control network,ECN)[10]、視覺網絡(visual network,VIN)[11]、感覺運動網絡(sensorimotor network,SEN)[12]、右側額頂網絡(right frontoparietal network,RFPN)[13]和左側額頂網絡(left frontoparietal network,LFPN)[14]等。這些異常的RSN 可作為潛在的無先兆偏頭痛影像學生物標志物用于個體化診斷。
當前,無先兆偏頭痛輔助診斷已經取得了一定的進展,但其智能化程度仍然較低。一些研究者從RsfMRI數據中提取特征,訓練支持向量機(support vector machines,SVM)等傳統機器學習分類器,取得了良好的診斷結果。根據Rs-fMRI數據特征提取方法的不同,可分為基于功能連接的特征提取方法以及基于體素和功能連接融合的特征提取方法。基于功能連接的特征提取方法主要有:Chong等人[15]以疼痛相關的33個大腦區域為種子點,計算該區域的功能連接特征,特征降維后使用對角二次判別分析算法構建無先兆偏頭痛和健康對照的輔助診斷模型,診斷準確率為86.1%。Tu等人[16]利用Dosenbach模板將大腦分為160個腦區并計算所選腦區之間的功能連接特征,經過特征篩選后訓練SVM算法,最終診斷準確率達到91.4%。該研究同時還發現VIN、DMN 和SEN 等可作為診斷無先兆偏頭痛的神經標志物?;隗w素和功能連接融合的特征提取方法主要有:Zhang等人[17]利用自動解剖標記模板將被試的大腦劃分為一系列感興趣區域(region of interests,ROIs),計算ROIs 區域基于體素的低頻振幅、局部一致性和局部功能連接特征。該研究還從多模態角度計算結構MRI的白質特征,將上述特征融合后訓練SVM算法,診斷準確率為84.0%。
以上對功能連接特征的提取主要是基于模型驅動的分析方法,該方法對特征的提取依賴預先定義的腦圖譜模板,不同的腦圖譜模板對大腦的劃分精細度不同,且ROIs 腦區的選擇需要豐富的先驗知識,可能影響后續的分析結果。此外,不同的特征提取方法得到的特征仍存在維度過高或冗余等問題,在訓練傳統分類器之前仍需再次進行篩選。
在基于數據驅動的功能連接分析方法中,獨立成分分析(independent component analysis,ICA)方法無需任何先驗模型即可自動從Rs-fMRI 數據中分離出有意義的RSN 和頭動等噪聲,目前已經廣泛應用于RSN 分析,為研究病理機制提供新的視角[18]。然而將ICA應用于多被試Rs-fMRI數據分析時,其輸出成分具有無序性和不可預估性等缺點,難以在不同被試間建立對應性,因此對于多被試Rs-fMRI數據分析常用組ICA方法[19]。常用的組ICA 方法有Back-Reconstruction[20]方法、Dual Regression[21]方法和組信息指導的獨立成分分析(group information guided ICA,GIG-ICA)[22]方法等。使用Back-Reconstruction 方法和Dual Regression 方法恢復個體成分時,均無法保證個體被試不同成分間的獨立性,且Back-Reconstruction 方法無法使用已有RSN 模板對新數據進行分析。而GIG-ICA 方法不但可以生成具有更強個體間獨立性、組間對應性以及更高準確性的RSN,還可將已有RSN模板作為先驗信息指導新被試的RSN生成,新被試生成的RSN 和先驗信息具有較強的對應性。因此,該方法是一種提高無先兆偏頭痛臨床診斷準確率的有效技術方法。
近年來,各種深度學習技術在Rs-fMRI圖像處理及精神疾病分析研究中取得了重要進展[23-24],其中卷積神經網絡(convolutional neural networks,CNN)因能自動創建具有分層表示的多級非線性模型且能很好地捕捉圖像數據中的空間結構信息獲得廣泛關注及應用[25]。Yang等[26]首次使用2D-CNN技術構造了無先兆偏頭痛、健康對照以及有先兆偏頭痛的診斷模型,但該方法將生成的3D特征轉化為2D切片特征,由于2D-CNN僅在兩個維度提取特征,該方法丟失了豐富的第三維度信息。與2D-CNN相比,3D-CNN可以從數據中直接提取三維特征,并對時空信息進行建模,通過組合不同的通道信息獲得最終的高維特征表示[27-28]。近期在精神分裂癥和阿爾茲海默癥上的工作也證明了3D-CNN 處理RSN 數據的優勢[29-30]。
為提高無先兆偏頭痛的診斷準確率及智能化程度,基于設計的新型3D-CNN技術,本文提出了一種無先兆偏頭痛智能輔助診斷算法MwoA3D-Net(3D convolutional neural network based diagnosis of migraine without aura)。該算法將GIG-ICA方法引入到無先兆偏頭痛的Rs-fMRI 功能連接分析中,用于生成被試的RSN,避免因腦圖譜模板不同導致的結果差異。此外,該算法可直接學習RSN的3D空間結構特征,且在算法設計時加入一系列針對醫學影像小樣本過擬合問題的優化策略,進一步提高無先兆偏頭痛的診斷準確率。
本研究在山東中醫藥大學附屬醫院頭痛門診招募60名無先兆偏頭痛患者,所有患者均被腦病科醫生根據ICHD-III 標準確診。同時選擇與無先兆偏頭痛組患者性別、年齡相匹配的65 名健康被試作為對照組。被試的人口學信息如表1所示。

表1 實驗被試人口學統計表Table 1 Demographic of all subjects
無先兆偏頭痛組納入標準:(1)符合ICHD-III 臨床診斷標準;(2)年齡在20~30歲之間;(3)掃描前3天內沒有發作且未服用相關藥物;(4)無MRI 掃描禁忌癥;(5)無腦部器質性疾??;(6)右利手;(7)知情同意者。
無先兆偏頭痛組排除標準:(1)年齡不在20~30 歲之間;(2)患有嚴重身體疾病或除無先兆偏頭痛以外其他神經疾??;(3)有藥物濫用史;(4)妊娠或哺乳期婦女;(5)有MRI掃描禁忌癥。
正常對照組納入標準:(1)年齡在20~30 歲之間;(2)無偏頭痛家族史、慢性疼痛及其他神經或精神類疾??;(3)無糖尿病、心臟病、高血壓以及其他慢性全身性疾??;(4)無認知障礙;(5)無藥物濫用史;(6)右利手;(7)無MRI掃描禁忌癥。
所有被試的影像數據均在Philips Achieva 3.0T 掃描儀上進行采集。掃描前,所有被試頭部均用海綿墊固定以減少頭部位移。在掃描過程中,要求受試者仰臥在MRI 設備上,保持閉眼靜息狀態,均勻呼吸,避免思考,但需保持清醒。Rs-fMRI 數據采用回聲平面成像(echo-planar imaging,EPI)序列獲取,參數設置如下:TR=3 000 ms,TE=35 ms,FA=90°,matrix=128×128,FOV=230 mm×230 mm。T1WI 結構像掃描參數設置如下:TR=8.0 ms,TE=3.8 ms,FOV=230 mm×230 mm,matrix=512×512,FA=12°。
本研究獲得山東中醫藥大學附屬醫院倫理委員會批準,所有被試在數據采集前均簽署書面知情同意書。
為減少因Rs-fMRI數據采集、被試生理學特征以及個體化差異帶來的誤差,影響后續的分析結果,本文采用SPM12軟件(https://www.fil.ion.ucl.ac.uk/spm)和DPARSF工具箱(http://www.rfmri.org/DPARSF)完成Rs-fMRI 數據的預處理工作,預處理步驟如下。
步驟1時間校正:去除被試圖像前10 個時間點的數據,并以第29層圖像數據作為參考進行時間校正。
步驟2頭動校正:計算被試在掃描過程中的頭動位移和頭動旋轉角度,為減少因頭動對被試數據采集的影響,將頭動位移超過2.5 mm或者頭動旋轉角度超過2.5°的被試剔除。本研究中所有被試的頭動位移和頭動旋轉角度均未超過上述標準。
步驟3空間標準化:將每個被試頭動校正后的圖像標準化到每個被試各自的結構像,隨后將被試的功能像標準化到統一的模板,并重采樣體素大小為2 mm 的立方體。
步驟4空間平滑:為減少機器不穩定或者生理運動所產生的干擾信號,對空間標準化后的圖像進行4 mm半寬全高高斯空間平滑。
對預處理后的被試數據,本文首先采用靜息態腦網絡生成模塊生成與無先兆偏頭痛相關的RSN,然后將得到的RSN用于訓練MwoA3D-Net算法完成無先兆偏頭痛的輔助診斷,具體流程如圖1所示。

圖1 MwoA3D-Net算法框架示意圖Fig.1 Framework of MwoA3D-Net algorithm
從有噪聲的Rs-fMRI 數據中精準地生成所需的RSN,是提高無先兆偏頭痛診斷準確率的關鍵步驟之一。GIG-ICA方法的計算主要包括以下兩個步驟:
(1)對所有被試數據進行組水平ICA,得到一系列組獨立成分,此步驟可表示為:

其中,X=[X1;X2;…;XN]為N個被試的功能連接矩陣;S=[S1;S2;…;SM]是估計的M個組水平的獨立成分。
(2)將組獨立成分作為參考信息輸入到基于多目標函數帶參考信號的ICA算法(ICA with reference,ICA-R),計算被試的獨立成分以及對應的時間序列。其多目標優化函數可表示為[22]:

式中,表示第k個被試的第l個獨立成分;v是高斯隨機變量;G是非二次函數;是對應的解混合矩陣;J()是的負熵;F()代表和的相似性。
最后,求解式(2)可得到:

GIG-ICA方法由GIFT工具(https://trendscenter.org)實現。基于Smith 等人[31]以大量健康被試制作的RSN模板,本文選取了與無先兆偏頭痛密切相關的8個RSN作為先驗信息,用于指導新被試的RSN 生成。本文所選RSN如圖2所示,分別為VIN、DMN、SEN、ECN、LFPN和RFPN。將預處理后的Rs-fMRI 數據以及所選RSN模板同時輸入GIG-ICA 方法,生成被試的3D 空間成分以及對應的1D 時間序列。在上述計算過程中,未設置任何閾值,因此可以完整保留所有被試的3D 空間成分信息,3D空間成分代表被試的RSN。最終,每個被試生成的8個3D RSN用于訓練MwoA3D-Net算法。

圖2 本文選取的8個RSN模板示意圖(顯示閾值為3~10)Fig.2 Eight RSN templates selected in this paper(display threshold:3~10)
MwoA3D-Net算法由RSN-Net和全連接層構成,如圖1 所示。該算法的訓練分兩階段進行:在第一階段,針對每個被試的8個RSN,本文分別訓練了8個RSN-Net完成對每個被試RSN的特征提取,保存8個RSN-Net最優模型并統計每個RSN 的診斷結果。在第二階段,MwoA3D-Net 算法加載8 個RSN-Net 最優模型,以被試所有的RSN作為輸入,進行前向傳播并計算損失,最后使用反向傳播對MwoA3D-Net算法的全連接層進行訓練。
2.2.1 RSN-Net結構
8 個RSN-Net 的網絡結構均相同,由4 個卷積層、4個批正則化層、3 個池化層、1 個全連接層和1 個輸出層構成。模塊的主要結構和功能如下:
(1)卷積層:卷積層將低層次的特征逐層地提取和組合形成更高級復雜的抽象特征,避免了傳統機器學習中復雜的手工提取過程。4 個卷積層分別包含16、16、32和64個卷積核,每個卷積層的卷積核大小為3×3×3,步長為1,并由Kaiming 均勻分布初始化。本文在卷積層與非線性激活函數之間增加了3D 批正則化層,可對每個卷積層生成的特征圖進行批次歸一化操作,以加快模型訓練時的收斂速度。
由于RSN 的激活區域在整個大腦中占比相對較高,可用較大的感受野去學習其復雜的空間特征。雖然大尺寸的卷積核可以帶來更大的感受野,但直接構造大尺度的卷積核可導致模型的計算量呈指數級增長,不利于深層模型的訓練。因此本文使用2個3×3×3的卷積核串聯,使得卷積層具有5×5×5 的感受野,可有效提升模型對RSN空間結構特征的提取和表達能力。相比直接使用5×5×5 卷積核,此方法還可有效降低算法的參數量[32]。此外,兩層卷積之間多使用一個非線性激活函數,使得模型增加了特征的非線性表達能力。本文卷積的計算公式為:

式中,l代表卷積層數;mj代表輸入層的感受野;是該層的神經元偏置項;g(?) 為非線性激活函數;h代表卷積核;是該層神經元j與前一層神經元i的連接強度。
(2)池化層:池化層的目的主要是完成卷積層特征圖的降維。本文的3 個池化層使用最大池化函數且池化核尺寸及步長均為2×2×2。
(3)全連接層:在分類任務中,全連接層通常在池化層之后起到整合局部信息的作用。RSN-Net 的全連接層擁有128個神經元,每一個神經元都與前面的池化層連接,最終生成128維的特征向量。
(4)輸出層:RSN-Net的輸出層擁有2個神經元用于生成類別信息,采用Softmax 分類器對模型輸出的結果的類別屬性進行概率量化。Softmax是邏輯回歸分類器在多類別領域的推廣,采用如下方式計算類別概率:

式中,gi代表第i個類別,C表示分類的類別個數。
2.2.2 全連接層結構
本文將RSN-Net輸出的診斷信息進行融合后,訓練MwoA3D-Net 算法的全連接層,用于完成輔助診斷任務。MwoA3D-Net 算法的全連接層擁有128 個神經元以及2 個輸出神經元,同樣采用Softmax 分類器對模型輸出結果的類別屬性進行概率判斷。
2.2.3 MwoA3D-Net算法優化策略
針對深度學習中醫學影像小數據的過擬合問題,本文采用如下措施解決:(1)對輸入的訓練數據進行3D數據增強使其產生更多訓練圖片供模型學習,主要包括隨機概率的水平反轉、高斯模糊、旋轉和縮放操作。對驗證和測試數據則不做數據增強。(2)使用“早停法”,即在算法訓練時監控驗證集的損失變化情況,如果在一定周期內損失值不下降則模型停止訓練。(3)使用L1 和L2正則化技術,即在損失函數的基礎上添加懲罰項,迫使算法學習較小的權重。(4)使用批正則化技術,即對所有卷積層生成的特征圖進行批次歸一化操作,以解決梯度消失和爆炸的問題,在一定程度上也能減緩過擬合。
MwoA3D-Net 算法采用PyTorch 1.4[33]深度學習框架編寫。算法在Intel 8268 CPU、NVIDIA Tesla V100 32 GB GPU 上進行訓練和測試。算法學習率設置為0.001,使用隨機梯度下降(stochastic gradient descent,SGD)作為優化器,動量參數為0.9,使用交叉熵作損失函數。
為驗證本文所提算法的有效性,將其與多個不同類型的分類器做性能對比實驗。首先與傳統機器學習方法構造的分類器對比。本文應用GIG-ICA 方法為每個被試生成8個RSN后,使用雙樣本T檢驗獲取RSN的個體差異,經PCA 算法降維后采用SVM 算法作為分類器完成輔助診斷。其次與深度學習方法中的2D-CNN 以及3D-CNN算法模型進行對比。本文對比的2D-CNN算法有AlexNet和CNN with Inception網絡[26],3D-CNN算法有MB-CNN[30]以及3D ResNet18[34]。其中,MB-CNN主要用于解決阿爾茨海默癥和輕度認知障礙的診斷問題,3D ResNet18 使用的3D 卷積可有效學習時空的三維特征信息。
為進一步驗證所提算法的設計優勢,本文做了不同的消融實驗。一是,研究驗證使用4種不同的過擬合優化策略對診斷結果的影響。二是,針對不同的池化函數做了對比實驗和分析。三是,對RSN 的診斷準確率做了比較研究。
為了測試算法的魯棒性及泛化性能,基于5折交叉驗證方法,本文將所有數據隨機平均分為5 份,每次取其中1 份數據作測試集,剩余4 份數據作訓練集。同時本文對4 份訓練集再隨機均分為5 份,取其中的1 份作驗證集,剩余4份作訓練集。最終,每折數據訓練集、驗證集、測試集的比例為16∶4∶5。
本文采用準確率(accuracy,ACC)、靈敏度(sensitivity,SEN)和特異性(specific,SPEC)指標評估本文所提算法的性能。上述指標計算公式分別為:

式中,TP為真陽性總數;FP為假陽性總數;TN為真陰性總數;FN為假陰性總數。
3.4.1 無先兆偏頭痛輔助診斷性能測試
MwoA3D-Net 算法以及其他對比算法在無先兆偏頭痛輔助診斷任務上的實驗結果,如表2所示。根據實驗結果可知,MwoA3D-Net 算法具有最高的診斷準確率,為98.40%,在靈敏度和特異性指標上也顯示出優越性。MwoA3D-Net算法的準確率比SVM算法高6.40個百分點,比AlexNet 和CNN with Inception 等2D-CNN分別高5.36 個百分點和2.24 個百分點,與Kam 等人提出的MB-CNN 相比高17.20 個百分點,與3D 算法中常用的3D ResNet18相比高18.40個百分點。

表2 算法模型性能比較結果統計表Table 2 Performance comparison of different methods %
由于直接將分類器應用于原始Rs-fMRI數據時,存在Rs-fMRI 數據噪音信息大以及Rs-fMRI 空間復雜性高等問題,當前的研究方法都是從Rs-fMRI數據中提取功能連接等特征,然后構造分類器以實現無先兆偏頭痛的輔助診斷。其中,計算功能連接特征多依賴于預先定義的腦圖譜模板,需要對大腦的區域進行手動選擇。這種人為選擇腦區的處理方法,易受主觀因素影響,后續研究結果常會出現較大偏差。為克服此問題,本文采用一種較為客觀的方式,即使用GIG-ICA方法計算被試的功能連接特征。
為克服訓練時由數據量少引發的過擬合現象,本文采用4 種過擬合處理技術以提高MwoA3D-Net 算法性能。根據實驗結果可知,MB-CNN雖能完成3D RSN的分類任務,但該算法未采取任何解決過擬合問題的措施,故其在本研究中得到的實驗結果很不理想。3D ResNet18 雖有殘差結構保證了網絡可以增加深度來改善圖像分類性能,但受限于RSN 的圖像尺寸小以及數據量少等因素影響,該算法的實驗結果也較差。將生成的3D RSN數據轉換為2D切片數據后,隨著訓練數據量的激增以及采用過擬合優化策略,2D-CNN,如AlexNet和CNN with Inception 網絡,取得了比MB-CNN 和3D ResNet18高的結果。
3.4.2 過擬合優化策略測試實驗
本文對4種過擬合優化策略做了實驗研究,實驗結果如表3 所示,分別是批正則化、3D 數據增強、L1 和L2正則化以及“早停法”。該實驗結果表明,4種不同策略均具有一定的有效性。其中,引入“早停法”的優化策略對實驗結果的提升較大?!霸缤7ā笔窃O置一定的次數閾值,如果模型的統計指標在閾值次數內仍然未改善,則直接停止模型的訓練,該策略優點是在較差情況下也可獲得較優的訓練模型。其他的過擬合優化策略,關注于改進模型的訓練過程,對實驗結果的提升較小。

表3 過擬合優化策略性能測試結果統計表Table 3 Performance comparison of different overfitting optimization methods%
池化函數對MwoA3D-Net 算法性能影響的實驗結果如圖3所示。從圖3中可以看出,基于最大池化函數的池化層在3個指標上,均較基于平均池化函數的池化層有所提高。分析其原因,是由于最大池化函數可保留特征中的最大值,特征中的最大值可能是異常的功能連接強度,因此使用最大池化函數會提升模型的性能。而平均池化函數對異常的功能連接特征進行了平均,雖保留了特征中的大量信息,但是噪聲也較大,導致這種大量的特征信息對模型性能的提升較弱。

圖3 池化層函數性能測試結果示意圖Fig.3 Performance comparison of different pooling operators
3.4.4 RSN輔助診斷無先兆偏頭痛貢獻率比較研究
本文同時對8 個RSN 的診斷準確率進行了對比分析,結果如圖4所示。從圖中可以得出,本文所選的RSN均具有較強的辨別性,其中診斷貢獻率最高的RSN 為LFPN。LFPN包括內額額葉和尾狀核等區域,主要負責語言的處理[35],但該網絡在無先兆偏頭痛患者中研究較少,尚未引起研究者的重視。優選后的LFPN等RSN可作為潛在的影像學生物標志物,用于無先兆偏頭痛的個體化診斷。

圖4 RSN診斷準確率比較結果示意圖Fig.4 Comparison of accuracy of different RSN
為提高無先兆偏頭痛的診斷準確率,本文設計了一個基于改進3D-CNN 的無先兆偏頭痛智能輔助診斷算法MwoA3D-Net。該算法采用GIG-ICA 方法生成被試的RSN,可有效避免因腦圖譜模板不同而導致的結果差異。此外,該算法可直接學習RSN 的3D 空間結構特征。針對醫學影像小樣本過擬合問題,本文也利用一系列的優化策略進行解決。研究結果表明本文算法的診斷準確率以及智能化水平較高,具有潛在的臨床應用價值。雖然MwoA3D-Net算法取得了較好的效果,但仍存在模型參數量較大,提取的功能連接特征缺乏對瞬時功能變化的捕捉等問題。在后續的工作中可以將動態腦網絡與RSN結合,設計更優的深度學習算法,進一步深入研究LFPN等RSN,提高無先兆偏頭痛的輔助診斷性能。