李 凌 金貞蘭 李 斌

1)(電子科技大學(xué)生命科學(xué)與技術(shù)學(xué)院，神經(jīng)信息教育部重點(diǎn)實(shí)驗(yàn)室，成都 610054)2)(電子科技大學(xué)物理電子學(xué)院，成都 610054)(2010年7月14收到;2010年8月6日收到修改稿)

基于因子分析方法的相位同步腦電源的時(shí)-空動(dòng)力學(xué)分析*

李 凌1)金貞蘭1)李 斌2)

1)(電子科技大學(xué)生命科學(xué)與技術(shù)學(xué)院，神經(jīng)信息教育部重點(diǎn)實(shí)驗(yàn)室，成都 610054)2)(電子科技大學(xué)物理電子學(xué)院，成都 610054)(2010年7月14收到;2010年8月6日收到修改稿)

頭皮腦電時(shí)間序列的相關(guān)性是大腦皮層源的相位同步性的一種體現(xiàn)，因此對(duì)相位同步源進(jìn)行定位，同時(shí)找到源對(duì)應(yīng)的時(shí)間序列在腦成像研究領(lǐng)域具有重要意義.基于 R ssler模型提出仿真相位同步偶極子源的時(shí)間序列的方法，利用時(shí)間序列進(jìn)行同心四層球頭模型正演，獲得仿真頭皮腦電數(shù)據(jù).提出了基于最大似然因子分析的相位同步腦電源的時(shí)-空動(dòng)力學(xué)分析方法，對(duì)仿真和真實(shí)頭皮腦電數(shù)據(jù)進(jìn)行了驗(yàn)證，并與主成分分析法進(jìn)行對(duì)比.仿真實(shí)驗(yàn)結(jié)果表明:最大似然因子分析法估計(jì)的時(shí)間序列與仿真源的時(shí)間序列具有更高的相關(guān)系數(shù)，同時(shí)估計(jì)源與仿真源的位置誤差更小，并且對(duì)空間分辨率和噪聲的魯棒性強(qiáng).空間選擇性注意實(shí)驗(yàn)的腦電數(shù)據(jù)分析結(jié)果也顯示因子分析法能獲得具有生理意義的時(shí)間序列和空間分布.

腦電圖，相位同步，因子分析，主成分分析

PACS:87.19.Le，05.90.+m，05.45.Pq

1.引 言

腦電圖(electroencephalogram，EEG)是神經(jīng)電活動(dòng)集群效應(yīng)的直接記錄，是研究人腦功能的最重要技術(shù)之一，具有高時(shí)間分辨率(ms)的優(yōu)勢(shì).近年來(lái)在腦電技術(shù)與方法上的不斷創(chuàng)新和進(jìn)步，以及腦電設(shè)備的易用性、無(wú)創(chuàng)性，價(jià)格相對(duì)較低，使得EEG成為心理物理學(xué)中最常用的研究手段之一［1］.由于空間分辨率較低的原因，EEG技術(shù)在基礎(chǔ)研究和臨床應(yīng)用中存在兩個(gè)有待解決的問(wèn)題:源或等效源的定位以及大腦動(dòng)力學(xué)特性的分析.對(duì)頭皮腦電進(jìn)行大腦皮層源的定位和求解，被稱(chēng)為腦電信號(hào)的反演、成像或逆問(wèn)題［2］.由于 EEG受到采集設(shè)備的空間分辨率的限制以及大腦容積導(dǎo)體對(duì)信號(hào)產(chǎn)生的非線(xiàn)性變換的影響，逆問(wèn)題求解具有相當(dāng)難度.

目前已發(fā)展了多種估計(jì)神經(jīng)活動(dòng)源的方法，包括頭表拉普拉斯地形圖技術(shù)(scalp Laplacian mapping，LM)［3，4］，低分辨率電磁斷層成像方法(low resolution electromagnetic tomography，LORETA)［5，6］，皮層表面電位成像技術(shù)(cortical surface potential mapping，CPM)［7］，多信號(hào)分類(lèi)法(multiple signal classification，MUSIC)［8］等.上述大部分方法是基于瞬時(shí)腦電電位的反演方法，只利用到了數(shù)據(jù)的空間信息，只有MUSIC方法是基于時(shí)間序列的反演方法，同時(shí)利用了數(shù)據(jù)的時(shí)間和空間信息.然而MUSIC算法卻對(duì)噪聲敏感，并且無(wú)法處理相關(guān)源的定位.

頭皮腦電時(shí)間序列的相關(guān)性是大腦皮層源的相位同步性的一種體現(xiàn).例如，在視覺(jué)刺激下會(huì)同時(shí)激活多個(gè)大腦區(qū)域，它們互相之間需要進(jìn)行信息的傳遞和交流，因此可能存在多個(gè)相位同步相關(guān)源.因而對(duì)相位同步源進(jìn)行定位、同時(shí)找到源對(duì)應(yīng)的時(shí)間序列(即大腦時(shí)-空動(dòng)力學(xué)分析)在腦成像研究領(lǐng)域具有重要意義.近年來(lái)人們已經(jīng)把各種線(xiàn)性與非線(xiàn)性動(dòng)力學(xué)的理論與方法應(yīng)用于腦電信號(hào)的時(shí)-空分析中，主要的方法包括線(xiàn)性技術(shù):分析信號(hào)之間的互相關(guān)或相干估計(jì);非線(xiàn)性技術(shù):相關(guān)維、互信息、廣義同步和復(fù)雜度［9，10，11］，以及相空間重構(gòu)［12］、相位同步分析方法［13］等.

基于數(shù)據(jù)驅(qū)動(dòng)的時(shí)-空分析方法:主成分分析法(principal components analysis，PCA)是降低 EEG 數(shù)據(jù)維數(shù)最常用的方法之一，被廣泛應(yīng)用于腦電信號(hào)處理領(lǐng)域.PCA可以把EEG信號(hào)分解成一組相互正交的時(shí)-空成分［14，15］，目前 PCA 被認(rèn)為是確定 EEG時(shí)間序列及其對(duì)應(yīng)的頭皮位置的一種客觀(guān)手段.然而，利用PCA分析EEG數(shù)據(jù)也存在一些爭(zhēng)論，例如Wood和McCarthy在一個(gè)仿真研究中證明了可能存在方差閾值選取不當(dāng)?shù)那闆r，以及處理相關(guān)源存在一定局限性［16］.

本文發(fā)展了一個(gè)相關(guān)源的動(dòng)力學(xué)分析新技術(shù)，首次提出基于最大似然因子分析(factor analysis，F(xiàn)A)的相位同步腦電源的時(shí)-空動(dòng)力學(xué)分析方法，對(duì)仿真和真實(shí)頭皮腦電數(shù)據(jù)進(jìn)行了驗(yàn)證，并與主成分分析法進(jìn)行結(jié)果的對(duì)比.目的不僅要對(duì)相位同步源進(jìn)行定位，同時(shí)還要找到這些等效源對(duì)應(yīng)的時(shí)間序列.為了驗(yàn)證方法的有效性，本文提出了一種基于R ssler模型仿真相位同步偶極子源的時(shí)間序列的方法，并利用偶極子時(shí)間序列進(jìn)行同心四層球頭模型正演，以獲得仿真頭皮腦電數(shù)據(jù).比較了隨機(jī)噪聲、空間相關(guān)噪聲和不同導(dǎo)聯(lián)數(shù)對(duì)最大似然因子分析和主成分分析方法的影響.最后，把因子分析法應(yīng)用到空間選擇性注意腦電數(shù)據(jù)上，獲得相應(yīng)的時(shí)間序列和空間分布.

2.腦電數(shù)據(jù)仿真和真實(shí)數(shù)據(jù)來(lái)源

2.2.腦電數(shù)據(jù)的仿真方法

2.1.1.相位同步偶極子源的時(shí)間序列仿真

為了描述成對(duì)神經(jīng)細(xì)胞群活動(dòng)的關(guān)系，利用兩對(duì)不同耦合強(qiáng)度的R ssler振子來(lái)產(chǎn)生四個(gè)偶極子源的時(shí)間序列.兩個(gè)耦合的 R ssler振蕩器(1和2)模型如下:

式中，參數(shù)ω1，2表示兩個(gè)振蕩器的振蕩頻率，C表示它們的耦合強(qiáng)度.

四個(gè)仿真偶極子源的時(shí)間過(guò)程由 x1，2的解確定.由于腦電波是以8—13 Hz的 α波為主，設(shè)兩個(gè)振蕩器的振蕩頻率分別為10+0.024 Hz和10—0.024 Hz，耦合強(qiáng)度 C分別為0.011和0.035.當(dāng) C為0.035時(shí)，兩個(gè)振蕩器達(dá)到相位同步的狀態(tài)［17］，即相位差小于某個(gè)常數(shù)，達(dá)到部分相關(guān);當(dāng) C為0.011時(shí)，表示兩個(gè)振蕩器是非常弱的耦合，即沒(méi)有相關(guān)性.假定仿真源 s1和 s2的 C 為0.011，s3和 s4的C為0.035，則仿真源 s3和s4有較強(qiáng)的相位同步性.圖1顯示了四個(gè)仿真源的時(shí)間序列，橫坐標(biāo)表示時(shí)間過(guò)程，縱坐標(biāo)表示幅度大小，但兩者都屬于無(wú)量綱變量.

圖1 四個(gè)仿真源的時(shí)間序列示意圖 (a)相位不同步的s1(實(shí)線(xiàn))和s2(虛線(xiàn))序列;(b)相位同步的 s3(實(shí)線(xiàn))和 s4(虛線(xiàn))序列(s1和 s2，s1和 s3，s1和 s4，s2和 s3，s2和 s4，s3和 s4的相位同步指數(shù)分別為 0.2943，0.1226，0.1120，0.2584，0.2796，0.9794)

2.1.2.正演頭皮數(shù)據(jù)方法

采用同心四層球頭模型，傳導(dǎo)特性和半徑參數(shù)取自的文獻(xiàn)［18］.四層球的半徑大小分別為7.9，8.1，8.5和 8.8 cm，表示內(nèi)部大腦半球的外部邊界、腦脊液層、顱骨層和頭皮層，傳導(dǎo)率分別為 0.461，1.39，0.0058 和 0.461 A(V m)－1，詳細(xì)的正問(wèn)題算法可以參考文獻(xiàn)［19］.正問(wèn)題算法中采用笛卡爾坐標(biāo)系，假設(shè)左、右耳垂連線(xiàn)中點(diǎn)為坐標(biāo)原點(diǎn)，原點(diǎn)指向右耳垂的方向?yàn)?x軸方向，原點(diǎn)指向鼻根的方向?yàn)?y軸方向，+z軸方向?yàn)榇怪庇趚y軸平面并且由原點(diǎn)指向頭頂?shù)姆较?

假設(shè)在頭球模型內(nèi)有四個(gè)+z方向的偶極子源(s1，s2，s3，s4)，它們的笛卡爾坐標(biāo)(x，y，z)分別為(－3.5，5，3.5)，(3.5， －5，3.5)，(3.5，5，3.5)和(－3.5，－5，3.5)cm，利用相位同步時(shí)間序列調(diào)制偶極子源的大小，則能仿真出相應(yīng)的頭皮腦電序列.四個(gè)單位偶極子的正演頭皮電位地形圖如圖2所示，電極導(dǎo)聯(lián)數(shù)為128.

圖2 四個(gè)單位仿真源的頭皮電位地形圖

2.2.真實(shí)腦電數(shù)據(jù)來(lái)源

實(shí)驗(yàn)范式為空間選擇性注意實(shí)驗(yàn)［20］，14名健康、右利手的被試(5名女性)參加了實(shí)驗(yàn)，年齡為20—23歲，實(shí)驗(yàn)前都簽署了同意書(shū)，參與者聽(tīng)覺(jué)和視覺(jué)正常，實(shí)驗(yàn)后給予了金錢(qián)的補(bǔ)償.刺激由兩種大小的圓形黑白棋盤(pán)格組成，大小分別為2.2°和1.4°，大的出現(xiàn)概率為80%，是標(biāo)準(zhǔn)刺激，小的出現(xiàn)概率為20%，是靶刺激.刺激出現(xiàn)在水平左或右視野，離中心點(diǎn)距離為5.2°，兩種刺激在兩個(gè)視野位置上隨機(jī)出現(xiàn)，刺激持續(xù)時(shí)間為50 ms，刺激間隔為350—550 ms.

實(shí)驗(yàn)過(guò)程中，被試坐在光線(xiàn)較暗、聲音和電磁屏蔽的房間里，要求被試雙眼注視屏幕中心的綠十字架(1.4°)，對(duì)注意視野中的靶刺激進(jìn)行按鍵反應(yīng).一共運(yùn)行16組，因而一共各有320個(gè)注意左或右視野，出現(xiàn)在左或右視野的標(biāo)準(zhǔn)刺激.采用128導(dǎo)(EGI System)腦電采集系統(tǒng)記錄EEG信號(hào)，參考電極為頭頂電極，采樣率為500 Hz，帶通濾波為0.1—48 Hz，分段時(shí)間范圍為刺激前200 ms—刺激后350 ms.幅度大于100 μV的眼動(dòng)、肌電等噪聲偽跡自動(dòng)進(jìn)行剔除，最后進(jìn)行平均疊加、帶通濾波(1—31 Hz).得到四組有待分析的腦電事件相關(guān)電位(event-related potential，ERP)數(shù)據(jù)，條件分別為:注意與刺激都在右視野(簡(jiǎn)稱(chēng)RA);注意與刺激都在左視野(LA);注意在右、刺激在左視野(RU);注意在左、刺激在右視野(LU).

3.因子分析方法和數(shù)據(jù)分析結(jié)果

3.1.最大似然因子分析方法

EEG信號(hào)一般少則32導(dǎo)(變量數(shù))，多則256導(dǎo)，且鄰近導(dǎo)聯(lián)信號(hào)一般具有強(qiáng)的相關(guān)性，因而具有相當(dāng)?shù)娜哂喽?因子分析與主成分分析具有類(lèi)似的功能，即對(duì)高維變量數(shù)據(jù)的簡(jiǎn)化，適合用于腦電數(shù)據(jù)分析中.然而它們也存在明顯的差異:主成分分析僅僅是原始變量線(xiàn)性組合產(chǎn)生的新綜合變量;而因子分析則是通過(guò)研究眾多變量之間的內(nèi)部依賴(lài)關(guān)系，并用少數(shù)幾個(gè)假想變量來(lái)表示其基本的數(shù)據(jù)結(jié)構(gòu)，因此需要構(gòu)造因子模型.即因子分析能通過(guò)EEG測(cè)量數(shù)據(jù)構(gòu)造出潛在因子(不可測(cè)量的腦內(nèi)源數(shù)據(jù))，與主成分分析法比較，因子分析模型更加適合腦電信號(hào)的時(shí)-空動(dòng)力學(xué)分析.

假設(shè)采集到的EEG信號(hào)表示為 X1，X2，…，Xp，其中Xi表示EEG數(shù)據(jù)中第i個(gè)電極導(dǎo)聯(lián)信號(hào)，因子分析模型為

式中的 F1，F(xiàn)2，…，F(xiàn)m，稱(chēng)為公共因子(或因子得分)，通常情況下因子的個(gè)數(shù) m小于變量個(gè)數(shù) p，εi稱(chēng)為第i個(gè)導(dǎo)聯(lián)的特殊因子.由aij組成的矩陣稱(chēng)為因子載荷矩陣A，因子模型的矩陣表示式為

為了估計(jì)出公共因子F(潛在變量)，首先需要估計(jì)出因子載荷矩陣A，本文采用最大似然法對(duì)A進(jìn)行估計(jì)，詳細(xì)算法步驟參考文獻(xiàn)［21，22］.

下一步則是估計(jì)公共因子F，本研究采用1939年Thomson提出的湯姆森回歸法進(jìn)行公共因子的估計(jì)，假設(shè)公共因子可對(duì)p個(gè)原始變量作回歸

即得公共因子的估計(jì)公式為

其中R是X的相關(guān)系數(shù)矩陣.

因子分析的目標(biāo)之一就是要對(duì)所提取的潛在變量F和因子載荷矩陣A的實(shí)際含義進(jìn)行合理解釋?zhuān)狙芯恐蠪代表的就是腦電源的時(shí)間序列，而A代表的是相應(yīng)的空間分布.直接計(jì)算的因子載荷矩陣往往不能突出公共因子，就需要進(jìn)行因子旋轉(zhuǎn)，使得每個(gè)變量在盡可能少的公共因子上有較大載荷，而在其余公共因子上的載荷較小或中等大小.因子旋轉(zhuǎn)采用了方差最大法，使得載荷矩陣中每列或每行元素的平方值向0和1兩極分化.還有一個(gè)關(guān)鍵問(wèn)題，即因子個(gè)數(shù)m的確定，與主成分分析法類(lèi)似，可以采用方差閾值法，取方差較大的成分.仿真計(jì)算中，方差閾值設(shè)為2，而PCA的主成分個(gè)數(shù)主動(dòng)選取前4個(gè)成分.實(shí)際數(shù)據(jù)中因子個(gè)數(shù)的選擇一方面要參考方差的大小，另一方面還要觀(guān)察對(duì)應(yīng)的公共因子是否能代表一定的生理、心理意義.

綜上，基于因子分析的時(shí)-空動(dòng)力學(xué)分析步驟如下:

1)獲得電極位置分布和頭皮腦電數(shù)據(jù);

2)對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，建立相關(guān)系數(shù)矩陣R;

3)利用最大似然估計(jì)因子載荷矩陣，得到因子模型;

4)利用方差最大法對(duì)因子載荷矩陣進(jìn)行因子旋轉(zhuǎn)，獲得m個(gè)具有生理意義的空間分布Ai;

5)估計(jì)m個(gè)公共因子 Fi，獲得m個(gè)估計(jì)的時(shí)間序列.

為了衡量仿真效果，用兩個(gè)參數(shù)評(píng)估方法的有效性:

1)估計(jì)的時(shí)間序列x與仿真偶極子時(shí)間序列y間的相關(guān)系數(shù)

其中E表示數(shù)學(xué)期望，var表示向量的方差.

2)距離(distances，Ds):大因子載荷的中心在頭皮表面的投影位置與對(duì)應(yīng)的仿真源的頭皮映射位置的距離.向量Pi表示第i個(gè)時(shí)間序列的大因子載荷的中心在頭皮表面的投影位置，向量Ps表示對(duì)應(yīng)的仿真源的頭皮映射位置.因子分析中載荷的閾值取為0.8，主成分分析中載荷的閾值取為0.3.

3.2.數(shù)據(jù)分析結(jié)果

3.2.1.仿真數(shù)據(jù)分析結(jié)果

由于實(shí)際腦電測(cè)量都會(huì)不同程度地引入各類(lèi)噪聲，因此在仿真中加入了不同噪聲信號(hào)比(noise signal ratio，NSR)的噪聲數(shù)據(jù)，NSR表示噪聲方差與信號(hào)方差的比值.我們一共仿真了三個(gè)因素的影響，分別為時(shí)-空隨機(jī)噪聲(NSR為 0.2，0.4，0.6，0.8，1，五種情況);在時(shí)-空隨機(jī)噪聲(NSR=0.2)的基礎(chǔ)上加上空間相關(guān)、時(shí)間隨機(jī)的噪聲(NSR 為 0.2，0.4，0.6，0.8，1，五種情況);空間分辨率不同的情況(電極導(dǎo)聯(lián)數(shù)分別為32，64，128三種情況).

以電極導(dǎo)聯(lián)為64導(dǎo)、含有時(shí)-空隨機(jī)噪聲(NSR為0.2)數(shù)據(jù)為例，有4個(gè)公共因子的方差大于2，因此公共因子m個(gè)數(shù)為4，圖3顯示了四個(gè)成分的因子載荷Ai，表示估計(jì)源的空間分布.圖4顯示了四個(gè)公共因子Fi，表示估計(jì)的時(shí)間序列.估計(jì)源與仿真源之間的相關(guān)系數(shù)和位置距離如表1所示，表1也相應(yīng)列出了主成分分析的結(jié)果.

圖3 因子載荷Ai(i=1，2，3，4)的空間分布圖，分別代表四個(gè)估計(jì)源的空間分布地形圖

表1 四個(gè)估計(jì)源與仿真源之間的相關(guān)系數(shù)和歸一化位置距離

圖3顯示對(duì)源的空間分布的估計(jì)與仿真源的空間分布非常一致，采用距離參數(shù)與主成分分析法比較發(fā)現(xiàn)兩種方法對(duì)位置的估計(jì)偏差大約在4.88%和6.95%(見(jiàn)表1的平均距離)，比較接近，偏差都較小.圖4顯示估計(jì)源的時(shí)間序列能很好地保持仿真源的相位同步與不同步特性，效果顯著.觀(guān)察表1的平均相關(guān)系數(shù)，發(fā)現(xiàn)估計(jì)的時(shí)間序列與仿真的時(shí)間序列相關(guān)性達(dá)到0.9682，而主成分分析法的相關(guān)系數(shù)只有0.8590，說(shuō)明因子分析方法對(duì)于相位同步源的估計(jì)效果遠(yuǎn)比主成分分析法要好.

圖4 四個(gè)估計(jì)的時(shí)間序列圖 (a)沒(méi)有相位同步特性的F1(實(shí)線(xiàn))和F2(虛線(xiàn))序列;(b)具有相位同步特性的F3(實(shí)線(xiàn))和F4(虛線(xiàn))序列(F1和 F2，F(xiàn)1和 F3，F(xiàn)1和 F4，F(xiàn)2和 F3，F(xiàn)2和 F4，F(xiàn)3和F4的 相 位 同 步 指 數(shù) 分 別 為 0.2657，0.1267，0.0674，0.1049，0.1882，0.9485)

進(jìn)一步對(duì)不同時(shí)-空隨機(jī)噪聲(NSR=0.4，0.6，0.8，1)、不同導(dǎo)聯(lián)數(shù)的仿真信號(hào)進(jìn)行計(jì)算，并用四個(gè)仿真源與估計(jì)源的相關(guān)系數(shù)和距離的平均值進(jìn)行衡量.圖5顯示了導(dǎo)聯(lián)數(shù)分別為32，64和128的平均相關(guān)系數(shù)與NSR的關(guān)系，同時(shí)給出PCA結(jié)果.圖5顯示FA比PCA方法估計(jì)的時(shí)間序列更加精確;無(wú)論是FA還是PCA方法對(duì)時(shí)-空隨機(jī)噪聲都不太敏感;FA方法對(duì)導(dǎo)聯(lián)數(shù)不敏感，即使在空間分辨率較低的32導(dǎo)聯(lián)條件下，總平均相關(guān)系數(shù)都能達(dá)到0.94;而PCA方法對(duì)導(dǎo)聯(lián)數(shù)則非常敏感，當(dāng)導(dǎo)聯(lián)數(shù)為64和128導(dǎo)時(shí)，總平均相關(guān)系數(shù)均達(dá)到0.84，但當(dāng)導(dǎo)聯(lián)數(shù)為32時(shí)，總平均相關(guān)系數(shù)只有0.69.

表2列出了不同導(dǎo)聯(lián)數(shù)和不同NSR對(duì)估計(jì)源與仿真源距離的影響，結(jié)果顯示當(dāng)導(dǎo)聯(lián)數(shù)為64和128時(shí)，兩種方法對(duì)位置的估計(jì)偏差為4.72%和8.08%，比較接近，對(duì)噪聲不敏感，偏差都較小;當(dāng)導(dǎo)聯(lián)數(shù)為32時(shí)，因子分析法的估計(jì)偏差(10.59%)遠(yuǎn)遠(yuǎn)小于主成分分析法(44.72%)，說(shuō)明主成分分析法對(duì)空間分辨率非常敏感，空間分辨率低會(huì)使得估計(jì)位置產(chǎn)生相當(dāng)大的偏差;當(dāng)導(dǎo)聯(lián)數(shù)為32時(shí)，因子分析法的估計(jì)偏差隨著噪聲的增加而有所提高.

圖5 FA和PCA處理含有時(shí)-空隨機(jī)噪聲數(shù)據(jù)，估計(jì)的時(shí)間序列與仿真時(shí)間序列的平均相關(guān)系數(shù)圖

表2 在不同時(shí)-空隨機(jī)噪聲和導(dǎo)聯(lián)數(shù)下，估計(jì)源與仿真源之間的平均歸一化位置距離

為了在頭皮EEG中加入一個(gè)空間上相關(guān)、時(shí)間上隨機(jī)的噪聲，我們假設(shè)在腦內(nèi)有一個(gè)時(shí)間序列是隨機(jī)的噪聲源，如果噪聲源的位置遠(yuǎn)離這四個(gè)仿真源，即使NSR為1，對(duì)估計(jì)的影響都是非常小的.因此，本文噪聲源的位置設(shè)在與第一個(gè)仿真源(s1)相同的坐標(biāo)位置上，在時(shí)-空隨機(jī)噪聲(NSR=0.2)的基礎(chǔ)上再加上空間相關(guān)、時(shí)間隨機(jī)的噪聲(NSR為0.2，0.4，0.6，0.8，1);這樣就構(gòu)造出了含有空間上相關(guān)、時(shí)間上隨機(jī)的噪聲的數(shù)據(jù).圖6顯示了導(dǎo)聯(lián)數(shù)分別為32，64和128的平均相關(guān)系數(shù)與空間相關(guān)NSR的關(guān)系，同時(shí)給出PCA結(jié)果.

圖6顯示FA比PCA方法估計(jì)的時(shí)間序列更加精確;無(wú)論是FA還是PCA方法均對(duì)空間相關(guān)噪聲比較敏感，隨著噪聲的增加，估計(jì)的效果越來(lái)越差;FA方法對(duì)導(dǎo)聯(lián)數(shù)不敏感，32，64，128導(dǎo)總平均相關(guān)系數(shù)分別為0.857，0.852，0.836;而 PCA方法對(duì)導(dǎo)聯(lián)數(shù)非常敏感，32，64，128導(dǎo)總平均相關(guān)系數(shù)分別為 0.688，0.777，0.803.

圖6 FA和PCA處理含有空間相關(guān)噪聲數(shù)據(jù)，估計(jì)的時(shí)間序列與仿真時(shí)間序列的平均相關(guān)系數(shù)圖

表3列出了不同導(dǎo)聯(lián)數(shù)和不同空間相關(guān)NSR對(duì)估計(jì)源與仿真源距離的影響，結(jié)果顯示當(dāng)導(dǎo)聯(lián)數(shù)為64和128時(shí)，兩種方法對(duì)位置的估計(jì)偏差為5.01%和8.93%，比較接近，對(duì)噪聲不敏感，偏差都較小;當(dāng)導(dǎo)聯(lián)數(shù)為32時(shí)，因子分析法的估計(jì)偏差(6.50%)遠(yuǎn)遠(yuǎn)小于主成分分析法(31.85%)，與前類(lèi)似，說(shuō)明主成分分析法對(duì)空間分辨率非常敏感.

表3 在不同空間相關(guān)噪聲和導(dǎo)聯(lián)數(shù)下，估計(jì)源與仿真源之間的平均歸一化位置距離

3.2.2.真實(shí)數(shù)據(jù)分析結(jié)果

四種條件下(RA，LA，RU，LU)都是有4個(gè)公共因子的方差大于2，因而m取4.以L(fǎng)A條件為例，四個(gè)空間分布如圖7所示.主要的活動(dòng)腦源位于刺激對(duì)側(cè)枕部(A1)，額-頂葉區(qū)域(A2)，刺激同側(cè)枕部(A3)，中線(xiàn)枕-額區(qū)域(A4).RA的空間分布與LA類(lèi)似，而非注意狀態(tài)(LU和RU)的第一個(gè)公共因子位于額-頂葉區(qū)域(與圖7的 A2分布類(lèi)似)，第二個(gè)公共因子位于刺激對(duì)側(cè)枕部(與圖7的A1分布類(lèi)似)，其他兩個(gè)成分與注意狀態(tài)分布對(duì)應(yīng)類(lèi)似.這些活動(dòng)源都具有一定的生理意義:視覺(jué)注意任務(wù)主要誘發(fā)的區(qū)域是視覺(jué)皮層(A1);額-頂葉區(qū)域(A2)與空間注意的轉(zhuǎn)移有關(guān);刺激同側(cè)視覺(jué)皮層的激活(A3)可能是刺激對(duì)側(cè)枕葉皮層信號(hào)的一個(gè)傳導(dǎo)結(jié)果;枕-額區(qū)域(A4)反映了一種遠(yuǎn)程連接，可能與注意控制有關(guān).

發(fā)現(xiàn)注意條件下的第一個(gè)公共因子與非注意條件下的第二個(gè)公共因子具有很高的相似性(LA與RU的相關(guān)系數(shù)達(dá)到0.8712，RA與LU的相關(guān)系數(shù)達(dá)到0.9649).圖8顯示了這兩對(duì)具有高相似性的時(shí)間序列.

圖8的時(shí)間序列顯示包含兩個(gè)成分:P1(第一個(gè)正波，80—150 ms)和 N1(第一個(gè)負(fù)波，150—200 ms)，與先前的研究結(jié)果一致［20］，是與外界視覺(jué)刺激相關(guān)的外源性成分.注意的波幅顯然大于非注意條件的波幅，清楚地顯示了一個(gè)眾所周知的現(xiàn)象:刺激出現(xiàn)在注意位置會(huì)在大腦對(duì)側(cè)后部皮層誘發(fā)出更大的 P1和 N1成分［23］.

圖7 LA條件下，四個(gè)因子載荷Ai空間分布圖(i=1，2，3，4)，分別代表四個(gè)估計(jì)源的空間分布地形圖

注意和非注意條件下的刺激同側(cè)視覺(jué)皮層(A3)對(duì)應(yīng)的公共因子之間也具有相似性，相關(guān)系數(shù)平均達(dá)到0.7884.而其他對(duì)公共因子之間的相關(guān)系數(shù)較低，說(shuō)明在注意和非注意狀態(tài)下，位于較高皮層的活動(dòng)源(包括頂葉、額葉區(qū)域)存在較大時(shí)間序列的差異，暗示兩種狀態(tài)下，額-頂網(wǎng)絡(luò)工作的機(jī)理可能是不同的.

圖8 注意與非注意條件下的估計(jì)源的時(shí)間過(guò)程 (a)刺激出現(xiàn)在左邊;(b)刺激出現(xiàn)在右邊

4.結(jié) 論

對(duì)比因子分析和主成分分析的相位同步腦電源的時(shí)-空動(dòng)力學(xué)分析結(jié)果，對(duì)仿真數(shù)據(jù)，無(wú)論是含有時(shí)-空隨機(jī)噪聲，還是空間相關(guān)噪聲，因子分析估計(jì)的時(shí)間序列準(zhǔn)確性都明顯高于主成分分析，說(shuō)明因子分析法具有更強(qiáng)的抗噪性.由于腦電測(cè)量設(shè)備導(dǎo)聯(lián)數(shù)不統(tǒng)一，科研機(jī)構(gòu)一般采用64導(dǎo)或128導(dǎo)腦電儀，具有較高的空間分辨率，兩種方法對(duì)源位置的估計(jì)誤差相似;然而在臨床醫(yī)院，一般采用32導(dǎo)腦電設(shè)備，具有較低空間分辨率，主成分分析的源位置估計(jì)具有40%左右的偏差，效果顯著比因子分析法差.因而因子分析對(duì)空間分辨率不敏感，具有更強(qiáng)的魯棒性和適應(yīng)性.把因子分析法應(yīng)用于真實(shí)腦電數(shù)據(jù)中，估計(jì)的時(shí)間序列和空間分布具有相應(yīng)的神經(jīng)生理學(xué)意義.

綜上所述，提出了一種基于R ssler模型仿真相位同步偶極子源時(shí)間序列的方法，并利用偶極子時(shí)間序列進(jìn)行同心四層球頭模型正演，獲得仿真頭皮腦電數(shù)據(jù).基于因子分析法的相位同步腦電源的時(shí)-空動(dòng)力學(xué)分析不僅能對(duì)相位同步源進(jìn)行定位，同時(shí)還能獲得這些等效源所對(duì)應(yīng)的時(shí)間序列.仿真和真實(shí)數(shù)據(jù)分析結(jié)果都表明了該方法的有效性和適用性，能較準(zhǔn)確地獲得時(shí)間序列和空間分布，為神經(jīng)動(dòng)力學(xué)分析提供了一個(gè)有效的工具.

［1］ Gevins A 1998 Electroenceph.Clin.Neurophysiol.106 165

［2］ Michel C M，Murray M M，Lantz G，Gonzalez S，Spinelli L，Rolando Grave de Peralta 2004 Clin.Neurophysiol.115 2195

［3］ Nunez P L，Silibertein R B，Cdush P J，Wijesinghe R S，Westdrop A F，Srinivasan R 1994 Electroenceph. Clin.Neurophysiol.90 40

［4］ Yao D 2002 Clin.Neurophysol.113 956

［5］ Pascual-Marqui R D，Michel C M，Lehmann D 1994 Int.J.Psychophysiol.18 49

［6］ Yao D，He B 2001 Ann.Biomed.Eng.29 1019

［7］ Yao D 2001 Phys.Med.Biol.46 3177

［8］ Mosher J C，Leahy R M 1998 IEEE Trans.Biomed.Eng.45 1342

［9］ Stam C J，van Dijk B W 2002 Phys.D 163 236

［10］ Fang X L，Jiang Z L 2007 Acta.Phys.Sin.56 7330(in Chinese)［方小玲、姜宗來(lái)2007物理學(xué)報(bào)56 7330］

［11］ Meng Q F，Zhou W D，Chen Y H，Peng Y H 2010 Acta.Phys.Sin.59 123(in Chinese)［孟慶芳、周衛(wèi)東、陳月輝、彭玉華2010物理學(xué)報(bào) 59 123］

［12］ Xie Y，Xu J X，Yang H J，Hu S J 2002 Acta.Phys.Sin.51 205［謝 勇、徐健學(xué)、楊紅軍、胡三覺(jué) 2002物理學(xué)報(bào) 51 205］

［13］ Benjamin L，Jurgen S，Andreas S 2006 J.Physiol.99 66

［14］ Van Boxtel G J M 1998 Behav.Res.Methods.Instrum.Comput.30 87

［15］ Kayser J，Tenke C E 2003 Clin.Neurophysiol.114 2307

［16］ Wood C C， McCarthy G 1984 Electroencephalogr. Clin.Neurophysiol.59 249

［17］ Rosenblum M G，Pikovsky A S，Kurths J 1996 Phys.Rev.Lett.76 1804

［18］ Cuffin B N，Cohen D 1997 IEEE Trans.Biomed.Eng.24 372

［19］ Yao D 2003 Phys.Med.Biol.48 1997

［20］ Li L，Yao D，Yin G 2009 Brain Res.1282 84

［21］ J reskog K G 1967 Psychometrika.32 443

［22］ Rubin D B，Thayer D T 1982 Psychometrika.47 69

［23］ Clark V，Hillyard S A 1996 J.Cogn.Neurosci.8 387

PACS:87.19.Le，05.90.+m，05.45.Pq

Spatiotemporal dynamic analysis of phase synchronized sources based on factor analysis*

Li Ling1)Jin Zhen-Lan1)Li Bin2)
1)(Key Laboratory for NeuroInformation of Ministry of Education，School of Life Science and Technology，University of Electronic Science and Technology of China，Chengdu 610054，China)2)(School of Physical Electronics，University of Electronic Science and Technology of China，Chengdu 610054，China)(Received 14 July 2010;revised manuscript received 6 August 2010)

Correlations in time courses of scalp electroencephalogram(EEG)may be represented by the phase synchronization in cerebral cortex sources to a certain some degree.Therefore，it is very important to localize the sources of phase synchronization and find corresponding time courses in the brain imaging study.Based on coupled R ssler oscillators with different coupling strengths，we propose a new method of simulating phase synchronized dipole sources and use a concentric 4-sphere head model to obtain simulation forward scalp EEG data.In addition，we propose the spatiotemporal dynamic analysis of phase synchronized sources based on the maximum likelihood factor analysis，verify the simulated and real scalp EEG data，and further compare the results with those of principal component analysis.Simulation results demonstrate that time courses estimated by maximum likelihood factor analysis have higher correlation with simulated sources，and less locational error between estimated sources and simulated sources.Factor analysis shows a better robust to the spatial resolution and the noise than principal component analysis.Furthermore，real data from spatial attention experiments show that factor analysis is capable of obtaining time courses and spatial distribution under the physiological base.

electroencephalogram，phase synchronization，factor analysis，principal component analysis

*國(guó)家自然科學(xué)基金青年科學(xué)基金(批準(zhǔn)號(hào):30800242)資助的課題.

.E-mail:libin@uestc.edu.cn

*Project supported by the Young Scientists Fund of the National Natural Science Foundation of China(Grant No.30800242). Corresponding author.E-mail:libin@uestc.edu.cn