基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法与流程

本发明本发明属于生物医学领域，具体涉及一种检测脑电信号间效应连通性的方法。

背景技术：

大脑不同区域信号之间的效应连通性对于确定病灶区域具有重要作用。目前对脑电信号间效应连通性的检测包括时域分析和频域分析。pdc(partialdirectedcoherence，偏直接相干函数)是一种常用的从频域上分析信号之间因果关系的方法，该方法可以辨别出两个信号之间直接或者间接的因果关系。但该方法是基于线性模型的，而真实的eeg(electroencephalography,脑电图)是具有非线性特性的，pdc算法无法完全检测到信号中非线性因果关系。npdc(nonlinearpdc，非线性pdc)是基于pdc所改进的适用于narx(nonlinearautoregressiveexogenous，非线性自回归模型)的方法，它可以检测到信号之间非线性的相互影响。目前常用的为二维npdc，二维npdc算法是一种基于二维pdc算法改进的适用于二维narx模型的大脑效应连通性算法，该算法能够检测信号之间的线性与非线性因果关系。但是二维npdc算法仅能处理二维信号，处理多维信号时无法区分信号之间的直接和间接因果关系。

技术实现要素：

发明目的：针对现有技术中的不足，本发明提供了一种基于三维npdc的脑电信号效应连通性检测方法，可以检测三维信号之间的因果关系。

技术方案：本发明采用如下技术方案：

基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，包括如下步骤：

(1)构造单输入多输出的非线性自回归模型，所述模型为：

其中为脑电信号在时刻k的采样值，n为模型的非线性阶数，nout为输出信号的个数，q、p分别是输入信号u(k-ki)和输出信号的非线性程度，且p+q≤n；ki为信号的延迟值，k为模型的延迟阶数，ki≤k，是u(k-ki)和的线性或非线性组合的模型项，为模型项的系数，为脑电信号的采样值和模型预测值之间的误差项；

(2)应用frols算法对步骤(1)构造的模型进行系数估计；

(3)对三维pds进行形式变换，得到用频率响应函数描述的信号yi对yj的pdc的定义式；

(4)应用volterra级数核函数的多维傅里叶变换对simonarx模型进行频域分析，计算出模型的非线性频率响应函数；

(5)将步骤(4)计算出的非线性频率响应函数代入步骤(3)中的pdc定义式，得到三维npdc，得出在同时考虑三维信号的情况下某一信号对另一信号的因果影响。

所述步骤(2)具体包括：

(2-1)将非线性自回归模型改写为线性参数形式：

其中pl(n)是脑电信号y(n-ki)和u(n-ki)的线性或非线性组合的模型项，l为候选项个数，θl为线性模型系数，e(n)为线性模型误差项；

(2-2)将步骤(2-1)中的线性参数形式模型转换为正交模型：

其中wl(n)是相互正交的，gl为正交模型的系数；

(2-3)令d＝{p1,p2,…,pl}为由l个候选基底所组成的初始字典，pl＝[pl(1),pl(2),…,pl(n)]^t，令ql＝pl且σ＝y^ty，l＝1,2,…,l，计算向量y为由脑电信号采样值构成的向量；

令即err[m1]＝max{err⁽¹⁾[l]:1≤l≤l}，则第一个重要的模型项被选出，第一个正交向量可以被选为

(2-4)假设在算法的第s-1步，选出了一个子集ds-1，由s-1个重要的模型项组成，所述s-1个模型项经过正交变换组成了新的正交向量q1,q2,…,qs-1；在算法的第s步，令l≠m1,l≠m2,…,l≠ms-1，对于l＝1,2,…,l，计算：

令那么第s个重要的模型项就可以被选出来，第s个正交向量可以经过正交变换得到

(2-5)重复步骤(2-4)，当被筛选出的所有模型项的err值的和达到预先设定的阈值，就停止筛选过程；

(2-6)设经过上述步骤从l个候选项中选出的l0个重要模型项的线性组合，l0<l，即正交模型为：其等价于其中参数向量可由计算得到。

所述步骤(3)中对三维pdc进行形式变换具体包括：

(3-1)利用三维线性自回归模型对时域中长度为n的三维脑电信号y1、y2和y3进行建模：

其中，yi(n)为脑电信号yi在时刻n的值，yi(n-k)称为模型项，k为模型项的延迟值，p为延迟阶数，k≤p；aij(k)为模型的系数，ei(n)为信号yi的采样值和由模型得出的预测值之间的误差项；

(3-2)将式(4)右边的yi(n-k)项移到左边，然后再对两边进行傅里叶变换，可得到：

其中，yi(f)为信号yi的频谱，ei(f)为误差项ei的频谱，a(f)为模型的频域系数矩阵，矩阵a(f)中的元素arl(f)可按式进行计算：

(3-3)式(5)中矩阵a(f)中的元素按照如下方式进行构造：

借助于式(6)，式(5)可以改写为另一种形式：

将式(7)展开：

(3-4)信号yi对yj的pdc的定义式，可以更改为用式中频率响应函数描述的表达式：

其中

所述步骤(4)具体包括：

(4-1)在某系统中存在m个信号，y1,y2,…,ym，利用simonarx模型对其进行建模，即将其中一个信号看作输入信号u，其余信号看作输出信号yjk，如式(1)所示；simonarx模型在频域中可以用volterra级数核函数的多维傅里叶变换描述这个模型中的输入输出关系：

其中，u为输入信号u的频谱，称为simonarx模型的n阶广义频响函数，由n阶volterra级数核函数的多维傅里叶变换定义；求取的过程如下：

其中

更高阶的gfrf由更低阶的gfrf递归计算得到，递归终止于一阶gfrf：

为了形成一个类似于线性系统的统一的频域表达式，式(10)可以改写为：

(4-2)对于三维信号y1、y2和y3，利用simonarx模型对其进行建模；以y1和y2为输出信号、y3为输入信号并表示成u为例，信号y1和y2根据式(1)可表示为：

式(11)具体化到信号y1和y2：

其中，和按照步骤(4-1)所述的推导过程计算；

式(12)中的函数h3→1和h3→2即为simonarx模型的非线性频率响应函数；同理，将y1看作输入u或y2看作输入u，剩余两信号看作输出，可得：

以y1和y2为输出信号y3为输入信号为例，信号y2的整体频谱也可表示为内在影响和因果影响之和，信号y2总的频谱结构，即y2可表示为：

其中，e2是e2的频谱，由式可计算出函数同理可得：

通过上述过程的计算，得到了h3→2、h1→2、h2→1、h3→1、h1→3、h2→3和

在同时考虑三维信号的情况下脑电信号yi对yj的因果影响为：

其中

有益效果：与现有技术相比，本发明公开的基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，将三维pdc扩展至非线性，提供了一种适用于单输入多输出非线性自回归模型的三维npdc算法，该算法不仅能够识别出三维脑电信号之间的线性和非线性因果关系，而且还可以区分信号之间直接和间接的因果关系。本发明公开的方法有助于在术前诊断中精确定位病灶区域。

附图说明

图1为实施例一中模型17和模型18中信号之间相互作用的关系图；

图2为实施例一中模型17和模型18中三个信号的频谱图；

图3为实施例一中模型17中三个信号之间二维npdc的结果示意图；

图4为实施例一中模型18中三个信号之间二维npdc的结果示意图；

图5为实施例一采用pdc检测模型18中信号y3对y2因果影响的结果示意图；

图6为实施例一中模型17中三个信号之间三维npdc的结果示意图；

图7为实施例一中模型18中三个信号之间三维npdc的结果示意图；

图8为实施例二中eeg信号模型中三个信号的频谱图；

图9为实施例二中对0～40hz频带取均值处理后的二维npdc检测结果示意图；

图10为实施例二中采用二维npdc检测eeg信号模型中信号y1、y2和y3之间的相互作用关系图；

图11为实施例二中对0～40hz频带取均值处理后的三维npdc检测结果示意图；

图12为实施例二中采用三维npdc检测eeg信号模型中信号y1、y2和y3之间的相互作用关系图；

图13为本发明公开的脑电信号间效应连通性检测方法的流程图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

如图13所示，基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，包括如下步骤：

步骤1、构造单输入多输出的非线性自回归模型(single-inputmulti-outputnonlinearautoregressiveexogenous,simonarx)，所述模型为：

其中为脑电信号在时刻k的采样值，n为模型的非线性阶数，nout为输出信号的个数，q、p分别是输入信号u(k-ki)和输出信号的非线性程度，且p+q≤n；ki为信号的延迟值，与输入项u(k-ki)和输出项相对应，k为模型的延迟阶数，ki≤k，是u(k-ki)和的线性或非线性组合的模型项，为模型项的系数，为脑电信号的采样值和模型预测值之间的误差项；

步骤2、应用frols(theforwardregressionorthogonalleastsquares)算法对步骤(1)构造的模型进行系数估计；具体包括：

(2-1)对于非线性自回归模型进行参数估计是较为复杂的，本发明应用frols方法对信号{y1(n),y2(n),y3(n)}的系数进行估计，核心思想是对式(1)中的模型定义引入一个额外的模型，将非线性自回归模型改写为线性参数形式：

其中pl(n)是式(1)中脑电信号y(n-ki)和u(n-ki)的线性或非线性组合的模型项，l为候选项个数，θl为线性模型系数，e(n)为线性模型误差项；

(2-2)式(2)是非正交化的，将步骤(2-1)中的线性参数形式模型转换为正交模型：

其中wl(n)是相互正交的，gl为正交模型的系数；

(2-3)令d＝{p1,p2,…,pl}为由l个候选基底所组成的初始字典，pl＝[pl(1),pl(2),…,pl(n)]^t，令ql＝pl且σ＝y^ty，l＝1,2,…,l，计算向量y为由脑电信号采样值构成的向量；

err(errorreductionratio，误差减少比率)提供了一个简单有效的标准度量每个模型候选项的重要程度，根据err值选择重要项丢弃非重要项。令即err[m1]＝max{err⁽¹⁾[l]:1≤l≤l}，则第一个重要的模型项被选出，第一个正交向量可以被选为

(2-4)假设在算法的第s-1步，选出了一个子集ds-1，由s-1个重要的模型项组成，这s-1个模型项经过正交变换组成了新的正交向量q1,q2,…,qs-1；在算法的第s步，令l≠m1,l≠m2,…,l≠ms-1，对于l＝1,2,…,l，计算：

令那么第s个重要的模型项就可以被选出来，第s个正交向量可以经过正交变换得到

(2-5)重复步骤(2-4)，当被筛选出的所有模型项的err值的和达到预先设定的阈值，就停止筛选过程；

(2-6)设经过上述步骤从l个候选项中选出的l0个重要模型项的线性组合，l0<l，即正交模型为：其等价于其中参数向量可由计算得到。

步骤3、对三维pds进行形式变换，得到用频率响应函数描述的信号yi对yj的pdc的定义式；

对三维pdc进行形式变换具体包括：

(3-1)利用三维线性自回归模型对时域中长度为n的三维脑电信号y1、y2和y3进行建模：

其中，yi(n)为脑电信号yi在时刻n的值，yi(n-k)称为模型项，k为模型项的延迟值，p为延迟阶数，k≤p；aij(k)为模型的系数，ei(n)为信号yi的采样值和由模型得出的预测值之间的误差项；

(3-2)将式(4)右边的yi(n-k)项移到左边，然后再对两边进行傅里叶变换，可得到：

其中，yi(f)为信号yi的频谱，ei(f)为误差项ei的频谱，a(f)为模型的频域系数矩阵，矩阵a(f)中的元素arl(f)可按式进行计算：

(3-3)式(5)中矩阵a(f)中的元素按照如下方式进行构造：

借助于式(6)，式(5)可以改写为另一种形式：

将式(7)展开：

(3-4)以式(4)中第一个式子为例，分别把信号y2和y3当作该模型的输入，其余信号当作输出，也就是说整个模型可以看作是单输入多输出的线性系统。那么，函数h2→1和h3→1实际上是模型的频率响应函数，描述系统的输入输出关系。在式(8)中，i＝j时称为内在影响，i≠j时称为因果影响。以信号y1为例，是内在影响部分，只涉及自身对自身的影响，在式(8)中表现为误差项的频谱e1(f)与函数的乘积；和是因果影响部分，涉及信号y2和y3对信号y1的影响，在式(8)中表现为信号y2和y3的频谱y2(f)和y3(f)与频率响应函数h2→1(f)和h3→1(f)的乘积。

信号yi对yj的pdc的定义式，可以更改为用式中频率响应函数描述的表达式：

其中

步骤4、应用volterra级数核函数的多维傅里叶变换对simonarx模型进行频域分析，计算出模型的非线性频率响应函数；具体包括：

(4-1)在某系统中存在m个信号，y1,y2,…,ym，利用simonarx模型对其进行建模，即将其中一个信号看作输入信号u，其余信号看作输出信号如式(1)所示；simonarx模型在频域中可以用volterra级数核函数的多维傅里叶变换描述这个模型中的输入输出关系：

其中，u(fi)为输入信号u(k-ki)的频谱，称为simonarx模型的n阶广义频响函数(generalizedfrequencyresponsefunction,gfrf)，由n阶volterra级数核函数的多维傅里叶变换定义；求取的过程如下：

其中

更高阶的gfrf由更低阶的gfrf递归计算得到，递归终止于一阶gfrf：

为了形成一个类似于线性系统的统一的频域上的表达式，式(10)可以改写为：

(4-2)对于三维信号y1、y2和y3，利用simonarx模型对其进行建模；以y1和y2为输出信号、y3为输入信号并表示成u为例，信号y1和y2根据式(1)可表示为：

式(11)具体化到信号y1和y2：

其中，和按照步骤(4-1)所述的推导过程计算；

式(12)中的函数h3→1和h3→2即为simonarx模型的非线性频率响应函数；同理，将y1看作输入u或y2看作输入u，剩余两信号看作输出，可得：

以y1和y2为输出信号y3为输入信号为例，信号y2的整体频谱也可表示为内在影响和因果影响之和，信号y2总的频谱结构，即y2可表示为：

其中，e2是e2的频谱，由式可计算出函数同理可得：

通过上述过程的计算，得到了h3→2、h1→2、h2→1、h3→1、h1→3、h2→3和

步骤5、将步骤(4)计算出的非线性频率响应函数代入步骤(3)中的pdc定义式，得到三维npdc，得出在同时考虑三维信号的情况下某一信号对另一信号的因果影响为：

其中

实施例一：

本实施例以自回归模型生成的模拟信号为例说明本发明的步骤。

首先利用线性和非线性自回归模型生成模拟信号数据，应用三维npdc算法检测信号之间的因果关系，以此来验证该算法的有效性。本实施例采用如下式所示的三维自回归模型生成模拟数据：

其中，e1(n)、e2(n)和e3(n)为均值为0方差为0.01的高斯白噪声。可以看出在式(17)表示的模型中(简称模型17)，信号y1对y2、信号y2对y3以及信号y3对y1存在直接影响，而信号y2对y1、信号y3对y2以及信号y1对y3存在间接影响。这些因果关系都是线性的。在式(18)表示的模型中(简称模型18)，除了存在与模型17相同的因果影响之外，还存在一个信号y3对y2直接的非线性的因果影响。模型17和18中三个信号之间的相互作用关系分别如图1-(a)和图1-(b)所示，图中箭头表示作用的方向，实线箭头表示存在线性影响，虚线箭头表示存在非线性影响。

模型17和模型18中三个信号的频谱分别如图2-(a)和图2-(b)所示，图中从上至下依次为y1、y2和y3的频谱图。从该图中可以看出：在模型17中信号y1、y2和y3主要在频率60hz左右有较明显的波峰；在模型18中信号y1和y3主要在频率60hz左右有较明显的波峰，而信号y2主要在频率60hz和115hz左右有较明显的波峰。

实验中，分别用式(17)和式(18)生成模拟数据，共生成100组模拟数据，每组模拟数据的长度为n＝1024，采样频率为256hz。对每组数据，首先采用frols算法来估计模型系数，本实验中frols算法估计系数所需的阈值取为0.9999；然后用取得的系数按照前述步骤，分别求得pdc、二维npdc和三维npdc的结果，图3到图7给出了100次实验结果的平均值。

模拟模型17中三个信号之间二维npdc的结果如图3所示，图中npdcij表示信号yi对yj的npdc。如图3所示，在模型17中，二维npdc算法识别出了信号y1对y2和信号y3对y1在60hz左右处的影响，识别出的信号y2对y3的影响相对较弱。这代表二维npdc算法识别出某一信号的频率分量产生对另一信号的直接因果影响。然而，二维npdc算法也识别出了信号y2对y1、信号y3对y2以及信号y1对y3在60hz左右处的影响。这代表二维npdc算法识别出某一信号的频率分量产生对另一信号的间接因果影响。

模拟模型18中三个信号之间二维npdc的结果如图4所示，图中npdcij表示信号yi对yj的npdc。

如图4所示，在模型18中，二维npdc算法识别出了信号y1对y2和信号y3对y1在60hz左右处的影响，而且识别出信号y3对y2在60hz左右和120hz左右处的影响，识别出的信号y2对y3的影响相对较弱。这代表二维npdc算法识别出某一信号的频率分量产生对另一信号的线性直接因果影响。然而，二维npdc算法也识别出了信号y1对y3在60hz左右处的影响，识别出的信号y2对y1的影响相对较弱。这代表二维npdc算法识别出某一信号的频率分量产生对另一信号的间接因果影响。实验结果表明，二维npdc算法不能区分信号之间的直接和间接因果影响。

采用pdc检测模型18中信号y3对y2因果影响，其pdc如图5所示。线性pdc只能识别出在60hz左右的影响，这代表信号y3的频率分量产生对信号y2的线性因果影响。而如图4所示的二维npdc算法则识别出60hz和115hz左右的影响，而115hz左右的频率分量在pdc中没有被识别出来，它揭示了信号y3对y2的非线性因果影响。这表示二维npdc算法识别出模型18中信号y3对y2的非线性因果影响。

三维npdc算法检测模型17和模型18的结果分别如图6和图7所示。从图6可以看到，npdc12、npdc23和npdc31在60hz左右都有一较大的波峰，这表明，三维npdc算法识别出了模型17中信号y1对y2、信号y2对y3和信号y3对y1的直接影响，而npdc21、npdc32和npdc13在整个频率段的值都很小，表明他们间不存在直接影响。从图7的结果中可以看出，三维npdc算法不仅检测出了在模型18中信号y1对y2、信号y2对y3和信号y3对y1的直接影响，而且也检测出了信号y3对y2的直接的非线性因果关系，此外npdc21和npdc13在整个频率段的值都几乎接近于零也表明了信号y2对y1，信号y1对y3没有直接的因果影响。需要注意的是，类似于二维npdc算法，三维npdc算法同样识别出模型18中信号y3对y2在60hz和120hz左右的影响。这表示二维npdc算法和三维npdc算法都能够识别信号y3对y2的非线性因果影响。

比较图3和图6的右上图，二维npdc21的结果揭示了信号y2对y1有因果影响，三维npdc21的结果表明了信号y2对y1没有直接因果影响，综合两个结果，可以得出的结论为：信号y2对y1有间接的因果影响，该结论与模型17中的实际情况一致。

由本实施例可以得出，三维npdc算法不仅能够识别出信号之间的线性和非线性因果关系，而且还可以区分信号之间直接和间接的因果关系。

实施例二：

本实施例以癫痫患者的eeg信号为例说明本发明的步骤。

为了探究癫痫患者发病时大脑不同区域间的效应连通性，本实施例对癫痫患者的eeg信号进行采样，选取了一组20通道的癫痫eeg信号进行研究。这20通道的信号名称分别为a2、a6、a11、b1、b6、b11、c1、c4、c9、d1、d5、f2、f8、h2、i2、p1、p4、p8、t1和t8。本实施例在这20通道的信号中选取信号p8、a2和i2组成eeg信号模型，分别对应信号y1、y2和y3，以此来探究eeg信号之间的相互作用关系。

eeg信号的采样数据的总长度为18432，采样频率为256hz，即采样总时长为72秒。其中，0～20秒为癫痫发作前期的eeg信号，21～52秒为癫痫发作期的eeg信号，53～72秒为癫痫发作后期的eeg信号。eeg信号模型中信号的频谱如图8所示，从上至下依次为y1、y2和y3的频谱。从该图中可以看出：在eeg信号模型中信号y1、y2和y3主要在0～40hz的频带范围内有较明显的波峰。

本实施例对eeg信号进行重叠加窗处理。选择一个窗长w与步长l，分别从eeg信号y1、y2和y3中提取长度为w的一部分窗信号。对每组窗信号，首先采用frols算法来估计模型系数，本实施例中frols算法估计系数所需的阈值取为0.9999；然后用取得的系数按照前述步骤，分别求得二维npdc和三维npdc的结果。然后，将窗向后滑动一段距离l，得到另一组窗信号，对该段信号进行上述处理。本实施例中，设定w＝1024，l＝512，可得到34组窗信号的检测结果。其中，1～9组、10～24组和25～34组的检测结果分别对应发作前期、发作期间和发作后期。

由于在eeg信号模型中信号y1、y2和y3主要在0～40hz频带范围内有较明显的波峰，本实施例中，二维npdc和三维npdc的检测结果分别对0～40hz频带取均值，来对比不同窗信号的实验结果。设定一个阈值，即npdcij不小于该阈值时可以认为信号yi对yj存在因果影响，小于该阈值时可以认为不存在因果影响。根据二维npdc和三维npdc的实验结果选取阈值为0.08。由此，绘制出eeg信号模型中信号y1、y2和y3分别在发作前期、发作期间和发作后期的相互作用关系图。图9到图12给出了对0～40hz频带取均值处理后的二维npdc和三维npdc的检测结果，以及信号之间的相互作用关系图。

图9为对0～40hz频带取均值处理后的二维npdc检测结果，图10为与图9对应的eeg信号模型中信号y1、y2和y3之间的相互作用关系图，图10-(a)至图10-(c)依次为癫痫发作前期、期间和后期eeg信号模型中信号y1、y2和y3之间的相互作用关系图。如图10所示，在发作前期，没有检测到信号之间的因果关系；在发作期间，检测到信号y3与y1之间相互的因果关系、信号y1对y2和信号y3对y2的因果影响；在发作后期，检测到信号y2与y1之间相互的因果关系、信号y1对y3和信号y2对y3的因果影响。

图11为对0～40hz频带取均值处理后的三维npdc检测结果，图12为与图11对应的eeg信号模型中信号y1、y2和y3之间的相互作用关系图，图12-(a)至图12-(c)依次为癫痫发作前期、期间和后期eeg信号模型中信号y1、y2和y3之间的相互作用关系图。如图12所示，在发作前期，没有检测到信号之间的因果关系；在发作期间，检测到信号y3与y1之间相互的因果关系和信号y1对y2的因果影响；在发作后期，检测到y1与y2之间相互的因果关系和信号y1对y3的因果影响。

比较图10和图12中的结果，首先，在癫痫发作前期(图10(a)和图12(a))，二维npdc和三维npdc给出了相同的结果，即三个信号之间还没有因果影响。其次，在癫痫发作期间(图10(b)和图12(b))，看信号y2与y3之间的因果联系，二维npdc的结果表明信号y3对y2有因果影响，可能是直接或者间接影响；而三维npdc的结果表明了信号y3对y2没有直接因果影响。综合两个结果，可以得出，信号y3对y2有间接(通过信号y1)的因果影响，对其他两个信号之间进行相似的分析，可以得出三个信号间的因果影响如图12(b)所示。最后，在癫痫发作后期(图10(c)和图12(c))，对信号y1与y2，y1与y3之间的因果联系，二维npdc和三维npdc都给出了相同的结果。再看信号y2与y3之间，二维npdc的结果揭示了信号y2对y3有因果影响(无法判断是直接还是间接影响)；三维npdc的结果表明了信号y2对y3没有直接因果影响。结合两个结果，可以判断：信号y2对y3有间接(通过信号y1)的因果影响，因此，可以得出在癫痫发作后期三个信号间的因果影响如图12(c)所示。

由本实施例可得，三维npdc算法不仅能够识别出癫痫eeg信号之间的线性和非线性因果关系，而且还可以区分信号之间直接和间接的因果关系。本发明公开的方法有助于在癫痫疾病的术前诊断中精确定位致痫区。