陰桂梅 王千山 姚 蓉 李海芳
(1.太原師范學(xué)院計(jì)算機(jī)系,山西晉中030619;2.太原理工大學(xué)信息與計(jì)算機(jī)學(xué)院,山西太原030024)
宇航員面臨著我們大多數(shù)人永遠(yuǎn)都不會(huì)經(jīng)歷的精神和心理挑戰(zhàn),自從美國女宇航員莉薩·諾瓦克事件后,NASA 一直致力于研究降低宇航員在太空飛行期間的精神健康問題風(fēng)險(xiǎn)。航天器的工程技術(shù)實(shí)現(xiàn)能力是實(shí)現(xiàn)載人航天任務(wù)的基礎(chǔ),但宇航員需面臨新的心理學(xué)和精神病學(xué)的挑戰(zhàn)是主要的限制性因素。為此,本文將從宇航員的視角,運(yùn)用拓?fù)鋽?shù)據(jù)分析中的持續(xù)同調(diào)理論,以精神分裂癥數(shù)據(jù)為例,分析其復(fù)雜腦網(wǎng)絡(luò)的持續(xù)拓?fù)涮卣?,以便能夠?yàn)樘崆邦A(yù)防和檢測(cè)宇航員精神疾病提供有效的腦網(wǎng)絡(luò)指標(biāo)。
拓?fù)鋽?shù)據(jù)分析(TDA)[1]是一個(gè)數(shù)據(jù)分析、代數(shù)拓?fù)洹⒂?jì)算幾何、計(jì)算機(jī)科學(xué)、統(tǒng)計(jì)等多領(lǐng)域相關(guān)的一個(gè)領(lǐng)域。TDA 的主要目標(biāo)是利用幾何學(xué)和拓?fù)鋵W(xué)的思想研究數(shù)據(jù)的定性特征。為了達(dá)到這個(gè)目標(biāo),TDA 需要精確的定義定性特性,還有在具體實(shí)踐應(yīng)用中的計(jì)算工具,以及保證這些特性穩(wěn)定性、健壯性的理論。解決這幾個(gè)問題的一種方法就是TDA 中的持久同源性(PH)理論。將PH理論運(yùn)用于大腦網(wǎng)絡(luò)分析是當(dāng)前正在興起的一個(gè)研究方向。
對(duì)大腦成像數(shù)據(jù)處理和分析時(shí),一般都是通過生成表示節(jié)點(diǎn)之間連接強(qiáng)度的矩陣,然后選擇合適的閾值對(duì)矩陣二值化最終生成鄰接矩陣來構(gòu)建腦網(wǎng)絡(luò)。閾值的選擇在網(wǎng)絡(luò)構(gòu)建中起著重要作用,因?yàn)樗绊懢W(wǎng)絡(luò)連接的密度和網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)。通常有三種網(wǎng)絡(luò)閾值化方法[2]。
1) 選取單個(gè)閾值。一般選某個(gè)連接密度作為單個(gè)閾值,當(dāng)閾值選取2lgN/N時(shí),網(wǎng)絡(luò)達(dá)到全連接,也就是網(wǎng)絡(luò)中不存在孤立點(diǎn)。該方法只適用于隨機(jī)網(wǎng)絡(luò),實(shí)際網(wǎng)絡(luò)中沒有意義;
2) 預(yù)定義閾值空間。該方法通常采用統(tǒng)計(jì)方法剔除偽連接或弱連接,間接達(dá)到選取閾值的目的。缺點(diǎn)是當(dāng)數(shù)據(jù)變化時(shí)確定閾值過程復(fù)雜且不具備普適性,另外在刪除的弱連接中可能會(huì)存在重要的信息傳遞路徑;
3) 條件限制下的閾值空間。該方法要求在選取的閾值空間下所構(gòu)建的腦網(wǎng)絡(luò)具有小世界屬性,若構(gòu)建隨機(jī)網(wǎng)絡(luò),則要求與原始網(wǎng)絡(luò)具有相同的節(jié)點(diǎn)數(shù)和節(jié)點(diǎn)度分布,且節(jié)點(diǎn)的平均度應(yīng)滿足>2lgN。
目前還有一些新的閾值化方法,如使用網(wǎng)絡(luò)的最小生成樹(MST)構(gòu)建無偏網(wǎng)絡(luò),MST 對(duì)閾值和密度值不敏感,被認(rèn)為是網(wǎng)絡(luò)二值化的良好技術(shù),但由于網(wǎng)絡(luò)的MST 非常稀疏,會(huì)造成許多重要的本地連接被忽略[3]。還有基于熱核高斯核的無窗口方法[4],該方法可減少系統(tǒng)中虛假的快速變化大腦連接的狀態(tài)空間,解決了滑動(dòng)窗口法運(yùn)用于動(dòng)態(tài)腦網(wǎng)絡(luò)分析時(shí)存在高頻噪聲問題。
雖然以上方法從不同的角度對(duì)于腦網(wǎng)絡(luò)構(gòu)建時(shí)閾值選擇問題提出了解決方法,但是網(wǎng)絡(luò)閾值的選擇依舊沒有金標(biāo)準(zhǔn),為此本文將拓?fù)鋽?shù)據(jù)分析(TDA)方法中持續(xù)同源性(PH)理論引入腦網(wǎng)絡(luò)分析中,該方法構(gòu)建腦網(wǎng)絡(luò)時(shí)無需閾值化,可以在全尺度范圍內(nèi)進(jìn)行分析,可實(shí)現(xiàn)跨多個(gè)尺度提取腦網(wǎng)絡(luò)中持續(xù)拓?fù)涮卣鳌?/p>
定義1 拓?fù)淇臻g[1,5]:設(shè)集合X上的一個(gè)拓?fù)淇臻gU是2X上的一個(gè)子集,即U?2X,如果滿足下列條件:(1)Φ,X?U;(2)u1,u2?U,u1∪u2?U;(3)u1,u2?U,u1∩u2?U;則稱(X,U) 為有限集X的拓?fù)淇臻g。
定義2 單純形[5]:設(shè)在實(shí)數(shù)域的n維向量空間Rn中,存在一組向量a0,a1,a2,…,an,使得{a1-a0,a2- a0,…,an - a0} 線性無關(guān)。設(shè)E ={θ0a0+θ1a1+…+ θnan}θ0+ θ1+…+ θn =1,θi >0} ,點(diǎn)集E就稱為一個(gè)n維單純形。
0 維單純形就是點(diǎn);1 維單純形就是線段;2 維單純形就是三角形;三維單純形就是立體三角形。
定義3 單純復(fù)形:設(shè)Κ是單純形的有限集合,若滿足如下條件[1]:(1)若σ∈K,則Κ中任意一個(gè)單純形的任意面仍屬于K;( 2) 對(duì)于σ1,σ2∈K,如果σ1∩σ2是空集,或者σ1∩σ2在σ1和σ2的公共面,那么稱Κ為單純復(fù)形。單純復(fù)形K中單純形維數(shù)的最大值稱為K的維數(shù),表示為
定義4 Rips 復(fù)形[6]:對(duì)于點(diǎn)云集合X,設(shè)d(,) 表示點(diǎn)云集合中兩點(diǎn)的距離,那么R(X,λ) 為Rips 復(fù)形當(dāng)且僅當(dāng)其k維單形[x0x1…xk] 滿足d(xi,xj) ≤λ,0 ≤i,j≤k。
根據(jù)持續(xù)同源數(shù)據(jù)分析方法[7],結(jié)合腦電信號(hào)處理的特點(diǎn),本文設(shè)計(jì)的全尺度腦網(wǎng)絡(luò)分析模型如圖1 所示。模型的輸入是腦電時(shí)間序列信號(hào),選取合適的度量空間后,腦電時(shí)序信號(hào)中的點(diǎn)就稱為點(diǎn)云,在該空間構(gòu)建點(diǎn)云的鄰接矩陣,然后通過計(jì)算持續(xù)同調(diào)性獲取持續(xù)拓?fù)涮卣?,最后通過持續(xù)特征的穩(wěn)定性分析確定網(wǎng)絡(luò)的持續(xù)不變特征。
將經(jīng)過預(yù)處理的EEG 時(shí)間序列信號(hào)輸入模型,選擇皮爾遜相關(guān)性度量空間,為EEG 信號(hào)各通道數(shù)據(jù)(即點(diǎn)云)構(gòu)造連接矩陣。針對(duì)腦電信號(hào)的特點(diǎn),實(shí)驗(yàn)構(gòu)造的是無向加權(quán)網(wǎng),電極通道為網(wǎng)絡(luò)節(jié)點(diǎn),也就是一維單純形。
圖1 基于PH 的全尺度腦網(wǎng)絡(luò)分析模型框架圖Fig.1 The framework of the whole-scale brain network analysis model based on PH
構(gòu)造大腦網(wǎng)絡(luò)嵌套復(fù)形的過程也就是計(jì)算持續(xù)同調(diào)的過程。持續(xù)同調(diào)分兩部分,即同源性和持續(xù)性。同源在群論稱作同調(diào),它是拓?fù)浼戏诸惖墓ぞ?,可以度量一個(gè)單純復(fù)形的特定結(jié)構(gòu);持續(xù)性是指給定一個(gè)ε,在ε所有可能值下計(jì)算哪些結(jié)構(gòu)是持續(xù)存在的,即獲得持續(xù)拓?fù)涮卣?。特征中能夠保持時(shí)間長(zhǎng)的就是有用的特征,而壽命較短的可能是噪聲,這個(gè)過程就稱為持續(xù)同調(diào)。構(gòu)造復(fù)形的關(guān)鍵步驟就是選擇合適的過濾算法和過濾閾值ε。
3.2.1 過濾閾值ε 的選擇
過濾閾值ε的選擇非常重要,一般方法是通過選擇不同ε構(gòu)造復(fù)形,找到有效結(jié)果所對(duì)應(yīng)的ε。如果ε過小,那么構(gòu)造的復(fù)形就可能是原始的點(diǎn)云,或者點(diǎn)云的幾條邊;如果ε過大,可能的結(jié)果就是原始點(diǎn)云構(gòu)成一個(gè)巨大的超維復(fù)形。
3.2.2 過濾算法的選擇
針對(duì)不同實(shí)際的應(yīng)用,需要構(gòu)造不同類型的單純復(fù)形,Vietoris-Rips complex 算法是基于圖的過濾,非常適合應(yīng)用于基于圖論的復(fù)雜腦網(wǎng)絡(luò)復(fù)形構(gòu)造,并且對(duì)處理高維數(shù)據(jù)有很好的性能。本實(shí)驗(yàn)采用Vietoris-Rips complex 算法過濾。
隨著過濾閾值ε的變化,Rips 復(fù)形的拓?fù)涮卣鲿?huì)發(fā)生變化。過濾過程中網(wǎng)絡(luò)的拓?fù)渥兓褂脳l形碼或者持續(xù)性圖可視化表示。在過濾過程中主要是計(jì)算p維貝蒂數(shù)間隔[εbirth,εdeath] ,其中εbirth就是單純復(fù)形中p維孔開始的時(shí)間,而εdeath是消失的時(shí)間,同時(shí)它們也是條形碼中一個(gè)條形碼的開始和結(jié)束點(diǎn),這些間隔用圖表示就是持續(xù)條形碼。與條形碼等效的是持續(xù)圖。條形碼中,橫坐標(biāo)表示持續(xù)特征出現(xiàn)的時(shí)間,即εbirth,縱坐標(biāo)是持續(xù)特征消失的時(shí)間εdeath,將過濾過程中求得的間隔集合[εbirth,εdeath] 作為持續(xù)圖中點(diǎn)的坐標(biāo),將集合中所有的間隔對(duì)表示到坐標(biāo)中就繪制了持續(xù)圖。在條形碼中,橫坐標(biāo)表示過濾閾值ε,用[εbirth,εdeath] 的長(zhǎng)度表示條形碼的長(zhǎng)度,條形碼長(zhǎng)度大的表示持續(xù)拓?fù)涮卣?,條形碼長(zhǎng)度很短或者只是一個(gè)點(diǎn)的表示噪聲,相對(duì)應(yīng)的在持續(xù)圖中,離對(duì)角線遠(yuǎn)的點(diǎn)表示持續(xù)特征,離對(duì)角線很近的點(diǎn)表示噪聲。
拓?fù)涮卣鞯姆€(wěn)定性分析也就是條形碼或持續(xù)圖的統(tǒng)計(jì)分析問題。本實(shí)驗(yàn)選用方法成熟且適合腦網(wǎng)絡(luò)分析的穩(wěn)定性度量標(biāo)準(zhǔn)。常用的穩(wěn)定性度量標(biāo)準(zhǔn)有Bottleneck 距離和Wasserstein 距離。如果對(duì)數(shù)據(jù)集的一個(gè)小擾動(dòng)只會(huì)在該指標(biāo)之前的持續(xù)圖中造成一個(gè)小變動(dòng),說明該指標(biāo)是穩(wěn)定的。
定義5:設(shè)p∈[1,∞),兩個(gè)圖X和Y之間的p階Wasserstein 距離定義為
式中:?:X→Y——從X到Y(jié)的映射。
當(dāng)p =∞時(shí),距離d是二維空間的度量,以上公式表示為
式中:W∞[d]∞——Bottleneck 距離。
Bottleneck 距離度量?jī)蓚€(gè)圖對(duì)應(yīng)匹配點(diǎn)之間的最大距離,可以捕獲持續(xù)圖大的變化。Wasserstein距離度量的是兩個(gè)圖對(duì)應(yīng)匹配點(diǎn)之間的總距離,可以提供持續(xù)圖之間相似性的總體變化,它對(duì)持續(xù)圖小的變化比較敏感。
實(shí)驗(yàn)數(shù)據(jù)采用來自北京市回龍觀醫(yī)院精神分裂癥工作記憶(WM)EEG 信號(hào),經(jīng)過重參考、分段、去除眼電、肌電等偽跡預(yù)處理,網(wǎng)絡(luò)節(jié)點(diǎn)規(guī)模為60個(gè),分5 個(gè)波段θ(4~7)Hz,α(7~14)Hz,β1(14~20)Hz,β2(20~30)Hz,γ(30 ~40)Hz。實(shí)驗(yàn)分三個(gè)階段,每個(gè)被試每個(gè)頻段每個(gè)階段下,選取20 個(gè)平穩(wěn)腦電信號(hào)拼接形成最終的EEG 時(shí)間序列。
本實(shí)驗(yàn)采用皮爾遜相關(guān)性度量空間構(gòu)造鄰接矩陣。將各節(jié)點(diǎn)間皮爾遜相關(guān)系數(shù)的倒數(shù)作為節(jié)點(diǎn)間連接的權(quán)重,得到60 ×60 無向加權(quán)網(wǎng)。精神分裂癥工作記憶編碼階段不同稀疏度下構(gòu)造的動(dòng)態(tài)鄰接矩陣如圖2 所示,可得到如下結(jié)論。
圖2 編碼階段不同連接密度下鄰接矩陣示意圖Fig.2 Adjacency matrixes at different connection densities in encoding stage
1)網(wǎng)絡(luò)連接密度較小時(shí)(20%左右),健康被試和精分病人腦網(wǎng)絡(luò)差異很大;
2)從網(wǎng)絡(luò)連接密度50%左右開始,健康被試和精分病人的連接矩陣變化逐漸變小,說明健康被試和精分病人在工作記憶編碼階段的連接矩陣是有“形狀”的,同時(shí)從圖3 無閾值條件下構(gòu)造的鄰接矩陣也可以看到同樣的結(jié)果。
實(shí)驗(yàn)采用斯坦福大學(xué)拓?fù)溆?jì)算小組基于PLEX庫[8]開發(fā)的軟件包JavaPlex[9]計(jì)算。
圖3 編碼階段不同頻段及全頻段鄰接矩陣圖Fig.3 Adjacency matrix at different bands and mean in encoding stage
構(gòu)造復(fù)形需要確定四個(gè)參數(shù)(1)由邊權(quán)矩陣構(gòu)造的點(diǎn)云坐標(biāo)文件(. txt);(2) 最大過濾值;(3)最大維數(shù);(4) 過濾步驟數(shù)(Filatration steps,簡(jiǎn)稱Fs)。以下根據(jù)實(shí)驗(yàn)情況分別討論如何確定這幾個(gè)參數(shù)可以達(dá)到最佳實(shí)驗(yàn)效果。
1)由邊權(quán)矩陣構(gòu)造點(diǎn)云坐標(biāo)文件
首先將鄰接矩陣轉(zhuǎn)換為邊權(quán)矩陣,矩陣中每一行為“ijωij”,然后用度量映射(Isometric Feature Mapping 簡(jiǎn)稱ISOMAP)算法[10],將高維鄰接矩陣中兩兩節(jié)點(diǎn)之間的距離在低維上(通常是2 維,也可以是任意維)找到一組新的樣本點(diǎn)(即點(diǎn)云),使降維后兩點(diǎn)間的距離與它們?cè)诟呔S上的距離相等。ISOMAP 算法既能夠保留非線性數(shù)據(jù)的本身的幾何結(jié)構(gòu),又能夠保持全局的結(jié)構(gòu)信息。
2)最大過濾值
構(gòu)造邊權(quán)矩陣后,分別取每個(gè)階段各個(gè)節(jié)點(diǎn)之間距離的最大值作為實(shí)驗(yàn)最大過濾值,全頻段及5個(gè)頻段下的最大過濾值見表1。
表1 最大過濾值Tab.1 Max filtration value
3)最大維數(shù)和過濾步驟數(shù)
最大維數(shù)初值設(shè)定為3,也就是提取dim0、dim1、dim2 和dim3 四個(gè)維度下的持續(xù)拓?fù)涮卣鳌?/p>
過濾步驟數(shù)即過濾步長(zhǎng),根據(jù)文獻(xiàn)[7]一般Fs=20,本實(shí)驗(yàn)取Fs=20,100,1 000 三種情況分別提取持續(xù)拓?fù)涮卣?,來確定模型中最優(yōu)Fs。實(shí)驗(yàn)結(jié)果見表2。
表2 三種過濾步驟數(shù)下實(shí)驗(yàn)結(jié)果比較Tab.2 Experiment results in three filtration steps
表2 中,運(yùn)行時(shí)間是在計(jì)算機(jī)配置為CPU:Interl(R)core(TM)i7-6700;內(nèi)存32G;操作系統(tǒng)WindowsX64 位機(jī)上運(yùn)行得到的。從表2 可知,三種情況下構(gòu)造的復(fù)形總數(shù)不變,第二種情況運(yùn)行時(shí)間比第一種情況多消耗21.14%,但是特征數(shù)變化不大,第三種情況運(yùn)行時(shí)間比第二種情況多消耗3.49%,特征數(shù)變化較大,因此數(shù)據(jù)量較大時(shí),權(quán)衡時(shí)間效率和特征數(shù),過濾步驟數(shù)可以選擇20 或者1000。以下通過持續(xù)圖可視化確定最終的過濾步驟數(shù)。Fs=20,100,1000 三種情況下健康被試持續(xù)拓?fù)涮卣鞯某掷m(xù)圖如圖4 至圖6 所示。
圖4 健康被試編碼階段持續(xù)圖(Fs=20)Fig.4 Healthy subjects' persistence diagrams(Fs=20)in encoding stage
圖5 健康被試編碼階段持續(xù)圖(Fs=100)Fig.5 Healthy subjects' persistence diagrams(Fs=100)in encoding stag
圖6 健康被試編碼階段持續(xù)圖(Fs=1000)Fig.6 Healthy subjects' persistence diagrams(Fs=1000)in encoding stage
從圖6 可知,dim0 特征基本相同,dim1 和dim2兩種情況,雖然Fs=100 和Fs=1000 時(shí)特征數(shù)變多,但這些特征多分布在對(duì)角線附近,也就是持續(xù)時(shí)間很短的噪聲,只有間隔為[3.474274,3.709818]始終存在也就是持續(xù)拓?fù)涮卣?。由此可以確定過濾步驟選Fs=20 效果最佳。同時(shí)從表2 可知,dim3 特征數(shù)始終是0,所以過濾的最大維數(shù)選擇2。
根據(jù)以上實(shí)驗(yàn)確定的過濾參數(shù)(1) 最大維數(shù)2;(2) 各頻段健康被試和精分病人最大過濾值對(duì)應(yīng)取表1 中值;(3) 過濾步驟Fs=20,計(jì)算全頻段、α、β1、β2、γ、θ 腦網(wǎng)絡(luò)的持續(xù)特征,分別用條形碼和持續(xù)圖可視化。編碼階段精分病人全頻段持續(xù)特征的條形碼如圖7 所示。
圖7 精分病人編碼階段三個(gè)維度的條形碼示意圖Fig.7 Patients persistence barcodes in encoding stage
本實(shí)驗(yàn)采用Bottleneck 距離和Wasserstein 距離兩個(gè)度量標(biāo)準(zhǔn)對(duì)持續(xù)圖進(jìn)行比較測(cè)度持續(xù)特征的穩(wěn)定性。實(shí)驗(yàn)采用Python 環(huán)境下Ripser 調(diào)用GUDHI 包計(jì)算Bottleneck 距離和Wasserstein 距離。
1) Bottleneck 距離
計(jì)算Bottleneck 距離[10]的重要參數(shù)是精確度e。本實(shí)驗(yàn)e取兩個(gè)值作比較,一個(gè)e=0.01,計(jì)算近似值,另一個(gè)e 取默認(rèn)值,計(jì)算真實(shí)值,計(jì)算結(jié)果見表3。
表3 持續(xù)圖間的Bottleneck 距離Tab.3 Bottleneck distances between diagrams
從表中可知,除dim0 的α頻段和dim1 的γ頻段外近似值和真實(shí)值之間的誤差很小,可能這兩個(gè)頻段的持續(xù)拓?fù)涮卣鞔嬖谄娈愔怠?/p>
2)Wasserstein 距離
健康被試和精分病人全頻段和各個(gè)頻段下兩個(gè)維度的Wasserstein 距離[11]計(jì)算結(jié)果見表4。
表4 持續(xù)圖間的Wasserstein 距離Tab.4 Wasserstein distances between diagrams
本文提出全尺度復(fù)雜腦網(wǎng)絡(luò)模型并應(yīng)用于精神分裂癥病人WM 數(shù)據(jù)分析中,經(jīng)實(shí)驗(yàn)分析確定模型中的重要參數(shù)和算法,在邊權(quán)矩陣和點(diǎn)云文件互換時(shí),運(yùn)用ISOMAP 算法降維后,矩陣由60 維基本上降到了34~38 維之間,為后續(xù)高效率處理數(shù)據(jù)提供了良好的基礎(chǔ)。
通過全方位的實(shí)驗(yàn)效果比較確定幾個(gè)重要參數(shù):1) 最大維數(shù)為2;2) 實(shí)驗(yàn)表明健康被試β2,θ和γ三個(gè)頻段在dim2,精分病人β2、和γ兩個(gè)頻段在dim2 都無持續(xù)拓?fù)涮卣?,所以在穩(wěn)定性分析時(shí),實(shí)驗(yàn)只分析dim0 和dim1 兩個(gè)維度;3) 最大過濾步驟為Fs=20 既不丟失重要特征,又可剔除噪聲,同時(shí)還可提高時(shí)間效率。
在穩(wěn)定性分析中,持續(xù)圖間的Bottleneck 距離結(jié)果顯示,精確度參數(shù)e采用默認(rèn)值實(shí)驗(yàn)效果更接近真實(shí)值,但當(dāng)持續(xù)圖沒有奇異值時(shí)真實(shí)值和差異值差別不大。另外,在α和θ頻段Bottleneck 距離較小,也就是健康被試和精分病人的持續(xù)圖大的變化較少,可以選取這兩個(gè)頻段的持續(xù)拓?fù)涮卣髯鳛槟P偷妮敵觥?/p>
本文提出基于持續(xù)同調(diào)性的全尺度腦網(wǎng)絡(luò)分析模型,并對(duì)模型中的每一數(shù)據(jù)處理步驟涉及到的算法和參數(shù)做了分析,研究了全尺度腦網(wǎng)絡(luò)構(gòu)造中的節(jié)點(diǎn)、邊權(quán)矩陣構(gòu)造、過濾閾值的選擇等關(guān)鍵性問題。并將模型應(yīng)用于精神分裂癥WM 持續(xù)拓?fù)涮卣鞣治鲋?,獲取持續(xù)拓?fù)涮卣?,并分析其穩(wěn)定性。實(shí)驗(yàn)結(jié)果表明本文提出的模型具有穩(wěn)定性和抗噪性,可以為精神類疾病醫(yī)學(xué)影像分析提供穩(wěn)定的生物參考指標(biāo)。