低復雜度的fMRI腦激活區定位的盲分離算法

陳安瑩，吳海鋒，李 棟

低復雜度的fMRI腦激活區定位的盲分離算法

陳安瑩，吳海鋒，李 棟

(云南民族大學電氣信息工程學院，云南 昆明 650500)

功能磁共振成像(FMRI)是一種醫學影像技術，由于具有非侵入性和較高的時空分辨率等優點現已被廣泛應用于腦區定位。然而傳統的FMRI信號分離算法復雜度太高，運行時間長，不利于FMRI技術更有效地應用于腦功能的研究。針對傳統FMRI腦區分離算法的計算復雜度問題，提出了一種基于二階哈達碼變換的盲分離算法。先計算fMRI數據中血氧水平依賴(BOLD)信號的相關函數，然后對其進行特征值分解得到解混矩陣，以此實現激活腦區定位。由于哈達碼只由1或-1構成，因此可減少BOLD信號相關矩陣計算的復雜度。仿真結果表明，相比高階統計量的獨立分量分析(ICA)和二階統計量的傅里葉變換盲分離算法，該算法的計算時間分別只有其25%和50%，而定位誤差卻較為接近。

功能磁共振成像；盲分離；獨立分量分析；二階統計量的盲辨識；腦激活區

人類的大腦皮層可以根據功能劃分成不同的區域，如視覺區、運動區等，若該區域遭受損傷，可通過醫學影像技術對其進行定位，從而確定相應的治療方案。由于功能磁共振成像(functional magnetic resonance imaging, fMRI)具有非侵入性和高空間分辨率等優點，近年來已成為醫學成像中的常用腦區定位技術。

在傳統fMRI激活腦區定位算法中，統計參數圖(statistical parametric mapping，SPM)[1]是一種常用的方法，其利用統計檢驗方法對受試者的不同成像結果進行比較來尋找具有統計學意義的激活區。然而，該方法需要預先知道設計矩陣，當設計矩陣不確定時往往會影響最終的定位效果。獨立分量分析(independent components analysis, ICA)算法[2]采用盲分離的方式，將腦區信號作為源信號分離出，相比SPM方法，ICA并不需要設計矩陣已知[3]。然而，ICA類盲算法通常需要計算高階統計量，且需要較多的迭代次數以保證收斂，因此計算復雜度較高。二階統計量的盲辨識(second order blind identifiability，SOBI)算法[4]僅需計算相關函數就能實現信號分離，不需要計算高階統計量，復雜度相對較低。然而，由于fMRI的腦區信號具有稀疏性，其相關函數信息量不夠豐富，因此影響了分離性能。頻域的SOBI(frequency SOBI, f-SOBI)算法[5]對稀疏信號進行反傅里葉變換，然后再求相關函數，因此使得原信號經過變換后不再稀疏，保證了變換后的相關函數包含了較多源信息，從而提高了分離性能。

針對fMRI腦區定位盲分離的計算復雜度問題，本文提出了一種哈達碼變換的二階盲分離算法(h-SOBI)，由于哈達碼變換矩陣只由1或-1構成，因此可降低相關函數計算復雜度，從而減少整個信號分離算法的計算量。實驗中采用SimTB (simulation toolbox)[6]工具箱產生fMRI仿真數據和實測數據，分別給出了fastICA，f-SOBI和本文算法的分離結果。實驗結果表明，本文算法的計算時間得到了有效的減少，而激活區的定位誤差卻與傳統分離算法相似。

1 相關工作

其中，為掃描時間點數；為設計矩陣；為解釋變量；為殘差。式(1)中，利用最小二乘準則估計，再對其做統計假設檢驗，即可推斷腦區的激活情況。然而該方法需要預先知道設計矩陣的信息，否則就難以被估計。

SOBI算法首先計算信號的相關函數，令和-為的第列和第-列，那么R()的相關函數為

其中，?{1,2,···,}為時延數。然后，通過聯合對角化(joint diagonalization，JD)[10]獲取解混矩陣。然而，由于fMRI數據具有稀疏性，其相關函數R()信息不夠豐富，導致相關函數經聯合對角化后難以獲取較準確的解混矩陣。f-SOBI算法將觀測信號進行反傅里葉變換后再求相關矩陣，變換后的相關函數第行和第列可表示為

2 問題提出

在fMRI數據分離算法中，第個被試者的BOLD信號矩陣(m)是一混合模型，省略字母后可表示為

由圖1可知，觀測數據其實是時間進程(time courses，TCs)與空間腦區映射(spatial maps，SMs)的乘積，而中的每一個腦區信號其實均為一稀疏的信號，經矩陣混合后得到的信號仍有可能為稀疏信號，若直接利用SOBI算法計算其相關函數[11]，會導致其包含的源信號信息不夠豐富，從而難以較準確分離源信號。f-SOBI利用反傅里葉變換對fMRI數據進行反稀疏變換[11]，解決了稀疏的fMRI數據相關函數信息不夠豐富的問題。然而，傅里葉變換的母函數為復數，與其相乘需要分別與實部和虛部各相乘一次。本文將從變換復雜度的問題出發，考慮另外一種變換方法，以減少乘法次數來降低計算復雜度。

圖1 fMRI盲分離算線性混合模型

3 哈達碼的二階盲分離算法

3.1 相關函數

h-SOBI利用哈達瑪變換來獲取fMRI數據的相關函數，由于哈達碼并不含有復數，且只由1或-1構成，因此在變換時只需與其實部相乘，而沒有虛部，從而可減少乘法次數，其相關函數計算步驟如下。

x和y分別是經過哈達瑪變換后的第行和第行中的第列元素，表示為

其中，u,j為×的哈達瑪變換矩陣中第行第列元素，該哈達瑪變換矩陣可表示為

其中，當=2時，哈達碼矩陣為

由式(8)和(9)可知，哈達瑪變換矩陣的元素是由1和-1構成的正交矩陣。

原相關函數R()經哈達瑪變換后的第行和第列的相關函數可表示為

將式(6)和(7)代入式(10)中，可得

由于哈達瑪變換矩陣是正交矩陣且元素之間滿足以下關系

因此，式(11)可變為

其為利用哈達瑪變換矩陣計算相關函數的表達式，由于式(13)中的u,k只有1或-1，因此可以減少相關函數求解過程的復雜度。

3.2 解混矩陣

從式(13)中獲得解混矩陣，先要對相關矩陣進行聯合對角化

以獲取酉矩陣，其中，為聯合對角化，則最后解混矩陣可通過

計算。其中，0為對信號進行白化的矩陣，(·)-1表示求逆。

3.3 h-SOBI的算法步驟

該算法先對數據進行PCA組降維，再求解降維后的數據的相關矩陣，然后通過聯合對角化相關矩陣來獲取解混矩陣，最后分離fMRI信號實現激活腦區的定位，具體算法如下：

輸入：個被試的fMRI信號(m)，=1,2,···,。

步驟1.對fMRI信號進行PCA組降維和白化處理處理后得到；

步驟2.由式(13)得到的相關矩陣()；

步驟3.由式(14)聯合對角化相關矩陣得到酉矩陣；

3.4 算法復雜度

4 仿真實驗與分析

4.1 數據來源

本文的fMRI實驗數據分別采用了仿真數據和實測數據。其中，仿真數據是由SimTB(Simulation Toolbox)工具箱仿真得到，該工具箱可在http:// mialab.mrn.org/software/simtb/index.html網站上下載，SimTB模擬創建30個腦區的fMRI數據，如圖2(a)所示。實驗中選取8個腦區作為源信號，如圖2(b)所示。在創建fMRI數據的空間源過程中，將體素個數設置為=148×148，被試者數設置為=5，對比度噪聲比(contrast to noise ratio，CNR)設置為模擬實際數據之間的空間差異性，每個被試者的腦區在正常范圍內進行水平或垂直方向上的平移、旋轉和擴展，其參數值分布見表1。其余參數，包括塊、頭部運動等模塊的設置可由SimTB中文件“experiment_params_aod.m”獲得。fMRI實測數據來自TReNDS實驗室的公開數據，可從http:// trendscenter.org/trends/software/gift/Index.html網頁獲得，該數據是通過讓被試者觀察一個棋盤圖案時所獲得的fMRI數據，其參數設置見表2。

圖2 fMRI數據的腦區模擬圖

表1 構造腦區空間差異性的參數值

表2 實測fMRI數據集的參數設置

4.2 實驗設置

實驗中將本文算法h-SOBI和傳統ICA算法和f-SOBI算法進行了比較，具體設置如下：

(1) fastICA，采用GIFT軟件中的fastICA算法對仿真數據和實測數據進行信號分離，該軟件下載地址為http://mialab.mrn.org/software/gift/index.html。將SimTB產生的仿真數據和實測數據分別輸入GIFT中，然后選取FastICA 算法對數據進行處理。其中，仿真數據和實測數據的腦區獨立分量個數分別為8和16。

(2) f-SOBI，相關函數中取值為4，其余參數采用文獻[3]的設置，將該算法嵌入GIFT軟件中，對SimTb產生的數據進行處理，分別獲取8個激活腦區的時間、空間分布圖。

(3) h-SOBI，相關函數中取值為4，其余步驟可參見表1，將該算法嵌入GIFT[12]軟件中，對SimTb產生的數據進行處理，分別獲取8個激活腦區的時間、空間分布圖。

4.3 結果分析

4.3.1 仿真數據

本小節給出了仿真數據的分離結果，為了直觀地觀測空間圖像分離的效果，圖3分別顯示了在閾值為1.5的情況下，ICA算法、f-SOBI算法以及h-SOBI算法對fMRI數據進行分離后獲得的腦區圖像。該圖像分別選取第1個、第13個、第16個、第20個腦區，圖中白色亮點及其周圍黑色區域為激活腦區。

圖3 3種算法提取的與任務有關的腦區圖

從圖3中可以看出，fastICA算法獲得的激活區大小是3種算法中最大的，f-SOBI算法獲得的激活區大小次之，而h-SOBI算法獲得的激活區大小則是最小的，說明fastICA算法在3種算法中分離效果最好，能清楚地觀察到分離的腦激活區，而h-SOBI相比其他2種算法，其分離效果較差，激活區不夠明顯。從圖3中還可以看到，腦激活區域最亮的是fastICA算法獲得的腦區，其次是f-SOBI，而h-SOBI算法獲得的腦區較暗，表明fastICA分離信號的強度要高于其他2種算法。但是，3種算法所得到的激活腦區的位置是一致的，表明3種算法均能實現激活區的定位。

為了定量的分析3種算法的fMRI信號分離性能，本實驗使用腦區定位的相對誤差(brain position relative error, BPRE)作為評判準則，主要是計算差異性極顯著的混合前與分離后腦激活區幅值最大點所對應的位置點的相對誤差，即

圖4給出了fastICA算法、f-SOBI算法以及h-SOBI算法分離fMRI信號的相對誤差，從中可以看出，fastICA算法和f-SOBI算法的腦區定位誤差均低于10%。h-SOBI算法除第5腦區的定位誤差為15%，其余腦區的誤差均低于10%，與其他2種算法的誤差較為接近。

圖4 3種算法的分離性能

本文實驗通過運行時間來評判3種算法的計算復雜度。由于對算法計算復雜度影響較大是其解混矩陣的獲取，因此表3給出了fastICA算法、f-SOBI算法和h-SOBI算法獲取解混矩陣的運行時間，3種算法均在統一的MatLab平臺下進行。

由表3可以看出，f-SOBI和h-SOBI算法獲取解混矩陣的運行時間要小于fastICA,僅是fastICA算法的25%，其中，h-SOBI算法的運行時間更短，僅是f-SOBI算法的50%，因此，本文的h-SOBI算法的計算復雜度最低，fastICA算法的計算復雜度最高。

4.3.2 實測數據

本小節給出了實測數據的腦區信號分離結果，根據文獻[13]可知，通過分離該真實數據可估計得到16個腦區分量，其中與視覺任務最為相關的2個腦區域分別是左視覺皮層和右視覺皮層，分別用藍色和紅色表示，如圖5所示。

表3 3種算法分離仿真數據的運行時間(s)

圖5 3種算法分離實測fMRI信號的腦區結果圖

圖5給出了fastICA算法、f-SOBI算法和h-SOBI算法分離實測fMRI數據得到的與任務相關度高的2個腦區分量切片圖，且給出了與之對應的時間進程曲線圖。將圖5中3種算法分離得到的腦切片圖進行對比可以看出，f-SOBI算法和h-SOBI算法分離得到的腦激活空間可以與fastICA算法分離得到的腦激活空間一一對應。在相同的閾值(1.0)下，h-SOBI算法顯示的腦激活空間范圍要小于f-SOBI算法和fastICA算法顯示的腦激活空間圖范圍，其中h-SOBI顯示的腦激活空間圖范圍更小。實驗結果表明，對于實測數據，f-SOBI算法和h-SOBI算法的分離性能與fastICA算法的分離性能相差不大。

表4為3種算法分離實測數據的運行時間，由長到短依次為fastICA，f-SOBI和h-SOBI，分別是4.025 s，1.639 s以及0.885 s，表明f-SOBI和h-SOBI的運行時間要遠小于fastICA，約為fastICA的50%，而h-SOBI的運行時間更短，約為f-SOBI的50%，與表3中的仿真數據結果基本一致。

表4 3種算法分離實測數據的運行時間(s)

5 結 論

針對傳統fMRI腦區分離算法ICA計算復雜度高的問題，本文提出了一種基于哈達瑪變換矩陣的二階盲分離算法h-SOBI，通過該算法分離fMRI數據來定位腦激活區。該方法將稀疏的腦信號變換為非稀疏信號，然后根據信號的自相關矩陣來實現腦區信號的分離，使算法的計算復雜度大大降低，而腦區的定位誤差與傳統方法較為接近。在仿真實驗中，本文分別對SimTB工具產生的fMRI仿真數據和TReNDS實驗室公開的實測數據進行了測試。對于模擬的fMRI數據，本文算法的分離誤差比fastICA算法高出4%，而運行時間卻分別只有fastICA算法和f-SOBI算法的25%和50%。對于真實的fMRI數據，h-SOBI算法分離的腦激活空間能夠與fastICA算法分離的腦激活空間一一對應，但是激活空間范圍相比于fastICA更小。在運行時間上，f-SOBI和h-SOBI的運行時間分別約為fastICA的40%和20%。

[1] 唐煥文, 潘麗麗, 唐一源. SPM的數學基礎及其在腦功能成像研究中的應用[J]. 應用基礎與工程科學學報, 2005, 13(3): 6-14.TANG H W, PAN L L, TANG Y Y. The mathematical foundation of SPM and its application in brain functional imaging[J]. Journal of Applied Basic Science and Engineering, 2005, 13(3): 6-14 (in Chinese).

[2] CALHOUN V D, ADALI T, HANSEN L K. ICA of function MRI data: an overview[J]. Avsh-Alom Elyada, 2003, 33(5): 281-288.

[3] 程明. 獨立成分分析算法應用于功能磁共振成像數據中[D]. 長春:吉林大學, 2013. CHENG M. Application of independent component analysis algorithm to functional magnetic resonance imaging data[D]. Changchun:Jilin University, 2013 (in Chinese).

[4] BELOUCHRANI A, ABED-MERAIM K, CARDOSO J-F. A blind source separation technique using second-order statistics[J]. IEEE Transactions on Signal Processing, 1997, 45(2): 434-444.

[5] NUZILLARD D, NUZILLARD J-M. Second-order blind source separation in the Fourier space of data[J]. Signal Processing, 2003, 83(3): 627-631.

[6] ALLEN E A, ERHARD E B, Yonghua WEI Y H, et al. A simulation toolbox for fMRI data: SimTB[EB/OL]. (2011-03-29) [2020-03-10]. http://mialab.mrn.org.

[7] COMON P. General linear model (GLM) applied to fMRI data analysis[J]. IEEE Transactions on Signal Processing, 2004, 36(2): 287-314.

[8] 彭堯, 熊馨. 一種改進的FastICA算法及其在fMRI數據中的仿真應用[J]. 軟件導刊, 2018, 17(4): 67-70. PENG Y, XIONG X. An improved FastICA algorithm and its simulation application in fMRI data[J]. Software Guide, 2018, 17(4): 67-70 (in Chinese).

[9] CALHOUN V D, ADALI T, PEARISON G D, et al. A method for making group inferences from functional MRI data using independent component analysis[J]. Human Brain Mapping, 2001, 14(3): 140-151.

[10] HU W W, XU G Z. DOA estimation with double L-shaped array based on Hadamard product and joint diagonalization in the presence of sensor gain-phase errors[J]. Multidimensional Systems and Signal Processing, 2019, 30: 465-491.

[11] NUZILLARD D. Adaptation de SOBI à des données fréquentielles[EB/OL]. [2020-01-19]. http://documents. irevues.inist.fr/handle/2042/13128.

[12] The GIFT Documentation Team. Group ICA/IVA of fMRI toolbox (GIFT) manual[EB/OL]. [2020-01-15]. http://mialab.mrn.org.

[13] ERHARDT E B, RACHAKONDA S, BEDRICK E J, et al. Comparison of multi‐subject ICA methods for analysis of fMRI data[J]. Human brain mapping, 2011, 32(12): 2075-2095.

A blind separation algorithm with low complexity for fMRI brain activation

CHEN An-ying, WU Hai-feng, LI Dong

(School of Electrical and Information, Yunnan Minzu University, Kunming Yunnan 650500, China)

Functional magnetic resonance imaging (FMRI) is a medical imaging technology widely employed in brain region positioning for its non-invasiveness and high spatiotemporal resolution. However, the traditional FMRI signal separation algorithm was too complex and time-consuming to effectively apply the FMRI technology to brain function research. Aiming at the computational complexity of traditional FMRI brain separation algorithms, a blind separation algorithm was proposed based on the second-order Hadamard transform. This algorithm first calculated the correlation function of the blood oxygen level dependent (BOLD) signal in the fMRI data, and then performed eigenvalue decomposition to obtain the unmixing matrix, thereby realizing the activation of brain regions. Given the composition of the Hadamard being only 1 or-1, the complexity can be reduced for the BOLD signal correlation matrix calculation. The simulation results show that compared with the independent component analysis (ICA) of high-order statistics and the Fourier transform blind separation algorithm of second-order statistics, the calculation time of this algorithm was only 25% and 50% of theirs, respectively, while the positioning error was close.

functional magnetic resonance imaging; blind separation; independent components analysis; second order blind identifiability; brain activation area

TN 911.73

10.11996/JG.j.2095-302X.2020060947

A

2095-302X(2020)06-0947-07

2020-05-09；

2020-07-16

9 May，2020；

16 July，2020

國家自然科學基金項目(61762093)；云南省應用基礎研究重點項目(2018FA036)；云南省高校智能傳感網絡及信息系統科技創新團隊；2018年云南民族大學研究生創新基金項目(2018YJCXS176)

National Natural Science Foundation of China (61762093); Yunnan Provincial Applied Fundamental Research Key Project (2018FA036); Yunnan University Intelligent Sensor Network and Information System Technology Innovation Team; 2018 Yunnan Nationalities University Graduate Innovation Fund Project (2018YJCXS176)

陳安瑩(1994-)，女，河南信陽人，碩士研究生。主要研究方向為fMRI信號處理。E-mail：857978569@qq.com

CHEN An-ying (1994-), female, master student. Her main research interest covers fMRI signal processing. E-mail：857978569@qq.com

吳海鋒(1977-)，男，云南昆明人，教授，博士。主要研究方向為信號處理、機器學習和RFID。E-mail：whf5469@gmail.com

WU Hai-feng (1977-), male, professor, Ph.D. His main research interest cover signal processing, machine learning and RFID. E-mail：whf5469@gmail.com