石玉虎 曾衛(wèi)明 鄧 金 王倪傳

1(上海海事大學信息工程學院，上海 201306) 2(淮海工學院計算機工程學院, 江蘇 連云港 222023)

引言

利用功能磁共振成像(functional magnetic resonance imaging, fMRI)技術(shù)來研究人腦的功能連通性及相關(guān)的認知活動特性，可以更好地揭示人腦相關(guān)的認知活動規(guī)律，促進fMRI技術(shù)在神經(jīng)信息處理以及神經(jīng)認知應用等領(lǐng)域發(fā)揮重要作用。目前，關(guān)于fMRI腦功能連通性的研究方法主要分為基于模型驅(qū)動的方法和基于數(shù)據(jù)驅(qū)動的方法。其中，基于模型驅(qū)動的方法主要有互相關(guān)分析方法[1]、一致性分析方法[2]以及統(tǒng)計參數(shù)圖分析方法[3]等，但這些方法均受限于預先給定的先驗信息。為了克服此種局限性，便產(chǎn)生了諸如聚類分析[4-7]、稀疏成分分析[8-10]及獨立成分分析(independent component analysis, ICA)[11-12]等數(shù)據(jù)驅(qū)動方法。與模型驅(qū)動方法相比，數(shù)據(jù)驅(qū)動方法因不需要任何的模型假設條件而被廣泛應用于fMRI信號分析的研究。這些方法在fMRI信號處理上展現(xiàn)出較高的準確性與實用性，在具體應用中有較大的優(yōu)勢，但同時也存在不足之處。

ICA是一種基于高階統(tǒng)計量的盲源信號分離方法，其目標是利用高階統(tǒng)計量的概率統(tǒng)計特性，在假定源信號相互獨立非高斯信號的基礎上，將采集到的混合信號分解為若干相互獨立的成分[13]。自從1998年McKeown等首次將該方法應用于fMRI數(shù)據(jù)分析以來[14]，ICA已經(jīng)成為該領(lǐng)域最為流行的方法之一，被廣泛應用于fMRI腦功能連通性檢測研究[15-21]。作為一種數(shù)據(jù)驅(qū)動分離技術(shù)，ICA無需依賴于任何先驗信息，因而具有獨特的優(yōu)勢。但同時研究也發(fā)現(xiàn)，ICA在fMRI分析過程中具有其自身的局限性，如成分個數(shù)不確定、成分輸出無序的問題等，由此可能導致大量計算時間和計算資源的浪費，特別是當待分析數(shù)據(jù)量比較大時問題就更加明顯。

針對上述問題，研究人員提出通過在ICA模型中加入可利用先驗信息來加以克服，而研究結(jié)果表明：在ICA的計算過程中引入先驗知識，確實能夠提高ICA的fMRI數(shù)據(jù)分析性能[22-26]。迄今為止，學者們已經(jīng)提出了許多方法，以解決關(guān)于ICA中包含先驗信息的問題。例如：Lu和Rajapakse在ICA的基礎上，率先提出了一種約束獨立成分分析(constrained ICA, CICA))來解決ICA融合先驗信息的問題[27]，并進一步提出了基于增廣拉格朗日方法和牛頓迭代算法的ICA with reference (ICA-R)方法[28]。隨后，Lin等在CICA和ICA-R方法的基礎上，選用不同的度量函數(shù)將空間先驗信息引入ICA模型，并采用快速不動點算法進行求解，提出了所謂的半盲空間ICA方法，并在實際中表現(xiàn)出良好的分析能力[29]。與ICA方法相比，加入先驗信息的CICA方法可以只提取所需要的源信號信息，這樣不僅有利于后續(xù)的應用分析，而且還能減少算法計算過程中所需要的時間和存儲容量[30]。另外，加入先驗信息能夠提高ICA算法的計算性能，使它更加適應于實時數(shù)據(jù)的分析環(huán)境[31-32]，并提高感興趣源信號檢測的質(zhì)量和精度[33]。

然而，目前的CICA方法都需要一個閾值參數(shù)來約束輸出信號與先驗參考信號之間的相似性，使得算法的輸出結(jié)果是與參考信號最為接近的唯一解。而源信號的不可知性，使得確定一個合適的閾值參數(shù)來約束輸出信號與先驗參考信號之間的近似度十分不容易，不合適的閾值往往會導致截然不同的結(jié)果，因此發(fā)展一種不依賴于閾值參數(shù)選擇的包含先驗信息的ICA方法顯得十分必要。最近，Du等提出了一種多目標的優(yōu)化策略來實現(xiàn)傳統(tǒng)CICA方法的求解問題，通過將傳統(tǒng)CICA方法中包含先驗信息的不等式約束條件轉(zhuǎn)變?yōu)橐粋€優(yōu)化目標函數(shù)，與原始ICA的獨立性目標度量函數(shù)構(gòu)成多目標優(yōu)化策略，這樣就規(guī)避了傳統(tǒng)CICA中閾值參數(shù)選擇不當導致分解結(jié)果不理想的問題[34]。

此外，目前能夠嵌入的先驗信息通常是一些具體且可直接利用的先驗知識，如實驗中激勵任務的刺激模式、研究比較成熟的腦網(wǎng)絡模板(如視覺網(wǎng)絡)、源信號的某些特定屬性等。但是，這些先驗信息并不總是存在的，或在某些情況下可能存在卻不易獲得，比如在靜息狀態(tài)或復雜認知任務狀態(tài)下的fMRI數(shù)據(jù)就難以獲取這樣的先驗信息。正是由于這些先驗信息難以獲取，使得一些研究不得不放棄CICA方法，轉(zhuǎn)而使用傳統(tǒng)的ICA方法恢復所有的獨立源信號，然后進行后期選擇。由此看來，挖掘數(shù)據(jù)隱含的先驗信息，不僅是為了克服CICA中先驗信息難以獲取的問題，而且對于CICA方法的廣泛應用也十分有益。因此，如何從已有條件中獲取可利用的隱含先驗信息就顯得尤為重要。

無論是從實驗中直接獲取的先驗信息，還是通過數(shù)據(jù)驅(qū)動方法挖掘的隱含先驗信息，它們都只是提供了一些有關(guān)源信號的粗糙信息，并不可能百分之百正確，而不同精度的先驗信息對提高ICA算法的分析能力也是不同。因此，如何減少CICA方法對參考信號的精度要求，用盡量少的先驗信息來最大程度提高ICA的分析能力，是一個值得研究且非常有意義的問題。此外，在有些可利用的先驗信息里面，可能還會或多或少含有一些錯誤信息，從而會誤導算法優(yōu)化過程中的收斂方向，使得輸出的結(jié)果不是所期望的結(jié)果。因此，如何提高包含先驗信息ICA算法的信息容錯能力，減少先驗信息中錯誤信息對算法的不利影響，是一個亟待解決的關(guān)鍵問題。Wang等研究發(fā)現(xiàn)，同時利用時間和空間先驗信息，要比單獨使用其中某種先驗信息能夠使CICA算法更具有魯棒性，減少了CICA方法對單一先驗信息的準確度要求[35]。

因此，本研究在多目標優(yōu)化框架的基礎上，提出了一種同時融合時空先驗信息的多目標優(yōu)化CICA (multi-objective optimization based CICA, MOPCICA)方法，有效規(guī)避了CICA方法中閾值參數(shù)難以選擇的問題。通過引入了時間和空間先驗信息，不僅降低了CICA方法對先驗信息的精度要求，而且提高了對先驗信息中錯誤信息的容錯能力。實驗結(jié)果表明，MOPCICA方法在時間和空間源信號恢復能力方面均展現(xiàn)出了良好性能，取得了不錯的腦功能連通性檢測效果。同時，提出了一種自適應挖掘本真參考信號的方法，通過從多被試fMRI數(shù)據(jù)中挖掘可利用的先驗信息，指導fMRI數(shù)據(jù)組被試水平上的ICA分析。實驗結(jié)果表明，該方式在fMRI組分析中具有優(yōu)越性能?？傊?，MOPCICA方法不僅克服和解決了現(xiàn)有CICA中存在的一些關(guān)鍵問題，同時也使得研究人員能夠更加方便地使用CICA方法，使CICA方法能夠更加廣泛地應用于fMRI的各項研究，進而提高ICA在fMRI數(shù)據(jù)分析方面的能力。

1 材料和方法

1.1 材料

1.1.1模擬數(shù)據(jù)

本研究中的模擬數(shù)據(jù)由SimTB工具包(http://mialab.mrn.org/spftware)生成，該工具包由Vince D. Calhoun博士領(lǐng)導下的醫(yī)學影像實驗室發(fā)布，能夠提供模擬信號的空間圖譜與相應的時間過程或時間腦動力模式，方便不同算法進行性能對比評測。具體來說，采用SimTB生成10例模擬數(shù)據(jù)，每個被試的切片大小為100體素×100體素，共模擬了20個源信號(這些源信號直接從SimTB提供的信號模板中選取)，120個時間點，TR為2 s。該模擬數(shù)據(jù)中基準幅值設定為800，基準圖譜如圖1(a)所示，20個源信號的空間分布如圖1(b)所示。為了模擬不同被試之間的空間差異，每個被試的源信號都通過正態(tài)偏差加入了少量的平移、旋轉(zhuǎn)和伸縮變換。其中，平移和旋轉(zhuǎn)變換分別服從正態(tài)分布N(0, 0.5)和N(0, 1)，而伸縮變換則服從正態(tài)分布N(1, 0.1)。

圖1 模擬數(shù)據(jù)集。(a)空間基準圖譜；(b)20個模擬源信號的空間分布；(c)所有模擬被試對應源信號s6的空間圖譜和相應的時間過程(其中空間圖譜中的橫坐標和縱坐標均表示模擬的體素個數(shù)，而時間過程圖中的橫坐標和縱坐標則分別代表時間點數(shù)和去中心化后的幅度大小)Fig.1 Simulated dataset. (a) Spatial baseline map；(b)Spatial distribution of the 20 simulated sources；(c)Spatial map and corresponding time course of source 6 in all simulated subjects，where the abscissa and ordinate in the spatial map represent the number of simulated voxels, while the abscissa and ordinate in the time course subgraph represent the number of time points and the amplitude after decentralization, respectively

源信號s6除了混合特定的事件刺激之外，其時間過程還共享了特定區(qū)塊事件刺激模式的影響，且區(qū)塊事件相關(guān)的信號強度(定義為區(qū)塊事件幅值與特定事件幅值之比)為2。區(qū)塊事件刺激模式為(OFF, ON)×3，其中OFF或ON區(qū)塊持續(xù)的時間長度均為40 s。通過將該刺激模式與典型的血氧動力學響應函數(shù)進行卷積，生成源信號s6的時間過程。剩下的19個源信號均未受到區(qū)塊事件調(diào)制模式的影響，其活動模式僅由特定的血氧動力學波動信號(即通過特定事件與血氧動力學響應函數(shù)進行卷積生成)來描述。

所有模擬信號均受特定事件的影響，特定事件在每個TR內(nèi)隨機發(fā)生的概率為0.2。然而，對于受區(qū)塊事件調(diào)制影響的源信號s6而言，特定事件的活動幅值略小，為0.4；對于不受區(qū)塊事件調(diào)制影響的源信號而言，特定事件的幅值為1。對于所有的模擬信號，信號在空間域的變化強度服從正態(tài)分布N(0.03, 0.25)；與此同時，模擬信號中引入了加性噪聲，其信噪比(signal-noise-ratio, SNR)服從0.2～1.2的均勻分布。最后，所有模擬被試對應源信號s6的空間圖譜和時間過程如圖1(c)所示。

1.1.2任務態(tài)數(shù)據(jù)

任務態(tài)數(shù)據(jù)共包括5名被試在視覺刺激任務下掃描獲取的BOLD-fMRI數(shù)據(jù)。在進行數(shù)據(jù)采集之前，所有被試均被告知此次采集實驗數(shù)據(jù)的目的，并獲得了所有被試的數(shù)據(jù)采集許可，本次任務實驗中的視覺刺激為具有放射狀藍黃團的旋轉(zhuǎn)棋盤，其相應的視覺刺激為(OFF, ON)×3，每個區(qū)塊(OFF或者ON狀態(tài))持續(xù)時間均為40 s；當刺激模式為ON時，刺激呈現(xiàn)屏中旋轉(zhuǎn)棋盤；當刺激模式為OFF時，則要求被試眼睛注視刺激呈現(xiàn)屏的十字中心。該數(shù)據(jù)集中所有被試數(shù)據(jù)均在上海市磁共振重點實驗室采集完成，采用磁場強度為3.0 T的西門子(Siemen)磁共振儀以及梯度回波平面成像技術(shù)進行覆蓋全腦掃描獲得，共掃描125個時間點。數(shù)據(jù)獲取的其他參數(shù)為：TR=2 s，切片數(shù)=36片，掃描分辨率=64×64，片內(nèi)分辨率=3.75 mm×3.75 mm，片厚=4 mm，切片間隔=1 mm。

1.1.3靜息態(tài)fMRI數(shù)據(jù)

實驗中使用的靜息態(tài)fMRI數(shù)據(jù)集來自國際神經(jīng)影像公共數(shù)據(jù)庫(http://www.nitrc.org/projects/fcon_1000/)，包括23名認知能力正常、身體健康被試的靜息態(tài)BOLD-fMRI掃描數(shù)據(jù)，被試年齡在20～40歲之間。所有被試數(shù)據(jù)通過磁場強度3.0T的飛利浦磁共振儀，以及單次激發(fā)敏感梯度回波平面成像技術(shù)覆蓋全腦掃描獲得，共掃描123個時間點。數(shù)據(jù)采集的其他參數(shù)為：TR=2.5 s，切片數(shù)=47，敏感加速因子=2.0，掃描分辨率=96×96，片內(nèi)分辨率=2.67 mm×2.67 mm，片厚度=3 mm。

1.2 方法

1.2.1數(shù)據(jù)預處理

采用DPARSF(http://rfmri.org/DPARSF)軟件對實驗中的任務態(tài)和靜息態(tài)fMRI數(shù)據(jù)進行預處理，包括時間層校正、頭動校正、空間標準化和空間平滑，其中高斯平滑核設定為4 mm。特別地，在模擬數(shù)據(jù)和任務態(tài)數(shù)據(jù)實驗中，僅實現(xiàn)了單被試水平上的ICA、包含時間先驗信息的CICA (CICA with temporal reference, CICA-tR)[27]、包含空間先驗信息的CICA (CICA with spatial reference, CICA-sR)[27]，以及本研究提出的MOPCICA方法；而在靜息態(tài)數(shù)據(jù)中，則同時實現(xiàn)了組被試水平上的ICA、包含空間先驗信息的基于牛頓迭代法的CICA(Newton′s method based CICA，CICA-nR)[28]和基于不動點迭代法的CICA(fixed point iteration based CICA，CICA-fR)[29]，以及MOPCICA方法和單被試水平上的ICA方法，并且在組分析中使用時空雙回歸方式來獲取每個被試對應的時空信息。另外，采用最小描述長度方法[36]來對真實fMRI數(shù)據(jù)中的源信號個數(shù)進行估計。同時，為了在分析過程中獲得更加穩(wěn)定的輸出結(jié)果，在ICA計算過程中使用ICASSO[37]來降低輸出結(jié)果的隨機性，其參數(shù)設置為運行100次，采用隨機初始化和自助法相結(jié)合的采樣方式。而對于其他方法則在計算過程取100次運行結(jié)果的平均值作為最終結(jié)果。最后，使用MRIcro軟件，對各種方法輸出結(jié)果中的空間網(wǎng)絡和位置進行展示。

1.2.2MOPCICA方法優(yōu)化

為了克服經(jīng)典CICA方法中閾值參數(shù)難以選擇的問題，同時考慮到時空先驗信息的問題，本研究在多目標優(yōu)化框架上建立CICA模型如下：

(1)

本研究采用G(v)=ln(cosh(v))，v是一個零均值單位方差的高斯隨機變量，因此

(2)

在多目標優(yōu)化問題式(1)中，每個解都對應一個最優(yōu)分離向量wi，并且滿足‖wi‖2=1。

在諸多求解多目標優(yōu)化問題的方法中，線性加權(quán)求和方法是一種最為簡單和有效的方法，只要權(quán)重參數(shù)嚴格為正且和為1，就能通過對目標函數(shù)的線性加權(quán)來進行單目標優(yōu)化求解。因此，下面將采用這種方法，對所研究的多目標優(yōu)化問題進行求解。只要在加權(quán)和函數(shù)中使用不同的權(quán)重值，就能找到不同解。為了避免優(yōu)化過程被具有較大幅值的成本函數(shù)所控制，使用反正切函數(shù)來對式(1)中的目標函數(shù)J(wi)進行標準化，即

(3)

式中，ci是一個可以自動確定的常數(shù)，這樣f1(wi)就和ε1(wi)、ε2(wi)一樣，其可能的取值范圍從0～1變化。

重新整理后，線性加權(quán)目標函數(shù)為

f(wi)=a1×f1(wi)+a2×ε1(wi)+a3×ε2(wi)

(4)

式中，ai(i=1,2,3)表示權(quán)重參數(shù)，且a1+a2+a3=1，通常是根據(jù)經(jīng)驗來決定的。

那么，優(yōu)化f(wi)的迭代算法推導如下：

(5)

(6)

(7)

a3×E[Z-1×RTi]

(8)

式中，g(·)表示G(·)的導數(shù)，因此g(v)=tanh(v)，E(·)可以用所有樣本的均值來進行估計，具體如下：

(9)

當目標函數(shù)式(4)的梯度被計算出來之后，就能建立最速上升的迭代公式，即

wi(k+1)=wi(k)+μ(k)×di(k)

(10)

對于P個參考信號RSp和RTp(p=1,2,…,P)，可以使用相同的方法逐個估計。最后，當獲得獨立成分Si后，其相應的時間過程可通過下列公式估計得到，即

Ti=Z-1×wi

(11)

1.2.3自適應先驗信息挖掘

為了克服先驗信息難以獲取的問題，提出一種從組被試fMRI數(shù)據(jù)中獲取空間先驗信息的方法。假設在一組fMRI數(shù)據(jù)中有K個被試，且所有被試在經(jīng)過標準化預處理后包含T個時間點和V個體素。

首先，進行單被試水平上的ICA分析，對于被試i，ICA可以表示為

Xi=Mi×Si(i=1,2,…,K)

(12)

式中：Xi表示T×V階的觀測數(shù)據(jù)；Mi表示T×Ni階的混合矩陣；Si=(Si1,Si2,…,SiNi)'表示Ni×V階的源信號矩陣，并且每行代表一個被試i在ICA分解過程中得到的獨立成分，Sij(j=1,…,Ni)是大小為V×1的列向量。

接下來，記Sini(i=1,2,…,K)表示與感興趣組獨立成分對應被試i的第ni個獨立成分，不同被試間的成分對應性可通過空間相關(guān)性進行度量[37]。

下面利用主成分分析(principal component analysis, PCA)，從由Sini(i=1,2,…,K)組成的K×v階矩陣R來提取隱含參考信號，有

R=[S1n1,S2n2,…,SKnK]T

(13)

1.2.4算法驗證

為了驗證MOPCICA在fMRI數(shù)據(jù)分析中的有效性能，接下來分別通過模擬數(shù)據(jù)、任務態(tài)和靜息態(tài)fMRI數(shù)據(jù)中的分析結(jié)果進行評估。其中，模擬數(shù)據(jù)和任務態(tài)fMRI數(shù)據(jù)主要用于評測MOPCICA的算法性能，而靜息態(tài)fMRI數(shù)據(jù)則同時用于評估本真先驗信息挖掘算法的有效性。雖然MOPCICA方法可以同時包含時間和空間先驗信息，但是為了不同類型數(shù)據(jù)中前后分析的一致性，在分析中僅考慮了MOPCICA包含空間先驗信息的情形。

圖2 ICA、CICA-tR、CICA-sR和MOPCICA在每個模擬被試上對s6的源信號恢復性能評估。(a)空間相關(guān)系數(shù)；(b)空間AUC；(c)時間相關(guān)系數(shù)Fig.2 The performance evaluation of source signal recovery of s6 by ICA, CICA-tR, ICA-sR and MOPCICA in each simulated subject. (a)Spatial correlation coefficients；(b)Spatial AUC；(c)Temporal correlation coefficients

同時，在數(shù)據(jù)分析過程中，MOPCICA分別實現(xiàn)了空間權(quán)重參數(shù)從0.1～0.9、以步長為0.1取值時對應的9種情況，而在結(jié)果部分僅呈現(xiàn)了參數(shù)取值為0.5時的結(jié)果，主要是因為該參數(shù)在3種類型數(shù)據(jù)中均表現(xiàn)出良好的分析性能。另外，為了比較目的，在模擬數(shù)據(jù)和任務態(tài)fMRI數(shù)據(jù)中，同時呈現(xiàn)了ICA、CICA-tR和CICA-sR等方法的分析結(jié)果，并且三者均通過牛頓迭代算法實現(xiàn)。而在靜息態(tài)fMRI數(shù)據(jù)中，則同時呈現(xiàn)了ICA、 CICA-nR和CICA-fR等方法的分析結(jié)果，并通過牛頓迭代算法實現(xiàn)ICA，還通過牛頓迭代算法實現(xiàn)CICA-nR和快速不動點算法實現(xiàn)CICA-fR。

最后，在3種類型的數(shù)據(jù)實驗中，分別采用多種評價指標來評測比較上述不同方法的數(shù)據(jù)分析性能，而不同方法所對應相同指標之間的對比結(jié)果均通過置信水平為95%的雙樣本t檢驗統(tǒng)計獲得。其中，在模擬數(shù)據(jù)實驗中，采用空間AUC和時空相關(guān)系數(shù)來評測ICA、CICA-tR、CICA-sR和MOPCICA的源信號恢復性能。在任務態(tài)fMRI數(shù)據(jù)實驗中，同樣采用時空相關(guān)系數(shù)來評測ICA、CICA-tR、CICA-sR和MOPCICA的源信號恢復性能。同時，使用峭度和負熵來評測它們所得源信號的空間獨立性。類似地，在靜息態(tài)fMRI數(shù)據(jù)實驗中，利用互信息、峭度和負熵來度量ICA、CICA-nR、CICA-fR和MOPCICA所得源信號的空間獨立性。進一步，根據(jù)它們計算得到的組空間成分，通過時空雙回歸方式獲得的組中每個被試對應空間成分之間的平均相關(guān)系數(shù)，評測其組分析性能。

2 結(jié)果

下面主要給出MOPCICA在模擬數(shù)據(jù)、任務態(tài)和靜息態(tài)fMRI數(shù)據(jù)中的分析結(jié)果，以及與其他方法結(jié)果之間的比較情況。

2.1 模擬數(shù)據(jù)結(jié)果

為了評估上述方法的空間源信號恢復性能，分別計算了它們根據(jù)每個模擬被試數(shù)據(jù)獲得的源信號成分與對應真實源信號之間的相關(guān)系數(shù)(即空間相關(guān)系數(shù))，如圖2(a)所示?？梢园l(fā)現(xiàn)，雖然與ICA相比無明顯差異，但是MOPCICA在大部分被試上計算得到的相關(guān)系數(shù)都要高于CICA-tR和CICA-sR計算得到的相關(guān)系數(shù)，除了CICA-tR中被試2、3和10以及CICA-sR中被試1～3和10之外，這些有可能是由于噪聲SNR不同而造成的結(jié)果。進一步，計算了它們在所有被試上的平均空間相關(guān)系數(shù)，對應ICA、CICA-tR、CICA-sR和MOPCICA分別為0.84±0.07、0.78±0.08、0.81±0.08和0.84±0.09。通過置信水平為95%的雙樣本t檢驗可知，MOPCICA顯著高于CICA-tR和CICA-sR(P<0.05)，從而說明MOPCICA的空間源信號恢復能力要優(yōu)于CICA-tR和CICA-sR的相關(guān)能力。

接著，ROC曲線中的能量分析也被用來評估這些方法的空間源信號恢復能力，表示為ROC曲線下的面積(AUC)。一般地，AUC越大，表示該方法對應的源信號恢復能力越強。圖2(b)分別給出了ICA、CICA-tR、CICA-sR和MOPCICA在不同模擬被試上計算得到與源信號s6相應成分的AUCs?？梢钥闯觯琈OPCICA在絕大部分被試上計算得到的AUCs都要高于ICA、CICA-tR和CICA-sR方法所得到的AUCs。它們在所有被試上的平均AUC分別為0.62±0.02、0.72±0.03、0.71±0.06和0.75±0.05，由雙樣本t檢驗可知，MOPCICA顯著高于ICA、CICA-tR和CICA-sR(P<0.05)，表明MOPCICA具有更強的空間源信號恢復能力。

圖3 ICA、CICA-tR、CICA-sR、MOPCICA和GLM計算得到的每個被試對應視覺相關(guān)網(wǎng)絡的腦空間圖譜及其時間過程(GLM除外)。其中，腦空間圖譜通過z-score變換在閾值為2的條件下計算得到，對應的MNI坐標為(-10, -77, 6)；而每個時間過程小圖中的橫坐標和縱坐標分別代表時間點數(shù)和去中心化后的幅度大小，紅線和藍線則分別表示先驗時間過程和實際計算得到的時間過程Fig.3 The spatial brain maps of the visual-related networks calculated by ICA, CICA-nR, CICA-fR, MOPCICA and GLM in each subject, as well as their time courses (except for GLM). The maps of brain networks are obtained with threshold |z| ≥ 2 after z-scored, and the corresponding MNI coordinates are (-10, -77, 6). The abscissa and ordinate of each time course subplot respectively represent the time points and the amplitude after decentralization, while the red line and the blue line respectively represent the prior time course and the time course obtained by actual calculation

此外，根據(jù)上述方法在每個模擬被試中計算得到的時間過程與相應真實時間過程之間的相關(guān)性，度量它們在時間域上的源信號恢復能力，如圖2(c)所示。分別給出了ICA、CICA-tR、CICA-sR和MOPCICA在每個模擬被試上計算得到的時間過程與該被試對應真實時間過程之間的相關(guān)系數(shù)，即時間相關(guān)系數(shù)。從圖中可以看出，除了少數(shù)被試之外，MOPCICA在大部分被試上的相關(guān)系數(shù)都要高于ICA、CICA-tR和CICA-sR計算得到的相關(guān)系數(shù)，這意味著MOPCICA在時間源信號恢復性能方面也要優(yōu)于ICA、CICA-tR和CICA-sR的相關(guān)能力。而且，它們在所有被試上的平均時間相關(guān)系數(shù)分別為0.67±0.04、0.74±0.09、0.77±0.13和0.81±0.13，由雙樣本t檢驗可知，MOPCICA顯著高于ICA、CICA-tR和CICA-sR(P<0.05)，從而進一步證明了MOPCICA在源信號檢測方面的優(yōu)越性能。

2.2 任務態(tài)數(shù)據(jù)結(jié)果

圖3分別給出了ICA、CICA-tR、CICA-sR、MOPCICA和GLM檢測得到每個被試對應視覺相關(guān)成分的腦空間圖譜和時間過程(GLM除外), 以及該任務實驗中的先驗時間過程?？紤]到大腦對外界刺激響應的延遲效應，這里的先驗時間過程由實驗中的Block刺激模式與延遲參數(shù)為2的伽馬函數(shù)卷積得到。從圖3中可以看到，MOPCICA檢測得到的視覺區(qū)域要明顯優(yōu)于CICA-tR和ICA檢測得到的視覺區(qū)域，但與CICA-sR檢測得到的視覺區(qū)域相比沒有明顯差異，需要進一步量化分析，這意味著MOPCICA的空間源信號恢復能力要優(yōu)于CICA-tR和ICA的空間源信號恢復能力。

圖4 ICA、CICA-tR、CICA-sR和MOPCICA在每個被試上對視覺相關(guān)成分的源信號恢復性能評估。(a)空間相關(guān)系數(shù)；(b)峭度；(c)負熵；(d)時間相關(guān)系數(shù)Fig.4 The performance evaluation of source signal recovery of visual-related component by ICA, CICA-tR, CICA-sR and MOPCICA in each subject. (a)Spatial correlation coefficients；(b)Kurtosis；(c)Negentropy；(d)Temporal correlation coefficients

為了進一步評估上述方法的空間源信號恢復能力，分別計算了ICA、CICA-tR、CICA-sR和MOPCICA從每個被試fMRI數(shù)據(jù)中計算得到的視覺相關(guān)成分與相應由GLM方法計算得到的視覺相關(guān)成分之間的相關(guān)系數(shù)，如圖4(a)所示?？梢园l(fā)現(xiàn)，除了少數(shù)被試差異不明顯之外，MOPCICA在所有被試上計算得到的相關(guān)系數(shù)都要顯著高于其他方法計算得到的相關(guān)系數(shù)。進一步計算了它們在所有被試上的平均相關(guān)系數(shù)，分別為0.73±0.11、0.28±0.05、0.72±0.03和0.84±0.04。由雙樣本t檢驗可知MOPCICA顯著高于ICA、CICA-tR和CICA-sR(P<0.05)，意味著MOPCICA的空間源信號恢復能力要優(yōu)于ICA、CICA-tR和CICA-sR的相應能力，與圖3中的結(jié)論相一致。

同時，給出了上述方法計算得到每個被試對應視覺相關(guān)成分的峭度和負熵，以度量它們所得源信號的空間獨立性，如圖4(b)和4(c)所示。從圖可以看出，除了ICA之外，MOPCICA在所有被試上計算獲得的峭度和負熵都要高于CICA-tR和CICA-sR的相應情況。同時，計算了它們在所有被試上的平均峭度和負熵，分別為93.97±50.39、17.60±13.22、36.71±13.43和69.20±23.36以及0.030 2±0.004 9、0.003 7±0.002 1、0.018 4±0.004 5和0.031 2±0.007 7。根據(jù)雙樣本t檢驗可知，MOPCICA顯著高于CICA-tR和CICA-sR，而平均峭度則低于ICA(P<0.05)，這說明MOPCICA估計所得源信號的獨立性要優(yōu)于CICA-tR和CICA-sR相應的獨立性。

與模擬數(shù)據(jù)的結(jié)果分析類似，由計算得到的每個被試對應視覺相關(guān)成分的時間過程與先驗時間過程之間的相關(guān)性被用來度量它們在時間域上的源信號恢復性能，圖4(d)分別給出了ICA、CICA-tR、CICA-sR和MOPCICA在每個被試上計算得到的時間過程與先驗時間過程之間的相關(guān)系數(shù)?？梢钥闯?，通過MOPCICA計算得到的相關(guān)系數(shù)要明顯高于ICA、CICA-tR和CICA-sR計算得到的相關(guān)系數(shù)，而且它們在所有被試上的平均時間相關(guān)系數(shù)分別為0.66±0.07、0.90±0.01、0.85±0.05和0.91±0.03。

圖5 ICA、CICA-nR、CICA-fR和MOPCICA計算得到的8種經(jīng)典靜息態(tài)腦網(wǎng)絡對應的腦空間圖譜及其相應的MNI坐標(腦網(wǎng)絡圖譜均通過z-score變換在閾值為2的條件下計算獲得)Fig.5 The spatial maps of nine resting-state brain networks and their MNI coordinates including DMN, VIN, LVN, AUN, SMN, ECN, RWMN, LWMN and CCN calculated by ICA、CICA-nR、CICA-fR and MOPCICA (The maps of brain networks are obtained with threshold |z| ≥ 2 after z-scored)

由雙樣本t檢驗可知，MOPCICA顯著高于ICA和CICA-sR(P<0.05)，說明MOPCICA在時間域上的分析性能要優(yōu)于ICA、CICA-tR和CICA-sR等方法的相應性能。

2.3 靜息態(tài)數(shù)據(jù)結(jié)果

圖5分別給出了ICA、CICA-nR、CICA-fR和MOPCICA檢測得到的空間獨立成分對應的腦網(wǎng)絡圖譜及其相應的MNI坐標，包括默認網(wǎng)絡(default mode network, DMN)、主視覺網(wǎng)絡(visual network, VIN)、兩側(cè)視覺網(wǎng)絡(bilateral visual network, BVN)、聽覺網(wǎng)絡(auditory network, AUN)、執(zhí)行控制網(wǎng)絡(executive control network, ECN)、突顯網(wǎng)絡(salience network, SAN)、右側(cè)工作記憶網(wǎng)絡(right working memory network, RWMN)和左側(cè)工作記憶網(wǎng)絡(left working memory network, LWMN)等8種經(jīng)典靜息態(tài)腦網(wǎng)絡，其中所有腦網(wǎng)絡圖譜均通過z-score變換在閾值為2的條件下計算獲得。由圖5可知，上述方法均獲得了明顯的靜息態(tài)腦功能網(wǎng)絡，但無法看出明顯差異，需要進一步定量分析。

圖6 ICA、CICA-tR、CICA-sR和MOPCICA在組水平上對8種靜息態(tài)腦網(wǎng)絡的源信號恢復性能評估。(a)互信息；(b)峭度；(c)負熵；(d)組成分與通過雙回歸方式計算的組中每個被試對應獨立成分之間的相關(guān)系數(shù)Fig.6 The performance evaluation of source signal recovery of eight resting-state brain networks by ICA, CICA-tR, ICA-sR and MOPCICA at group level. (a)Mutual information；(b)Kurtosis；(c)Negentropy；(d) The correlation coefficients between the group independent components and the corresponding independent components of each subject in the group obtained by dual-regression

圖6給出了ICA、CICA-nR、CICA-fR和MOPCICA檢測得到的8種靜息態(tài)腦網(wǎng)絡對應成分的空間獨立性比較結(jié)果，其中獨立性度量指標包括互信息(mutual information)、峭度(kurtosis)和負熵(negentropy)等3種情形分別如(a)～(c)所示，而(d)則分別給出了這4種方法計算得到的組空間成分與通過時空雙回歸方式獲得的組中每個被試對應空間獨立成分之間的相關(guān)系數(shù)。

由圖6(a)可知，除了AUN對應的空間成分之外，MOPCICA計算得到的其他7種網(wǎng)絡對應空間成分的互信息要低于其他方法計算得到的空間獨立成分的互信息；而圖6(b)和(c)的結(jié)果表明，除了少數(shù)幾個腦網(wǎng)絡之外，MOPCICA腦網(wǎng)絡對應空間成分的峭度和負熵均要高于其他方法的值。進一步，計算了它們在所有被試上對應8種靜息態(tài)腦網(wǎng)絡的平均互信息、平均峭度和平均負熵，分別為0.44±0.04、0.44±0.04、0.44±0.04和0.43±0.04，8.68±3.28、8.85±3.32、8.78±3.25和9.15±3.64，以及5.09±5.10、4.23±4.30、4.00±4.50和4.88±4.81，根據(jù)t檢驗可知，除了ICA負熵之外，MOPCICA均顯著優(yōu)于其他3種方法(P<0.05)，從而說明MOPCICA計算得到的空間成分具有更好的獨立性，更加符合ICA分析的前提假設條件，因而結(jié)果的可信性更高。

由圖6(d)可知，相比于其他方法，MOPCICA計算得到的組成分與組中每個被試相應獨立成分之間具有更高的相關(guān)性，而且它們在所有被試上對應8種靜息態(tài)腦網(wǎng)絡的平均相關(guān)系數(shù)為0.44±0.02、0.45±0.02、0.44±0.02和0.46±0.03。根據(jù)t檢驗可知，MOPCICA顯著高于ICA、CICA-nR和CICA-fR(P<0.05)，說明MOPCICA計算得到的組成分更能代表組中被試的共性，具有更好的fMRI腦功能連通性檢測性能。

3 討論

在經(jīng)典的CICA方法中，先驗參考信號是通過不等式約束條件g(y)=ε(y,r)-ξ≤0的方式引入ICA算法的估計過程，y表示輸出信號，r表示參考信號，ε(y,r)是距離測度，而ξ表示閾值參數(shù)，用于約束感興趣信號是唯一滿足不等式約束條件的輸出信號[27-28, 30]。然而，在實際應用中，往往很難預先確定閾值參數(shù)ξ，因為獨立成分是未知的，所以ξ的選擇完全依賴于CICA的應用經(jīng)驗。不合適的ξ會導致兩種可能的結(jié)果：當ξ大于可行域的上限時，輸出信號可能是不感興趣的信號；反之，當ξ小于可行域的下限時，可能無法輸出任何的信號[34]。因此，往往需要額外的努力來決定一個合適的參數(shù)ξ。本研究采用了多目標優(yōu)化策略，很好地規(guī)避了閾值參數(shù)ξ的選擇問題；并且Shi等也采用了相同的策略，實驗結(jié)果表明該策略的優(yōu)越性能[38]。

在利用線性加權(quán)求和方法來求解式(1)中的多目標優(yōu)化問題時，式(4)中的權(quán)重參數(shù)a1、a2和a3分別反映了和函數(shù)f(wi)中對應目標函數(shù)f1(wi)、ε1(wi)和ε2(wi)的重要性。加權(quán)求和方法的目標是在輸出信號獨立性和參考信號相似性之間尋求一個平衡，從而獲得一個與參考信號最接近且獨立性最大的源信號。根據(jù)Klamroth等提出的應用加權(quán)求和方法[39]，求解多目標優(yōu)化問題的理論，只要權(quán)重參數(shù)滿足它們嚴格為正且和為1的條件，那么每一組這樣的權(quán)重參數(shù)都能獲得Pareto最優(yōu)集中的一個解。

在MOPCICA算法驗證過程中，在模擬數(shù)據(jù)和任務態(tài)fMRI數(shù)據(jù)分析過程中，均只考慮一種源信號的情形。例如，在模擬數(shù)據(jù)中，僅考慮與特定區(qū)塊事件刺激模式相關(guān)的源信號s6，主要是因為它在所有模擬被試上均有一個類似于區(qū)塊事件刺激模式Block形狀的時間過程[40]。特別地，這里的Block刺激模式為 (OFF, ON)×3，其中OFF或ON區(qū)塊持續(xù)的時間長度均為40 s，這樣可方便分析結(jié)果的評測；并且在計算過程中，分別利用該Block刺激模式和相應源信號s6的真實空間模板，作為時間和空間參考信號。同時，在分析過程中，模擬數(shù)據(jù)成分個數(shù)設置為21(包含20個源信號和1個背景信號)[41]。

任務態(tài)fMRI數(shù)據(jù)主要考慮與視覺任務相關(guān)的源信號成分，其對應的空間參考信號利用wfu_pickatlas軟件(http://www.rad.wfubmc.edu/fmri)并根據(jù)布羅德曼模板產(chǎn)生，包括BA17和BA18兩個腦區(qū)[42]。同時，利用先驗的Block刺激模式作為時間參考信號。為了保證不同方法結(jié)果之間比較的一致性，結(jié)果分析中選取了GLM計算得到的視覺空間成分作為基準模板來進行比較驗證，其中 GLM通過SPM(http://www.fil.ion.ucl.ac.uk/spm/)軟件包實現(xiàn)，參數(shù)設定為使用FWE校正，p值取0.01，體素閾值個數(shù)設為10[3]。另外，在數(shù)據(jù)分析過程中，5名被試對應的源信號個數(shù)分別通過最小描述長度估計為50、52、54、54和44。

但是，這并不意味著MOPCICA僅適用于一種源信號的情形，因此在靜息態(tài)fMRI數(shù)據(jù)分析過程中，考慮了多種源信號(對應于多種腦網(wǎng)絡)的情形。由于靜息態(tài)fMRI數(shù)據(jù)無法獲得類似于模擬數(shù)據(jù)和任務態(tài)fMRI數(shù)據(jù)的時間先驗信息，因此在數(shù)據(jù)分析過程中，首先通過本真先驗信息挖掘算法獲得幾種經(jīng)典腦網(wǎng)絡(見圖5)對應的空間先驗信息，其中每個被試對應的源信號個數(shù)通過最小描述長度方法估計得到；然后利用ICA以及包含空間先驗信息的CICA-sR和MOPCICA方法進行組水平數(shù)據(jù)分析，其中組水平源信號個數(shù)通過最小描述長度估計為26(為所有被試源信號估計個數(shù)的平均值)。上述方法在組分析過程中均采用時間級聯(lián)方式，即將所有被試數(shù)據(jù)根據(jù)時間維度進行連接分析[43]，所以在實驗數(shù)據(jù)預處理過程中，需要進行單被試和組水平兩次主成分分析，且主成分個數(shù)分別為23和26。

在獲取空間信息過程中，由于個體自身的差異性以及每個被試在采集數(shù)據(jù)過程中所受到的噪聲影響不同，使得組中每個被試的成分不僅包含相同的有用信息，同時也包含不同的噪聲信息[44]。因此，選擇一種合適的方法，從這些成分中挖掘出最主要的信息十分關(guān)鍵。作為一種經(jīng)典的多元數(shù)據(jù)處理方法，PCA不僅能夠降低數(shù)據(jù)維度，提取出原數(shù)據(jù)中包含的主要信息，而且還能起到降噪的作用。所以，本研究采用PCA技術(shù)，從組被試的成分中挖掘出可利用的先驗信息。因為第一主成分包含有原數(shù)據(jù)中最大量的信息，所以選取它為相應的參考信號。特別地，當特征值λ1=λ2時，選擇它們對應主成分的平均值作為相應的參考信號。除此之外，其他方式(如平均等)也可被用來從組被試數(shù)據(jù)中求解隱含的參考信號，值得后續(xù)進一步研究。

4 結(jié)論

本研究在多目標優(yōu)化的基礎上，提出了一種同時包含時間和空間參考信號的MOPCICA方法，克服了經(jīng)典CICA方法中閾值參數(shù)難以選擇的問題，降低了CICA方法對參考信號的準確性要求，提高了它對錯誤信息的容錯能力。同時，利用從fMRI自身挖掘的空間先驗信息來指導組數(shù)據(jù)分析，使得MOPCICA獲得的空間源信號成分能更好地代表組中被試的共性。