拓撲數(shù)據(jù)分析在復(fù)雜腦網(wǎng)絡(luò)分析中的應(yīng)用

陰桂梅 王千山 姚 蓉 李海芳

（1.太原師范學院計算機系，山西晉中030619；2.太原理工大學信息與計算機學院，山西太原030024）

1 引 言

宇航員面臨著我們大多數(shù)人永遠都不會經(jīng)歷的精神和心理挑戰(zhàn)，自從美國女宇航員莉薩·諾瓦克事件后，NASA 一直致力于研究降低宇航員在太空飛行期間的精神健康問題風險。航天器的工程技術(shù)實現(xiàn)能力是實現(xiàn)載人航天任務(wù)的基礎(chǔ)，但宇航員需面臨新的心理學和精神病學的挑戰(zhàn)是主要的限制性因素。為此，本文將從宇航員的視角，運用拓撲數(shù)據(jù)分析中的持續(xù)同調(diào)理論，以精神分裂癥數(shù)據(jù)為例，分析其復(fù)雜腦網(wǎng)絡(luò)的持續(xù)拓撲特征，以便能夠為提前預(yù)防和檢測宇航員精神疾病提供有效的腦網(wǎng)絡(luò)指標。

拓撲數(shù)據(jù)分析（TDA）［1］是一個數(shù)據(jù)分析、代數(shù)拓撲、計算幾何、計算機科學、統(tǒng)計等多領(lǐng)域相關(guān)的一個領(lǐng)域。TDA 的主要目標是利用幾何學和拓撲學的思想研究數(shù)據(jù)的定性特征。為了達到這個目標，TDA 需要精確的定義定性特性，還有在具體實踐應(yīng)用中的計算工具，以及保證這些特性穩(wěn)定性、健壯性的理論。解決這幾個問題的一種方法就是TDA 中的持久同源性（PH）理論。將PH理論運用于大腦網(wǎng)絡(luò)分析是當前正在興起的一個研究方向。

對大腦成像數(shù)據(jù)處理和分析時，一般都是通過生成表示節(jié)點之間連接強度的矩陣，然后選擇合適的閾值對矩陣二值化最終生成鄰接矩陣來構(gòu)建腦網(wǎng)絡(luò)。閾值的選擇在網(wǎng)絡(luò)構(gòu)建中起著重要作用，因為它影響網(wǎng)絡(luò)連接的密度和網(wǎng)絡(luò)的拓撲結(jié)構(gòu)。通常有三種網(wǎng)絡(luò)閾值化方法［2］。

1） 選取單個閾值。一般選某個連接密度作為單個閾值，當閾值選取2lgN/N時，網(wǎng)絡(luò)達到全連接，也就是網(wǎng)絡(luò)中不存在孤立點。該方法只適用于隨機網(wǎng)絡(luò)，實際網(wǎng)絡(luò)中沒有意義；

2） 預(yù)定義閾值空間。該方法通常采用統(tǒng)計方法剔除偽連接或弱連接，間接達到選取閾值的目的。缺點是當數(shù)據(jù)變化時確定閾值過程復(fù)雜且不具備普適性，另外在刪除的弱連接中可能會存在重要的信息傳遞路徑；

3） 條件限制下的閾值空間。該方法要求在選取的閾值空間下所構(gòu)建的腦網(wǎng)絡(luò)具有小世界屬性，若構(gòu)建隨機網(wǎng)絡(luò)，則要求與原始網(wǎng)絡(luò)具有相同的節(jié)點數(shù)和節(jié)點度分布，且節(jié)點的平均度應(yīng)滿足＞2lgN。

目前還有一些新的閾值化方法，如使用網(wǎng)絡(luò)的最小生成樹（MST）構(gòu)建無偏網(wǎng)絡(luò)，MST 對閾值和密度值不敏感，被認為是網(wǎng)絡(luò)二值化的良好技術(shù)，但由于網(wǎng)絡(luò)的MST 非常稀疏，會造成許多重要的本地連接被忽略［3］。還有基于熱核高斯核的無窗口方法［4］，該方法可減少系統(tǒng)中虛假的快速變化大腦連接的狀態(tài)空間，解決了滑動窗口法運用于動態(tài)腦網(wǎng)絡(luò)分析時存在高頻噪聲問題。

雖然以上方法從不同的角度對于腦網(wǎng)絡(luò)構(gòu)建時閾值選擇問題提出了解決方法，但是網(wǎng)絡(luò)閾值的選擇依舊沒有金標準，為此本文將拓撲數(shù)據(jù)分析（TDA）方法中持續(xù)同源性（PH）理論引入腦網(wǎng)絡(luò)分析中，該方法構(gòu)建腦網(wǎng)絡(luò)時無需閾值化，可以在全尺度范圍內(nèi)進行分析，可實現(xiàn)跨多個尺度提取腦網(wǎng)絡(luò)中持續(xù)拓撲特征。

2 拓撲分析持續(xù)同調(diào)性理論

定義1 拓撲空間［1，5］:設(shè)集合X上的一個拓撲空間U是2X上的一個子集，即U?2X，如果滿足下列條件:（1）Φ，X?U；（2）u1，u2?U，u1∪u2?U；（3）u1，u2?U，u1∩u2?U；則稱（X，U） 為有限集X的拓撲空間。

定義2 單純形［5］:設(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｝ ，點集E就稱為一個n維單純形。

0 維單純形就是點；1 維單純形就是線段；2 維單純形就是三角形；三維單純形就是立體三角形。

定義3 單純復(fù)形:設(shè)Κ是單純形的有限集合，若滿足如下條件［1］:（1）若σ∈K，則Κ中任意一個單純形的任意面仍屬于K；（ 2） 對于σ1，σ2∈K，如果σ1∩σ2是空集，或者σ1∩σ2在σ1和σ2的公共面，那么稱Κ為單純復(fù)形。單純復(fù)形K中單純形維數(shù)的最大值稱為K的維數(shù)，表示為

定義4 Rips 復(fù)形［6］:對于點云集合X，設(shè)d（，） 表示點云集合中兩點的距離，那么R（X，λ） 為Rips 復(fù)形當且僅當其k維單形［x0x1…xk］ 滿足d（xi，xj） ≤λ，0 ≤i，j≤k。

3 基于PH 的全尺度腦網(wǎng)絡(luò)分析模型

根據(jù)持續(xù)同源數(shù)據(jù)分析方法［7］，結(jié)合腦電信號處理的特點，本文設(shè)計的全尺度腦網(wǎng)絡(luò)分析模型如圖1 所示。模型的輸入是腦電時間序列信號，選取合適的度量空間后，腦電時序信號中的點就稱為點云，在該空間構(gòu)建點云的鄰接矩陣，然后通過計算持續(xù)同調(diào)性獲取持續(xù)拓撲特征，最后通過持續(xù)特征的穩(wěn)定性分析確定網(wǎng)絡(luò)的持續(xù)不變特征。

3.1 構(gòu)造鄰接矩陣

將經(jīng)過預(yù)處理的EEG 時間序列信號輸入模型，選擇皮爾遜相關(guān)性度量空間，為EEG 信號各通道數(shù)據(jù)（即點云）構(gòu)造連接矩陣。針對腦電信號的特點，實驗構(gòu)造的是無向加權(quán)網(wǎng)，電極通道為網(wǎng)絡(luò)節(jié)點，也就是一維單純形。

圖1 基于PH 的全尺度腦網(wǎng)絡(luò)分析模型框架圖Fig.1 The framework of the whole-scale brain network analysis model based on PH

3.2 構(gòu)造大腦網(wǎng)絡(luò)復(fù)形

構(gòu)造大腦網(wǎng)絡(luò)嵌套復(fù)形的過程也就是計算持續(xù)同調(diào)的過程。持續(xù)同調(diào)分兩部分，即同源性和持續(xù)性。同源在群論稱作同調(diào)，它是拓撲集合分類的工具，可以度量一個單純復(fù)形的特定結(jié)構(gòu)；持續(xù)性是指給定一個ε，在ε所有可能值下計算哪些結(jié)構(gòu)是持續(xù)存在的，即獲得持續(xù)拓撲特征。特征中能夠保持時間長的就是有用的特征，而壽命較短的可能是噪聲，這個過程就稱為持續(xù)同調(diào)。構(gòu)造復(fù)形的關(guān)鍵步驟就是選擇合適的過濾算法和過濾閾值ε。

3.2.1 過濾閾值ε 的選擇

過濾閾值ε的選擇非常重要，一般方法是通過選擇不同ε構(gòu)造復(fù)形，找到有效結(jié)果所對應(yīng)的ε。如果ε過小，那么構(gòu)造的復(fù)形就可能是原始的點云，或者點云的幾條邊；如果ε過大，可能的結(jié)果就是原始點云構(gòu)成一個巨大的超維復(fù)形。

3.2.2 過濾算法的選擇

針對不同實際的應(yīng)用，需要構(gòu)造不同類型的單純復(fù)形，Vietoris-Rips complex 算法是基于圖的過濾，非常適合應(yīng)用于基于圖論的復(fù)雜腦網(wǎng)絡(luò)復(fù)形構(gòu)造，并且對處理高維數(shù)據(jù)有很好的性能。本實驗采用Vietoris-Rips complex 算法過濾。

3.3 持續(xù)拓撲特征的可視化

隨著過濾閾值ε的變化，Rips 復(fù)形的拓撲特征會發(fā)生變化。過濾過程中網(wǎng)絡(luò)的拓撲變化使用條形碼或者持續(xù)性圖可視化表示。在過濾過程中主要是計算p維貝蒂數(shù)間隔［εbirth，εdeath］ ，其中εbirth就是單純復(fù)形中p維孔開始的時間，而εdeath是消失的時間，同時它們也是條形碼中一個條形碼的開始和結(jié)束點，這些間隔用圖表示就是持續(xù)條形碼。與條形碼等效的是持續(xù)圖。條形碼中，橫坐標表示持續(xù)特征出現(xiàn)的時間，即εbirth，縱坐標是持續(xù)特征消失的時間εdeath，將過濾過程中求得的間隔集合［εbirth，εdeath］ 作為持續(xù)圖中點的坐標，將集合中所有的間隔對表示到坐標中就繪制了持續(xù)圖。在條形碼中，橫坐標表示過濾閾值ε，用［εbirth，εdeath］ 的長度表示條形碼的長度，條形碼長度大的表示持續(xù)拓撲特征，條形碼長度很短或者只是一個點的表示噪聲，相對應(yīng)的在持續(xù)圖中，離對角線遠的點表示持續(xù)特征，離對角線很近的點表示噪聲。

3.4 持續(xù)拓撲特征的穩(wěn)定性分析

拓撲特征的穩(wěn)定性分析也就是條形碼或持續(xù)圖的統(tǒng)計分析問題。本實驗選用方法成熟且適合腦網(wǎng)絡(luò)分析的穩(wěn)定性度量標準。常用的穩(wěn)定性度量標準有Bottleneck 距離和Wasserstein 距離。如果對數(shù)據(jù)集的一個小擾動只會在該指標之前的持續(xù)圖中造成一個小變動，說明該指標是穩(wěn)定的。

定義5:設(shè)p∈［1，∞），兩個圖X和Y之間的p階Wasserstein 距離定義為

式中:?:X→Y——從X到Y(jié)的映射。

當p =∞時，距離d是二維空間的度量，以上公式表示為

式中:W∞［d］∞——Bottleneck 距離。

Bottleneck 距離度量兩個圖對應(yīng)匹配點之間的最大距離，可以捕獲持續(xù)圖大的變化。Wasserstein距離度量的是兩個圖對應(yīng)匹配點之間的總距離，可以提供持續(xù)圖之間相似性的總體變化，它對持續(xù)圖小的變化比較敏感。

4 實驗及分析

4.1 實驗數(shù)據(jù)及預(yù)處理

實驗數(shù)據(jù)采用來自北京市回龍觀醫(yī)院精神分裂癥工作記憶（WM）EEG 信號，經(jīng)過重參考、分段、去除眼電、肌電等偽跡預(yù)處理，網(wǎng)絡(luò)節(jié)點規(guī)模為60個，分5 個波段θ（4～7）Hz，α（7～14）Hz，β1（14～20）Hz，β2（20～30）Hz，γ（30 ～40）Hz。實驗分三個階段，每個被試每個頻段每個階段下，選取20 個平穩(wěn)腦電信號拼接形成最終的EEG 時間序列。

4.2 基于PH 的鄰接矩陣構(gòu)建

本實驗采用皮爾遜相關(guān)性度量空間構(gòu)造鄰接矩陣。將各節(jié)點間皮爾遜相關(guān)系數(shù)的倒數(shù)作為節(jié)點間連接的權(quán)重，得到60 ×60 無向加權(quán)網(wǎng)。精神分裂癥工作記憶編碼階段不同稀疏度下構(gòu)造的動態(tài)鄰接矩陣如圖2 所示，可得到如下結(jié)論。

圖2 編碼階段不同連接密度下鄰接矩陣示意圖Fig.2 Adjacency matrixes at different connection densities in encoding stage

1）網(wǎng)絡(luò)連接密度較小時（20%左右），健康被試和精分病人腦網(wǎng)絡(luò)差異很大；

2）從網(wǎng)絡(luò)連接密度50%左右開始，健康被試和精分病人的連接矩陣變化逐漸變小，說明健康被試和精分病人在工作記憶編碼階段的連接矩陣是有“形狀”的，同時從圖3 無閾值條件下構(gòu)造的鄰接矩陣也可以看到同樣的結(jié)果。

4.3 構(gòu)造精神分裂癥WM 數(shù)據(jù)的單純復(fù)形

實驗采用斯坦福大學拓撲計算小組基于PLEX庫［8］開發(fā)的軟件包JavaPlex［9］計算。

圖3 編碼階段不同頻段及全頻段鄰接矩陣圖Fig.3 Adjacency matrix at different bands and mean in encoding stage

構(gòu)造復(fù)形需要確定四個參數(shù)（1）由邊權(quán)矩陣構(gòu)造的點云坐標文件（. txt）；（2） 最大過濾值；（3）最大維數(shù)；（4） 過濾步驟數(shù)（Filatration steps，簡稱Fs）。以下根據(jù)實驗情況分別討論如何確定這幾個參數(shù)可以達到最佳實驗效果。

1）由邊權(quán)矩陣構(gòu)造點云坐標文件

首先將鄰接矩陣轉(zhuǎn)換為邊權(quán)矩陣，矩陣中每一行為“ijωij”，然后用度量映射（Isometric Feature Mapping 簡稱ISOMAP）算法［10］，將高維鄰接矩陣中兩兩節(jié)點之間的距離在低維上（通常是2 維，也可以是任意維）找到一組新的樣本點（即點云），使降維后兩點間的距離與它們在高維上的距離相等。ISOMAP 算法既能夠保留非線性數(shù)據(jù)的本身的幾何結(jié)構(gòu)，又能夠保持全局的結(jié)構(gòu)信息。

2）最大過濾值

構(gòu)造邊權(quán)矩陣后，分別取每個階段各個節(jié)點之間距離的最大值作為實驗最大過濾值，全頻段及5個頻段下的最大過濾值見表1。

表1 最大過濾值Tab.1 Max filtration value

3）最大維數(shù)和過濾步驟數(shù)

最大維數(shù)初值設(shè)定為3，也就是提取dim0、dim1、dim2 和dim3 四個維度下的持續(xù)拓撲特征。

過濾步驟數(shù)即過濾步長，根據(jù)文獻［7］一般Fs=20，本實驗取Fs=20，100，1 000 三種情況分別提取持續(xù)拓撲特征，來確定模型中最優(yōu)Fs。實驗結(jié)果見表2。

表2 三種過濾步驟數(shù)下實驗結(jié)果比較Tab.2 Experiment results in three filtration steps

表2 中，運行時間是在計算機配置為CPU:Interl（R）core（TM）i7-6700；內(nèi)存32G；操作系統(tǒng)WindowsX64 位機上運行得到的。從表2 可知，三種情況下構(gòu)造的復(fù)形總數(shù)不變，第二種情況運行時間比第一種情況多消耗21.14%，但是特征數(shù)變化不大，第三種情況運行時間比第二種情況多消耗3.49%，特征數(shù)變化較大，因此數(shù)據(jù)量較大時，權(quán)衡時間效率和特征數(shù)，過濾步驟數(shù)可以選擇20 或者1000。以下通過持續(xù)圖可視化確定最終的過濾步驟數(shù)。Fs=20，100，1000 三種情況下健康被試持續(xù)拓撲特征的持續(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ù)變多，但這些特征多分布在對角線附近，也就是持續(xù)時間很短的噪聲，只有間隔為［3.474274，3.709818］始終存在也就是持續(xù)拓撲特征。由此可以確定過濾步驟選Fs=20 效果最佳。同時從表2 可知，dim3 特征數(shù)始終是0，所以過濾的最大維數(shù)選擇2。

4.4 精神分裂癥WM 持續(xù)拓撲特征

根據(jù)以上實驗確定的過濾參數(shù)（1） 最大維數(shù)2；（2） 各頻段健康被試和精分病人最大過濾值對應(yīng)取表1 中值；（3） 過濾步驟Fs=20，計算全頻段、α、β1、β2、γ、θ 腦網(wǎng)絡(luò)的持續(xù)特征，分別用條形碼和持續(xù)圖可視化。編碼階段精分病人全頻段持續(xù)特征的條形碼如圖7 所示。

圖7 精分病人編碼階段三個維度的條形碼示意圖Fig.7 Patients persistence barcodes in encoding stage

4.5 精神分裂癥WM 持續(xù)拓撲特征穩(wěn)定性分析

本實驗采用Bottleneck 距離和Wasserstein 距離兩個度量標準對持續(xù)圖進行比較測度持續(xù)特征的穩(wěn)定性。實驗采用Python 環(huán)境下Ripser 調(diào)用GUDHI 包計算Bottleneck 距離和Wasserstein 距離。

1） Bottleneck 距離

計算Bottleneck 距離［10］的重要參數(shù)是精確度e。本實驗e取兩個值作比較，一個e=0.01，計算近似值，另一個e 取默認值，計算真實值，計算結(jié)果見表3。

表3 持續(xù)圖間的Bottleneck 距離Tab.3 Bottleneck distances between diagrams

從表中可知，除dim0 的α頻段和dim1 的γ頻段外近似值和真實值之間的誤差很小，可能這兩個頻段的持續(xù)拓撲特征存在奇異值。

2）Wasserstein 距離

健康被試和精分病人全頻段和各個頻段下兩個維度的Wasserstein 距離［11］計算結(jié)果見表4。

表4 持續(xù)圖間的Wasserstein 距離Tab.4 Wasserstein distances between diagrams

4.6 結(jié)果分析

本文提出全尺度復(fù)雜腦網(wǎng)絡(luò)模型并應(yīng)用于精神分裂癥病人WM 數(shù)據(jù)分析中，經(jīng)實驗分析確定模型中的重要參數(shù)和算法，在邊權(quán)矩陣和點云文件互換時，運用ISOMAP 算法降維后，矩陣由60 維基本上降到了34～38 維之間，為后續(xù)高效率處理數(shù)據(jù)提供了良好的基礎(chǔ)。

通過全方位的實驗效果比較確定幾個重要參數(shù):1） 最大維數(shù)為2；2） 實驗表明健康被試β2，θ和γ三個頻段在dim2，精分病人β2、和γ兩個頻段在dim2 都無持續(xù)拓撲特征，所以在穩(wěn)定性分析時，實驗只分析dim0 和dim1 兩個維度；3） 最大過濾步驟為Fs=20 既不丟失重要特征，又可剔除噪聲，同時還可提高時間效率。

在穩(wěn)定性分析中，持續(xù)圖間的Bottleneck 距離結(jié)果顯示，精確度參數(shù)e采用默認值實驗效果更接近真實值，但當持續(xù)圖沒有奇異值時真實值和差異值差別不大。另外，在α和θ頻段Bottleneck 距離較小，也就是健康被試和精分病人的持續(xù)圖大的變化較少，可以選取這兩個頻段的持續(xù)拓撲特征作為模型的輸出。

5 結(jié)束語

本文提出基于持續(xù)同調(diào)性的全尺度腦網(wǎng)絡(luò)分析模型，并對模型中的每一數(shù)據(jù)處理步驟涉及到的算法和參數(shù)做了分析，研究了全尺度腦網(wǎng)絡(luò)構(gòu)造中的節(jié)點、邊權(quán)矩陣構(gòu)造、過濾閾值的選擇等關(guān)鍵性問題。并將模型應(yīng)用于精神分裂癥WM 持續(xù)拓撲特征分析中，獲取持續(xù)拓撲特征，并分析其穩(wěn)定性。實驗結(jié)果表明本文提出的模型具有穩(wěn)定性和抗噪性，可以為精神類疾病醫(yī)學影像分析提供穩(wěn)定的生物參考指標。