一种用于模式重建和预测的广义多元奇异谱分析方法

摘要

本发明公开了一种用于模式重建和预测的广义多元奇异谱分析方法，其中，包括：通过嵌入方式获得轨迹矩阵A

说明书

技术领域

本发明涉及奇异谱分析技术领域，尤其涉及一种用于模式重建和预测的广义多元奇异谱分析方法。

背景技术

众所周知，任何模式分类或预测系统的性能通常都会受到噪声的显著影响。模式重建和/或特征增强方法是在提取特征之前减少此性能下降的普遍方法。这种预处理的目标是重建关键拓扑结构，或产生与干净数据的特征近似的特征估计。过去几十年中已出现了许多用于图像模式重建的方法。Schonfeld和Goutsias提出了一种形态滤波器，可实现有噪声二进制图像的最优重建。他们的模式重建方法可得出最优解，方法是使用数学形态学和决策理论的基本工具。他们将差集距离函数定义为图像间的比较尺度。该函数用于引入平均差函数，作为对形态滤波造成的几何和拓扑失真程度的定量测度。Laine等人使用超完备多分辨率表示对乳腺x线图像特征重建进行了研究。他们通过测量已知乳腺x 线图像特征的局部对比，对三种不同种类的小波变换与传统的图像增强技术的效果进行了比较。他们证明了，从多分辨率表示中提取的特征可以提供完成局部对比增强的适应机制。Murray等人提出了调幅与调频(AM-FM)解调方法，用于重建和分析一般图像。基于这些技术，近期出现了一些新的混合方法或变体方法，可进一步改进图像模式表示和重建。

对于一元时间序列(工程学术语中的信号)，还可以通过基于小波变换或矩阵分解的类似技术来增强观测到的混合干净数据和噪声源的特征。奇异谱分析 (SSA)是用于时间序列分析和模式识别的先进技术，融合了经典时间序列分析、线性代数、多元统计、多元几何、模式识别和信号处理等元素。SSA旨在将时间序列分解成少量可解读分量的总和，例如缓慢变化的趋势、振荡分量和“无结构”噪声。基于这些分量，已经出现了一系列针对广泛领域进行模式重建、特征增强、去噪、特征提取、预测和分类的不同算法。以SSA作为预处理过程，可以显著减小神经网络或支持向量机模型的分类或预测误差。基本SSA算法仅适用于单点记录，即一元时间序列。到目前为止，各个领域中的大多数SSA应用都在进行此类一元时间序列分析。

当前，SSA已经发展到多元或多通道奇异谱分析(MSSA)，以适应多元时间序列。MSSA和SSA的主要区别在于其最初的嵌入步骤，该步骤将基于每个时间序列构造的多个轨迹矩阵以行或列的方式堆叠，形成一个大型轨迹矩阵。然后采用奇异值分解(SVD)，以与SSA相同的方式将该整体矩阵分解为特征三元组。Hasani等人将MSSA应用到8项英国工业生产指数的预测中，预测食品、化学品、基本金属和金属制品、机械、电机、车辆、电力和天然气行业的实际产出。他们发现MSSA的性能优于矢量自回归(VAR)，在随机性方面，预测误差较小。Ghodsi等人使用MSSA检测胎儿畸形，在联合6导联心电图中成功提取了胎儿信号。

传感器技术的创新和数字设备更强的可购性成就了今天的“大数据时代”，“大数据”是表示大量数字数据的数量和多样性激增的涵盖性术语。大数据的特点是维度高和样本量大。这两个特征带来了独特的挑战：(I)高维度带来了噪声累积、伪相关和偶然齐次性等问题；(II)维度高且样本量大会导致诸如计算成本高和算法不稳定等问题。MSSA在分析一些实际多元时间序列方面取得了一定的成功。然而，当所有变量的轨迹矩阵连接成一个矩阵时，这个矩阵的维数非常大。由于矩阵较大，而训练样本相对较少，因此难以准确地估计协方差矩阵。除了“高维度”困难和“样本量少”问题外，计算大矩阵的奇异向量还需要很大的计算量。此外，在堆叠单矩阵中，多元时间序列中的空间(结构) 信息会丢失或失真。还存在其他问题，如后续预测或模式分类中的数值不稳定性，以及如何降低计算复杂度和存储要求等。

因此，因此，如何通过奇异谱分析实现对高维度和样本量少的矩阵进行模式重建和预测成为本领域亟待解决的技术问题。

发明内容

本发明旨在至少解决现有技术中存在的技术问题之一，提供一种用于模式重建和预测的广义多元奇异谱分析方法，以解决现有技术中的问题。

作为本发明的第一个方面，提供一种用于模式重建和预测的广义多元奇异谱分析方法，其中，所述用于模式重建和预测的广义多元奇异谱分析方法包括：

通过嵌入方式获得轨迹矩阵A_i，

其中，表示给定的多元时间序列数据集，其中的每个时间序列样本为d是窗口长度，m＝n-d+1表示嵌入的延迟移位行向量的数量，每个时间序列T_i将变换为轨迹矩阵；

对所述轨迹矩阵A_i(i＝1，…，c)进行秩分解得到d个秩为1的矩阵

将所有所述秩为1的矩阵分成多组，并分别对每组内的矩阵进行求和，包括：

以L＝{l₁，l₂，…，l_p}表示一组指数l₁，l₂，…，l_p，

将与组L对应的矩阵定义为

将索引集{1，2，…，d}拆分为不相交的子集L₁，L₂，…，L_s，得到重组矩阵的求和结果

将中的每个重组矩阵转换为长度为n的初始时间序列。

优选地，所述对所述轨迹矩阵A_i(i＝1，…，c)进行秩分解得到d个秩为1的矩阵包括：

分别定义轨迹矩阵集中的平均水平协方差矩阵F和平均垂直协方差矩阵G为：和

将U作为F的特征向量矩阵以及将V作为G的特征向量矩阵得到：

和

定义M_i，j＝U^TA_iV(：，j)＝U^TA_iv_j(i＝1，…，c；j＝1，…，d)；

定义得到d个秩为1的矩阵

优选地，所述将中的每个重组矩阵转换为长度为n的初始时间序列包括：

采用对角平均法，对正斜率对角线分组求平均值，生成每个结式矩阵的重构元素，其中所述重构元素与原始序列中相应元素的位置相同，矩阵包含元素y_gh，1≤g≤m，1≤h≤d；

令d^*＝min(m，d)、m^*＝max(m，d)及n＝m+d-1；

若m＜d，则

根据对角线平均法将矩阵转换为序列其中，对应于对角线i+j＝t+1上矩阵元素的平均值，的公式为：

原始序列T_i被分解为S个序列之和：

优选地，所述将所有所述秩为1的矩阵分成多组，并分别对每组内的矩阵进行求和还包括：

当所述序列包含信号和噪声时，则采用两组索引集L1＝{1，2，…，k}和>2＝{k+1，…，d}；

分别将L₁＝{1，2，…，k}与信号分量关联以及将L₂＝{k+1，…，d}与噪声分量关联。

优选地，所述将所有所述秩为1的矩阵分成多组，并分别对每组内的矩阵进行求和还包括：

根据的样本熵划分每个的索引集{1，2，…，d}，包括：

定义矩阵A_i的样本熵等效于时间序列T_i的样本熵，其中A_i的第e个向量序列为：1≤e≤n-d+1；

向量及其邻近向量之间的距离定义为：

在公差r范围内，与其所有邻近向量之间的平均相似度定义为：

其中，Θ表示赫维赛德函数，且

两个向量序列与d点匹配的概率为：

通过T_i构建A_i的d+1维嵌入向量并按照与相同的方式定义平均相似度及概率

获得在n→∞下的矩阵A_i或时间序列T_i的样本熵估值：

其中，SampEn表示样本熵。

本发明提供的用于模式重建和预测的广义多元奇异谱分析方法，通过将多元时间序列同时分解成其组成分量，而无需将每个时间序列的轨迹矩阵重新排列成一个矩阵，能够用于多元时间序列模式重建、特征增强和去噪，并能够准确地估计协方差矩阵，且在多元时间序列中的空间信息不会发生失真。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：

图1为本发明提供的用于模式重建和预测的广义多元奇异谱分析方法的流程图。

图2为本发明提供的10％NL下有噪声的Lorenz序列的广义多元奇异谱的分解图。

图3为本发明提供的10％NL下有噪声的序列的广义多元奇异谱的分解图。

图4为本发明提供的以2％为间隔从2％到10％的各种噪声级别下有噪声的Lorenz序列的样本熵。

图5为本发明提供的以2％为间隔从2％到10％的各种噪声级别下有噪声的序列的样本熵。

图6为本发明提供的纯Lorenz的重建拓扑结构以及基于添加了10％噪声级别的各种分量之和得出的重构拓扑图。

图7为本发明提供的时间序列的4D模拟图，以及3个不同视角下序列的3D演变图。

图8为本发明提供的纯序列与10％高斯噪声下通过GMSSA获得的前3个(左侧)和全部10个(右侧)分量之和的模式重建图。

图9为本发明提供的10％NL下Lorenz和序列的多步预测均方误差图。

图10为本发明提供的10％NL下Lorenz和序列的多步预测相关系数图。

图11为本发明提供的上海证券交易所(SSE)综合指数2000个数据点的样本熵。

图12为本发明提供的上海证券交易所(SSE)综合指数的实际收盘价格，以及使用支持向量机与广义多元奇异谱结合和支持向量机的预测的价格图。

具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

作为本发明的一个方面，提供一种用于模式重建和预测的广义多元奇异谱分析方法，其中，如图1所示，所述用于模式重建和预测的广义多元奇异谱分析方法包括：

S110、通过嵌入方式获得轨迹矩阵A_i，

其中，表示给定的多元时间序列数据集，其中的每个时间序列样本为d是窗口长度，m＝n-d+1表示嵌入的延迟移位行向量的数量，每个时间序列T_i将变换为轨迹矩阵；

S120、对所述轨迹矩阵A_i(i＝1，…，c)进行秩分解得到d个秩为1的矩阵

S130、将所有所述秩为1的矩阵分成多组，并分别对每组内的矩阵进行求和，包括：

以L＝{l₁，l₂，…，l_p}表示一组指数l₁，l₂，…，l_p，

将与组L对应的矩阵定义为

S140、将索引集{1，2，…，d}拆分为不相交的子集L₁，L₂，…，L_s，得到重组矩阵的求和结果

S150、将中的每个重组矩阵转换为长度为n的初始时间序列。

本发明提供的用于模式重建和预测的广义多元奇异谱分析方法，通过将多元时间序列同时分解成其组成分量，而无需将每个时间序列的轨迹矩阵重新排列成一个矩阵，能够用于多元时间序列模式重建、特征增强和去噪，并能够准确地估计协方差矩阵，且在多元时间序列中的空间信息不会发生失真。

作为一种具体地实施方式，提供了秩分解的具体实施步骤，所述对所述轨迹矩阵A_i(i＝1，…，c)进行秩分解得到d个秩为1的矩阵包括：

分别定义轨迹矩阵集中的平均水平协方差矩阵F和平均垂直协方差矩阵G为：和

将U作为F的特征向量矩阵以及将V作为G的特征向量矩阵得到：

和

定义M_i，j＝U^TA_iV(：，j)＝U^TA_iv_j(i＝1，…，c；j＝1，…，d)；

定义得到d个秩为1的矩阵

作为另一种具体地实施方式，提供了时间序列重建的具体实施过程，所述将中的每个重组矩阵转换为长度为n的初始时间序列包括：

采用对角平均法，对正斜率对角线分组求平均值，生成每个结式矩阵的重构元素，其中所述重构元素与原始序列中相应元素的位置相同，矩阵包含元素y_gh，1≤g≤m，1≤h≤d；

令d^*＝min(m，d)、m^*＝max(m，d)及n＝m+d-1；

若m＜d，则

根据对角线平均法将矩阵转换为序列其中，对应于对角线i+j＝t+1上矩阵元素的平均值，的公式为：

原始序列T_i被分解为S个序列之和：

进一步地，当所述序列包括信号和噪声时，提供了将每个所述秩为1的矩阵分成多组的具体实施过程，所述将所有所述秩为1的矩阵分成多组，并分别对每组内的矩阵进行求和还包括：

当所述序列包含信号和噪声时，则采用两组索引集L₁＝{1，2，…，k}和>2＝{k+1，…，d}；

分别将L₁＝{1，2，…，k}与信号分量关联以及将L₂＝{k+1，…，d}与噪声分量关联。

所述将所有所述秩为1的矩阵分成多组，并分别对每组内的矩阵进行求和还包括：

根据的样本熵划分每个的索引集{1，2，…，d}，包括：

定义矩阵A_i的样本熵等效于时间序列T_i的样本熵，其中A_i的第e个向量序列为：1≤e≤n-d+1；

向量及其邻近向量之间的距离定义为：

在公差r范围内，与其所有邻近向量之间的平均相似度定义为：

其中，Θ表示赫维赛德函数，且

两个向量序列与d点匹配的概率为：

通过T_i构建A_i的d+1维嵌入向量并按照与相同的方式定义平均相似度及概率

获得在n→∞下的矩阵A_i或时间序列T_i的样本熵估值：

其中，SampEn表示样本熵。

通过采用修正版的样本熵，对于噪声可实现更好的相对一致性和更强的稳健性。奇异谱分析旨在将原始序列分解成少量可解读的独立分量的总和，例如缓慢变化的趋势、振荡分量和无结构噪声。样本熵值越大，光谱带越宽，而样本熵越小，光谱带越窄。因此，将样本熵值较大的分组到噪声中是很正常的，而小的样本熵值对应于趋势，中等值对应于振荡分量。对于实际多元时间序列，每个原始时间序列的组数取决于生成单个时间序列的子系统的物理特性。应当理解的是，使用样本熵来区分嵌入矩阵复杂度时，每个的组数可能有所不同，即各个时间序列的L₁，L₂，…，L_s会有所不同。与现有技术中的总方差或最大特征值差距标准相比不同的是，多元时间序列利用相同的截止特征值。

从广义低秩逼近的角度来看，广义多元奇异谱分析方法将矩阵分解为多个秩为一的分量之和，相当于下式的最小二乘法最优解：

其中，v_j是协方差矩阵G的第j个特征向量。为进行简化，假设>1，v₂，…，v_k)和M_i≡(M_i1，M_i2，…，M_ik)，因此上式的最优解可表示为：

J的近似值J^opt的全局最优解可以由下式得出：

R＝V(：，1：k)＝(v₁，v₂，…，v_k)，M_i＝A_iV_k，

由于||A||²＝Tr(A^TA)且Tr(AB)＝Tr(BA)，因此，

当斜率为零时，此二次函数的M_i可以达到最小值：因此，可以得出M_i＝A_iR，进而得出：

因而，J的最简化形式为：

从线性代数角度看，R的最优解由R＝(v₁，v₂，…，v_k)Π给出，其中Π是任意k>由于V是满秩正交矩阵，因此，V^TV＝I，因而轨迹可以表示为：

因此，最优解可以表示为：

因此，广义多元奇异谱为基于原始多元时间序列构造的矩阵提供了低秩逼近的全局最优解。

前面对广义多元奇异谱分析方法进行了详细描述，并证明它是低秩逼近的全局最优解。使用广义多元奇异谱分析方法可以在广泛的应用中以更有效的方式开发用于多元时间序列模式重建、特征增强、特征提取、分类、预测和压缩的各种算法。模式重建和预测是两个基本的重要问题，这不仅限于模式识别和学习系统，在其他领域中也是如此。下面详细说明使用广义多元奇异谱分析方法进行拓扑结构重建以及多元时间序列预测的结果。

1、关于多元时间序列的模式重建。

基于观测到的有噪声时间序列重建拓扑结构通常是模式识别的第一步。然后可采用适当的统计特征或其他特征提取方法来训练分类器或预测器。为了证明所提出的MGSSA方法的可用性，我们利用它来重建嵌入到噪声中的两个基准数据集：三元Lorenz混沌时间序列和四元超混沌时间序列。

Lorenz时间序列具有三个变量x、y、z，而时间序列具有四元变量x、y、z、w。通过在各种噪声级别(NL)下添加独立同分布(i.i.d.)高斯白噪声，将噪声叠加到时间序列上。基于伪邻点法，数据集的嵌入维度设置为10。图2和图3说明了10％NL下有噪声的Lorenz序列和序列的广义多元奇异谱的分解。图2和图3两幅图中，最上方一排均为原始的有噪声数据，下面的10排为状态变量的10个。

其中，图2(A)10％NL下有噪声的Lorenz序列的3个状态变量(x、 y、z)。图2(B)使用广义多元奇异谱时每个变量的10个重构分量为T1到T10。第一个分量显示主要趋势，类似于两个数据集的原始序列。这些结果表明，广义多元奇异谱可以同时将多个原始时间序列分解成其组成分量，突出显示了其在各种多元分析中的应用，包括高维数据模式重建、特征增强、分类、预测和压缩。

图3(A)10％NL下有噪声的序列的4个状态变量(x、y、z、 w)。图3(B)使用广义多元奇异谱时每个变量的10个重构分量为T1到T10。

如图4和图5所示，图4和图5说明了以2％为间隔从2％到10％的各种噪声级别下，相同的Lorenz序列和序列的样本熵。样本熵估值中，公差参数r设置为0.2。图4和图5的两图中，从左到右所有图像中的样本熵对应于从主特征值到次特征值的顺序。对于每个图像，样本熵最初会增大，对应于低频带到高频带的主要分量。然后，样本熵会因各种噪声级别而严重波动，这表示每个变量中的噪声分量。蒙特卡罗模拟表明，当公差r为0.2时，高斯噪声和均匀分布噪声的样本熵都会超过1.3。因此，我们将每个曲线上升阶段中样本熵小于1.3的分量选作主分量。在发明中基于样本熵标准，我们可以准确地得出前三个分量包含有噪声的Lorenz和时间序列中最多的信息。

图6显示了纯Lorenz的重建拓扑结构以及基于添加了10％噪声级别的各种分量之和得出的重构拓扑，其中的第一个图像显示第一个分量的拓扑。经目视检查发现，纯拓扑包络超出了图像A中第一个分量的拓扑包络，表示仅使用第一个分量重建模式时会发生信息丢失。图6的图像C和D分别显示由前5 个和10个分量之和形成的重构拓扑结构。这两个图像的包络超出了干净数据的包络。因而前三个分量包含多元Lorenz时间序列最必要的信息(如图4所示)。因此，可以通过选择使用广义多元奇异谱得出的适当数量的分量来重建有噪声的Lorenz时间序列的拓扑轨迹。

四元序列具有更复杂的动态特性。图7中的图像A显示了时间序列的独特4D结构，而其余3个图像从3个不同的视角展示其3D演变。

图8为纯序列与10％高斯噪声下通过GMSSA获得的前3个 (左侧)和全部10个(右侧)分量之和的模式重建图。左侧(图像A1-A4) 显示纯序列的3D拓扑结构，以及基于增加了10％噪声级别的前三个分量之和重建的3D拓扑结构。该图右侧(图像B1-B4)显示干净序列的3D 拓扑结构，以及有噪声的序列所有分量之和的3D拓扑结构。与Lorenz序列类似，在表示4个子系统(即(x，y，z)、(x，y，w)、(y，z，w)和(x，z，w))的图像 A1-A4中，序列干净数据的拓扑结构与前三个分量之和完全匹配。该示例进一步说明广义多元奇异谱可以使用样本熵标准有效且高效地同时分解多元时间序列，以及基于有噪声的时间序列重建模式。

2、关于多元时间序列的预测。

通过采用将支持向量机与广义多元奇异谱结合，实现对多元时间序列的预测，具体可以包括如下步骤：

(a1)将多元数据的每个时间序列嵌入具有近似维度的矩阵中；

(a2)使用前文所述的广义多元奇异谱分解矩阵，并选择最富信息的分量，从而使用样本熵标准重构矩阵

(a3)构建预测模型输入输出对。输入集包括同一行中的向量相应的输出由时间序列T_i的组成，其中h表示所需预测时域。

(a4)训练支持向量机的预测模型。

(a5)按所需提前期h预测时间序列。

下面通过实验对支持向量机与广义多元奇异谱结合的方法进行评估。从前三个基本矩阵(即))之和获得的行向量用作构建预测模型的输入。Lorenz和模型的输入向量维度分别为30和40。这两个模型的输出都是Lorenz和序列的x变量。训练集由1500个点组成，测试集由500个点组成，用于评估预测效果。Matlab工具箱LIVSVM用于训练支持向量机模型。根据均方误差(MSE)和相关系数(CC)评估预测效果。图9 显示10％NL下，Lorenz和序列1到20步预测的MSE。还评估了不含GMSSA预增强的SVM模型相应的MSE，以便进行比较。图10显示同一预测模型1到20步预测的相关系数。随着预测提前期增长，可以很明显地看出MSE增长，CC降低。图9和图10所示的结果清楚地表明，支持向量机与广义多元奇异谱结合的这一方法可以有效提高多元时间序列的预测精确度。

下面将广义多元奇异谱应用于金融数据集-上海证券交易所(SSE)综合指数，以测试广义多元奇异谱对实际数据的适用性。

由于证券市场是非常复杂的非线性系统，其价格波动的统计特征不稳定。将收盘价格、开盘价格、最高价格和最低价格时间序列合并为一个多元时间序列，预测收盘价格的演变。这种多源输入已应用于使用ANN预测收盘价格。与上述评估支持向量机与广义多元奇异谱结合的方法类似的方式预测收盘价格。数据长度是自2005/09/26开始的2000个每日数据点。前1500个数据点用于训练，其余数据点用于预测测试。图11显示4个序列的样本熵分布，基于此分布将主分量的数目确定为2个。图12显示实际收盘价和使用这两种方法提前5天预测的收盘价。表1总结了预测时间从1到5天的MSE和CC，进一步说明了混合广义多元奇异谱方法对于实际数据预测的优越性。

表1上海证券交易所综合指数预测的RMS和CC

本发明提供的用于模式重建和预测的广义多元奇异谱分析方法，可以将多个时间序列同时分解成各自的组成分量，包括缓慢变化的趋势、振荡分量和无结构噪声。可以使用通过广义多元奇异谱获得的这些分量来开发各种模式重建、特征增强、去噪、特征提取、分类和预测方法。由于广义多元奇异谱处理每个时间序列的轨迹矩阵，而不是所有轨迹矩阵的高维矩阵堆叠，因此在广义多元奇异谱中不再需要解决维度难题。此外，由于要分析的输入特征向量实际上是训练样本的行向量，所以特征集明显扩大。因此，广义多元奇异谱中不存在样本量小的问题。因此，广义多元奇异谱中协方差矩阵的估值比多通道奇异谱分析更加稳健，更加准确。此外，通过使用原始矩阵，二维空间信息能够完好地保留下来。确定分量的最佳数量是分量学习或建模的根本问题。本发明中提供的采用样本熵标准估计主要分量，可以用于任何种类的奇异谱分析、多元奇异谱分析或其他分量分析技术。

可以理解的是，以上实施方式仅仅是为了说明本发明的原理而采用的示例性实施方式，然而本发明并不局限于此。对于本领域内的普通技术人员而言，在不脱离本发明的精神和实质的情况下，可以做出各种变型和改进，这些变型和改进也视为本发明的保护范围。