一种非线性数据全局互相关性的并行定量计算方法
【技术领域】
[0001 ]本发明属于数据分析技术领域，涉及一种数据相关性分析方法，尤其涉及一种同 时揭示相关强度与相关方向的、全局的、非线性数据全局互相关性的并行定量计算方法。
【背景技术】
[0002] 科学与实际工程应用中，复杂系统的监测数据(如大脑神经信号)具有明显的非线 性特性，复杂系统多个组成部分之间同步现象的研究是当今复杂系统科学领域中最活跃的 课题之一。以神经信号为例，同步测量是对大脑活动时产生的两个或更多连续时间序列的 同步估计，独立的时间序列同步值低，相关的时间序列同步值高。利用多通道脑电信号分析 方法研究神经信号的同步特征有助于深入理解生理病理情况下大脑各区域间的协作和调 制机制，为脑部神经疾病的诊断提供依，从而进行正确的治疗，通过研究理解导致脑部神经 疾病的机理，也有助于对疾病的提前预防。（[文献1 ])
[0003] 现有的一些二元信号同步性分析方法，比如交叉相关性（[文献1])，频谱相干性 ([文献2])，相位同步性（[文献3])和非线性互相关性(Non-linear Interdependency，NLI) ([文献4])等方法。非线性互相关性(NLI)是一种非对称性的相关性测量方法，用来描述信 号的相位同步强度，可用于衡量广义同步。NLI算法不进能够计算成对二元信号的强度，还 能计算二元非线性数据同步的方向（[文献5])，并且具有很好的抗噪声干扰能力。
[0004]而多通道信号比二元信号所包含的信息更多。为了进一步分析多通道信号之间的 同步性关系，（如，为了探索不同大脑区域的工作机理我们需要分析其不同通道的神经信号 的相关性），多通道信号分析需要一个计算全局相关性的方法。因此，相比较于二元信号相 关性分析，由于多通道信号的普适性，多通道信号分析收到了越来越多的重视。现有的一些 多通道信号互相关性分析方法，S估计器（[文献6])，关联性矩阵分析（[文献7])，互信息，局 部定向相关性（[文献8])。文献9对这些方法做了一个比较，非线性互相关性方法(NLI)对噪 声的鲁棒性明显强于其它的方法。因此，相比较而言，NLI算法更适合于多通道数据全局互 相关性的定量计算。
[0005] 然而，在多通道信号分析上应用NLI算法是高度计算密集型的，NLI算法主要应用 于研究两通道间信号关系，虽然该算法在处理脑电信号上有其独特的优势，但当数据量较 大或通道数较多时，该算法的执行效率明显下降，从而制约了该算法的发展和使用。因此亟 待发明一种快速有效的多通道非线性数据相关性分析算法。
[0006] [文献l]C.Carmeli，M.G.Knyazeva，G.M. Innocenti，and O.D.Feo, "Assessment of EEG synchronization based on state-space analysis，''Neuroimage，vol?25，no?2， pp.339-354,2005.
[0007] [文南犬2]A.Kraskov，"Synchronization and interdependence measures and theirapp1icat ions to the electroencephalogram of epilepsy patients and clustering of data,NIC-Directors Ph.D. ,Jillich,Germany,2004.
[0008] [文献 3 ] F ? Mormann，K ? Lehner t z，P ? Dav i d，and C ? E ? E 1 ger，"Mean phasecoherence as a measure for phase synchronization and its applicationto the EEG of epilepsy patients，"Physica D，vol.l44，no.3-4，pp.358_369,2000
[0009] [文献4]M.Le Van Quyen，J.Soss，V.Navarro，R.Robertson，M.Chavez，M.Baulac， and J.Martinerie，"Preictal state identification by synchronization changes in longterm intracranial EEG recordings，'，Clin.Neurophysiol ?，vol ? 116，no.3，pp.559_ 568,2005.
[0010] [文献5] M.Breakspear，J.R. Terry，K.J.Fris ton，A. W.F.Harr is，L.M. Williams， K.Brown?J.Brennan?and E.Gordon?UA disturbance ofnonlinear interdependence in scalp EEG of subjects with first episodeschizophrenia，'，NeuroImage，vol?20， no.1，pp.466-478,2003.
[0011] [文献6]M.G.Knyazeva，G.M. Innocenti，C? Carmeli，and 0? D.Feo，"Assessment of EEG synchronization based on state-space analysis，'，Neuroimage，vol ? 25，no ? 2， pp.339-354,2005.
[0012] [文献7])(.1^;[，0.。11；[，？.<1；[1'1181^，<1.￡{01，}(.￥&0，&11(1 <1.6.1?.知亡亡6^8， "Synchronization measurement of multiple neuronal populations，'， J.Neurophysiol.?vol.98?pp.3341-3348?2007.
[0013] [文南犬8]K.Sameshima and L.A.Baccala，"Using partial directed coherence todescribe neuronal ensemble interactions，'，J.Neurosci.Methods，vol.94，pp.93-103,1999.
[0014] [文南犬 9]T.Kreuz，F. Mormann，R ? Andrze jak，A ? Kraskov，K ? Lehnertz，and P?Grassberger，"Measuring synchronization in coupled model systems:A comparison of different approaches/^hysica D?vol.225?no.1?pp.29-42?2007.
【发明内容】
[0015] 为了解决上述技术问题，本发明提供了一种同时掲示相关强度与相关方向的、全 局的、非线性数据全局互相关性的并行定量计算方法。
[0016] 本发明所采用的技术方案是：一种非线性数据全局互相关性的并行定量计算方 法，其特征在于，包括以下步骤：
[0017]步骤1:对M通道原始脑电数据shuru进行直接分解，M>2,以两两通道数据最为原 始数据进入下层循环;设定初始通道i = l，j = l;
[0018]步骤2:取第i通道原始数据shurui ;
[0019]步骤3:取第j通道原始数据shuruj ;
[0020] 步骤4:设置数据时间窗Epoch和滑动距离Over lap;
[0021] 步骤5:设定初始循环次数为t = l，计算循环次数N:
[0022] N=length(shuru-Epoch)/(Epoch-Overlap)
[0023 ]步骤4:按照时间窗从数据输入shuru i和shuru j里取数据，计算当前窗口数据偏移 量(1：-1)*(￡?〇(：11-〇61']^?)，得到当前数据父1；(1〇和￥1；(1〇;
[0024] 步骤5:运用NLI算法计算两个数据的相关性，得到独立性测量S、H、N、M;
[0025] 步骤6:判断t>N是否成立；
[0026]若是，则顺序执行下述步骤6;
[0027]若否，则t = t+l，并回转执行所述步骤4;
[0028] 步骤7:判断j>M是否成立；
[0029]若是，则顺序执行下述步骤8;
[0030]若否，则j = j+l，并回转执行所述步骤3;
[0031] 步骤8:判断i>M是否成立；
[0032]若是，则顺序执行下述步骤9;
[0033]若否，则i = i+1，并回转执行所述步骤2;
[0034] 步骤9:对独立性测量S、H、N、M组成的相关性矩阵进行特征值特征向量分解；
[0035] 对数据矩阵Xmxt= {11(1〇}，1 = 1，~，]/[，1^=1，~，1'的协方差矩阵(：进行特征值分 解，其中M表示通道数量，T表示时间窗内的数据点数，获得特征值h，并对其进行归一化，得 到归一化特征值：
[0037]式中tr(C)为协方差矩阵的迹；
[0038]步骤10:基于特征值特征向量计算S估计器，S估计器的定义为：
[0040] 作为优选，步骤5中所述NLI算法，其具体实现过程包括以下子步骤：
[0041] 步骤5.1:给定两个输入的信号数据X和Y，重构两个输入信号的延迟向量^ = (xn，"_，xn-(m-1)1)和yn= (yn,…，yn-(m-1)1)，其中n=l，…，N，m为嵌入维数，T为预定义的延迟 时间；N = T-(m-l)T，1 < n < N;所有延迟向量组成了新的矩阵X = (xi，X2,…，xn)，Y = (yi, y2,…，yN);
[0042] 步骤5.2:分别计算X，Y中两两延迟向量之间的距离distXX，distXY，即distxxn = |xn-xj | |，distxyn= | |yn-yj | |，其中 j = l，2,…，n;对计算结果 distXX=(distxxi， distxx2,…，distxxn)，distXY= (distxyi，distxy2，…，distxyn)按由小到大的顺序进行排 序，取距离最小的k个计算结果所对应的时间索引序列r，s;
[0043] 步骤5.3:令ry和sy分别表示与X4Pyn两个信号序列的最近的延迟向量的时间索 引，j = l，…，k，对于每个&，其最近的k个延迟向量的均方欧几里德距离定义为：
⑴，
[0045]同理，在时