一种多变量时序信号因果连接识别方法、系统及存储介质

本发明涉及多变量时间序列处理、数据驱动系统辨识和条件因果推断，尤其涉及一种多变量时序信号因果连接识别方法、系统及存储介质，在工业应用中用于生成工业过程变量之间的故障传播路径，在脑电应用中用于生成脑电通道或脑区之间的有向有效连接网络及脑电连接拓扑特征。

背景技术：

1、连续工业过程、能源系统和复杂机电装备通常布设流量、压力、液位、温度、电流和振动等传感器。上游扰动、设备故障或控制参数改变会沿工艺链传播，使多个测点近似同时发生变化。仅根据相关性或两变量预测改善建立因果边，容易把由中间变量传播产生的间接相关误判为直接作用，也可能在工况变化后继续输出已经失效的旧因果图。

2、线性多变量格兰杰方法通过比较全模型和删除候选变量后的约简模型残差判断方向性预测作用，但线性滞后模型难以充分拟合状态依赖的非线性耦合；残差方差比较本身也没有提供与候选原因变量逐时对应的多尺度信号对象。

3、现有多变量时序信号因果关系分析方法通常侧重于单一尺度的预测关系、静态或阶段性因果结构识别，或者依赖特定模型及预设数据条件进行连接分析。在存在非线性耦合、多尺度动态变化和结构切换的情况下，仍可能出现间接作用被识别为直接作用、双向关系难以确定方向，以及系统状态变化后继续使用已失效因果连接结果等问题。同时，部分现有方法缺少能够与候选变量作用逐时对应的多尺度分析对象以及对因果连接结果有效性的在线控制机制。

4、因此，需要建立一条不依赖特定工业变量或脑电通道名称、可重复实施的统一技术链，来完成对多变量时序信号之间的直接有向预测作用进行在线识别，从而在复杂非线性和多尺度条件下识别直接有向预测作用，抑制间接路径和方向未决关系，进而在系统结构发生变化时对因果连接结果的有效性进行在线判断和控制。

技术实现思路

1、本发明的目的在于克服现有技术的不足，提供一种多变量时序信号因果连接识别方法、系统及存储介质，利用非线性全模型与约简模型比较、条件预测残差差信号、集合经验模态分解、多尺度图判决和在线状态控制，对多变量时序信号之间的直接有向预测作用进行在线识别。

2、本发明的目的是通过以下技术方案来实现的：

3、处理器通过多通道信号采集设备100和数据采集与同步接口110获取m个同步观测变量。所述多通道信号采集设备可以是工业过程传感器组或脑电信号采集设备。离散时刻t的观测向量表示为：

4、

5、其中，m为变量数，为第i个变量在时刻t的数字样本。窗口管理模块120维护长度为的预测回归窗口和长度为的因果分析窗口。预测回归窗口严格截止于时刻t－1，用于预测尚未参与模型估计的当前样本x(t)；因果分析窗口包含当前样本x(t)，用于在候选图更新时刻重新估计候选边评价模型。

6、

7、

8、和均为预设正整数，可以相同或不同，二者之间不存在必须满足的大小关系；各窗口长度应不小于完成对应的滞后模型估计、残差计算和信号分解所需的最小样本数。样本不足固定长度时使用已到达样本形成扩展窗，达到固定长度后使用滑动窗。冻结实施例取、，但该组数值不构成对窗口长度的限定。

9、预处理包括按统一时钟进行时间戳对齐、标记缺失样本、对明显越界值或伪迹样本进行标记或裁剪，并分别在预测回归窗口和因果分析窗口内记录变量尺度。适用数据为具有统一采样时基、预设有限滞后阶并在相应窗口内具有足够有效变化的多变量时序数据。冻结实施方式不使用完整序列末段或未来样本确定收敛带。

10、具体地，在线状态预测模型与候选边评价模型如下：

11、对于任一有序变量对，将称为候选原因变量，将称为待解释结果变量。该称谓表示本轮检验中的变量角色；在完成后续图判决之前，并不预先认定二者已经构成因果关系。同一过程变量在不同有序变量对中可以分别作为候选原因变量或待解释结果变量。

12、给定最大滞后阶p和最高多项式阶q，定义不含截距的原始非线性滞后特征向量。该向量由各变量第1阶至第p阶滞后值、各变量滞后值的第2阶至第q阶幂以及同一滞后阶不同变量的成对乘积组成；各模型的截距均作为独立的常数项设置：

13、

14、其中，d为原始非截距特征数，q取1至m，取1至p，r取2至q，a和b满足1≤a

15、在线状态预测模型仅用于评价当前预测误差和控制在线状态。在每一采样时刻t，处理器利用严格截止于t－1的预测回归窗口构建设计矩阵。以该窗口内每一非截距特征的样本均值和样本标准差为标准化参数，将原始特征向量φ转换为预测窗口标准化特征向量，并独立估计包含全部变量特征的在线状态预测模型；尚未参与本次模型估计的当前样本的一步预测值为：

16、

17、其中，上标(on)表示online，即在线状态预测；和分别为在线状态预测模型针对结果变量独立估计的截距和非截距系数向量， 为使用预测回归窗口统计量标准化后的当前时刻特征向量。仅用于后述瞬时预测误差和滚动误差统计量，不作为条件预测残差差信号中的候选边评价全模型预测值；上标(on)与表示最大滞后阶的p相互独立。

18、如图3所示，候选边评价全模型和候选边评价约简模型仅在候选图更新或内部图检查时利用包含当前样本x(t)的因果分析窗口重新估计。令为从该窗口扣除形成前p阶滞后特征所需起始样本后得到的时间对齐索引集合；以该索引集合内每一非截距特征的样本均值和样本标准差为标准化参数，得到因果窗口标准化特征向量。待解释结果变量的候选边评价全模型 保留全部d个非截距特征，其窗口内拟合值和残差分别为：

19、

20、

21、其中，上标(f)表示full；和分别为候选边评价全模型的估计截距和d维非截距系数向量， 为同一因果分析窗口内索引τ处的拟合值， 为对应全模型残差。该模型用于生成候选边比较所需的窗口内残差，不用于预测窗口之外的当前样本。

22、针对候选原因变量，从因果窗口标准化特征向量中删除变量集合包含j的全部分量，得到，其中小于d。候选边评价约简模型使用同一因果分析窗口、同一时间对齐索引和删除后的设计矩阵独立重新估计，不沿用候选边评价全模型中任何保留特征的系数；其窗口内拟合值和残差分别为：

23、

24、

25、其中，上标(r,j)表示reduced-by-j，即针对候选原因变量形成的约简模型；和分别为该约简模型独立估计的截距和维非截距系数向量，为对应约简模型残差。全模型残差和约简模型残差均在同一上计算并逐时对齐；前者只需对每个结果变量估计一次，后者针对每一候选原因变量分别删除特征并重新估计。

26、在线状态预测模型、候选边评价全模型和候选边评价约简模型均采用截距不受惩罚的岭回归。对任一模型m，令，将全1截距列与该模型的标准化非截距特征矩阵合并为增广设计矩阵，将对应结果变量在模型有效索引上的样本组成目标向量，并将截距和非截距系数组成增广系数向量，则估计式统一写为：

27、

28、其中，若模型m具有个有效回归样本和个非截距特征，则；λ为非负岭参数，，其首个对角元素为0以使截距不受惩罚。在线模型m=on使用预测回归窗口且不包含x(t)；全模型m=f与约简模型m=(r,j)使用同一包含x(t)的因果分析窗口和同一时间对齐索引，区别仅在于约简模型删除相关特征并独立重估。

29、条件预测误差比指标与条件预测残差差信号

30、候选原因变量指向待解释结果变量的条件预测误差比指标定义为当前因果分析窗口经滞后对齐后，候选边评价约简模型残差方差与候选边评价全模型残差方差之比的对数：

31、

32、其中，为防止分母为零的正常数。是无量纲的残差方差对数比：其数值为正表示删除相关特征后约简模型残差方差增大，数值越大表示对的一步预测越重要。该指标直接由残差方差比计算，不是互信息、传递熵或其他香农信息量；本说明仅用于界定其数学性质，后文不再重复该区分。

33、为了将候选变量删除操作转换为可进行逐时多尺度分解的信号对象，构造条件预测残差差信号：

34、

35、条件预测残差差信号表示删除相关特征后，约简模型相对于全模型在时刻τ新增的残差部分。该信号具有与相同的物理量纲，可直接作为逐时多尺度分解对象。

36、双对象集合经验模态分解

37、如图4所示，将因果分析窗口经滞后对齐后的候选原因变量信号和条件预测残差差信号作为两个分别处理的分解对象。为保证可重复性并避免两个对象直接复用同一已生成噪声数组，处理器根据确定性基础种子、候选边以及信号类型标识1或2，分别生成第一伪随机噪声种子和第二伪随机噪声种子：

38、

39、其中，h为确定性散列映射，数字1表示候选原因变量信号，数字2表示条件预测残差差信号。冻结软件把基础种子、候选边标识和信号类型标识组合为散列输入，采用sha-256摘要的前4个字节构造32位无符号整数种子。相同的基础种子、候选边和信号类型产生相同的分解结果；两种信号类型分别初始化对应的集合经验模态分解过程并分别生成可复现的辅助噪声序列。

40、对候选原因变量信号和条件预测残差差信号分别执行集合经验模态分解；每个对象在每次集合试验中加入由其对应种子产生的白噪声，并对同阶分量求集合均值。两组本征模态函数分别按照分解输出阶次排列；若模态数量分别为和，则保留个同阶模态对。冻结实施方式仅按阶次配对，不执行中心频率或能量重心匹配。

41、多尺度相位相干、平均模态方差与方差加权

42、对第m个候选原因变量信号本征模态函数和第m个条件预测残差差信号本征模态函数进行hilbert变换，获得瞬时相位和。模态相位相干值为：

43、

44、其中，n为参与相干计算的窗口样本数，abs表示复数模函数，位于闭区间。为同时反映候选原因变量模态和残差差模态的幅值变化，首先计算两侧方差的算术平均：

45、

46、全部平均模态方差之和为  ，其中为模态求和索引。当该方差之和大于正常数时，对各平均模态方差进行归一化；当该方差之和不大于时，将全部权重和聚合分数置零，而不使用均匀权重：

47、

48、

49、当上述方差之和不大于时，令全部且。为候选边原始多尺度因果分数并位于闭区间。所述原始分数用于后续跨待解释结果变量的路径比较和方向比较。

50、七、双阈值联合筛选、间接路径批量筛除和双向边方向优势判定

51、对于同一待解释结果变量，将各候选原因变量的原始多尺度因果分数归一化为入边权重：

52、

53、其中，q遍历除i以外的全部候选原因变量，表示候选边的原始多尺度因果分数。当分母不大于时，令。归一化入边权重位于闭区间，且仅用于同一待解释结果变量内部的候选来源比较。当不小于入边阈值且不小于误差比阈值时，将保留至初始候选有向图；该处理称为双阈值联合筛选。冻结实施方式采用“不小于”比较，代码和公式均包含等号边界。

54、在双阈值联合筛选后得到的固定候选有向图上完成全部间接路径判定。令表示从j到i、长度为2至m-1、节点不重复并且不包含待比较直接边的全部间接有向路径的集合。对于其中任一路径，路径瓶颈分数定义为该路径全部候选边原始多尺度因果分数中的最小值：

55、

56、候选边的最大间接解释分数定义为全部合格间接路径的路径瓶颈分数中的最大值：

57、

58、若至少存在一条所述间接路径且，则把直接候选边加入待删除集合，其中为正的间接路径抑制系数。处理器在固定候选有向图上完成全部候选边和全部间接路径的判定后，一次性删除待删除集合中的全部候选边；该处理称为间接路径批量筛除，删除结果不依赖候选原因变量、待解释结果变量或路径的遍历顺序。对于三变量情形，所述间接路径退化为。

59、完成间接路径批量筛除后，对仍同时存在和的变量对执行双向边方向优势判定，其中为大于1的方向优势系数。若，仅保留；若，仅保留；否则将该变量对标记为方向未决并删除两个方向。由此得到有向邻接矩阵：

60、

61、其中，表示在当前有效状态下识别出候选原因变量对待解释结果变量的直接有向预测作用；表示未通过双阈值联合筛选、经间接路径批量筛除被删除、经双向边方向优势判定被标记为方向未决，或者不存在该直接边。所述直接有向预测作用是由本发明的条件预测和多尺度图判决规则确定的连接关系；在脑电应用中，将其称为有向有效连接时，不将该计算结果直接等同于解剖连接。

62、滚动平均误差稳定判据与变化候选图二取三表决

63、在每一采样时刻t，处理器使用严格截至t－1的预测回归窗口建立在线状态预测模型，获得尚未参与本次模型估计的当前样本的一步预测值。以所述预测回归窗口内各变量标准差和正的尺度下限归一化预测误差，瞬时多变量归一化误差为：

64、

65、长度为l的滚动平均误差为：

66、

67、令最近h个滚动平均误差构成集合：

68、

69、当且仅当以下三个条件同时满足时，认为滚动误差统计量满足滚动平均误差稳定判据并且已积累完整因果分析窗口：

70、

71、其中，为当前辨识起点，为因果分析窗口长度；在初始辨识期间，为初始辨识起点；在重新辨识期间，为由变化候选图二取三表决确定的变化证据时刻。第三个条件表示自起累计的样本数已达到，从而用于因果辨识的整个因果分析窗口均由当前辨识阶段的数据构成。冻结实施方式不从完整序列末段估计所谓最终稳定带。

72、如图5所示，系统首次使用三个状态时，将其定义为辨识状态（identify）、有效状态（valid）和重新辨识状态（reidentify）。在辨识状态或重新辨识状态下满足上述三个条件时，进入有效状态。

73、在有效状态（valid）下，处理器在进入该状态的时刻计算一次有向邻接矩阵，随后每隔u个样本执行一次内部候选图更新；相邻图更新时刻之间保持上一有效有向邻接矩阵。所述保持行为不等于用变化候选图替换有效图。

74、当连续d个时刻不小于误差报警阈值时，产生预测误差报警并立即促发一次内部候选图检查，但预测误差报警本身不直接宣布有效因果图失效。若内部候选图与当前有效图一致，则清除变化候选记录并保持有效状态。

75、当内部候选图与当前有效图不同时，记录该变化候选图。冻结实施方式连续保存最近三个变化候选图并执行变化候选图二取三表决：对于任一候选边，当该边在三个变化候选图中至少出现两次时，将其纳入表决图。若所述表决图仍与当前有效图不同，则状态进入重新辨识状态（reidentify），并将参与表决的最早变化候选图时刻记为变化证据时刻。周期性内部图检查和预测误差报警促发的内部检查采用同一表决判据。

76、重新辨识状态（reidentify）期间继续执行一步预测和模型辨识，但不把变化候选图作为有效输出。自变化证据时刻起累计的样本数达到且再次满足滚动平均误差稳定判据后，状态恢复为有效状态（valid），并将该时刻记为。若物理切换时刻在真值实验中已知，则检测延迟和恢复延迟分别为：

77、

78、在线真实部署中物理切换时刻通常未知，输出数据包因而记录、、状态和原因码；在具有已知切换时刻的验证实验中报告和。

79、九、系统和输出接口

80、如图2所示，多变量时序信号在线因果连接识别系统包括多通道信号采集设备100、数据采集与同步接口110、处理器、存储器和结果输出接口180。处理器执行窗口管理模块120、状态预测与全模型及约简模型辨识模块130、条件预测量计算模块140、集合经验模态分解与多尺度相干聚合模块150、图判决模块160以及收敛判定与状态控制模块170。

81、存储器保存算法版本、确定性基础种子、窗口长度、模型阶次、阈值、当前状态、上一有效有向邻接矩阵和状态转换日志。结果输出接口180输出含序号、时间戳、算法版本、状态枚举、有效标记和原因码的结果数据包。在辨识状态（identify）和重新辨识状态（reidentify）下，有效标记为0，结果输出接口180不提供变化候选图，并禁止相应的下游分析模块采用变化候选图。在有效状态（valid）下，有效标记为1，结果输出接口180在结果数据包中提供有效有向邻接矩阵。当多通道信号采集设备100为工业过程传感器组时，结果输出接口180将所述有效有向邻接矩阵和有效标记提供给工业故障传播监测模块；所述工业故障传播监测模块将矩阵中取值为1的元素所对应的有向连接确定为工业过程变量之间的直接故障传播边，并按照有向边的方向及首尾节点连接关系生成工业故障传播路径。当多通道信号采集设备100为脑电信号采集设备时，结果输出接口180将所述有效有向邻接矩阵和有效标记提供给脑电意图状态分析模块；所述脑电意图状态分析模块将所述有效有向邻接矩阵确定为脑电通道或脑区之间的有向有效连接网络，并根据该网络生成脑电连接拓扑特征，作为意图状态分析的输入数据。本发明还能够以存储于存储介质中的程序实现，处理器执行该程序时完成数据窗口更新、在线状态预测、候选边评价模型估计、条件预测量计算、双对象eemd、多尺度聚合、确定性图判决、状态控制和有效结果输出。

82、本发明的有益效果是：与仅使用线性残差方差比较的方案相比，本发明通过多项式滞后项和成对耦合项提高非线性一步预测能力，并通过独立重估的约简模型保证残差比较对应于候选变量删除操作。

83、与分别分解全模型残差和约简模型残差的方案相比，本发明直接分解约简模型残差减去全模型残差所形成的条件预测残差差信号，使多尺度分析对象聚焦于删除候选原因变量相关特征后新增的逐时残差部分。

84、通过根据确定性基础种子、候选边标识和信号类型标识，分别为候选原因变量信号和条件预测残差差信号初始化集合经验模态分解过程并生成可复现的辅助噪声序列，能够避免两个分解对象直接复用同一已生成噪声数组，同时保留实验复现能力。

85、通过在固定候选有向图上执行间接路径批量筛除，即完成节点不重复的间接路径判定、计算路径瓶颈分数并一次性删除待删除集合，以及执行双向边方向优势判定，图判决结果不依赖遍历顺序，并能够在多于三个变量时处理长度超过两跳的间接解释路径。

86、通过严格仅使用历史数据的滚动平均误差稳定判据、预测误差报警与变化候选图二取三表决的双证据、完整因果窗口约束和无效状态输出门控，能够在结构变化后暂停变化候选图，并在重新辨识完成后恢复有效有向图输出；同一核心处理链可在不改变候选边评价、双对象eemd和状态控制规则的情况下，分别处理工业过程变量和多通道脑电信号，并将仅在有效状态下输出的有向邻接矩阵分别转换为工业故障传播路径或者脑电连接拓扑特征，避免下游模块基于尚未稳定或者已经变化的候选图进行故障传播分析或者意图状态分析。