本发明公开了一种多尺度时间序列生理信号的因果网络分析方法，包括以下步骤：S1.利用NA‑MEMD算法对待分析的生理信号u1,u2,...,um进行分解；S2.对于两个不同的生理信号ui,uj，进行因果分析，i＝1,2,...m,j＝1,2...,m,且i≠j，得到两个信号之间的因果关系；S3.对于u1,u2,...,um中任意两个信号，重复执行步骤S2，直至得到u1,u2,...,um中各个信号两两之间的因果关系，形成因果网络。本发明能够有效分析出生理信号的因果网络，为生理信号的应用提供了便利。
1.一种多尺度时间序列生理信号的因果网络分析方法，其特征在于：包括以下步骤：S1.输入待分析的生理信号：
u1＝{u1,1,u1,2,...,u1,t}u2＝{u2,1,u2,2,...,u2,t}......
um＝{um,1,um,2,...,um,t}利用NA‑MEMD算法对待分析的生理信号u1,u2,...,um进行分解：......
其中， 表示利用NA‑MEMD算法对信号进行分解，NA‑MEMD算法是指基于噪声辅助的+ +
多元经验模式分解算法，m表示输入的生理信号数目，并且满足m≥2，t∈N ，N 表示正整数；
表示NA‑MEMD算法选择的辅助噪声， 互不相关且均为随机高斯噪声，为选择的辅助噪声数目，n表示每个输入的生理信号分解后得到的IMF信号个数；
S2.对于两个不同的生理信号ui,uj，进行因果分析，i＝1,2,...m,j＝1,2...,m,且i≠j，分析过程如下：
S201.将生理信号ui分解得到的信号{IMFi,1,IMFi,2,...,IMFi,n}与生理信号uj分解得到的信号{IMFj,1,IMFj,2,...,IMFj,n}进行配对，得到n个IMF对：(IMFi,1,IMFj,1),(IMFi,2,IMFj,2),....,(IMFi,n,IMFj,n)；
S202.计算每一个IMF对的平均瞬时相位差，并将计算得到的平均瞬时相位差与预设的阈值进行比较，筛选出平均瞬时相位差小于预设阈值的IMF对，生成ICC集：其中， 中k1的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第k1个信号， 中k1的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第k1个信号；
中k2的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第k2个信号， 中k2的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第k2个信号；
同理， 中 的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第 个信号， 中的表示 为{IMFj，1,IMFj，2,...,IMFj，n}中的第 个信号；
S203 .分别计算ICC集中每一个IMF对的相位相关性；
其中， T表示IMFi,k和IMFj,k的时间长度，φi,k(t)表示IMFi,k在时刻t的瞬时相位，φj,k(t)表示IMFj,k在时刻t的瞬时相位；
从ICC集中各序号对应的IMF对中，筛选出频率最大的IMF对：由于NA‑MEMD算法分解出来的IMF信号频率是由大到小排列的，故筛选出频率最大的IMF对为利用原始信号uj减去 得到u′j，利用u′j替换掉输入信号组u1,u2,...,um中uj，然后对替换后的输入信号组进行NA‑MEMD分解；
分解后取出u′j对应的分解信号{IMF′j,1,IMF′j,2,...,IMF′j,n}利用原始信号ui减去 得到u′i，利用u′i替换掉输入信号组u1,u2,...,um中ui，然后对替换后的输入信号组进行NA‑MEMD分解；
分解后取出u′i对应的分解信号{IMF′i,1,IMF′i,2,...,IMF′i,n}；
S205.计算ui对uj的因果 以及uj对ui的因果式中， 是ui分解得到的第k个IMF的方差， 表示uj分解得到的第k个IMF的方差；Wk为中间变量；
S206.根据绝对因果强度ACS，计算如下比值：若该比值大于1，则ui为因，uj为果；
若该比值等于1，则ui、uj无因果关系或者互为因果；
S3.对于u1,u2,...,um中任意两个信号，重复执行步骤S2，直至得到u1,u2,...,um中各个信号两两之间的因果关系，形成因果网络。
2.根据权利要求1所述的一种多尺度时间序列生理信号的因果网络分析方法，其特征在于：所述步骤S202包括：
S2021.设定各个IMF对的平均瞬时相位差阈值δ1,δ2,...δn；
S2022.计算第h个IMF对(IMFi,h,IMFj,h)的平均瞬时相位差：设mean(φi,h)表示IMFi,h的平均瞬时相位，即IMFi,h在时间长度内的瞬时相位平均值；
mean(φj,h)表示IMFj,h的平均瞬时相位，即IMFj,h时间长度内的瞬时相位平均值；
则第h个IMF对(IMFi,h,IMFj,h)的平均瞬时相位差为：|mean(φi,h)‑mean(φj,h)|；
将|mean(φi,h)‑mean(φj,h)|与对应阈值δh进行比较，判断是否满足：|mean(φi,h)‑mean(φj,h)|＜δh；
如果满足，则将对应的(IMFi,h,IMFj,h)加入ICC集中；
S2023.在h＝1,2,...n时，分别重复执行步骤S2022，最终得到ICC集为：
一种多尺度时间序列生理信号的因果网络分析方法
技术领域
[0001] 本发明涉及生理信号处理，特别是涉及一种多尺度时间序列生理信号的因果网络分析方法。
背景技术
[0002] 多尺度时间序列生理信号的因果网络分析方法可在脑机接口技术、脑功能和结构机制研究、脑器交互技术(大脑与眼、心、肺、肌等人体器官构成的生理网络)中得以应用和
推广。
[0003] 基于Gall提出的神经科学功能分离和功能整合表明，不同的大脑功能在专门的皮质区域是存在局部局限性的，因此通过功能连接(FC)来评估大脑感知、认知和行为的科学
标准一直由格兰杰因果关系分析主导。FC被用作为一种视觉工具，以显示大脑区域之间空
间距离的时间相关性。
[0004] 在大多数真实情况下，因果关系可能是依赖时间的、同时的和相互的。然而，现有的时间序列因果关系分析方法大多基于预测的概念，可能无法描述瞬时事件之间的相互因
果关系运动。在这种情况下，格兰杰因果关系是基于以下假设：因果关系是可分离的，这在
许多线性随机系统中很有用，但可能不适用于复杂的动态过程(例如与大脑相关的生理网
络)。此外，还开发了收敛交叉映射(CCM)方法，以适应因果效应的不可分离性。在大脑网络
的因果建模中，动态因果建模和转移熵也很普遍，但基本上基于贝叶斯预测。杨等人提出了
一种通过希尔伯特‑黄变换(HHT)评估因果相互作用的方法。然而，它的理论框架不能直接
在大脑相关的复杂信息系统上进行推广应用，而且在先验因果相互作用方面也缺乏充分必
要的实践。
[0005] 协变量方法所依据的是伽利略和休谟的因果归纳，由于因果关系以某些方式从感官输入中编码：即，候选因果与结果之间的共变量可以定义为在候选原因存在时效果概率
与在候选原因不存在的情况下，该概率之间的差值。Kantian的因果强度表示解释了因果信
息的先验知识的存在。与通用的基于HHT的因果分解方法相比，经验模式分解(EMD)不能保
证跨多个时间序列的IMF数相等，并导致模式(比例)对齐(跨通道相互依赖)，从而导致多元
时间序列中时间尺度的延伸。例如，在大脑‑生理网络中，从信号采集设备收集到的振荡分
量通常以不同的频段采样。另外，EMD方法还容易受到导致模式混合问题(单通道独立性)的
噪声影响，导致不足以识别分解的固有因果分量(ICCs)。
[0006] 到目前为止，多时间尺度的因果关系分析已基本应用于时间序列的ICC评估。多尺度格兰杰因果关系通过过滤和下采样步骤来解决，以获得双变量过程的再缩放表示。基于
转移熵的多尺度分析主要基于傅里叶和小波变换得到的固定的基函数。但是，当应用于自
然继承的非线性和非固定属性的复杂过程时，这些分解过程不会捕获不同时间尺度上的重
要特征。
发明内容
[0007] 本发明的目的在于克服现有技术的不足，提供一种多尺度时间序列生理信号的因果网络分析方法，能够有效分析出生理信号的因果网络，为生理信号的应用提供了便利。
[0008] 本发明的目的是通过以下技术方案来实现的：一种多尺度时间序列生理信号的因果网络分析方法，包括以下步骤：
[0009] S1.输入待分析的生理信号：
[0010] u1＝{u1,1,u1,2,...,u1,t}
[0011] u2＝{u2,1,u2,2,...,u2,t}
[0012] ......
[0013] um＝{um,1,um,2,...,um,t}
[0014] 利用NA‑MEMD算法对待分析的生理信号u1,u2,...,um进行分解：
[0015]
[0016] 其中， 表示利用NA‑MEMD算法对信号进行分解，NA‑MEMD算法是指基于噪声辅+ +
助的多元经验模式分解算法，m表示输入的生理信号数目，并且满足m≥2，t∈N ，N 表示正整
数； 表示NA‑MEMD算法选择的辅助噪声， 互不相关且均为随机高斯噪
声， 选择的辅助噪声数目，n表示每个输入的生理信号分解后得到的IMF信号个数；
[0017] S2.对于两个不同的生理信号ui,uj，进行因果分析，i＝1,2,...m,j＝1,2...,m,且i≠j，分析过程如下：
[0018] S201.将生理信号ui分解得到的信号{IMFi,1,IMFi,2,...,IMFi,n}与生理信号uj分解得到的信号{IMFj,1,IMFj,2,...,IMFj,n}进行配对，得到n个IMF对：
[0019] (IMFi,1,IMFj,1),(IMFi,2,IMFj,2),....,(IMFi,n,IMFj,n)；
[0020] 每个IMF对中，两个IMF信号的时间长度相等；
[0021] S202.计算每一个IMF对的平均瞬时相位差，并将计算得到的平均瞬时相位差与预设的阈值进行比较，筛选出平均瞬时相位差小于预设阈值的IMF对，生成ICC集：
[0022]
[0023] 其中， 中k1的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第k1个信号，中k1的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第k1个信号；
[0024] 中k2的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第k2个信号， 中k2的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第k2个信号；
[0025] 同理， 中 的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第 个信号，中 的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第 个信号；
[0026] 表示ICC集中的IMF对个数；
[0027] S204.分别计算ICC集中每一个IMF对的相位相关性；
[0028]
[0029] 其中， T表示IMFi,k和IMFj,k的时间长度，φi,k(t)表示IMFi,k在时刻t的瞬时相位，φj,k(t)表示IMFj,k在时刻t的瞬时相位；
[0030] S204.信号重分解：
[0031] 从ICC集中各序号对应的IMF对中，筛选出频率最大的IMF对：由于NA‑MEMD算法分解出来的IMF信号频率是由大到小排列的，故筛选出频率最大的IMF对为
[0032] 利用原始信号uj减去 得到u′j，利用u′j替换掉输入信号组u1,u2,...,um中uj，然后对替换后的输入信号组进行NA‑MEMD分解；
[0033] 分解后取出u′j对应的分解信号{IMF′j,1,IMF′j,2,...,IMF′j,n}
[0034] 利用原始信号ui减去 得到u′i，利用u′i替换掉输入信号组u1,u2,...,um中ui，然后对替换后的输入信号组进行NA‑MEMD分解；
[0035] 分解后取出u′i对应的分解信号{IMF′i,1,IMF′i,2,...,IMF′i,n}；
[0036] S205.计算ui对uj的因果 以及uj对ui的因果
[0037]
[0038] 式中， 是ui分解得到的第k个IMF的方差， 表示uj分解得到的第k个IMF的方差；Wk为中间变量；
[0039] 获取绝对因果强度ACS：
[0040]
[0041] S206.根据绝对因果强度ACS，计算如下比值：
[0042]
[0043] 若该比值大于1，则ui为因，uj为果；
[0044] 若该比值小于1，则ui为果，uj为因；
[0045] 若该比值等于1，则ui、uj无因果关系或者互为因果；
[0046] 由此得到ui、uj的因果分析结果；
[0047] S3.对于u1,u2,...,um中任意两个信号，重复执行步骤S2，直至得到u1,u2,...,um中各个信号两两之间的因果关系，形成因果网络。
[0048] 其中，所述步骤S202包括：
[0049] S2021.设定各个IMF对的平均瞬时相位差阈值δ1,δ2,...δn；
[0050] S2022.计算第h个IMF对(IMFi,h,IMFj,h)的平均瞬时相位差：
[0051] 设mean(φi,h)表示IMFi,h的平均瞬时相位，即IMFi,h在时间长度内的瞬时相位平均值；mean(φj,h)表示IMFj,h的平均瞬时相位，即IMFj,h时间长度内的瞬时相位平均值；
[0052] 则第h个IMF对(IMFi,h,IMFj,h)的平均瞬时相位差为：
[0053] |mean(φi,h)‑mean(φj,h)|；
[0054] 将|mean(φi,h)‑mean(φj,h)|与对应阈值δh进行比较，判断是否满足：
[0055] |mean(φi,h)‑mean(φj,h)|＜δh；
[0056] 如果满足，则将对应的(IMFi,h,IMFj,h)加入ICC集中；
[0057] 若不满足，则抛弃(IMFi,h,IMFj,h)；
[0058] S2023.在h＝1,2,...n时，分别重复执行步骤S2022，最终得到ICC集为：
[0059]
[0060] 本发明的有益效果是：本发明提供了一种多尺度时间序列生理信号的因果网络分析方法，有效分析出生理信号的因果网络，相较于传统方法，本申请的分析结果能够更加完
整地体现对时间序列的因果关系，为生理信号在脑机接口技术、脑功能和结构机制研究、脑
器交互技术上的应用提供了技术条件。
附图说明
[0061] 图1为本发明的方法流程图。
具体实施方式
[0062] 下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。
[0063] 如图1所示，一种多尺度时间序列生理信号的因果网络分析方法，包括以下步骤：
[0064] S1.输入待分析的生理信号：
[0065] u1＝{u1,1,u1,2,...,u1,t}
[0066] u2＝{u2,1,u2,2,...,u2,t}
[0067] ......
[0068] um＝{um,1,um,2,...,um,t}
[0069] 利用NA‑MEMD算法对待分析的生理信号u1,u2,...,um进行分解：
[0070]
[0071] 其中， 表示利用NA‑MEMD算法对信号进行分解，NA‑MEMD算法是指基于噪声辅+ +
助的多元经验模式分解算法，m表示输入的生理信号数目，并且满足m≥2，t∈N ，N 表示正整
数； 表示NA‑MEMD算法选择的辅助噪声， 互不相关且均为随机高斯噪
声， 为选择的辅助噪声数目，n表示每个输入的生理信号分解后得到的IMF信号个数；
[0072] S2.对于两个不同的生理信号ui,uj，进行因果分析，i＝1,2,...m,j＝1,2...,m,且i≠j，分析过程如下：
[0073] S201.将生理信号ui分解得到的信号{IMFi,1,IMFi,2,...,IMFi,n}与生理信号uj分解得到的信号{IMFj,1,IMFj,2,...,IMFj,n}进行配对，得到n个IMF对：
[0074] (IMFi,1,IMFj,1),(IMFi,2,IMFj,2),....,(IMFi,n,IMFj,n)；
[0075] 每个IMF对中，两个IMF信号的时间长度相等；
[0076] S202.计算每一个IMF对的平均瞬时相位差，并将计算得到的平均瞬时相位差与预设的阈值进行比较，筛选出平均瞬时相位差小于预设阈值的IMF对，生成ICC集：
[0077]
[0078] 其中， 中k1的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第k1个信号，中k1的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第k1个信号；
[0079] 中k2的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第k2个信号， 中k2的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第k2个信号；
[0080] 同理， 中 的表示 为{IMFi,1,IMFi,2,...,IMFi,n}中的第 个信号，中 的表示 为{IMFj,1,IMFj,2,...,IMFj,n}中的第 个信号；
[0081] 表示ICC集中的IMF对个数；
[0082] S204.分别计算ICC集中每一个IMF对的相位相关性；
[0083]
[0084] 其中， T表示IMFi,k和IMFj,k的时间长度，φi,k(t)表示IMFi,k在时刻t的瞬时相位，φj,k(t)表示IMFj,k在时刻t的瞬时相位；
[0085] S204.信号重分解：
[0086] 从ICC集中各序号对应的IMF对中，筛选出频率最大的IMF对：由于NA‑MEMD算法分解出来的IMF信号频率是由大到小排列的，故筛选出频率最大的IMF对为
[0087] 利用原始信号uj减去 得到u′j，利用u′j替换掉输入信号组u1,u2,...,um中uj，然后对替换后的输入信号组进行NA‑MEMD分解；
[0088] 分解后取出u′j对应的分解信号{IMF′j,1,IMF′j,2,...,IMF′j,n}
[0089] 利用原始信号ui减去 得到u′i，利用u′i替换掉输入信号组u1,u2,...,um中ui，然后对替换后的输入信号组进行NA‑MEMD分解；
[0090] 分解后取出ui′对应的分解信号{IMF′i,1,IMF′i,2,...,IMF′i,n}；
[0091] S205.计算ui对uj的因果 以及uj对ui的因果
[0092]
[0093] 式中， 是ui分解得到的第k个IMF的方差， 表示uj分解得到的第k个IMF的方差；Wk为中间变量；
[0094] 获取绝对因果强度ACS：
[0095]
[0096] S206.根据绝对因果强度ACS，计算如下比值：
[0097]
[0098] 若该比值大于1，则ui为因，uj为果；
[0099] 若该比值小于1，则ui为果，uj为因；
[0100] 若该比值等于1，则ui、uj无因果关系或者互为因果；
[0101] 由此得到ui、uj的因果分析结果；
[0102] S3.对于u1,u2,...,um中任意两个信号，重复执行步骤S2，直至得到u1,u2,...,um中各个信号两两之间的因果关系，形成因果网络。
[0103] 其中，所述步骤S202包括：
[0104] S2021.设定各个IMF对的平均瞬时相位差阈值δ1,δ2,...δn；
[0105] S2022.计算第h个IMF对(IMFi,h,IMFj,h)的平均瞬时相位差：
[0106] 设mean(φi,h)表示IMFi,h的平均瞬时相位，即IMFi,h在时间长度内的瞬时相位平均值；mean(φj,h)表示IMFj,h的平均瞬时相位，即IMFj,h时间长度内的瞬时相位平均值；
[0107] 则第h个IMF对(IMFi,h,IMFj,h)的平均瞬时相位差为：
[0108] |mean(φi,h)‑mean(φj,h)|；
[0109] 将|mean(φi,h)‑mean(φj,h)|与对应阈值δh进行比较，判断是否满足：
[0110] |mean(φi,h)‑mean(φj,h)|＜δh；
[0111] 如果满足，则将对应的(IMFi,h,IMFj,h)加入ICC集中；
[0112] 若不满足，则抛弃(IMFi,h,IMFj,h)；
[0113] S2023.在h＝1,2,...n时，分别重复执行步骤S2022，最终得到ICC集为：
[0114]
[0115] 在本申请的实施例中，两个生理信号A、B之间的因果关系定义为：如果A和B中分解的IMF处于相同的时间尺度，当B的IMF从B本身删除时，若A和B中的IMF之间的瞬时相位依赖
关系消失，则变量A会导致变量B，即A是B的因(A为因，B为果)；若当B的IMF从B本身删除时，A
和B中的IMF之间的瞬时相位依赖关系没有消失，则变量A不会导致B，即A不是B的因。
[0116] 因果关系需要满足如下条件：
[0117] (1)任何因果关系都基于ICC在多个时间序列中的瞬时相位相关性；
[0118] (2)基于相位的因果本身是能够分离的。
[0119] 以上虽然描述了本发明的具体实施方法，但是本领域的技术人员应当理解，这些仅是举例说明，在不背离本发明原理和实现的前提下，可以对这些实施方案做出多种变更
或修改，因此，本发明的保护范围由所附权利要求书限定。
