本發(fā)明屬于腦電分析領(lǐng)域，涉及一種新的癲癇灶點(diǎn)粗定位方法，尤其涉及基于eeg源成像的癲癇灶點(diǎn)粗定位方法，對(duì)于頭皮腦電有良好的效果，該方法基于腦電信號(hào)的時(shí)頻特征和形態(tài)特征，具有計(jì)算復(fù)雜度低的優(yōu)點(diǎn)。

背景技術(shù)：

癲癇是大腦神經(jīng)元突發(fā)性異常放電，導(dǎo)致短暫大腦功能障礙的一種慢性疾病，由于其突發(fā)性和反復(fù)性，給廣大患者帶來極大痛苦。由于現(xiàn)有藥物和手術(shù)等手段的局限性，一些頑固性癲癇患者得難以得到有效治療。及早對(duì)癲癇發(fā)作進(jìn)行預(yù)警，對(duì)于治療方法的建立和患者生活質(zhì)量的提高有重要意義。約有30％的難治性癲癇患者表現(xiàn)出抗藥性，需要手術(shù)切除癲癇灶點(diǎn)或者使用電刺激抑制癲癇發(fā)作。

目前，臨床診斷中癲癇灶點(diǎn)定位通常需集中神經(jīng)內(nèi)科和神經(jīng)外科癲癇專家，結(jié)合患者的發(fā)作狀態(tài)和多模態(tài)醫(yī)學(xué)成像技術(shù)進(jìn)行綜合判斷。人工癲癇灶點(diǎn)定位耗時(shí)耗力，依賴醫(yī)生的經(jīng)驗(yàn)?；谛畔⒂?jì)算的eeg源成像技術(shù)可以推測(cè)頭皮eeg活動(dòng)的源信號(hào)，有望自動(dòng)確定引發(fā)癲癇異常腦電的皮層可疑源灶點(diǎn)，為后續(xù)的精確定位縮小范圍，輔助醫(yī)生做出診斷。同時(shí)，seeg為非植入式信號(hào)采集方法，也為病人減輕了痛苦和風(fēng)險(xiǎn)。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的是提供一種基于eeg源成像的癲癇灶點(diǎn)粗定位方法，該方法通過參數(shù)化癲癇動(dòng)態(tài)信息進(jìn)行源成像，能夠?qū)Πd癇灶點(diǎn)位置做粗定位。

本發(fā)明方法通過以下步驟實(shí)現(xiàn)：

(1)腦電數(shù)據(jù)預(yù)處理：首先，將原始腦電數(shù)據(jù)通過帶通濾波去除偽跡，選取有效頻段進(jìn)行分析和處理；其次，將數(shù)據(jù)進(jìn)行分段，以便后續(xù)數(shù)據(jù)處理；

步驟(1)采用濾波器為2階巴特沃茲濾波器(butterworth)，濾波參數(shù)為1hz-70hz；數(shù)據(jù)分段采用滑動(dòng)時(shí)間窗的方法，時(shí)間窗長(zhǎng)度1秒，步長(zhǎng)0.2秒。

(2)盲源信號(hào)分離：根據(jù)腦電信號(hào)的在時(shí)域和頻域的特性，把步驟(1)處理得到的腦電信號(hào)數(shù)據(jù)段進(jìn)行參數(shù)化，分離出獨(dú)立成分；

步驟(2)中盲源信號(hào)分解采用獨(dú)立成分分析(independentcomponentsanalysis，簡(jiǎn)稱ica)的方法，采集到的sensor空間的腦電信號(hào)通過獨(dú)立成分分析方法，被轉(zhuǎn)化為sensor空間的統(tǒng)計(jì)獨(dú)立的信號(hào)源的線性組合。

盲源信號(hào)分離的主要操作如下：

第一步，對(duì)數(shù)據(jù)進(jìn)行中心化，減少幅度不確定性，并保持?jǐn)?shù)據(jù)區(qū)分度不變。中心化的計(jì)算表達(dá)式如下：

第二步，對(duì)中心化后的數(shù)據(jù)進(jìn)行白化，目的是消除特征的直接相關(guān)性，降低輸入的冗余性，具體做法是將原始向量乘以一個(gè)矩陣變成使得的協(xié)方差矩陣是i,即滿足

具體做法是，首先進(jìn)行特征值分解

e(xx^t)＝ede^t

得到特征向量矩陣e和特征值矩陣d,然后計(jì)算計(jì)算表達(dá)式如下，

第三步，將白化后的數(shù)據(jù)做主成分分析(principlecomponentanalysis),目的是降低數(shù)據(jù)特征維度，并去除噪聲。做法是把數(shù)據(jù)的n個(gè)特征用數(shù)目更少的m個(gè)特征取代，新特征是舊特征的線性組合，這些線性組合最大化樣本方差，盡量使這m個(gè)特征互不相關(guān)。

第四步，求解出獨(dú)立成分，即

有n個(gè)信號(hào)源s(s1,s2,…,sn)^t,s∈rⁿ,a是一個(gè)mixingmatrix,x是觀測(cè)到的信號(hào)，即

x＝as

在a，s均未知的情況下推斷出s

為了解決這個(gè)問題，我們做出三個(gè)假設(shè)

1.觀測(cè)信號(hào)數(shù)目m不小于源信號(hào)數(shù)目n，為了方便，取m＝n

2.源信號(hào)s的各個(gè)分量之間相互統(tǒng)計(jì)獨(dú)立

3.源信號(hào)的各個(gè)分量最多只有一個(gè)高斯分布

設(shè)s＝a^-1x＝wx,每個(gè)si有概率密度ps,且各自獨(dú)立，那么源信號(hào)的聯(lián)合分布為

似然函數(shù)為

用梯度下降法等優(yōu)化方法就可以求得w,進(jìn)而得到s。

(3)癲癇成分選擇：根據(jù)信號(hào)的時(shí)頻特性選擇與癲癇相關(guān)的成分

步驟(3)中采用了短時(shí)傅里葉變換(short-timefouriertransformation)和相關(guān)性，主要操作如下：

第一步，用短時(shí)傅里葉變換計(jì)算已知癲癇通道信號(hào)的時(shí)頻特征(timefrequencyrepresentation,tfr)均值eeg-tfr，

第二步，計(jì)算各獨(dú)立成分的時(shí)頻特征ic-tfr，

第三步，計(jì)算每個(gè)ic-tfr與eeg-tfr的相關(guān)性，計(jì)算公式如下：

第四步，將與eeg-tfr相關(guān)性超過閾值p＝0.1的獨(dú)立成分作為癲癇成分；

(4)源空間映射：根據(jù)腦電信號(hào)的在時(shí)域和頻域的特性，把步驟(3)處理得到的腦電信號(hào)數(shù)據(jù)段映射為源空間的信號(hào)；

步驟(4)中源空間映射采用邊界元法(boundaryelementmethod，簡(jiǎn)稱bem)的方法，采集到的sensor空間的腦電信號(hào)通過邊界元法，被轉(zhuǎn)化為源空間的腦電信號(hào)。

源空間映射的主要操作如下：

第一步，對(duì)數(shù)據(jù)進(jìn)行平滑處理，通過滑動(dòng)窗的方式，對(duì)步驟(3)中的數(shù)據(jù)段進(jìn)行均值濾波，參數(shù)為：時(shí)間窗長(zhǎng)度0.02秒，步長(zhǎng)0.02秒。

第二步，邊界元法求解；

要解決的問題：

ax＝c

其中

是n*t的sensorspace信號(hào)矩陣，n為電極數(shù)

是m*t的源信號(hào)矩陣，m為源空間的維度

a是n*m的正向傳播矩陣

bem方法是求解在t時(shí)刻下的上述方程，常用的求解方法為廣義最小余量法(gmres,generalizedminimalresidualmethod)，該方法是利用krylov向量子空間的迭代算法。

n次krylov向量子空間為

kn＝kn(a,c)＝span{c,ac,a²c,…,aⁿc}

gmres的基本求解步驟：

選擇初始向量x0∈rⁿ，計(jì)算r0＝c-ax0,β＝||r0||,v1＝r0/β；

設(shè)收斂殘差為∈,對(duì)于k＝1,2,...，直到收斂

對(duì)于j＝1,2,3,...,k,用arnoldi方法計(jì)算

求解最小二乘法問題

得到y(tǒng),進(jìn)一步得到xk＝x0+vky其中ei＝{1,0,0,…}^t,vk＝[v1,v2,...,vk]，rk＝c-axk,

如果||rk||/β＜∈,迭代收斂,rk就是解，否則返回步驟2

上述算法中，

||.||表示范數(shù)

表示修正的(k+1)*k階上三角hessenberg矩陣，矩陣元素是hij

vk是krylov矢量子空間。

本發(fā)明方法通過參數(shù)化癲癇動(dòng)態(tài)信息進(jìn)行源成像，能夠?qū)Πd癇灶點(diǎn)位置做粗定位。本發(fā)明方法基于腦電信號(hào)的時(shí)頻特征和形態(tài)特征，具有特征數(shù)量少、計(jì)算復(fù)雜度低、實(shí)時(shí)性好等優(yōu)點(diǎn)。本發(fā)明基于信息計(jì)算的eeg源成像技術(shù)可以推測(cè)頭皮eeg活動(dòng)的源信號(hào)，實(shí)現(xiàn)自動(dòng)確定引發(fā)癲癇異常腦電的皮層可疑源灶點(diǎn)，為后續(xù)的精確定位縮小范圍，輔助醫(yī)生做出診斷。本發(fā)明為非植入式信號(hào)采集方法，也為病人減輕了痛苦和風(fēng)險(xiǎn)，對(duì)于頭皮腦電有良好的效果，該方法基于腦電信號(hào)的時(shí)頻特征和形態(tài)特征，具有計(jì)算復(fù)雜度低的優(yōu)點(diǎn)。

附圖說明

圖1為癲癇灶點(diǎn)粗定位流程圖。

圖2為時(shí)間窗方法示意圖。

圖3為癲癇成分結(jié)果圖。

圖4為源空間映射結(jié)果圖。

具體實(shí)施方式

本發(fā)明結(jié)合附圖和實(shí)施例作進(jìn)一步的說明。

實(shí)施例1

腦電數(shù)據(jù)預(yù)處理

腦電信號(hào)預(yù)處理包括濾波和分段兩個(gè)步驟。

(1)原始腦電數(shù)據(jù)帶通濾波

通過帶通濾波可以有效去除原始腦電數(shù)據(jù)的偽跡，選取有效頻段進(jìn)行分析和處理。采用濾波器為2階巴特沃茲濾波器(butterworth)，濾波參數(shù)為1hz-70hz；

(2)數(shù)據(jù)分段

數(shù)據(jù)分段的目的在于，確立數(shù)據(jù)處理的最小單元，以及在實(shí)時(shí)數(shù)據(jù)處理中的計(jì)算頻率。這里采用滑動(dòng)時(shí)間窗的方法，如所示，具體參數(shù)為時(shí)間窗長(zhǎng)度1秒，步長(zhǎng)0.2秒。

具體方法為：令原始數(shù)據(jù)為x(t)，其中t＝1,2,3…n是數(shù)據(jù)點(diǎn)序數(shù)，采樣率為1000hz，則每個(gè)時(shí)間窗內(nèi)包含1000個(gè)點(diǎn)，下一個(gè)時(shí)間窗的起點(diǎn)為當(dāng)前時(shí)間窗起點(diǎn)后的0.2秒即200個(gè)數(shù)據(jù)點(diǎn)。

通過滑動(dòng)時(shí)間窗的方法，可以將數(shù)據(jù)分段，并且可以通過調(diào)整步長(zhǎng)和時(shí)間窗長(zhǎng)度控制定位效果?；瑒?dòng)時(shí)間窗進(jìn)行數(shù)據(jù)分段見圖2。

盲源信號(hào)分離

盲源信號(hào)分離采用獨(dú)立成分分析的方法，將數(shù)據(jù)分離獨(dú)立成分。

(1)數(shù)據(jù)中心化

中心化的步驟如下：

1.找出每個(gè)通道的信號(hào)的均值；

2.數(shù)據(jù)減去各自通道的均值；

(2)對(duì)中心化后的數(shù)據(jù)進(jìn)行白化

白化的步驟如下：

1.進(jìn)行特征值分解

e(xx^t)＝ede^t

2.得到特征向量矩陣e和特征值矩陣d,然后計(jì)算計(jì)算表達(dá)式如下，

其中x為原始信號(hào)

(3)對(duì)白化后的數(shù)據(jù)做主成分分析

主成分分析的步驟如下：

1.計(jì)算協(xié)方差矩陣

2.計(jì)算協(xié)方差矩陣的特征矢量和特征值

3.選取主成分

(3)求解出獨(dú)立成分

步驟如下：

1.用梯度下降法優(yōu)化

2.計(jì)算獨(dú)立成分s＝wx

癲癇成分選擇

癲癇成分選擇主要采用了短時(shí)傅里葉變換和相關(guān)性，具體步驟如下：

(1)tfr計(jì)算

1.用短時(shí)傅里葉變換計(jì)算已知癲癇通道信號(hào)的時(shí)頻特征(timefrequencyrepresentation,tfr)均值eeg-tfr

2.計(jì)算各獨(dú)立成分的時(shí)頻特征ic-tfr

(2)計(jì)算每個(gè)ic-tfr與eeg-tfr的相關(guān)性

(3)將與eeg-tfr相關(guān)性超過閾值p＝0.1的獨(dú)立成分作為癲癇成分選擇出的癲癇成分如圖3所示。

源空間映射

源空間映射采用邊界元法(boundaryelementmethod，簡(jiǎn)稱bem)的方法，采集到的sensor空間的腦電信號(hào)通過邊界元法，被轉(zhuǎn)化為源空間的腦電信號(hào)。

(1)均值濾波

均值濾波也稱為線性濾波，其采用的主要方法為時(shí)間窗平均法。實(shí)際計(jì)算中，采用步長(zhǎng)0.02秒，時(shí)間窗長(zhǎng)度為0.02秒的滑動(dòng)時(shí)間窗，將時(shí)間窗中點(diǎn)的數(shù)據(jù)值賦值為時(shí)間窗均值，再通過插值方法得到平滑過的數(shù)據(jù)段；這里均值濾波的時(shí)間窗參數(shù)選擇應(yīng)當(dāng)由數(shù)據(jù)實(shí)際屬性確定。

(2)邊界元法求解

1.對(duì)每一時(shí)刻t的信號(hào)c，運(yùn)用廣義最小余量法(gmres,generalizedminimalresidualmethod)求解t時(shí)刻源映射方程:

ax＝c

其中:

是n*t的sensorspace信號(hào)矩陣，n為電極數(shù)

是m*t的源信號(hào)矩陣，m為源空間的維度

a是n*m的正向傳播矩陣

2.根據(jù)得到的正向傳播矩陣a，計(jì)算原信號(hào)向量x

癲癇成分映射到源空間的結(jié)果如圖4所示。

實(shí)驗(yàn)結(jié)果

我們通過癲癇病人頭皮腦電數(shù)據(jù)測(cè)試本算法。測(cè)試所用頭皮腦電信號(hào)包括3位病人7次癲癇發(fā)作。數(shù)據(jù)通道數(shù)為32，采樣率1000hz。

能夠粗定位p1的2個(gè)灶點(diǎn)，p2的1個(gè)灶點(diǎn)，p3的4個(gè)灶點(diǎn)中的3個(gè)，有一定的準(zhǔn)確性和可行性，能夠?yàn)楹罄m(xù)植入式精確定位縮小范圍。

如圖3所示，大部分頭皮區(qū)域的信號(hào)基本不發(fā)生變化保持為背景信號(hào)，只有與灶點(diǎn)有關(guān)的頭皮位置信號(hào)才會(huì)發(fā)生顯著變化。背景區(qū)域的細(xì)微變化為癲癇成分選擇中誤選的少量非癲癇成分引起。

如圖4所示，大部分皮層區(qū)域的信號(hào)基本不發(fā)生變化保持為背景信號(hào)，只有與灶點(diǎn)有關(guān)的皮層位置信號(hào)才會(huì)發(fā)生顯著變化。