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

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

Info

Publication number
CN109124623B
CN109124623B CN201810554759.2A CN201810554759A CN109124623B CN 109124623 B CN109124623 B CN 109124623B CN 201810554759 A CN201810554759 A CN 201810554759A CN 109124623 B CN109124623 B CN 109124623B
Authority
CN
China
Prior art keywords
model
signal
dimensional
signals
nonlinear
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810554759.2A
Other languages
English (en)
Other versions
CN109124623A (zh
Inventor
杨淳沨
杨文琪
刘彦超
伍家松
孔佑勇
姜龙玉
杨冠羽
舒华忠
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southeast University
Original Assignee
Southeast University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southeast University filed Critical Southeast University
Priority to CN201810554759.2A priority Critical patent/CN109124623B/zh
Publication of CN109124623A publication Critical patent/CN109124623A/zh
Application granted granted Critical
Publication of CN109124623B publication Critical patent/CN109124623B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/369—Electroencephalography [EEG]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Surgery (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Signal Processing (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Psychiatry (AREA)
  • Animal Behavior & Ethology (AREA)
  • Physics & Mathematics (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Power Engineering (AREA)
  • Psychology (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Complex Calculations (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明公开了一种基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，包括如下步骤：(1)构造单输入多输出的非线性自回归模型；(2)应用FROLS算法对步骤(1)构造的模型进行系数估计；(3)对三维PDS进行形式变换，得到用频率响应函数描述的信号yi对yj的PDC的定义式；(4)应用Volterra级数核函数的多维傅里叶变换对SIMO NARX模型进行频域分析，计算出模型的非线性频率响应函数；(5)将步骤(4)计算出的非线性频率响应函数代入步骤(3)中的PDC定义式，得到三维NPDC，得出在同时考虑三维信号的情况下某一信号对另一信号的因果影响。该方法可以检测三维脑电信号之间的因果关系。

Description

基于三维非线性偏直接相干函数的脑电信号间效应连通性检 测方法
技术领域
本发明本发明属于生物医学领域，具体涉及一种检测脑电信号间效应连通性的方法。
背景技术
大脑不同区域信号之间的效应连通性对于确定病灶区域具有重要作用。目前对脑电信号间效应连通性的检测包括时域分析和频域分析。PDC(Partial DirectedCoherence，偏直接相干函数)是一种常用的从频域上分析信号之间因果关系的方法，该方法可以辨别出两个信号之间直接或者间接的因果关系。但该方法是基于线性模型的，而真实的EEG(Electro Encephalo Graphy,脑电图)是具有非线性特性的，PDC算法无法完全检测到信号中非线性因果关系。NPDC(Nonlinear PDC，非线性PDC)是基于PDC所改进的适用于NARX(Nonlinear AutoRegressive eXogenous，非线性自回归模型)的方法，它可以检测到信号之间非线性的相互影响。目前常用的为二维NPDC，二维NPDC算法是一种基于二维PDC算法改进的适用于二维NARX模型的大脑效应连通性算法，该算法能够检测信号之间的线性与非线性因果关系。但是二维NPDC算法仅能处理二维信号，处理多维信号时无法区分信号之间的直接和间接因果关系。
发明内容
发明目的：针对现有技术中的不足，本发明提供了一种基于三维NPDC的脑电信号效应连通性检测方法，可以检测三维信号之间的因果关系。
技术方案：本发明采用如下技术方案：
基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，包括如下步骤：
(1)构造单输入多输出的非线性自回归模型，所述模型为：
Figure GDA0002908620750000011
其中
Figure GDA0002908620750000012
为脑电信号
Figure GDA0002908620750000013
在时刻k的采样值，N为模型的非线性阶数，Nout为输出信号的个数，q、p分别是输入信号u(k-ki)和输出信号
Figure GDA0002908620750000014
的非线性程度，且p+q≤N；ki为信号的延迟值，K为模型的延迟阶数，ki≤K，
Figure GDA0002908620750000021
Figure GDA0002908620750000022
是u(k-ki)和
Figure GDA0002908620750000023
的线性或非线性组合的模型项，
Figure GDA0002908620750000024
为模型项的系数，
Figure GDA0002908620750000025
为脑电信号
Figure GDA0002908620750000026
的采样值和模型预测值之间的误差项；
(2)应用FROLS算法对步骤(1)构造的模型进行系数估计；
(3)对三维PDS进行形式变换，得到用频率响应函数描述的信号yi对yj的PDC的定义式；
(4)应用Volterra级数核函数的多维傅里叶变换对SIMO NARX模型进行频域分析，计算出模型的非线性频率响应函数；
(5)将步骤(4)计算出的非线性频率响应函数代入步骤(3)中的PDC定义式，得到三维NPDC，得出在同时考虑三维信号的情况下某一信号对另一信号的因果影响。
所述步骤(2)具体包括：
(2-1)将非线性自回归模型改写为线性参数形式：
Figure GDA0002908620750000027
其中pl(n)是脑电信号y(n-ki)和u(n-ki)的线性或非线性组合的模型项，L为候选项个数，θl为线性模型系数，e(n)为线性模型误差项；
(2-2)将步骤(2-1)中的线性参数形式模型转换为正交模型：
Figure GDA0002908620750000028
其中wl(n)是相互正交的，gl为正交模型的系数；
(2-3)令D＝{p1,p2,…,pL}为由L个候选基底所组成的初始字典，pl＝[pl(1),pl(2),…,pl(N)]T，令ql＝pl且σ＝yTy，l＝1,2,…,L，计算
Figure GDA0002908620750000029
Figure GDA00029086207500000210
向量y为由脑电信号采样值构成的向量；
令
Figure GDA00029086207500000211
即ERR[m1]＝max{ERR(1)[l]:1≤l≤L}，则第一个重要的模型项
Figure GDA0002908620750000031
被选出，第一个正交向量可以被选为
Figure GDA0002908620750000032
(2-4)假设在算法的第s-1步，选出了一个子集Ds-1，由s-1个重要的模型项组成，
Figure GDA0002908620750000033
所述s-1个模型项经过正交变换组成了新的正交向量q1,q2,…,qs-1；在算法的第s步，令l≠m1,l≠m2,…,l≠ms-1，对于l＝1,2,…,L，计算：
Figure GDA0002908620750000034
Figure GDA0002908620750000035
Figure GDA0002908620750000036
令
Figure GDA0002908620750000037
那么第s个重要的模型项
Figure GDA0002908620750000038
就可以被选出来，第s个正交向量可以经过正交变换得到
Figure GDA0002908620750000039
(2-5)重复步骤(2-4)，当被筛选出的所有模型项的ERR值的和达到预先设定的阈值，就停止筛选过程；
(2-6)设经过上述步骤从L个候选项
Figure GDA00029086207500000310
中选出的L0个重要模型项的线性组合，L0<L，即正交模型为：
Figure GDA00029086207500000311
其等价于
Figure GDA00029086207500000312
其中参数向量
Figure GDA00029086207500000313
可由
Figure GDA00029086207500000314
计算得到。
所述步骤(3)中对三维PDC进行形式变换具体包括：
(3-1)利用三维线性自回归模型对时域中长度为N的三维脑电信号y1、y2和y3进行建模：
Figure GDA00029086207500000315
其中，yi(n)为脑电信号yi在时刻n的值，yi(n-k)称为模型项，k为模型项的延迟值，p为延迟阶数，k≤p；aij(k)为模型的系数，ei(n)为信号yi的采样值和由模型得出的预测值之间的误差项；
Figure GDA0002908620750000045
(3-2)将式(4)右边的yi(n-k)项移到左边，然后再对两边进行傅里叶变换，可得到：
Figure GDA0002908620750000041
其中，Yi(f)为信号yi的频谱，Ei(f)为误差项ei的频谱，A(f)为模型的频域系数矩阵，矩阵A(f)中的元素Arl(f)可按式进行计算：
Figure GDA0002908620750000042
(3-3)式(5)中矩阵A(f)中的元素按照如下方式进行构造：
Figure GDA0002908620750000043
借助于式(6)，式(5)可以改写为另一种形式：
Figure GDA0002908620750000044
将式(7)展开：
Figure GDA0002908620750000051
(3-4)信号yi对yj的PDC的定义式，可以更改为用式中频率响应函数描述的表达式：
Figure GDA0002908620750000052
其中
Figure GDA0002908620750000058
所述步骤(4)具体包括：
(4-1)在某系统中存在M个信号，y1,y2,…,yM，利用SIMO NARX模型对其进行建模，即将其中一个信号看作输入信号u，其余信号看作输出信号
Figure GDA0002908620750000059
如式(1)所示；SIMO NARX模型在频域中可以用Volterra级数核函数的多维傅里叶变换描述这个模型中的输入输出关系：
Figure GDA0002908620750000053
其中，U为输入信号u的频谱，
Figure GDA0002908620750000054
称为SIMO NARX模型的n阶广义频响函数，由n阶Volterra级数核函数的多维傅里叶变换定义；求取
Figure GDA0002908620750000055
的过程如下：
Figure GDA0002908620750000056
其中
Figure GDA0002908620750000057
Figure GDA0002908620750000061
更高阶的GFRF由更低阶的GFRF递归计算得到，递归终止于一阶GFRF：
Figure GDA0002908620750000062
Figure GDA0002908620750000063
为了形成统一的频域表达式，式(10)改写为：
Figure GDA0002908620750000064
(4-2)对于三维信号y1、y2和y3，利用SIMO NARX模型对其进行建模；将y1和y2作为输出信号、y3作为输入信号并表示成u，信号y1和y2根据式(1)可表示为：
Figure GDA0002908620750000065
式(11)具体化到信号y1和y2：
Figure GDA0002908620750000066
其中，
Figure GDA0002908620750000067
和
Figure GDA0002908620750000068
按照步骤(4-1)所述求取
Figure GDA0002908620750000069
的过程计算；
式(12)中的函数
Figure GDA00029086207500000610
和
Figure GDA00029086207500000611
即为SIMO NARX模型的非线性频率响应函数；同理，将y1看作输入u或y2看作输入u，剩余两信号看作输出，可得：
Figure GDA0002908620750000071
将y1和y2作为输出信号，y3作为输入信号，信号y2的整体频谱也可表示为内在影响和因果影响之和，信号y2总的频谱结构，即Y2可表示为：
Figure GDA0002908620750000072
其中，E2是e2的频谱，由上式可计算出函数
Figure GDA0002908620750000073
同理可得：
Figure GDA0002908620750000074
通过上述过程的计算，得到了SIMO NARX模型的非线性频率响应函数
Figure GDA0002908620750000075
Figure GDA0002908620750000076
和
Figure GDA0002908620750000077
在同时考虑三维信号的情况下脑电信号yi对yj的因果影响为：
Figure GDA0002908620750000078
其中
Figure GDA0002908620750000079
有益效果：与现有技术相比，本发明公开的基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，将三维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-Input Multi-OutputNonlinear AutoRegressive eXogenous, SIMO NARX)，所述模型为：
Figure GDA0002908620750000081
其中
Figure GDA0002908620750000082
为脑电信号
Figure GDA00029086207500000811
在时刻k的采样值，N为模型的非线性阶数，Nout为输出信号的个数，q、p分别是输入信号u(k-ki)和输出信号
Figure GDA0002908620750000083
的非线性程度，且p+q≤N；ki为信号的延迟值，与输入项u(k-ki)和输出项
Figure GDA0002908620750000084
相对应，K为模型的延迟阶数，ki≤K，
Figure GDA0002908620750000085
Figure GDA0002908620750000086
是u(k-ki)和
Figure GDA0002908620750000087
的线性或非线性组合的模型项，
Figure GDA0002908620750000088
为模型项的系数，
Figure GDA0002908620750000089
为脑电信号
Figure GDA00029086207500000810
的采样值和模型预测值之间的误差项；
步骤2、应用FROLS(The Forward Regression Orthogonal Least Squares)算法对步骤(1)构造的模型进行系数估计；具体包括：
(2-1)对于非线性自回归模型进行参数估计是较为复杂的，本发明应用FROLS方法对信号{y1(n),y2(n),y3(n)}的系数进行估计，核心思想是对式(1)中的模型定义引入一个额外的模型，将非线性自回归模型改写为线性参数形式：
Figure GDA0002908620750000091
其中pl(n)是式(1)中脑电信号y(n-ki)和u(n-ki)的线性或非线性组合的模型项，L为候选项个数，θl为线性模型系数，e(n)为线性模型误差项；
(2-2)式(2)是非正交化的，将步骤(2-1)中的线性参数形式模型转换为正交模型：
Figure GDA0002908620750000092
其中wl(n)是相互正交的，gl为正交模型的系数；
(2-3)令D＝{p1,p2,…,pL}为由L个候选基底所组成的初始字典，pl＝[pl(1),pl(2),…,pl(N)]T，令ql＝pl且σ＝yTy，l＝1,2,…,L，计算
Figure GDA0002908620750000093
Figure GDA0002908620750000094
向量y为由脑电信号采样值构成的向量；
ERR(Error Reduction Ratio，误差减少比率)提供了一个简单有效的标准度量每个模型候选项的重要程度，根据ERR值选择重要项丢弃非重要项。令
Figure GDA0002908620750000095
即ERR[m1]＝max{ERR(1)[l]:1≤l≤L}，则第一个重要的模型项
Figure GDA0002908620750000096
被选出，第一个正交向量可以被选为
Figure GDA0002908620750000097
(2-4)假设在算法的第s-1步，选出了一个子集Ds-1，由s-1个重要的模型项组成，
Figure GDA0002908620750000098
这s-1个模型项经过正交变换组成了新的正交向量q1,q2,…,qs-1；在算法的第s步，令l≠m1,l≠m2,…,l≠ms-1，对于l＝1,2,…,L，计算：
Figure GDA0002908620750000099
Figure GDA00029086207500000910
Figure GDA0002908620750000101
令
Figure GDA0002908620750000102
那么第s个重要的模型项
Figure GDA0002908620750000103
就可以被选出来，第s个正交向量可以经过正交变换得到
Figure GDA0002908620750000104
(2-5)重复步骤(2-4)，当被筛选出的所有模型项的ERR值的和达到预先设定的阈值，就停止筛选过程；
(2-6)设经过上述步骤从L个候选项
Figure GDA0002908620750000105
中选出的L0个重要模型项的线性组合，L0<L，即正交模型为：
Figure GDA0002908620750000106
其等价于
Figure GDA0002908620750000107
其中参数向量
Figure GDA0002908620750000108
可由
Figure GDA0002908620750000109
计算得到。
步骤3、对三维PDS进行形式变换，得到用频率响应函数描述的信号yi对yj的PDC的定义式；
对三维PDC进行形式变换具体包括：
(3-1)利用三维线性自回归模型对时域中长度为N的三维脑电信号y1、y2和y3进行建模：
Figure GDA00029086207500001010
其中，yi(n)为脑电信号yi在时刻n的值，yi(n-k)称为模型项，k为模型项的延迟值，p为延迟阶数，k≤p；aij(k)为模型的系数，ei(n)为信号yi的采样值和由模型得出的预测值之间的误差项；
Figure GDA00029086207500001011
(3-2)将式(4)右边的yi(n-k)项移到左边，然后再对两边进行傅里叶变换，可得到：
Figure GDA0002908620750000111
其中，Yi(f)为信号yi的频谱，Ei(f)为误差项ei的频谱，A(f)为模型的频域系数矩阵，矩阵A(f)中的元素Arl(f)可按式进行计算：
Figure GDA0002908620750000112
(3-3)式(5)中矩阵A(f)中的元素按照如下方式进行构造：
Figure GDA0002908620750000113
借助于式(6)，式(5)可以改写为另一种形式：
Figure GDA0002908620750000114
将式(7)展开：
Figure GDA0002908620750000115
(3-4)以式(4)中第一个式子为例，分别把信号y2和y3当作该模型的输入，其余信号当作输出，也就是说整个模型可以看作是单输入多输出的线性系统。那么，函数H2→1和H3→1实际上是模型的频率响应函数，描述系统的输入输出关系。在式(8)中，
Figure GDA0002908620750000126
i＝j时称为内在影响，i≠j时称为因果影响。以信号y1为例，
Figure GDA00029086207500001210
是内在影响部分，只涉及自身对自身的影响，在式(8)中表现为误差项的频谱E1(f)与函数
Figure GDA0002908620750000127
的乘积；
Figure GDA00029086207500001211
和
Figure GDA00029086207500001212
是因果影响部分，涉及信号y2和y3对信号y1的影响，在式(8)中表现为信号y2和y3的频谱Y2(f)和Y3(f)与频率响应函数H2→1(f)和H3→1(f)的乘积。
信号yi对yj的PDC的定义式，可以更改为用式中频率响应函数描述的表达式：
Figure GDA0002908620750000121
其中
Figure GDA0002908620750000128
步骤4、应用Volterra级数核函数的多维傅里叶变换对SIMO NARX模型进行频域分析，计算出模型的非线性频率响应函数；具体包括：
(4-1)在某系统中存在M个信号，y1,y2,…,yM，利用SIMO NARX模型对其进行建模，即将其中一个信号看作输入信号u，其余信号看作输出信号
Figure GDA0002908620750000129
如式(1)所示；SIMO NARX模型在频域中可以用Volterra级数核函数的多维傅里叶变换描述这个模型中的输入输出关系：
Figure GDA0002908620750000122
其中，U(fi)为输入信号u(k-ki)的频谱，
Figure GDA0002908620750000123
称为SIMO NARX模型的n阶广义频响函数(Generalized Frequency Response Function，GFRF)，由n阶Volterra级数核函数的多维傅里叶变换定义；求取
Figure GDA0002908620750000124
的过程如下：
Figure GDA0002908620750000125
其中
Figure GDA0002908620750000131
Figure GDA0002908620750000132
更高阶的GFRF由更低阶的GFRF递归计算得到，递归终止于一阶GFRF：
Figure GDA0002908620750000133
Figure GDA0002908620750000134
为了形成统一的频域上的表达式，式(10)改写为：
Figure GDA0002908620750000135
(4-2)对于三维信号y1、y2和y3，利用SIMO NARX模型对其进行建模；将y1和y2作为输出信号、y3作为输入信号并表示成u，信号y1和y2根据式(1)可表示为：
Figure GDA0002908620750000136
式(11)具体化到信号y1和y2：
Figure GDA0002908620750000137
其中，
Figure GDA0002908620750000138
和
Figure GDA0002908620750000139
按照步骤(4-1)所述求取
Figure GDA00029086207500001310
的过程计算；
式(12)中的函数
Figure GDA0002908620750000141
和
Figure GDA0002908620750000142
即为SIMO NARX模型的非线性频率响应函数；同理，将y1看作输入u或y2看作输入u，剩余两信号看作输出，可得：
Figure GDA0002908620750000143
将y1和y2作为输出信号y3作为输入信号，信号y2的整体频谱也可表示为内在影响和因果影响之和，信号y2总的频谱结构，即Y2可表示为：
Figure GDA0002908620750000144
其中，E2是e2的频谱，由上式可计算出函数
Figure GDA0002908620750000145
同理可得：
Figure GDA0002908620750000146
通过上述过程的计算，得到了SIMO NARX模型的非线性频率响应函数
Figure GDA0002908620750000147
Figure GDA0002908620750000148
和
Figure GDA0002908620750000149
步骤5、将步骤(4)计算出的非线性频率响应函数代入步骤(3)中的PDC定义式，得到三维NPDC，得出在同时考虑三维信号的情况下某一信号对另一信号的因果影响为：
Figure GDA00029086207500001410
其中
Figure GDA00029086207500001412
实施例一：
本实施例以自回归模型生成的模拟信号为例说明本发明的步骤。
首先利用线性和非线性自回归模型生成模拟信号数据，应用三维NPDC算法检测信号之间的因果关系，以此来验证该算法的有效性。本实施例采用如下式所示的三维自回归模型生成模拟数据：
Figure GDA00029086207500001411
Figure GDA0002908620750000151
其中，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信号之间的线性和非线性因果关系，而且还可以区分信号之间直接和间接的因果关系。本发明公开的方法有助于在癫痫疾病的术前诊断中精确定位致痫区。

Claims (5)

1.基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，其特征在于，包括如下步骤：
(1)构造单输入多输出的非线性自回归模型，所述模型为：
Figure FDA0002908620740000011
其中
Figure FDA00029086207400000110
为脑电信号
Figure FDA0002908620740000019
在时刻k的采样值，N为模型的非线性阶数，Nout为输出信号的个数，q、p分别是输入信号u(k-ki)和输出信号
Figure FDA00029086207400000111
的非线性程度，且p+q≤N；ki为信号的延迟值，K为模型的延迟阶数，ki≤K，
Figure FDA0002908620740000012
Figure FDA0002908620740000013
是u(k-ki)和
Figure FDA0002908620740000014
的线性或非线性组合的模型项，
Figure FDA0002908620740000015
为模型项的系数，
Figure FDA0002908620740000016
为脑电信号
Figure FDA0002908620740000017
的采样值和模型预测值之间的误差项；
(2)应用FROLS算法对步骤(1)构造的模型进行系数估计；
(3)对三维PDS进行形式变换，得到用频率响应函数描述的信号yi对yj的PDC的定义式；
(4)应用Volterra级数核函数的多维傅里叶变换对SIMO NARX模型进行频域分析，计算出模型的非线性频率响应函数；
(5)将步骤(4)计算出的非线性频率响应函数代入步骤(3)中的PDC定义式，得到三维NPDC，得出在同时考虑三维信号的情况下某一信号对另一信号的因果影响。
2.根据权利要求1所述的基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，其特征在于，所述步骤(2)具体包括：
(2-1)将非线性自回归模型改写为线性参数形式：
Figure FDA0002908620740000018
其中pl(n)是脑电信号y(n-ki)和u(n-ki)的线性或非线性组合的模型项，L为候选项个数，θl为线性模型系数，e(n)为线性模型误差项；
(2-2)将步骤(2-1)中的线性参数形式模型转换为正交模型：
Figure FDA0002908620740000021
其中wl(n)是相互正交的，gl为正交模型的系数；
(2-3)令D＝{p1,p2,…,pL}为由L个候选基底所组成的初始字典，pl＝[pl(1),pl(2),…,pl(N)]T，令ql＝pl且σ＝yTy，l＝1,2,…,L，计算
Figure FDA0002908620740000022
Figure FDA0002908620740000023
向量y为由脑电信号采样值构成的向量；
令
Figure FDA0002908620740000024
即ERR[m1]＝max{ERR(1)[l]:1≤l≤L}，则第一个重要的模型项
Figure FDA0002908620740000025
被选出，第一个正交向量可以被选为
Figure FDA0002908620740000026
(2-4)假设在算法的第s-1步，选出了一个子集Ds-1，由s-1个重要的模型项组成，
Figure FDA0002908620740000027
所述s-1个模型项经过正交变换组成了新的正交向量q1,q2,…,qs-1；在算法的第s步，令l≠m1,l≠m2,…,l≠ms-1，对于l＝1,2,…,L，计算：
Figure FDA0002908620740000028
Figure FDA0002908620740000029
Figure FDA00029086207400000210
令
Figure FDA00029086207400000211
那么第s个重要的模型项
Figure FDA00029086207400000218
就可以被选出来，第s个正交向量可以经过正交变换得到
Figure FDA00029086207400000212
(2-5)重复步骤(2-4)，当被筛选出的所有模型项的ERR值的和达到预先设定的阈值，就停止筛选过程；
(2-6)设经过上述步骤从L个候选项
Figure FDA00029086207400000213
中选出的L0个重要模型项的线性组合，L0<L，即正交模型为：
Figure FDA00029086207400000214
其等价于
Figure FDA00029086207400000215
其中参数向量
Figure FDA00029086207400000216
可由
Figure FDA00029086207400000217
计算得到。
3.根据权利要求1所述的基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，其特征在于，所述步骤(3)中对三维PDC进行形式变换具体包括：
(3-1)利用三维线性自回归模型对时域中长度为N的三维脑电信号y1、y2和y3进行建模：
Figure FDA0002908620740000031
其中，yi(n)为脑电信号yi在时刻n的值，yi(n-k)称为模型项，k为模型项的延迟值，p为延迟阶数，k≤p；aij(k)为模型的系数，ei(n)为信号yi的采样值和由模型得出的预测值之间的误差项；
Figure FDA0002908620740000035
(3-2)将式(4)右边的yi(n-k)项移到左边，然后再对两边进行傅里叶变换，可得到：
Figure FDA0002908620740000032
其中，Yi(f)为信号yi的频谱，Ei(f)为误差项ei的频谱，A(f)为模型的频域系数矩阵，矩阵A(f)中的元素Arl(f)可按式进行计算：
Figure FDA0002908620740000033
(3-3)式(5)中矩阵A(f)中的元素按照如下方式进行构造：
Figure FDA0002908620740000034
借助于式(6)，式(5)可以改写为另一种形式：
Figure FDA0002908620740000041
将式(7)展开：
Figure FDA0002908620740000042
(3-4)信号yi对yj的PDC的定义式，为用频率响应函数描述的表达式：
Figure FDA0002908620740000043
其中
Figure FDA0002908620740000044
4.根据权利要求1所述的基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，其特征在于，所述步骤(4)具体包括：
(4-1)在某系统中存在M个信号，y1,y2,…,yM，利用SIMO NARX模型对其进行建模，即将其中一个信号看作输入信号u，其余信号看作输出信号
Figure FDA0002908620740000045
如式(1)所示；SIMO NARX模型在频域中可以用Volterra级数核函数的多维傅里叶变换描述这个模型中的输入输出关系：
Figure FDA0002908620740000046
其中，U为输入信号u的频谱，
Figure FDA0002908620740000047
称为SIMO NARX模型的n阶广义频响函数，由n阶Volterra级数核函数的多维傅里叶变换定义；求取
Figure FDA0002908620740000048
的过程如下：
Figure FDA0002908620740000051
其中
Figure FDA0002908620740000052
Figure FDA0002908620740000053
更高阶的GFRF由更低阶的GFRF递归计算得到，递归终止于一阶GFRF：
Figure FDA0002908620740000054
Figure FDA0002908620740000055
为了形成统一的频域表达式，式(10)改写为：
Figure FDA0002908620740000056
(4-2)对于三维信号y1、y2和y3，利用SIMO NARX模型对其进行建模；将y1和y2作为输出信号、y3作为输入信号并表示成u，信号y1和y2根据式(1)可表示为：
Figure FDA0002908620740000057
式(11)具体化到信号y1和y2：
Figure FDA0002908620740000061
其中，
Figure FDA0002908620740000062
和
Figure FDA0002908620740000063
按照步骤(4-1)所述求取
Figure FDA0002908620740000064
的过程计算；
式(12)中的函数
Figure FDA0002908620740000065
和
Figure FDA0002908620740000066
即为SIMO NARX模型的非线性频率响应函数；同理，将y1看作输入u或y2看作输入u，剩余两信号看作输出，可得：
Figure FDA0002908620740000067
将y1和y2作为输出信号，y3作为输入信号，信号y2总的频谱结构，即Y2可表示为：
Figure FDA0002908620740000068
其中，E2是e2的频谱，由上式可计算出函数
Figure FDA0002908620740000069
同理可得：
Figure FDA00029086207400000610
通过上述过程的计算，得到了SIMO NARX模型的非线性频率响应函数
Figure FDA00029086207400000611
Figure FDA00029086207400000612
和
Figure FDA00029086207400000613
5.根据权利要求4所述的基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法，其特征在于，所述步骤(5)中在同时考虑三维信号的情况下脑电信号yi对yj的因果影响为：
Figure FDA00029086207400000614
其中
Figure FDA00029086207400000615
CN201810554759.2A 2018-06-01 2018-06-01 基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法 Active CN109124623B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810554759.2A CN109124623B (zh) 2018-06-01 2018-06-01 基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810554759.2A CN109124623B (zh) 2018-06-01 2018-06-01 基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法

Publications (2)

Publication Number Publication Date
CN109124623A CN109124623A (zh) 2019-01-04
CN109124623B true CN109124623B (zh) 2021-03-19

Family

ID=64801933

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810554759.2A Active CN109124623B (zh) 2018-06-01 2018-06-01 基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法

Country Status (1)

Country Link
CN (1) CN109124623B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109893126A (zh) * 2019-03-21 2019-06-18 杭州电子科技大学 基于脑功能网络特征的癫痫发作预测方法
CN111973180B (zh) * 2020-09-03 2021-09-17 北京航空航天大学 一种基于meg和eeg融合的脑结构成像系统和方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104305993A (zh) * 2014-11-12 2015-01-28 中国医学科学院生物医学工程研究所 一种基于格兰杰因果性的脑电源定位方法
CN104473636A (zh) * 2014-12-30 2015-04-01 天津大学 一种基于偏定向相干的疲劳脑网络分析方法
CN104720796A (zh) * 2015-02-12 2015-06-24 西安交通大学 一种用于癫痫发作时间段的自动检测系统及方法
CN105844111A (zh) * 2016-04-07 2016-08-10 杭州电子科技大学 一种基于脑效应网络的新型脑电卒中评估方法
CN107126193A (zh) * 2017-04-20 2017-09-05 杭州电子科技大学 基于滞后阶数自适应选择的多变量因果关系分析方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104305993A (zh) * 2014-11-12 2015-01-28 中国医学科学院生物医学工程研究所 一种基于格兰杰因果性的脑电源定位方法
CN104473636A (zh) * 2014-12-30 2015-04-01 天津大学 一种基于偏定向相干的疲劳脑网络分析方法
CN104720796A (zh) * 2015-02-12 2015-06-24 西安交通大学 一种用于癫痫发作时间段的自动检测系统及方法
CN105844111A (zh) * 2016-04-07 2016-08-10 杭州电子科技大学 一种基于脑效应网络的新型脑电卒中评估方法
CN107126193A (zh) * 2017-04-20 2017-09-05 杭州电子科技大学 基于滞后阶数自适应选择的多变量因果关系分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
A nonlinear causality measure in the frequency domain: Nonlinear partial directed coherence with applications to EEG;Fei He et al.;《Journal of Neuroscience Methods》;20141231;第71-80页 *

Also Published As

Publication number Publication date
CN109124623A (zh) 2019-01-04

Similar Documents

Publication Publication Date Title
Nazari et al. Successive variational mode decomposition
Mohebbian et al. Single channel high noise level ECG deconvolution using optimized blind adaptive filtering and fixed-point convolution kernel compensation
Kang et al. A method of denoising multi-channel EEG signals fast based on PCA and DEBSS algorithm
CN109124623B (zh) 基于三维非线性偏直接相干函数的脑电信号间效应连通性检测方法
Keshavamurthy et al. Review paper on denoising of ECG signal
Meng et al. Gaussian mixture models of ECoG signal features for improved detection of epileptic seizures
Krak et al. Electrocardiogram classification using wavelet transformations
Luengo et al. Blind analysis of atrial fibrillation electrograms: a sparsity-aware formulation
Lerga et al. Number of EEG signal components estimated using the short-term Rényi entropy
Tsui et al. Modified maternal ECG cancellation for portable fetal heart rate monitor
Ghosh et al. Boosting Algorithms based Cuff-less Blood Pressure Estimation from Clinically Relevant ECG and PPG Morphological Features
Martinez et al. Methods to evaluate the performance of fetal electrocardiogram extraction algorithms
Logesparan et al. Assessing the impact of signal normalization: Preliminary results on epileptic seizure detection
Chen et al. A fast ECG diagnosis using frequency-based compressive neural network
Khandait et al. ECG signal processing using classifier to analyses cardiovascular disease
Pham et al. Unsupervised discrimination of motor unit action potentials using spectrograms
Grilo Jr et al. Artifact removal for emotion recognition using mutual information and Epanechnikov kernel
Raghav et al. Fractal feature based ECG arrhythmia classification
Daqrouq et al. Arrhythmia detection using wavelet transform
Bansod et al. A new approach for removal of baseline wande in ECG signal using Empirical Mode Decomposition & Hurst Exponent
Jayasudha Comparision of preprocess techniques for brain image using machine learning
Zhou et al. ECG data enhancement method using generate adversarial networks based on Bi-LSTM and CBAM
Ciocoiu Single channel fetal ECG recovery using sparse redundant representations
Thuy-Duong et al. Separation of nonstationary EEG epileptic seizures using time-frequency-based blind signal processing techniques
Yang et al. Two Simple Multi-lead ECG Quality Assessment Algorithms

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant