當前復雜網絡的研究已經成為腦電信號研究的熱點。因為腦電網絡產生的時間序列保存了網絡結點信息，因此研究網絡產生的時間序列同樣能夠達到研究癲癇腦電信號的目的。基于此，本文提出了一種基于改進的k-最近鄰網絡產生時間序列來分析癲癇腦電信號的方法。研究結果表明，研究網絡產生的時間序列的功率譜比直接研究原始癲癇腦電信號的功率譜更容易區分正常人和癲癇患者。此外，研究改進的k-最近鄰網絡的聚類系數也能區分正常人和癲癇患者。通過本文研究結果，期望能夠為癲癇研究及其今后的臨床診斷提供相關參考依據。
引用本文： 余嫻, 劉程程, 戴加飛, 李錦, 王俊, 侯鳳貞. 基于改進的k-最近鄰網絡的癲癇腦電信號分析. 生物醫學工程學雜志, 2016, 33(6): 1039-1045. doi: 10.7507/1001-5515.20160167 復制
引言
癲癇是一種常見的由于神經功能異常而引起的腦部疾病，全世界的癲癇患者數量約有5 000萬[1-2]。癲癇發作是由大腦的異常放電所引起的，這種放電活動會對癲癇患者的身體和心理造成巨大傷害。為了能讓癲癇患者盡早確診病情，及時進行針對性的治療，對癲癇的正確診斷就顯得尤為重要。早期的癲癇診斷主要通過醫生觀察腦電圖(electroencephalogram，EEG)波形是否存在異常的方式進行判斷，這種診斷過程往往伴有醫生的主觀意識，診斷結果準確率也受到影響[3]。利用線性方法對腦電信號進行癲癇預測是較為成熟的方法，但是靈敏度不高，很難分析腦電信號的動力學結構，因此無法反映大腦活動的本質特征。隨著非線性信號處理理論的發展，其在癲癇疾病的研究已得到了廣泛應用，目前利用關聯維、多尺度熵、近似熵和Lyapunov指數等方法可以有效地分析癲癇患者的腦電信號[3-6]。
近年來，復雜網絡理論得到高速發展，利用復雜網絡理論來研究癲癇腦電信號已經成為人體生理信號研究的熱點[7-8]。在基于復雜網絡理論的癲癇腦電信號研究過程中，通常將癲癇腦電信號構建成復雜網絡，然后利用復雜網絡的靜態特性，如網絡的度、聚類系數等來研究癲癇腦電信號[9]。
將癲癇腦電信號轉換成復雜網絡已經成為研究癲癇腦電信號的重要方法，但卻忽略了逆向研究，即將構建的復雜網絡再轉換成時間序列分析。近年來，已有將復雜網絡轉換成時間序列的方法開展相應研究[10]。Yutaka等[11]提出了特征向量法將網絡轉換成時間序列；Weng等[12]提出了隨機游走方法將網絡轉換成時間序列。本文對癲癇腦電信號研究采取逆向研究思路，利用改進的k-最近鄰網絡算法將癲癇腦電信號構建成復雜網絡，然后利用特征向量法將復雜網絡轉換成時間序列，最后對轉換后的時間序列進行功率譜分析，還對前面改進的k-最近鄰網絡算法構建成的復雜網絡進行了網絡聚類系數分析，以期為癲癇腦電信號的動力學研究及臨床診斷提供重要的參考依據。
1 復雜網絡構建及網絡轉換成時間序列
1.1 相空間重構
混沌時間序列分析方法在很多科研和工程領域中已得到廣泛應用[13-14]。相空間重構是混沌時間序列分析的基礎。Takens等[15]提出了用延遲坐標法對混沌時間序列x=｛xi|i=1，2，…，N｝進行重構:
| $X={{X}_{i}}{{X}_{~}}i=[{{x}_{i}},{{x}_{i+t}},\ldots ,{{x}_{i+\left( m-1 \right)t}}],i=1,2,\ldots ,M$ |
其中，m為嵌入維，t為延遲，M=N-(m-1)t為相空間中的點數。Takens定理證明了如果嵌入維m≥2d+1，d為系統動力學階數，則重構的動力系統與原動力學系統在拓撲意義上等價。
本文將通過相關積分-相關積分(correlation integral-correlation integral,C-C)算法求延遲時間t和最佳嵌入維m[16]。在C-C算法的計算中將會用到相關積分：
| ${{C}^{m}}_{s}\left( m,N,r,t \right)=\frac{2}{M\left( M-1 \right)}{{\Sigma }_{1\le i\le j＜＜M}}\theta \{r-\|{{X}_{i}}-{{X}_{j}}\|\}$ |
| $\theta \left( a \right)=\left\{ \begin{array}{*{35}{l}} 1 & a\ge 0 \\ 0 & a0 \\ \end{array} \right.$ |
其中N為數據集合的長度，t是時間尺度，M=N-(m-1)t是m維相空間重構嵌入后的點數。C-C算法的研究與函數S(m,N,r,t)=C(m,N,r,t)－Cm(1,N,r,t)的值有關。將時間序列分為t個子序列，對S(m,N,r,t)的計算采用分塊平均的策略。S(m,N,r,t)可以表示為
| $\begin{align} & S\left( m,N,r,t \right)=\frac{1}{t}\sum\limits_{s=1}^{t}{\{{{C}_{s}}\left( m,N/t,r,t \right)-} \\ & {{C}^{m}}_{s}\left( 1,N/t,r,t \right)\} \\ \end{align}$ |
最后當N→∞,S(m,N,r,t)式將變為
| $S\left( m,r,t \right)=\frac{1}{t}\sum\limits_{s=1}^{t}{\{{{C}_{s}}\left( m,r,t \right)-}{{C}^{m}}_{s}\left( 1,r,t \right)\}$ |
現在定義差量ΔS(m,t)，其用表達式描述為
| $\begin{align} & \Delta S\left( m,t \right)=\max \left\{ S\left( m,{{r}_{j}},t,N \right) \right\}- \\ & \min \{S(m,{{r}_{j}},t,N)\} \\ \end{align}$ |
根據C-C算法統計理論，一般取m=2,3,4,5，${{r}_{j}}=\frac{i\sigma }{2}$ ，其中i=1,2,3,4[16]。根據m和rj取值分別求S(t)、ΔS(t)和Scor(t)。其表達式分別如下：
| $S\left( t \right)=\frac{1}{16}\sum\limits_{m-2}^{5}{\sum\limits_{j=1}^{4}{S(m,{{r}_{j}},t)}}$ |
| $\Delta S\left( t \right)=\frac{1}{4}\sum\limits_{m-2}^{5}{\Delta S\left( m,t \right)}$ |
| ${{S}_{cor}}\left( t \right)=\Delta S\left( t \right)+S\left( t \right)$ |
通過計算S(t)或者ΔS(t)的第一個極小值來確定重構相空間的延遲時間t，通過計算Scor的最小值確定延遲窗口τw，最后可以根據嵌入維和延遲時間的關系τw=(m-1)t來確定相空間重構的最佳嵌入維數m。
1.2 k-最近鄰網絡
k-最近鄰網絡構建思想是對于網絡中的任意一個結點尋找周圍最近的k個結點與之相連接[17]。每個點都以這個規律和其他點相連接。具體步驟如下：
(1) 首先要計算相空間中的歐式距離矩陣Dij,在距離矩陣的第i行包含其他點與第i個點之間的距離。
(2) 對于每一行i，根據距離矩陣選出距離i的k個最鄰近點，可以表示成Mi=｛x(j1)…x(jk)｝，其中Mi是對于第i行的k個最鄰近點的集合，x(j1)為第一個鄰近點，x(jk)為第k個鄰近點。
(3) 構建鄰接矩陣R，矩陣R中的第i行第j列表示為Rij。對于鄰接矩陣的第i行，如果第j個點屬于集合Mi，則Rij=1，否則Rij=0。按照這種方法完成對每一行的鄰接矩陣構建。
通過以上3個步驟，可以得到一個有向網絡的鄰接矩陣。該網絡即為k-最近鄰網絡。該網絡的度分布為固定值，即P(kout)≡δ(k)，邊的總數為N×k個，其中N為相空間點的總數，k為設置的k個鄰近點數目。
1.3 改進的k-最近鄰網絡
k-最近鄰網絡是一個有向網絡，并且擁有固定度分布的良好特性。而利用求特征值法將網絡轉換成時間序列是基于無向網絡。為了利用k-最近鄰網絡內在的良好特性，必須把有向網絡變換成無向網絡。以下是構建改進的k-最近鄰網絡的過程。
(1) 首先要計算相空間中的歐式距離矩陣Dij，在距離矩陣的第i行包含其他點與第i個點之間的距離。
(2) 對于每一行i，根據距離矩陣對其他點進行從小到大排序，排除與自身距離的比較，假設相空間有N個點，可以表示成為Mi=｛x(j1)…x(jN-1)｝，其中Mi是對于第i行的N-1個最鄰近點的集合，x(j1)為第一個鄰近點，x(jN-1)為最后一個鄰近點。
(3) 構建鄰接矩陣R,矩陣R中的第i行第j列表示為Rij。對于鄰接矩陣的每一行i，首先判斷已經有多少個點與第i個點連接，即求Rij1的個數，其中j=1,…,i-1,i+1,…,N，假設為s個點與i連接。若k-s＞0，則從Mi中選取前k-s個點與i相連接，并置Rij=1，并且Rij=1，在Mi中選取的點必須符合Rij=0，其中j是屬于Mi中的點。若k-s≤0，則跳過該行，對R矩陣的下一行進行構建。
1.4 特征值法
特征值法將網絡轉換成時間序列，主要是利用網絡鄰接矩陣的特征值[11]。此方法的主要思路是求鄰接矩陣的最大特征值對應的特征向量，該特征向量即為網絡轉換時間后的時間序列。本文將利用求鄰接矩陣特征值法將網絡轉換成時間序列。
假設A(={aij})是一個N×N網絡的鄰接矩陣。定義頂點vi和vj之間距離為dij。如果aij=0(i≠j)，則dij=w(＞1)，否則dij=aij。平方距離矩陣D定義為D=｛dij2｝,接下來將D代入表達式
| $G=-\frac{1}{2}{{J}_{N}}D{{J}^{T}}_{N}$ |
其中${{J}_{N}}=E-\frac{1}{N}{{1}_{N}}{{1}^{T}}_{N}$，E是N×N的單位矩陣，1N是一個包含N個1的列向量。將G分解為特征值和特征向量表示為：
| $G=P\Lambda {{P}^{T}}=\left( P{{\Lambda }^{1/2}} \right)\times {{\left( P{{\Lambda }^{1/2}} \right)}^{T}}\equiv X{{X}^{T}}$ |
其中$\Lambda =\text{diag}({{\lambda }_{1}},{{\lambda }_{2}},\ldots ,{{\lambda }_{h}}),{{\Lambda }^{(1/2)}}=\text{diag}(\sqrt{{{\lambda }_{1}}},\sqrt{{{\lambda }_{2}}},\ldots ,\sqrt{{{\lambda }_{h}}}),P=({{p}_{1}}{{p}_{2}}\ldots {{p}_{h}}){{p}_{m}}={{\left( {{p}_{m1}}{{p}_{m2}}\ldots {{p}_{mh}} \right)}^{T}}$，h是矩陣G非零特征值的個數。矩陣X=(x1，x2，…，xN)T，其中xm=(xm1,xm2,…,xmN)T是頂點Vm在h維空間中相應的值。定義時間序列為
| ${{s}_{m}}\left( t \right)=\sqrt{{{\lambda }_{m}}}{{p}_{m1}}\left( 1\le m\le h,1\le t\le N \right)$ |
最大特征值對應主特征向量，因此在本文研究中，都是選取最大特征值對應的時間序列。
2 實驗結果及其分析
2.1 實驗數據預處理及統計方法
本文實驗數據來源于南京軍區總醫院從臨床診斷中采集的數據，包含癲癇患者腦電信號組和正常健康者腦電信號。實驗分成兩組進行，分別為正常組和癲癇患者組。正常組和癲癇患者組的人數為各15人，年齡為20~60歲。腦電信號記錄采用標準10~20系統，使用放置有16個電極(FP1,FP2,F3,F4,C3,C4,P3,P4,O1,O2,F7,F8,T3,T4,T5,T6)的電極帽采集，腦電信號記錄數據采樣頻率為512 Hz，采集數據時受試者都處于清醒狀態。實驗選取了C3導聯腦電信號作為分析對象，并對該導聯數據利用帶陷濾波器濾除50 Hz工頻干擾。
本文采用的統計方法是t檢驗，t檢驗常用于檢測樣本含量少的兩樣本差異是否具有統計學意義[7-8]。當P值小于0.05時，表示兩樣本差異具有統計學意義，能有效進行區分；否則表示兩樣本差異無統計學意義，不能進行區分。本文采用SPSS軟件對實驗數據進行t檢驗。
2.2 相空間重構及網絡構建
復雜網絡構建前必須進行相空間的重構，進行相空間重構前要選擇合適的延遲時間和嵌入維。本文利用C-C算法來獲取數據相空間重構的延遲時間和嵌入維。在C-C算法中，S(t)或者ΔS(t)首次出現極小值處可以確定相空間延遲時間，Scor(t)首次出現極小值處用來確定延遲窗口。S(t)、ΔS(t)和Scor(t)隨延遲時間t變化的結果如圖 1所示。縱軸表示的含義是C-C算法中用到的計算公式在不同延遲時間下的均值，橫軸表示延遲時間，其變化范圍是0~50 s。
圖1
C-C算法實驗結果
Figure1.
Experimental result of the C-C algorithm
如圖 1所示，S(t)、ΔS(t)和Scor(t)極小值首次出現都是在第10 s處，因此最佳延遲時間t=10 s，最佳延遲窗口τw為10 s。根據延遲時間、嵌入維和延遲窗口三者關系τw=(m-1)t可以得到最佳嵌入維m為2，本文將以延遲時間t=10 s，嵌入維m=2重構相空間。
本文運用改進的k-最近鄰網絡的方法構建規模為500個結點的網絡。為驗證改進k-最近鄰網絡的可行性，將k=20和k=300時的改進k-最近鄰網絡度分布圖分別作圖，如圖 2所示，橫軸表示網絡結點度大小，縱軸表示度分布的概率大小。
圖2
不同網絡下的度分布
Figure2.
Degree distribution under different networks
改進的k-最近鄰網絡將有向網絡變換成無向網絡，并且該網絡度為k的概率接近為1，保留了k-最近鄰網絡的特征。說明改進的k-最近鄰網絡的方法是可行的。
2.3 時間序列的功率譜分析
首先，運用改進的k-最近鄰網絡方法構建網絡，構建過程中取k值為20，然后以網絡轉換成時間序列中的w為自變量，w的變量范圍取值是1~2，間隔取0.2，計算相應w下產生的時間序列的功率譜的最大值，再對每個w下對應的功率譜最大值全體求平均。
對正常組和癲癇患者組樣本分別構建480×480大小的網絡，然后將網絡根據不同w值得到不同時間序列，對不同時間序列求功率譜最大值，再將不同w下對應的功率譜最大值求平均。實驗結果如表 1、表 2所示。
表1
(正常組)在不同時間序列下功率譜最大值的平均值(以w因子為變量)
Table1.
(for the normal group) Mean values of maximum power spectrum during different time series versus w factor
表2
(癲癇患者組)不同時間序列下功率譜最大值的平均值(以w因子為變量)
Table2.
(for the group of patients with epilepsy) Mean values of maximum power spectrum during different time series versus w factor
對表 1和表 2中數據進行t檢驗，P=0.015。該實驗結果中正常人組腦電的總體均值大于癲癇患者組腦電的總體均值，這與直接分析腦電信號功率譜最大值的結果相反，所以該研究只能進行正常人和癲癇患者的區分，網絡轉換成時間序列后退化了其本質意義。
如果不對腦電數據進行網絡變換和時間序列變化，即直接對原始信號進行功率譜分析，實驗結果如表 3、表 4所示。
表3
(正常組)直接對原始腦電信號進行求功率譜最大值
Table3.
(for the normal group )Maximum value of the power spectrum of the original electroencephalogram signal
表4
(癲癇患者組) 直接對原始腦電信號進行求功率譜最大值
Table4.
(for the group of patients with epilepsy) Maximum value of the power spectrum of original electroencephalogram signal
對表 3和表 4的實驗結果進行t檢驗,P=0.081，兩組間差異沒有統計學意義，不能對兩組受試者進行區分。而前述網絡轉換時間序列的t檢驗結果有統計學意義，說明對經過變換后的信號求功率譜比直接對原始腦電信號求功率譜更容易區分正常人和癲癇患者。
以上是以w為變量求得不同時間序列下功率譜的最大值。以下將w設置一個固定值，不妨令w=1.6，然后以k為變量，根據k值變化得到不同的網絡，再求網絡的功率譜的最大值，其中k值的取值范圍是20~300，取值間隔為20。實驗結果如表 5、表 6所示。
表5
(正常組)不同時間序列下功率譜最大值的平均值(以k為變量)
Table5.
(for the normal group) Mean values of maximum power spectrum during different time series versus k
表6
(癲癇患者組)不同時間序列下功率譜最大值的平均值(以k為變量)
Table6.
(for the group of patients with epilepsy) Mean values of maximum power spectrum during different time series versus k
根據表 5、表 6發現，隨著k值不斷增大得到的網絡轉換的時間序列的最大功率譜的總體均值中，正常人的功率譜最大值的樣本總體均值小于癲癇患者。對表 5和表 6的數據進行t檢驗，P=0.004，差異具有統計學意義，因此也可用于區分正常人和癲癇患者。
表7
(正常組)構建改進的k-最近鄰網絡的聚類系數
Table7.
(for the normal group) Clustering coefficient of network based on improved k-nearest network
表8
(癲癇患者組) 構建改進的k-最近鄰網絡的聚類系數
Table8.
(for the group of patients with epilepsy) Clustering coefficient of network based on improved k-nearest network
由上述結果可知，無論以w因子還是k為變量，對經過網絡轉換生成的時間序列求功率譜，均可區分正常人和癲癇患者，而直接對原始腦電時間序列求功率譜則兩組間差異無統計學意義，足見本文方法具有明顯的優越性。
2.4 改進的k-最近鄰網絡的聚類系數分析
實驗以k為變量，根據k值的變化求得相應網絡的聚類系數，其中k值的變化范圍是20~300，變化的間隔是40。對實驗結果進行30次平均，然后取其均值作為實驗的結果。實驗結果如表 7、表 8所示。
數據表明正常組的網絡聚類系數總體均值小于癲癇患者組。對表 7和表 8的數據進行t檢驗，P=0.007，表明對正常組和癲癇患者組構建的網絡其聚類系數具有明顯差異，這種差異性也可能用于區分正常組和癲癇患者組。
3 結論
本文主要是基于改進的k-最近鄰網絡轉換成時間序列的癲癇研究。改進的k-最近鄰網絡是一個無向網絡，并且該網絡度為k的概率接近1，具有k-最近鄰網絡的特性。將網絡轉換成時間序列是運用求鄰接矩陣的特征值特征向量方法，實驗中是選取最大特征值對應的時間序列進行研究。
通過實驗結果表明，基于改進的k-最近鄰網絡轉換而來的時間序列功率譜分析比直接對原始腦電信號數據作功率分析更容易區分正常組和癲癇患者組。實驗結果還發現， 正常組和癲癇患者組間改進的k-最近鄰網絡的聚類系數也具有明顯差異。本研究將網絡轉換成時間序列的思想應用于癲癇研究，網絡構建采用改進的k-最近鄰網絡算法，經過實驗，結果證明本文方法所得特征值在正常組和癲癇患者組間有明顯差異，這不僅能夠為癲癇領域研究提供新的研究思路，而且能夠為癲癇臨床診斷提供相關參考依據。但由于本文僅僅是對網絡轉換成時間序列的功率譜進行研究，將來還可從其他角度來分析網絡轉換后的時間序列。
引言
癲癇是一種常見的由于神經功能異常而引起的腦部疾病，全世界的癲癇患者數量約有5 000萬[1-2]。癲癇發作是由大腦的異常放電所引起的，這種放電活動會對癲癇患者的身體和心理造成巨大傷害。為了能讓癲癇患者盡早確診病情，及時進行針對性的治療，對癲癇的正確診斷就顯得尤為重要。早期的癲癇診斷主要通過醫生觀察腦電圖(electroencephalogram，EEG)波形是否存在異常的方式進行判斷，這種診斷過程往往伴有醫生的主觀意識，診斷結果準確率也受到影響[3]。利用線性方法對腦電信號進行癲癇預測是較為成熟的方法，但是靈敏度不高，很難分析腦電信號的動力學結構，因此無法反映大腦活動的本質特征。隨著非線性信號處理理論的發展，其在癲癇疾病的研究已得到了廣泛應用，目前利用關聯維、多尺度熵、近似熵和Lyapunov指數等方法可以有效地分析癲癇患者的腦電信號[3-6]。
近年來，復雜網絡理論得到高速發展，利用復雜網絡理論來研究癲癇腦電信號已經成為人體生理信號研究的熱點[7-8]。在基于復雜網絡理論的癲癇腦電信號研究過程中，通常將癲癇腦電信號構建成復雜網絡，然后利用復雜網絡的靜態特性，如網絡的度、聚類系數等來研究癲癇腦電信號[9]。
將癲癇腦電信號轉換成復雜網絡已經成為研究癲癇腦電信號的重要方法，但卻忽略了逆向研究，即將構建的復雜網絡再轉換成時間序列分析。近年來，已有將復雜網絡轉換成時間序列的方法開展相應研究[10]。Yutaka等[11]提出了特征向量法將網絡轉換成時間序列；Weng等[12]提出了隨機游走方法將網絡轉換成時間序列。本文對癲癇腦電信號研究采取逆向研究思路，利用改進的k-最近鄰網絡算法將癲癇腦電信號構建成復雜網絡，然后利用特征向量法將復雜網絡轉換成時間序列，最后對轉換后的時間序列進行功率譜分析，還對前面改進的k-最近鄰網絡算法構建成的復雜網絡進行了網絡聚類系數分析，以期為癲癇腦電信號的動力學研究及臨床診斷提供重要的參考依據。
1 復雜網絡構建及網絡轉換成時間序列
1.1 相空間重構
混沌時間序列分析方法在很多科研和工程領域中已得到廣泛應用[13-14]。相空間重構是混沌時間序列分析的基礎。Takens等[15]提出了用延遲坐標法對混沌時間序列x=｛xi|i=1，2，…，N｝進行重構:
| $X={{X}_{i}}{{X}_{~}}i=[{{x}_{i}},{{x}_{i+t}},\ldots ,{{x}_{i+\left( m-1 \right)t}}],i=1,2,\ldots ,M$ |
其中，m為嵌入維，t為延遲，M=N-(m-1)t為相空間中的點數。Takens定理證明了如果嵌入維m≥2d+1，d為系統動力學階數，則重構的動力系統與原動力學系統在拓撲意義上等價。
本文將通過相關積分-相關積分(correlation integral-correlation integral,C-C)算法求延遲時間t和最佳嵌入維m[16]。在C-C算法的計算中將會用到相關積分：
| ${{C}^{m}}_{s}\left( m,N,r,t \right)=\frac{2}{M\left( M-1 \right)}{{\Sigma }_{1\le i\le j＜＜M}}\theta \{r-\|{{X}_{i}}-{{X}_{j}}\|\}$ |
| $\theta \left( a \right)=\left\{ \begin{array}{*{35}{l}} 1 & a\ge 0 \\ 0 & a0 \\ \end{array} \right.$ |
其中N為數據集合的長度，t是時間尺度，M=N-(m-1)t是m維相空間重構嵌入后的點數。C-C算法的研究與函數S(m,N,r,t)=C(m,N,r,t)－Cm(1,N,r,t)的值有關。將時間序列分為t個子序列，對S(m,N,r,t)的計算采用分塊平均的策略。S(m,N,r,t)可以表示為
| $\begin{align} & S\left( m,N,r,t \right)=\frac{1}{t}\sum\limits_{s=1}^{t}{\{{{C}_{s}}\left( m,N/t,r,t \right)-} \\ & {{C}^{m}}_{s}\left( 1,N/t,r,t \right)\} \\ \end{align}$ |
最后當N→∞,S(m,N,r,t)式將變為
| $S\left( m,r,t \right)=\frac{1}{t}\sum\limits_{s=1}^{t}{\{{{C}_{s}}\left( m,r,t \right)-}{{C}^{m}}_{s}\left( 1,r,t \right)\}$ |
現在定義差量ΔS(m,t)，其用表達式描述為
| $\begin{align} & \Delta S\left( m,t \right)=\max \left\{ S\left( m,{{r}_{j}},t,N \right) \right\}- \\ & \min \{S(m,{{r}_{j}},t,N)\} \\ \end{align}$ |
根據C-C算法統計理論，一般取m=2,3,4,5，${{r}_{j}}=\frac{i\sigma }{2}$ ，其中i=1,2,3,4[16]。根據m和rj取值分別求S(t)、ΔS(t)和Scor(t)。其表達式分別如下：
| $S\left( t \right)=\frac{1}{16}\sum\limits_{m-2}^{5}{\sum\limits_{j=1}^{4}{S(m,{{r}_{j}},t)}}$ |
| $\Delta S\left( t \right)=\frac{1}{4}\sum\limits_{m-2}^{5}{\Delta S\left( m,t \right)}$ |
| ${{S}_{cor}}\left( t \right)=\Delta S\left( t \right)+S\left( t \right)$ |
通過計算S(t)或者ΔS(t)的第一個極小值來確定重構相空間的延遲時間t，通過計算Scor的最小值確定延遲窗口τw，最后可以根據嵌入維和延遲時間的關系τw=(m-1)t來確定相空間重構的最佳嵌入維數m。
1.2 k-最近鄰網絡
k-最近鄰網絡構建思想是對于網絡中的任意一個結點尋找周圍最近的k個結點與之相連接[17]。每個點都以這個規律和其他點相連接。具體步驟如下：
(1) 首先要計算相空間中的歐式距離矩陣Dij,在距離矩陣的第i行包含其他點與第i個點之間的距離。
(2) 對于每一行i，根據距離矩陣選出距離i的k個最鄰近點，可以表示成Mi=｛x(j1)…x(jk)｝，其中Mi是對于第i行的k個最鄰近點的集合，x(j1)為第一個鄰近點，x(jk)為第k個鄰近點。
(3) 構建鄰接矩陣R，矩陣R中的第i行第j列表示為Rij。對于鄰接矩陣的第i行，如果第j個點屬于集合Mi，則Rij=1，否則Rij=0。按照這種方法完成對每一行的鄰接矩陣構建。
通過以上3個步驟，可以得到一個有向網絡的鄰接矩陣。該網絡即為k-最近鄰網絡。該網絡的度分布為固定值，即P(kout)≡δ(k)，邊的總數為N×k個，其中N為相空間點的總數，k為設置的k個鄰近點數目。
1.3 改進的k-最近鄰網絡
k-最近鄰網絡是一個有向網絡，并且擁有固定度分布的良好特性。而利用求特征值法將網絡轉換成時間序列是基于無向網絡。為了利用k-最近鄰網絡內在的良好特性，必須把有向網絡變換成無向網絡。以下是構建改進的k-最近鄰網絡的過程。
(1) 首先要計算相空間中的歐式距離矩陣Dij，在距離矩陣的第i行包含其他點與第i個點之間的距離。
(2) 對于每一行i，根據距離矩陣對其他點進行從小到大排序，排除與自身距離的比較，假設相空間有N個點，可以表示成為Mi=｛x(j1)…x(jN-1)｝，其中Mi是對于第i行的N-1個最鄰近點的集合，x(j1)為第一個鄰近點，x(jN-1)為最后一個鄰近點。
(3) 構建鄰接矩陣R,矩陣R中的第i行第j列表示為Rij。對于鄰接矩陣的每一行i，首先判斷已經有多少個點與第i個點連接，即求Rij1的個數，其中j=1,…,i-1,i+1,…,N，假設為s個點與i連接。若k-s＞0，則從Mi中選取前k-s個點與i相連接，并置Rij=1，并且Rij=1，在Mi中選取的點必須符合Rij=0，其中j是屬于Mi中的點。若k-s≤0，則跳過該行，對R矩陣的下一行進行構建。
1.4 特征值法
特征值法將網絡轉換成時間序列，主要是利用網絡鄰接矩陣的特征值[11]。此方法的主要思路是求鄰接矩陣的最大特征值對應的特征向量，該特征向量即為網絡轉換時間后的時間序列。本文將利用求鄰接矩陣特征值法將網絡轉換成時間序列。
假設A(={aij})是一個N×N網絡的鄰接矩陣。定義頂點vi和vj之間距離為dij。如果aij=0(i≠j)，則dij=w(＞1)，否則dij=aij。平方距離矩陣D定義為D=｛dij2｝,接下來將D代入表達式
| $G=-\frac{1}{2}{{J}_{N}}D{{J}^{T}}_{N}$ |
其中${{J}_{N}}=E-\frac{1}{N}{{1}_{N}}{{1}^{T}}_{N}$，E是N×N的單位矩陣，1N是一個包含N個1的列向量。將G分解為特征值和特征向量表示為：
| $G=P\Lambda {{P}^{T}}=\left( P{{\Lambda }^{1/2}} \right)\times {{\left( P{{\Lambda }^{1/2}} \right)}^{T}}\equiv X{{X}^{T}}$ |
其中$\Lambda =\text{diag}({{\lambda }_{1}},{{\lambda }_{2}},\ldots ,{{\lambda }_{h}}),{{\Lambda }^{(1/2)}}=\text{diag}(\sqrt{{{\lambda }_{1}}},\sqrt{{{\lambda }_{2}}},\ldots ,\sqrt{{{\lambda }_{h}}}),P=({{p}_{1}}{{p}_{2}}\ldots {{p}_{h}}){{p}_{m}}={{\left( {{p}_{m1}}{{p}_{m2}}\ldots {{p}_{mh}} \right)}^{T}}$，h是矩陣G非零特征值的個數。矩陣X=(x1，x2，…，xN)T，其中xm=(xm1,xm2,…,xmN)T是頂點Vm在h維空間中相應的值。定義時間序列為
| ${{s}_{m}}\left( t \right)=\sqrt{{{\lambda }_{m}}}{{p}_{m1}}\left( 1\le m\le h,1\le t\le N \right)$ |
最大特征值對應主特征向量，因此在本文研究中，都是選取最大特征值對應的時間序列。
2 實驗結果及其分析
2.1 實驗數據預處理及統計方法
本文實驗數據來源于南京軍區總醫院從臨床診斷中采集的數據，包含癲癇患者腦電信號組和正常健康者腦電信號。實驗分成兩組進行，分別為正常組和癲癇患者組。正常組和癲癇患者組的人數為各15人，年齡為20~60歲。腦電信號記錄采用標準10~20系統，使用放置有16個電極(FP1,FP2,F3,F4,C3,C4,P3,P4,O1,O2,F7,F8,T3,T4,T5,T6)的電極帽采集，腦電信號記錄數據采樣頻率為512 Hz，采集數據時受試者都處于清醒狀態。實驗選取了C3導聯腦電信號作為分析對象，并對該導聯數據利用帶陷濾波器濾除50 Hz工頻干擾。
本文采用的統計方法是t檢驗，t檢驗常用于檢測樣本含量少的兩樣本差異是否具有統計學意義[7-8]。當P值小于0.05時，表示兩樣本差異具有統計學意義，能有效進行區分；否則表示兩樣本差異無統計學意義，不能進行區分。本文采用SPSS軟件對實驗數據進行t檢驗。
2.2 相空間重構及網絡構建
復雜網絡構建前必須進行相空間的重構，進行相空間重構前要選擇合適的延遲時間和嵌入維。本文利用C-C算法來獲取數據相空間重構的延遲時間和嵌入維。在C-C算法中，S(t)或者ΔS(t)首次出現極小值處可以確定相空間延遲時間，Scor(t)首次出現極小值處用來確定延遲窗口。S(t)、ΔS(t)和Scor(t)隨延遲時間t變化的結果如圖 1所示。縱軸表示的含義是C-C算法中用到的計算公式在不同延遲時間下的均值，橫軸表示延遲時間，其變化范圍是0~50 s。
圖1
C-C算法實驗結果
Figure1.
Experimental result of the C-C algorithm
如圖 1所示，S(t)、ΔS(t)和Scor(t)極小值首次出現都是在第10 s處，因此最佳延遲時間t=10 s，最佳延遲窗口τw為10 s。根據延遲時間、嵌入維和延遲窗口三者關系τw=(m-1)t可以得到最佳嵌入維m為2，本文將以延遲時間t=10 s，嵌入維m=2重構相空間。
本文運用改進的k-最近鄰網絡的方法構建規模為500個結點的網絡。為驗證改進k-最近鄰網絡的可行性，將k=20和k=300時的改進k-最近鄰網絡度分布圖分別作圖，如圖 2所示，橫軸表示網絡結點度大小，縱軸表示度分布的概率大小。
圖2
不同網絡下的度分布
Figure2.
Degree distribution under different networks
改進的k-最近鄰網絡將有向網絡變換成無向網絡，并且該網絡度為k的概率接近為1，保留了k-最近鄰網絡的特征。說明改進的k-最近鄰網絡的方法是可行的。
2.3 時間序列的功率譜分析
首先，運用改進的k-最近鄰網絡方法構建網絡，構建過程中取k值為20，然后以網絡轉換成時間序列中的w為自變量，w的變量范圍取值是1~2，間隔取0.2，計算相應w下產生的時間序列的功率譜的最大值，再對每個w下對應的功率譜最大值全體求平均。
對正常組和癲癇患者組樣本分別構建480×480大小的網絡，然后將網絡根據不同w值得到不同時間序列，對不同時間序列求功率譜最大值，再將不同w下對應的功率譜最大值求平均。實驗結果如表 1、表 2所示。
表1
(正常組)在不同時間序列下功率譜最大值的平均值(以w因子為變量)
Table1.
(for the normal group) Mean values of maximum power spectrum during different time series versus w factor
表2
(癲癇患者組)不同時間序列下功率譜最大值的平均值(以w因子為變量)
Table2.
(for the group of patients with epilepsy) Mean values of maximum power spectrum during different time series versus w factor
對表 1和表 2中數據進行t檢驗，P=0.015。該實驗結果中正常人組腦電的總體均值大于癲癇患者組腦電的總體均值，這與直接分析腦電信號功率譜最大值的結果相反，所以該研究只能進行正常人和癲癇患者的區分，網絡轉換成時間序列后退化了其本質意義。
如果不對腦電數據進行網絡變換和時間序列變化，即直接對原始信號進行功率譜分析，實驗結果如表 3、表 4所示。
表3
(正常組)直接對原始腦電信號進行求功率譜最大值
Table3.
(for the normal group )Maximum value of the power spectrum of the original electroencephalogram signal
表4
(癲癇患者組) 直接對原始腦電信號進行求功率譜最大值
Table4.
(for the group of patients with epilepsy) Maximum value of the power spectrum of original electroencephalogram signal
對表 3和表 4的實驗結果進行t檢驗,P=0.081，兩組間差異沒有統計學意義，不能對兩組受試者進行區分。而前述網絡轉換時間序列的t檢驗結果有統計學意義，說明對經過變換后的信號求功率譜比直接對原始腦電信號求功率譜更容易區分正常人和癲癇患者。
以上是以w為變量求得不同時間序列下功率譜的最大值。以下將w設置一個固定值，不妨令w=1.6，然后以k為變量，根據k值變化得到不同的網絡，再求網絡的功率譜的最大值，其中k值的取值范圍是20~300，取值間隔為20。實驗結果如表 5、表 6所示。
表5
(正常組)不同時間序列下功率譜最大值的平均值(以k為變量)
Table5.
(for the normal group) Mean values of maximum power spectrum during different time series versus k
表6
(癲癇患者組)不同時間序列下功率譜最大值的平均值(以k為變量)
Table6.
(for the group of patients with epilepsy) Mean values of maximum power spectrum during different time series versus k
根據表 5、表 6發現，隨著k值不斷增大得到的網絡轉換的時間序列的最大功率譜的總體均值中，正常人的功率譜最大值的樣本總體均值小于癲癇患者。對表 5和表 6的數據進行t檢驗，P=0.004，差異具有統計學意義，因此也可用于區分正常人和癲癇患者。
表7
(正常組)構建改進的k-最近鄰網絡的聚類系數
Table7.
(for the normal group) Clustering coefficient of network based on improved k-nearest network
表8
(癲癇患者組) 構建改進的k-最近鄰網絡的聚類系數
Table8.
(for the group of patients with epilepsy) Clustering coefficient of network based on improved k-nearest network
由上述結果可知，無論以w因子還是k為變量，對經過網絡轉換生成的時間序列求功率譜，均可區分正常人和癲癇患者，而直接對原始腦電時間序列求功率譜則兩組間差異無統計學意義，足見本文方法具有明顯的優越性。
2.4 改進的k-最近鄰網絡的聚類系數分析
實驗以k為變量，根據k值的變化求得相應網絡的聚類系數，其中k值的變化范圍是20~300，變化的間隔是40。對實驗結果進行30次平均，然后取其均值作為實驗的結果。實驗結果如表 7、表 8所示。
數據表明正常組的網絡聚類系數總體均值小于癲癇患者組。對表 7和表 8的數據進行t檢驗，P=0.007，表明對正常組和癲癇患者組構建的網絡其聚類系數具有明顯差異，這種差異性也可能用于區分正常組和癲癇患者組。
3 結論
本文主要是基于改進的k-最近鄰網絡轉換成時間序列的癲癇研究。改進的k-最近鄰網絡是一個無向網絡，并且該網絡度為k的概率接近1，具有k-最近鄰網絡的特性。將網絡轉換成時間序列是運用求鄰接矩陣的特征值特征向量方法，實驗中是選取最大特征值對應的時間序列進行研究。
通過實驗結果表明，基于改進的k-最近鄰網絡轉換而來的時間序列功率譜分析比直接對原始腦電信號數據作功率分析更容易區分正常組和癲癇患者組。實驗結果還發現， 正常組和癲癇患者組間改進的k-最近鄰網絡的聚類系數也具有明顯差異。本研究將網絡轉換成時間序列的思想應用于癲癇研究，網絡構建采用改進的k-最近鄰網絡算法，經過實驗，結果證明本文方法所得特征值在正常組和癲癇患者組間有明顯差異，這不僅能夠為癲癇領域研究提供新的研究思路，而且能夠為癲癇臨床診斷提供相關參考依據。但由于本文僅僅是對網絡轉換成時間序列的功率譜進行研究，將來還可從其他角度來分析網絡轉換后的時間序列。