CN113704688A - 基于变分贝叶斯平行因子分解的缺失振动信号的恢复方法 - Google Patents

Abstract

本发明公开了一种基于变分贝叶斯平行因子分解的缺失振动信号的恢复方法，通过分割采样点将采集的时域振动信号构造成三维张量，将分解后的三维张量结合贝叶斯方法，引入似然模型，有效利用因子矩阵的先验信息，引入有效精度的后验分布，采用贝叶斯方法处理该模型，推断出包括因子矩阵和超参数在内的所有未知数的参数的后验分布，采用变分贝叶斯算法推导出因子矩阵和超参数的后验分布，从而进一步推断出缺失信号的分布预测。利用均方根误差对该方法的性能进行评估，变分贝叶斯平行因子分解算法相较于传统的低秩张量补全算法，误差更小，能够更加有效的恢复缺失的信号，有效地解决了振动信号分析中因传感器失效而引起的信号缺失的问题。

Description

基于变分贝叶斯平行因子分解的缺失振动信号的恢复方法

技术领域

本发明涉及信号处理技术，特别涉及一种基于变分贝叶斯平行因子分解的缺失振动信号恢复方法。

背景技术

机械设备中经常使用的部件在高温、重载和长时间在线使用等恶劣的工作环境下，很容易发生潜在故障。而设备故障会影响生产效率，导致严重的生产、经济损失甚至是人身伤害。因此，对产生故障的部件做出精准的分析，实现机械故障的准确诊断和识别具有很高的实用价值。目前在振动信号分析领域，各种振动信号分析方法因其动态信息分析能力而被提出。例如：小波变换(WT)[1]、经验模式分解(EMD)[2]、变分模态分解(VMD)[3]、同步提取变换[4]、包络分析、流形学习和短时傅里叶变换(STFT)[5]等等。屈海清等[6]等通过小波阈值去噪和平滑处理的对比分析，发现小波变换对机械信号的噪声消除具有良好的效果。魏永合等[7]采取了集合经验模态分解和支持向量机相结合进行滚动轴承的退化状态识别方法，建立退化状态识别模型，准确识别滚动轴承退化状态。岳应娟等[8]将Rihaczek复能量密度分布和变分模态分解两种方法有效结合，解决了内燃机中弱故障特征难以提取的问题，振动信号强耦合的问题，并且得到了时频聚集性较好的振动谱图像。李继猛等[9]提出了一种基于同步压缩-交叉小波变换的滚动轴承故障特征增强方法，可有效提取轴承在时频域内的细节特征,使轴承特征频率在时频域上的可读性增强，进而实现轴承故障的精确可靠诊断。杜伟等[10]提出基于独立特征选择(Individual FeatureSelection,IFS)与流形学习的故障诊断方法，将得到的低维特征输入多分类故障诊断模型进行识别，液压泵故障诊断实验表明，所提方法具备较高的诊断准确率。李恒等[11]将短时傅里叶变换和卷积神经网络相结合，通过轴承测试实验，表明新的方法对不同类型的故障有的更高的识别能力。

虽然以上信号处理方法可以有效地实现机械振动信号分析，但是众多学者的实验研究都是建立在采集到完整的振动数据的基础上[12]，机械设备的实际工作环境复杂多变，在实际测量和收集数据的过程中，由于不同的原因，无论是硬件故障，传感器故障和传输中断，都有可能会出现数据丢失(缺失)的情况[13]，并且这种情况是经常发生的，一旦信号丢失，上述各种振动信号处理方法将具有一定的局限性。如果不能够解决数据丢失的问题，在单次数据不方便重新采样的情况下，许多重要的数据将会被丢弃或者不正确的分析，对后续振动信号分析和故障诊都会造成很大的影响。

目前，缺失的振动信号恢复的方法几乎被忽略了，余路等[14]提出了一种基于过完备字典的缺失振动数据压缩感知重构算法，通过k奇异值分解算法得到过完备字典，解决了航空发动机的数据丢失问题。陈琳升等[15]将数据分组处理算法(GMDH)与经粒子群优化的支持向量机(PSO-SVM)算法相结合，用于轴承故障诊断领域，预测出缺失的数据。然而，这些振动信号缺失的恢复方法都需要大量的训练样本，而实际采集数据量少，这就会出现样本不足的情况，容易造成过拟合。马云飞等[16]将经验模态分解和贝叶斯压缩感知相结合，用来解决信号随机缺失问题，最后通过实验成功验证了所提方法的有效性。但是此方法在随机缺失的修复过程中，复杂度高，计算效率较低。

为了解决上述方法的局限性，同时考虑到张量分解算法在信号处理中的广泛应用[17-18]，面对数据丢失的情况，张量分解[19]也是数据恢复的一种方法，并且取得了很好的效果。文献[20]采用低秩张量补全方法(LRTC)对缺失的信号进行恢复，这种方法虽然很容易实现，但是在实施的过程中收敛速度很慢，每次迭代过程需要计算数个大规模的奇异值分解。

发明内容

为了解决工程中由于人为因素和自然界不可抗拒的因素，有时会造成传感器失效，从而造成信号采集的缺失的问题，本发明提供一种基于变分贝叶斯平行因子分解的缺失振动信号的恢复方法，推断出缺失元素的分布预测，解决振动分析领域中一直存在的信号缺失值，更好的解决信号丢失问题。

本发明采用以下技术方案实现上述目的。基于变分贝叶斯平行因子分解的缺失振动信号的恢复方法，其特征在于，具体步骤如下：

1)通过便捷式动态信号采集器采集机械零部件的时域振动信号S；

2)通过分割采样点将采集的时域振动信号S构造成三维张量Y，步骤如下；

(1)设置时域振动信号S包含K个采样点，然后将信号分成不重叠的数据段，数据段的数量为Q，则每一个数据段包含了P＝K/Q个数据点；且分割成了列数为Q的矩阵X，如下式：

将矩阵X进行转置，得到转置矩阵X^T表示为：

(2)将转置矩阵X^T沿水平进行分割切片构造成三维张量Y；

Y(:,P,:)＝S(:,(1+(P-1)*Q):(P*Q)) (3)

式中，P＝1、2、3…N，而数值Q不变；

3)将三维张量Y进行平行因子分解，其过程如下：

(1)三维张量Y是一个N阶不完全张量，丢失的数据大小为I₁×I₂×...×I_N，Y_Ω表示被观察的数据(i₁,i₂...i_N∈Ω)，其中Ω表示元素的指标集；定义张量O如下式子：

(2)假设三维张量Y为一个真实潜在张量X的噪声观测值，即Y＝X+ε，ε为高斯白噪声，三维张量

的平行因子分解为：

式中，

表示向量的外积，

为张量X对应的第n个模式因子矩阵，ε服从独立同分布的高斯分布，即：

由上式可知，对于第n个模式因子矩阵A⁽ⁿ⁾的表示形式如下式：

4)将分解后的三维张量Y结合贝叶斯方法，引入似然模型，表达式为：

式中：参数τ表示噪声精度，<…>表示N(N≥3)个向量的内积，

表示A⁽ⁿ⁾的第i_n行的行向量；

5)有效利用因子矩阵的先验信息，引入有效精度的后验分布：

式中：Λ^-1＝diag(λ)表示矩阵方差的逆且依赖所有模式的因子矩阵；假设超参数λ是独立的，则先验概率函数如下：

式中：Ga(x|a,b)表示Gamma分布；同时设置噪声精度τ的先验概率函数：

P(τ)＝Ga(τ|a₀,b₀) (10)

6)采用贝叶斯方法处理该模型，推断出包括因子矩阵和超参数在内的所有未知数的参数Θ的后验分布为：

推出缺失信号Y_/Ω的分布预测：

p(Y_/Ω|Y_Ω)＝∫p(Y_/Ω|Θ)p(Θ|Y_Ω)dΘ (12)

7)采用变分贝叶斯算法推导出因子矩阵和超参数的后验分布，从而进一步推断出缺失信号Y_/Ω的分布预测，其过程如下：

采用变分贝叶斯(VB)进行近似求解，通过极小化KL散度来寻找q(Θ)分布近似逼近真实的后验分布p(Θ|Y_Ω)，KL散度定义如下：

式中：lnp(Y_Ω)表示模型的证据因子为常数；

被定义为模型因子的下界，记为L(q)；由于lnp(Y_Ω)为常数，当KL散度为0时，即q(Θ)＝p(Θ|Y_Ω)，取到式(13)中模型因子的下界的最大值；采用平均场近似理论，假设变分分布被分解为：

式中：q(Θ)是总体分布，则各个因子q_j(Θ_j)的函数形式可以依次推导出来；基于下界L(q)的最大化给出q(Θ)的第j个因子的优化形式为：

式中：

表示除了Θ_j所有变量联合分布的期望，C为常数；所有变量的先验分布取自指数族且是共轭分布的，即得到未知数的后验分布；其中：

(1)因子矩阵A⁽ⁿ⁾的后验分布：

第n个模式因子矩阵A⁽ⁿ⁾的后验分布由观测数据信息和其他的第n个模式因子矩阵A^(k)(k≠n)的先验分布、超参τ的先验分布推出；由式(11)第n个模式因子矩阵A⁽ⁿ⁾(n＝1,…,N)的各行独立的

服从高斯分布；因此，对于任意n∈{n＝1,…N}的A⁽ⁿ⁾的后验分布可分解为：

其中，后验参数通过下式更新：

式中：Y_(n)代表张量数据Y的第n个模式因子矩阵，V⁽ⁿ⁾为辅助矩阵，E_q[·]代表包含所有变量的后验期望；式中V⁽ⁿ⁾首先利用因子先验Ε_q[Λ]和其它因子矩阵的协方差更新，而因子先验和其它因子矩阵的协方差的权重通过Ε_q[τ]来调节；然后

通过

旋转，并根据Ε_q[τ]进行尺度变换；

(2)超参数λ的后验分布：

超参数λ的后验分布可由N个因子矩阵的信息和超参数λ的先验信息推导出；通过式(15)可知，每个λ_r(r∈{1,…,R})独立的服从Gamma分布，即

其中R为张量的k秩，其中

表示从M个观测数据得到的后验参数，且其更新公式如下：

(3)噪声精度τ的后验分布：

噪声精度τ的后验分布可由观测数据和其本身的超先验信息推导出；结合式(15)可知，噪声精度τ的变分后验分布服从Gamma分布，即q_τ(τ)＝Ga(τ|a_M,b_M)，后验参数a_M、b_M的更新公式如下：

本发明将平行因子分解模型和变分贝叶斯算法相结合，针对振动信号采集可能造成的信号缺失问题，提出了一种基于变分贝叶斯平行因子分解缺失振动信号恢复算法。首先将一维的振动信号构造成三维张量，然后建立了变分贝叶斯平行因子分解(VBPF)的预测模型，采用变分贝叶斯(VB)推断信道参数以及因子矩阵的后验信息，更精准的恢复缺失信号，通过实验结果验证了提出的方法的有效性，解决了工程中无论是硬件故障，传感器故障和传输中断而引起的信号缺失的问题，有效地解决了振动信号分析中因传感器失效而引起的信号缺失的问题，具有重要的理论价值和工程应用价值，应用前景广阔。

附图说明

图1是贝叶斯平行因子分解的概率模型图；

图2是源信号时域波形图；

图3是源信号频谱图；

图4是不同缺失比例的均方根误差变化对比图；

图5a是缺失比例为10％的频谱图；

图5b是缺失比例为30％的频谱图；

图5c是缺失比例为50％的频谱图；

图6a是缺失比例为10％的恢复频谱图；

图6b是缺失比例为30％的恢复频谱图；

图6c是缺失比例为50％的恢复频谱图。

具体实施方式

以下结合附图和实施例对本发明作进一步说明。基于变分贝叶斯平行因子分解的缺失信号的恢复方法，步骤如下：

1)通过便捷式动态信号采集器采集机械零部件的时域振动信号S；

2)通过分割采样点将采集的时域振动信号S构造成三维张量Y，步骤如下；

(1)设置时域振动信号S包含K个采样点，然后将信号分成不重叠的数据段，数据段的数量为Q，则每一个数据段包含了P＝K/Q个数据点；且分割成了列数为Q的矩阵X，如下式：

将矩阵X进行转置，得到转置矩阵X^T表示为：

(2)将转置矩阵X^T沿水平进行分割切片构造成三维张量Y；

Y(:,P,:)＝S(:,(1+(P-1)*Q):(P*Q)) (3)

式中，P＝1、2、3…N，而数值Q不变；

3)将三维张量Y进行平行因子分解，其过程如下：

(1)三维张量Y是一个N阶不完全张量，丢失的数据大小为I₁×I₂×...×I_N，Y_Ω表示被观察的数据(i₁,i₂...i_N∈Ω)，其中Ω表示元素的指标集；定义张量O如下式子：

(2)假设三维张量Y为一个真实潜在张量X的噪声观测值，即Y＝X+ε，ε为高斯白噪声，三维张量

的平行因子分解为：

式中，

表示向量的外积，

为张量X对应的第n个模式因子矩阵，ε服从独立同分布的高斯分布，即：

由上式可知，对于第n个模式因子矩阵A⁽ⁿ⁾的表示形式如下式：

4)将分解后的三维张量Y结合贝叶斯方法，引入似然模型，表达式为：

式中：参数τ表示噪声精度，<…>表示N(N≥3)个向量的内积，

表示A⁽ⁿ⁾的第i_n行的行向量；

5)有效利用因子矩阵的先验信息，引入有效精度的后验分布：

式中：Λ^-1＝diag(λ)表示矩阵方差的逆且依赖所有模式的因子矩阵；假设超参数λ是独立的，则先验概率函数如下：

式中：Ga(x|a,b)表示Gamma分布；同时设置噪声精度τ的先验概率函数：

P(τ)＝Ga(τ|a₀,b₀) (10)

根据上述步骤1)到步骤5)，可以得到贝叶斯平行因子分解的概率图模型(如图1所示)。

6)采用贝叶斯方法处理该模型，推断出包括因子矩阵和超参数在内的所有未知数的参数Θ的后验分布为：

推出缺失信号Y_/Ω的分布预测：

p(Y_/Ω|Y_Ω)＝∫p(Y_/Ω|Θ)p(Θ|Y_Ω)dΘ (12)

7)采用变分贝叶斯算法推导出因子矩阵和超参数的后验分布，从而进一步推断出缺失信号Y_/Ω的分布预测，其过程如下：

采用变分贝叶斯(VB)进行近似求解，通过极小化KL散度来寻找q(Θ)分布近似逼近真实的后验分布p(Θ|Y_Ω)，KL散度定义如下：

式中：lnp(Y_Ω)表示模型的证据因子为常数；

被定义为模型因子的下界，记为L(q)；由于lnp(Y_Ω)为常数，当KL散度为0时，即q(Θ)＝p(Θ|Y_Ω)，取到式(13)中模型因子的下界的最大值；采用平均场近似理论，假设变分分布被分解为：

式中：q(Θ)是总体分布，则各个因子q_j(Θ_j)的函数形式可以依次推导出来；基于下界L(q)的最大化给出q(Θ)的第j个因子的优化形式为：

式中：

表示除了Θ_j所有变量联合分布的期望，C为常数；所有变量的先验分布取自指数族且是共轭分布的，即得到未知数的后验分布；其中：

(1)因子矩阵A⁽ⁿ⁾的后验分布：

第n个模式因子矩阵A⁽ⁿ⁾的后验分布由观测数据信息和其他的第n个模式因子矩阵A^(k)(k≠n)的先验分布、超参τ的先验分布推出；由式(11)第n个模式因子矩阵A⁽ⁿ⁾(n＝1,…,N)的各行独立的

服从高斯分布；因此，对于任意n∈{n＝1,…N}的A⁽ⁿ⁾的后验分布可分解为：

其中，后验参数通过下式更新：

式中：Y_(n)代表张量数据Y的第n个模式因子矩阵，V⁽ⁿ⁾为辅助矩阵，E_q[·]代表包含所有变量的后验期望；式中V⁽ⁿ⁾首先利用因子先验Ε_q[Λ]和其它因子矩阵的协方差更新，而因子先验和其它因子矩阵的协方差的权重通过Ε_q[τ]来调节；然后

通过

旋转，并根据Ε_q[τ]进行尺度变换；

(2)超参数λ的后验分布：

超参数λ的后验分布可由N个因子矩阵的信息和超参数λ的先验信息推导出；通过式(15)可知，每个λ_r(r∈{1,…,R})独立的服从Gamma分布，即

其中R为张量的k秩，其中

表示从M个观测数据得到的后验参数，且其更新公式如下：

(3)噪声精度τ的后验分布：

噪声精度τ的后验分布可由观测数据和其本身的超先验信息推导出；结合式(15)可知，噪声精度τ的变分后验分布服从Gamma分布，即q_τ(τ)＝Ga(τ|a_M,b_M)，后验参数a_M、b_M的更新公式如下：

实验验证：

为了验证基于变分贝叶斯平行因子分解的缺失振动信号恢复方法的有效性，实验滚动轴承(型号：LDK UER204)，内圈滚道直径为29.30mm，外圈滚道直径为39.80mm，轴承中径为34.55mm，滚珠直径为7.92mm，滚珠个数为8。实验中通过使用DT9837便捷式动态信号采集器采集振动信号，其中电机转速为2250r/min，相应的转频为f_r＝37.5Hz，采样频率为f_s＝25600Hz，采样间隔为1min，每次采样时长为1.28s。从中选取10000个数据为采样点，观测数据的时域波形和频谱(如图2和3所示)。从图3中可以清晰的看出观测信号的一倍频为37.5Hz、二倍频为73.4Hz、三倍频为115.6Hz、五倍频为184.4Hz和六倍频为231.3Hz。(补充图中所表示的数值)

对上述观测信号进行随机缺失，不同缺失比例情况下两种方法的均方根误差变化对比(如图4所示)。从图4中可以看出，当缺失比例为10％时，VBPF算法的均方根误差为0.1402，LRTC算法的均方根误差为0.1414，当缺失比例为20％时，VBPF算法的均方根误差为0.1724，LRTC算法的均方根误差为0.1984，两种算法的差距不大，随着缺失比例的增大，均方根误差逐渐增大，相较于VBPF算法，LRTC算法对缺失比例的变化更为敏感，变化的幅度更大，当缺失比例为30％时，VBPF算法的均方根误差为0.2035，LRTC算法的均方根误差为0.2450，当缺失比例为40％时，VBPF算法的均方根误差为0.2394，LRTC算法的均方根误差为0.2811，最后当缺失比例达到50％时，VBPF算法的均方根误差为0.2685，LRTC算法的均方根误差为0.3160，两种方法对信号的恢复能力存在明显的差距。

同样的，将随机缺失的数据用0来代替，选取三种缺失比例，得到的缺失信号的频谱(如图5a、5b和5c所示)。在图5a中，当缺失比例为10％时，缺失信号的频谱与原频谱差距不大，在图5b中，当缺失比例为30％时，信号的频谱特征相对于原始的不是很明显，在图5c中，当缺失比例为50％时，,信号的频谱特征已经变得混乱，不能够清晰的判断出原始信号的频谱特征，为后续信号的诊断造成了障碍。为了更直观的对比本发明所提的算法，同时将选中的三种缺失比例的信号用两种不同的方法来进行恢复，恢复的频域结果如图6所示。

在图6a中，缺失比例为10％时，无论是VBPF算法还是LRTC算法，两种方法都能很好的恢复出观测信号的频谱，对比观测信号的频谱，恢复信号的一倍频，二倍频，三倍频、五倍频和六倍频都可以在图6a中显示，并且两种方法的效果差距不大。图6b中，当缺失比例为30％时，LRTC算法并没有将三倍频和五倍频清晰地恢复出来，同时频谱对应的幅值也低于VBPF算法恢复的幅值。图6c中,当缺失比例为50％时，LRTC算法只恢复出了一倍频，其余的倍频已经非常混杂，频率对应的幅值已经不能够清晰的显现，而VBPF算法依旧能够很好的体现三倍频、五倍频和六倍频，且频率对应的幅值也很明显。综上对比实验的频谱图，可以进一步得出VBPF算法针对实际工况的轴承信号，恢复的效果更佳，实际应用更强，且明显优于传统的LRTC算法。

由以上结果可知，本发明建立了变分贝叶斯平行因子分解(VBPF)的预测模型，采用变分贝叶斯(VB)推断信道参数以及因子矩阵的后验信息，更精准的恢复缺失信号。实验结果表明，在相同的缺失比例情况下，本发明所提的方法的均方根误差更小，恢复的性能更好，进一步对比不同恢复方法得到的频谱图，随着缺失比例的增大，VBPF算法的优势体现的更明显，能够有效的恢复信号的频谱特征，更接近原始的信号的频谱，再次验证了所提方法的有效性，解决了振动信号分析中因各种突发情况而引起的振动信号缺失的问题，为缺失振动信号恢复提供了一种崭新的思路。

参考文献

[1]程铁栋,吴义文,罗小燕,等.基于EWT_Hankel_SVD的矿山微震信号特征提取及分类方法[J].仪器仪表学报,2019,40(06):184-194.

[2]Yong L,Yuan R,Wang T,et al.Health Degradation Monitoring and EarlyFault Diagnosis of a Rolling Bearing Based on CEEMDAN and Improved MMSE[J].Materials,2018,11(6):1009.

[3]李志农,朱明.基于变分模态分解的机械故障诊断方法研究[J].兵工学报,2017,38(03):593-599.

[4]胡志峰,李志农,朱彤,王成军.基于同步提取变换的滚动轴承微弱特征增强与提取方法[J].兵器装备工程学报,2021,42(02):234-238.

[5]王友仁,王俊,黄海安.基于非线性短时傅里叶变换阶次跟踪的变速行星齿轮箱故障诊断[J].中国机械工程,2018,29(14):54-61.

[6]屈海清,段腾龙,程汉列.基于小波变换的机械振动信号预处理研究[J].石油和化工设备,2017(2).

[7]魏永合,王明华.基于EEMD和SVM的滚动轴承退化状态识别[J].计算机集成制造系统,2015.

[8]岳应娟,王旭,蔡艳平.内燃机变分模态Rihaczek谱纹理特征识别诊断[J].仪器仪表学报,2017,038(010):2437-2445.

[9]李继猛,黄梦君,谢平,等.同步压缩-交叉小波变换及滚动轴承故障特征增强[J].计量学报,2018,39(02):95-99.

[10]杜伟,房立清,齐子元.基于独立特征选择与流形学习的故障诊断[J].振动与冲击,2018(16):77-82.

[11]李恒,张氢,秦仙蓉,等.基于短时傅里叶变换和卷积神经网络的轴承故障诊断方法[J].振动与冲击,2018,37(19):132-139.

[12]周奇才,沈鹤鸿,赵炯,等.基于改进堆叠式循环神经网络的轴承故障诊断[J].同济大学学报(自然科学版),2019,47(10):122-129.

[13]刘胜,吴迪,李芃.机床主轴承多源信息融合故障诊断[J].航天制造技术,2018,000(002):57-62.

[14]余路.基于过完备字典的缺失振动数据压缩感知重构算法[J].系统工程与电子技术,2017,39(08):209-215.

[15]陈琳升,吴永明,李少波.基于GMDH与SVM的信号缺失情景下轴承故障诊断[J].组合机床与自动化加工技术,2020,No.555(05):87-90.

[16]马云飞,贾希胜,胡起伟,等.基于EMD和BCS的振动信号数据修复方法[J].仪器仪表学报,2019,040(003):154-162.

[17]陈向俊,毛晓松,李黎苹,等.改进张量分解算法及在机械故障诊断中的应用[J].机械设计与制造,2019,000(004):47-50.

[18]Pham T D,Yan H.Tensor Decomposition of Gait Dynamics inParkinson's Disease[J].IEEE Transactions on Biomedical Engineering,2017:1820-1827.

[19]Akmal M,Zubair S,Jochumsen M,et al.A Tensor-Based Method forCompletion of Missing Electromyography Data[J].IEEE Access,2019,(99):1-1.

杨立东,王晶,谢湘,等.基于低秩张量补全的多声道音频信号恢复方法[J].电子与信息学报,2016,38(002):394-399。

Claims (1)

1.基于变分贝叶斯平行因子分解的振动信号缺失的恢复方法，其特征在于，具体步骤如下：

1)通过便捷式动态信号采集器采集机械零部件的时域振动信号S；

2)通过分割采样点将采集的时域振动信号S构造成三维张量Y，步骤如下；

(1)设置时域振动信号S包含K个采样点，然后将信号分成不重叠的数据段，数据段的数量为Q，则每一个数据段包含了P＝K/Q个数据点；且分割成了列数为Q的矩阵X，如下式：

将矩阵X进行转置，得到转置矩阵X^T表示为：

(2)将转置矩阵X^T沿水平进行分割切片构造成三维张量Y；

Y(:,P,:)＝S(:,(1+(P-1)*Q):(P*Q)) (3)

式中，P＝1、2、3…N，而数值Q不变；

3)将三维张量Y进行平行因子分解，其过程如下：

(1)三维张量Y是一个N阶不完全张量，丢失的数据大小为I₁×I₂×...×I_N，Y_Ω表示被观察的数据(i₁,i₂...i_N∈Ω)，其中Ω表示元素的指标集；定义张量O如下式子：

(2)假设三维张量Y为一个真实潜在张量X的噪声观测值，即Y＝X+ε，ε为高斯白噪声，三维张量

的平行因子分解为：

式中，

表示向量的外积，

为张量X对应的第n个模式因子矩阵，ε服从独立同分布的高斯分布，即：

由上式可知，对于第n个模式因子矩阵A⁽ⁿ⁾的表示形式如下式：

4)将分解后的三维张量Y结合贝叶斯方法，引入似然模型，表达式为：

式中：参数τ表示噪声精度，<…>表示N(N≥3)个向量的内积，

表示A⁽ⁿ⁾的第i_n行的行向量；

5)有效利用因子矩阵的先验信息，引入有效精度的后验分布：

式中：Λ^-1＝diag(λ)表示矩阵方差的逆且依赖所有模式的因子矩阵；假设超参数λ是独立的，则先验概率函数如下：

式中：Ga(x|a,b)表示Gamma分布；同时设置噪声精度τ的先验概率函数：

P(τ)＝Ga(τ|a₀,b₀) (10)

6)采用贝叶斯方法处理该模型，推断出包括因子矩阵和超参数在内的所有未知数的参数Θ的后验分布为：

推出缺失信号Y_/Ω的分布预测：

p(Y_/Ω|Y_Ω)＝∫p(Y_/Ω|Θ)p(Θ|Y_Ω)dΘ (12)

7)采用变分贝叶斯算法推导出因子矩阵和超参数的后验分布，从而进一步推断出缺失信号Y_/Ω的分布预测，其过程如下：

采用变分贝叶斯(VB)进行近似求解，通过极小化KL散度来寻找q(Θ)分布近似逼近真实的后验分布p(Θ|Y_Ω)，KL散度定义如下：

式中：lnp(Y_Ω)表示模型的证据因子为常数；

被定义为模型因子的下界，记为L(q)；由于ln p(Y_Ω)为常数，当KL散度为0时，即q(Θ)＝p(Θ|Y_Ω)，取到式(13)中模型因子的下界的最大值；采用平均场近似理论，假设变分分布被分解为：

式中：q(Θ)是总体分布，则各个因子q_j(Θ_j)的函数形式可以依次推导出来；基于下界L(q)的最大化给出q(Θ)的第j个因子的优化形式为：

式中：

表示除了Θ_j所有变量联合分布的期望，C为常数；所有变量的先验分布取自指数族且是共轭分布的，即得到未知数的后验分布；其中：

(1)因子矩阵A⁽ⁿ⁾的后验分布：

第n个模式因子矩阵A⁽ⁿ⁾的后验分布由观测数据信息和其他的第n个模式因子矩阵A^(k)(k≠n)的先验分布、超参τ的先验分布推出；由式(11)第n个模式因子矩阵A⁽ⁿ⁾(n＝1,…,N)的各行独立的

服从高斯分布；因此，对于任意n∈{n＝1,…N}的A⁽ⁿ⁾的后验分布可分解为：

其中，后验参数通过下式更新：

式中：Y_(n)代表张量数据Y的第n个模式因子矩阵，V⁽ⁿ⁾为辅助矩阵，E_q[·]代表包含所有变量的后验期望；式中V⁽ⁿ⁾首先利用因子先验Ε_q[Λ]和其它因子矩阵的协方差更新，而因子先验和其它因子矩阵的协方差的权重通过Ε_q[τ]来调节；然后

通过

旋转，并根据Ε_q[τ]进行尺度变换；

(2)超参数λ的后验分布：

超参数λ的后验分布可由N个因子矩阵的信息和超参数λ的先验信息推导出；通过式(15)可知，每个λ_r(r∈{1,…,R})独立的服从Gamma分布，即

其中R为张量的k秩，其中

表示从M个观测数据得到的后验参数，且其更新公式如下：

(3)噪声精度τ的后验分布：

噪声精度τ的后验分布可由观测数据和其本身的超先验信息推导出；结合式(15)可知，噪声精度τ的变分后验分布服从Gamma分布，即q_τ(τ)＝Ga(τ|a_M,b_M)，后验参数a_M、b_M的更新公式如下：

Images