CN108021871A - 一种基于主成分分析的特征频率提取方法 - Google Patents

Abstract

本发明公开了一种基于主成分分析的特征频率提取方法，包括以下步骤：首先，两个位移传感器以相互垂直的方式采用一定的采样频率采集试验台的转子位移数据；第二步，将采集到的信号x(t)，利用消直流算法滤出原始信号的直流分量，然后构造Hankel矩阵X，第三步，求矩阵X的协方差矩阵C，并构造协方差矩阵特征值分布图，根据待提取频率在幅值谱中的排列顺序，在分布图中选出对应的两个特征值，并选取对应的特征向量进行信号重构，然后，将重构后的矩阵叠加原始信号的均值并获得矩阵从矩阵中恢复出待提出的信号。本发明的方法对指定的特征频率(单个或多个特征频率)的提取或消除上效果十分显著。

Description

一种基于主成分分析的特征频率提取方法

技术领域

本发明涉及属于机械故障诊断与模态分析领域，特别涉及一种基于主成分分析的特征频率提取方法。

背景技术

从旋转机械的动态信号中提取出特征信息是故障诊断和模态分析的基础，近年来，学者们提出了多种用于轴承、齿轮的故障诊断方法，主要包括：小波分析、经验模态分解(EMD)、盲源分离和稀疏分解等算法，这些方法在特征频率的提取上，取得了非常好的效果，但也存在一些不足。例如稀疏分解算法在提取特征频率的时候需要了解故障信号的特点，构造符合信号特征的核函数，才能取得良好的提取效果，而且此算法的计算量较大。小波分析属于一种带通滤波器，存在着一定的带宽和过渡带，易受频带内的噪声干扰。

近年来，主成分分(Principal Component Analysis，PCA)在图像处理、数据压缩、故障诊断、神经网络、模式识别、小波变换等方面获得了广泛的应用。上述的这些应用均是利用PCA算法进行信号降噪或数据压缩。其中，通过累积贡献率L_k确定有效主成分个数，如下式：

其中，L_k为累积贡献率；λ_i为协方差矩阵C的特征值；k为有效特征值的个数；p为所有特征值的个数；i为1到k正整数，j为1到p正整数。

虽然这种通过采用累积贡献率的方法在信号降噪、数据降维等方面取得了良好的效果，但本发明提出一种基于PCA算法提取指定的特征频谱(单个频率)，或消除指定的特征频率的算法。在此之前PCA算法主要是用来信号降噪，将降噪后的算法结合其它算法进行故障特征的提取。单独的采用PCA算法提取特征频率的方法尚未见报道。

发明内容

针对现有技术中存在的技术问题，本发明的目的是：提供一种预制梁柱节点。

本发明的目的通过下述技术方案实现：一种基于主成分分析的特征频率提取方法，包括以下步骤：

首先，两个位移传感器以相互垂直的方式采用一定的采样频率采集试验台的转子位移数据；

第二步，将采集到的信号x(t)，利用快速傅里叶变换滤出原始信号的直流分量，然后将消除直流分量后的信号构造Hankel矩阵X，具体如下；

式中，X为m×n阶的矩阵，x_ij为X向量中的第i行，第j的元素，其中1≤i≤m，1≤j≤n，m为矩阵的行数，n为矩阵的列数；

同样地，将X表示成列向量的形式，即X＝[x₁，x₂…x_i…x_m]^T，x_i表示矩阵X中第i行，包含n个元素的行向量，T表示向量的转置；

当信号x(t)的数据长度为L时，L为偶数时，m＝L/2，n＝L/2+1；L为奇数时，m＝(L+1)/2，n＝(L+1)/2；

第三步，求矩阵X的协方差矩阵C，如下式，并求其特征值λ_i，按从大到小排列为λ₁,λ₂…λ_m，同时，构造协方差矩阵特征值分布图，并求出对应特征值的特征向量为α₁,α₂…α_m；

式中：cov(x_i,x_j)＝E[(x_i-E(x_i))(x_j-E(x_j))]；E(x_i)为求向量x_i的均值，cov()求协方差运算；C为协方差矩阵；x_i表示矩阵X中第i行，包含n个元素的行向量，1≤i≤m；

根据特征值λ_i在协方差矩阵特征值分布图中的排序，选择与某一具体频率对应的两个特征值及相应的特征向量进行重构，若某一频率在原始信号幅值谱中的幅值大小排序为k，则选择协方差矩阵特征值分布图中的第2k-1与第2k个特征值对应的特征向量α_i进行信号重构，得到一个新的矩阵如下所示；

然后，将重构后的矩阵叠加原始信号的均值并获得矩阵

一般，从矩阵中恢复出信号有两种方法；

方法1：采用平均法恢复信号

式中，为近似重构矩阵中第i行，第j列的元素；m为原始矩阵X的行数，n为列数；i，j，m，n，k均为正整数；

方法2：采用伪逆矩阵的方式恢复信号

式中，为恢复后的矩阵；为恢复的一维信号；M⁺称为M的伪逆；M为m个阶数为n的单位矩阵构成；其中，M⁺＝(M^TM)^-1M^T∈R^mn×l，

l＝m+n-1，m，n，l为正整数；

信号就是从原始信号x(t)中提取的特征频率成分。

优选的，选取工频干扰50Hz在特征值分布图中对应的两个特征值及对应的特征向量进行重构，然后从原始信号中减去重构后的特征信号，便滤除了50Hz工频干扰的信号。

优选的，完成特征频率的提取后，进行轴心轨迹的合成，将第一个加速度传感器采集到的信号进行特征频率的提取，将处理后的信号作为x轴；然后，将处理后的第二个加速度传感器采集的信号作为y轴，进行二维轴心轨迹的合成。

优选的，提取后的信号为基频、二倍频、三倍频或1/2基频的组合。

优选的，合成三维轴心轨迹，即将t作为z轴，合成位移量与时间的三维轴心轨迹。

优选的，滤除原始信号中的直流分量的方法采用平均值法或高通滤波器法。

优选的，设有m个随机向量x₁，x₂…x_m，每个向量x_i含有n个样本，即x_i＝(x_i1，x_i2…x_in)，则可构造一个m×n阶的矩阵X：

式2中，X为m×n阶的矩阵；x_ij为X向量中的第i行，第j的元素，其中1≤i≤m，1≤j≤n；m为矩阵的行数，n为矩阵的列数；

X可表示为X＝[x₁，x₂…x_m]^T，对矩阵X进行主成分分析，可得l个新的变量指标y_i，i＝1,2…l(l≤m)表达式如下：

式3中，y_i为l个新的变量指标，y_i∈R^1×n，其中1≤i≤l≤n；x_i为矩阵X的行向量；X为m×n阶的矩阵；α_i为X的协方差矩阵第i个特征值由大到小对应的特征向量，α_ij为特征向量α_i内的第j个元素，1≤j≤n；T为转置标识符；i，j为正整数；

其中y_i∈R^1×n，由主成分分析的定义可知α_i＝(α_i1,α_i2,…,α_im)^T为X的协方差矩阵第i个特征值由大到小对应的特征向量，且α_i满足：

式4中，α_ij为特征向量α_i内的第j个元素，1≤i≤n，1≤j≤n，i，j为正整数；

式2所示矩阵X的协方差矩阵为：

式5中：cov(x_i,x_j)＝E[(x_i-E(x_i))(x_j-E(x_j))]，E(x_i)为求向量的均值，cov()求协方差运算；C为协方差矩阵；x_i为第i个变量的n个样本组成的行向量，1≤i≤m；

由主成分分析原理可知，协方差矩阵C的特征方程如下：

Cα_i＝λ_iα_i 式6

式6中，C为协方差矩阵；α_i为协方差矩阵C的第i个征向量；λ_i是矩阵C的特征值；

式6中，λ_i是矩阵C的特征值，α_i是与λ_i对应的特征向量，把特征值按由大到小的顺序排列λ₁＞λ₂＞…＞λ_l。根据公式1和公式3便可以实现数据降维，由原来的m个变量转换成新的l个变量；若要进行信号处理，则根据公式3进行信号的重构，两边左乘α_i，并进行求和，得：

式7中，I为单位矩阵；y_i为降维后的矩阵；m为X向量的行数；α_i为协方差矩阵C的第i个特征向量，T为转置标识符；

若根据累积贡献率L_k选择前k个主成分进行重构，可得到一个近似的矩阵：

式8中，为恢复后的矩阵；y_i为降维后的矩阵；k为有效特征值的个数；α_i为协方差矩阵C的第i个特征向量；T为转置标识符；

与原矩阵X相比，近似重构的矩阵保留了原始矩阵的大部分信息，却消除了原始矩阵的冗余信息；

获得重构矩阵后按矩阵构成方式恢复信号即可，通过矩阵的逆变换从重构矩阵得到恢复后的信号即：

式9中，为恢复后的矩阵；为恢复的一维信号；M⁺称为M的伪逆；M为m个阶数为n的单位矩阵构成。其中，M⁺＝(M^TM)^-1M^T∈R^mn×l，l＝m+n-1，m，n，l为正整数。

式中，M⁺＝(M^TM)^-1M^T∈R^mn×l，M⁺称为M的伪逆。通过求伪逆矩阵恢复出信号实际就是对矩阵各副对角线上的元素求平均，这与采用平均法恢复信号的结果一致，由于矩阵M为稀疏矩阵，尤其是当m与n的取值较大时，计算伪逆矩阵M⁺时所需要的内存和时间会成倍的增加，因而本文采用平均法从矩阵中恢复信号表达式如下：

式10中，为近似重构矩阵中第i行，第j列的元素；m为原始矩阵X的行数，n为列数；i，j，m，n，k均为正整数；根据公式十便可以恢复出信号

本发明相对于现有技术具有如下的优点及效果：

1、本发明的方法对指定的特征频率(单个或多个特征频率)的提取或消除上效果十分显著，而一般的带通滤波器是很难做到的，滤波器总是会存在一定的带宽，并且还有过渡带的影响，除非两个频率相隔较远，否则滤波器在分离指定的特征频率时总是会受到临近频率的影响。而采用本发明提出的PCA信号分离方法，即使两个频率相差1Hz也能够不受影响地将它们提取出来。

附图说明

图1a是含有1种频率成分信号1的协方差矩阵C的特征值分布图；

图1b是含有1种频率成分信号2的协方差矩阵C的特征值分布图；

图1c是含有1种频率成分信号3的协方差矩阵C的特征值分布图；

图2a是含有2种频率成分信号1的协方差矩阵C的特征值分布图；

图2b是含有2种频率成分信号2的协方差矩阵C的特征值分布图；

图2c是含有2种频率成分信号3的协方差矩阵C的特征值分布图；

图3a是含有3种频率成分信号1的协方差矩阵C的特征值分布图；

图3b是含有3种频率成分信号2的协方差矩阵C的特征值分布图；

图3c是含有3种频率成分信号3的协方差矩阵C的特征值分布图；

图4a是仿真信号的时域图；

图4b是仿真信号的频域图；

图4c是仿真信号的特征值分布图；

图4d是第一组特征值重构的幅值谱；

图4e是第二组特征值重构的幅值谱；

图4f是第三组特征值重构的幅值谱；

图4g是第1、2和3组特征值重构与原始无噪信号对比；

图5是在试验台转子A侧的轴颈处的电涡流位移传感器安装示意图；

图6是在试验台转子B侧的轴颈处的电涡流位移传感器安装示意图；

图7a是D11信号隔离直流分量后的时域波形；

图7b是D12信号隔离直流分量后的时域波形；

图7c是D11信号(x轴)的幅值谱；

图7d是D12信号(y轴)的幅值谱；

图8是A端面D11信号与D1信号2原始轴心轨迹图；

图9a是D11信号的协方差矩阵的特征值分布图；

图9b是D12信号的协方差矩阵的特征值分布图；

图10a是D11信号利用本发明的方法提取前两倍的频域图；

图10b是D12信号利用本发明的方法提取前两倍的频域图；

图11是利用本发明的方法提取到的基频和二倍频信号合成的轴心轨迹图；

图12a是D11提取前5倍频的频域图；

图12b是D12提取前5倍频的频域图；

图13a是D11提取前5倍频的频谱图；

图13b是D12提取前5倍频的频谱图；

图14是D11信号与D12信号前5倍频合成的轴心轨迹；

图15a是小波包提取D11的第1倍信号的频谱图；

图15b是小波包提取D11的第2倍信号的频谱图；

图15c是小波包提取D11的第3倍信号的频谱图；

图15d是小波包提取D11的第4倍信号的频谱图；

图16a是谐波小波提取D11一倍频信号的频谱图；

图16b是谐波小波提取D11二倍频信号的频谱图；

图16c是D11前2倍频谱图；

图16d是谐波小波提取D12一倍频信号的频谱图；

图16e是谐波小波提取D12二倍频信号的频谱图；

图16f是D12前2倍频谱图；

图17是谐波小波提取D11与D12的前2倍合成的轴向轨迹图；

图18是D11滤除50Hz后的频谱图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例一：

一种基于主成分分析的特征频率提取方法，包括以下步骤：

首先，采用两个位移传感器以采用相互垂直的安装方式(如图5和图6)采用一定的采样频率采集试验台的转子位移数据；

第二步，将采集到的信号x(t)，利用快速傅里叶变换滤出原始信号的直流分量，然后将消除直流分量后的信号构造Hankel矩阵X，具体如下；

式中，X为m×n阶的矩阵，x_ij为X向量中的第i行，第j的元素，其中1≤i≤m，1≤j≤n，m为矩阵的行数，n为矩阵的列数；

同样地，将X表示成列向量的形式，即X＝[x₁，x₂…x_i…x_m]^T，x_i表示矩阵X中第i行，包含n个元素的行向量，T表示向量的转置；

当信号x(t)的数据长度为L时，L为偶数时，m＝L/2，n＝L/2+1；L为奇数时，m＝(L+1)/2，n＝(L+1)/2；

第三步，求矩阵X的协方差矩阵C，如下式，并求其特征值λ_i，并按从大到小排列为λ₁,λ₂…λ_m，同时，构造协方差矩阵特征值分布图，并求出对应特征值的特征向量为α₁,α₂…α_m；

式中：cov(x_i,x_j)＝E[(x_i-E(x_i))(x_j-E(x_j))]；E(x_i)为求向量的均值，cov()求协方差运算；C为协方差矩阵；x_i表示矩阵X中第i行，包含n个元素的行向量，1≤i≤m；

根据特征值的λ_i在协方差矩阵特征值分布图中的排序，选择与某一具体频率对应的两个特征值及相应的特征向量进行重构，若某一频率在原始信号幅值谱中的幅值大小排序为k，则选择协方差矩阵特征值分布图中的第2k-1与第2k个特征值对应的特征向量α_i进行信号重构，得到一个新的矩阵

将重构后的矩阵叠加原矩阵的均值获得矩阵

从矩阵中恢复出信号有两种方法；

方法1：采用平均法恢复信号

式中，为近似重构矩阵中第i行，第j列的元素；m为原始矩阵X的行数，n为列数；i，j，m，n，k均为正整数；

方法2：采用伪逆矩阵的方式恢复信号

式中，为恢复后的矩阵；为恢复的一维信号；M⁺称为M的伪逆；M为m个阶数为n的单位矩阵构成；其中，M⁺＝(M^TM)^-1M^T∈R^mn×l，

l＝m+n-1，m，n，l为正整数；

信号就是从原始信号x(t)中提取的特征频率成分。

优选的，选取工频干扰50Hz在特征值分布图中对应的两个特征值及对应的特征向量进行重构，然后从原始信号中减去重构后的特征信号，便滤除了50Hz工频干扰的信号。

优选的，完成特征频率的提取后，进行轴心轨迹的合成，将第一个加速度传感器采集到的信号进行特征频率的提取，将提取后的信号作为x轴；然后，将处理后的第二个加速度传感器采集的信号作为y轴，进行二维轴心轨迹的合成。

优选的，提取后的信号为基频、二倍频或三倍频的组合。

优选的，合成三维轴心轨迹，即将t作为z轴，合成位移量与时间的三维轴心轨迹。

优选的，滤除原始信号中的直流分量的方法采用平均值法或高通滤波器法。

实施例二：

一种基于主成分分析的特征频率提取方法，设有m个随机向量x₁，x₂…x_m，每个向量x_i含有n个样本，即x_i＝(x_i1，x_i2…x_in)，则可构造一个m×n阶的矩阵X：

式(2)中，X为m×n阶的矩阵；x_ij为X向量中的第i行，第j的元素，其中1≤i≤m，1≤j≤n；m为矩阵的行数，n为矩阵的列数。

X可表示为X＝[x₁，x₂…x_m]^T，对矩阵X进行主成分分析，可得l个新的变量指标y_i，i＝1,2…l(l≤m)表达式如下：

式(3)中，y_i为l个新的变量指标，y_i∈R^1×n，其中1≤i≤l≤n；x_i为矩阵X的行向量；X为m×n阶的矩阵；α_i为X的协方差矩阵第i个特征值(由大到小)对应的特征向量，α_ij为特征向量α_i内的第j个元素，1≤j≤n；T为转置标识符；i，j为正整数。

其中y_i∈R^1×n，由主成分分析的定义可知α_i＝(α_i1,α_i2,…,α_im)^T为X的协方差矩阵第i个特征值(由大到小)对应的特征向量，且α_i满足：

式(4)中，α_ij为特征向量α_i内的第j个元素，1≤i≤n，1≤j≤n，i，j为正整数。式(2)所示矩阵X的协方差矩阵为：

式(5)中：cov(x_i,x_j)＝E[(x_i-E(x_i))(x_j-E(x_j))]，E(x_i)为求向量的均值，cov()求协方差运算；C为协方差矩阵；x_i为第i个变量的n个样本组成的行向量，1≤i≤m。由主成分分析原理可知，协方差矩阵C的特征方程如下：

Cα_i＝λ_iα_i 式(6)

式中，C为协方差矩阵；α_i为协方差矩阵C的第i个特征向量；λ_i是矩阵C的特征值。

式中，λ_i是矩阵C的特征值，α_i是与λ_i对应的特征向量。把特征值按由大到小的顺序排列λ₁＞λ₂＞…＞λ_l。根据公式(1)和公式(3)便可以实现数据降维，由原来的m个变量转换成新的l个变量。若要进行信号处理，则根据公式(3)进行信号的重构，两边左乘α_i，并进行求和，得：

式中，I为单位矩阵；y_i为降维后的矩阵；m为X向量的行数；α_i为协方差矩阵C的第i个特征向量，T为转置标识符。

若根据累积贡献率L_k选择前k个主成分进行重构，可得到一个近似的矩阵：

式中，为恢复后的矩阵；y_i为降维后的矩阵；k为有效特征值的个数；α_i为协方差矩阵C的第i个特征向量；T为转置标识符。

与原矩阵X相比，近似重构的矩阵保留了原始矩阵的大部分信息，却消除了原始矩阵的冗余信息(例如噪声与工频干扰)。

获得重构矩阵后只需按矩阵构成方式恢复信号即可。一般通过矩阵的逆变换从重构矩阵得到恢复后的信号即：

式(9)中，为恢复后的矩阵；为恢复的一维信号；M⁺称为M的伪逆；M为m个阶数为n的单位矩阵构成。其中，M⁺＝(M^TM)^-1M^T∈R^mn×l，l＝m+n-1，m，n，l为正整数。

式中，M⁺＝(M^TM)^-1M^T∈R^mn×l，M⁺称为M的伪逆。通过求伪逆矩阵恢复出信号实际就是对矩阵各副对角线上的元素求平均，这与文献中报道的采用平均法恢复信号的结果一致。由于矩阵M为稀疏矩阵，尤其是当m与n的取值较大时，计算伪逆矩阵M⁺时所需要的内存和时间会成倍的增加，因而此处采用平均法从矩阵中恢复信号表达式如下：

式(10)中，为近似重构矩阵中第i行，第j列的元素；m为原始矩阵X的行数，n为列数；i，j，m，n，k均为正整数。根据公式(10)便可以恢复出信号

实施例三：

为了探讨信号有效特征值的数量与频率个数的内在关系，构造不同幅值、频率和相位的信号，如下式：

式中，k为有效频率的个数；A_i为信号的幅值；f_i为信号的频率；为信号的相位。

对x(t)信号构造Hankel矩阵，所构造矩阵行数为m，列数为n，采样频率1024Hz，采样点数为4096，采用上述的PCA分解算法，进行信号分解与重构，此时m＝2048，n＝2049。

(1)首先，探讨当k＝1，A_i＝1，构造3组信号，每组信号仅含有1个有效频率成分，f_i取10Hz，20Hz，30Hz、相位取10，20，30时，主成分特征值分布如图1所示，特征值个数q取范围[0,2048]，而本例仅显示前50个特征值的分布情况。

从图1可以看出，当k＝1时，所构造的每组信号仅含有一个有效频率成分，幅值A_i相同，频率f_i和相位不同，由三组信号的协方差矩阵C的特征值分布图中可看出，每个信号产生两个非零的特征值λ₁与λ₂，一前一后紧密排列，而且三组信号所产的两个特征值对应大小相同。

(2)其次，探讨当k＝2，A_i取1和0.8时，同样构造3组信号，每组含有2个有效频率成分，f_i取10Hz和20Hz，30Hz和40Hz，50Hz和60Hz，而取10，20，30，得到的特征值分布如图2所示；

从图2可以看出，当k＝2，A_i取1和0.8时，所构造的三组信号中，每组信号都含有两个有效频率成分，三组信号中相同的频率成分的幅值A_i大小相同，但对应的频率f_i和相位不同。由三组信号协方差矩阵C的特征值分布图可看出，每个信号产生两组非零的特征值，每组特征值包含两个特征值λ_i与λ_i+1，一前一后紧密排列。同时还可以看出三组信号所产生的两组特征值大小相同。

图2与图1比较可以发现，两幅图中第一组特征值的大小基本相同，而图2是在图1幅值为1的频率成分的基础上增加了一个幅值为0.8的频率成分。所以可以确定图2中的第一组非零的特征值λ₁与λ₂是由幅值A₁为1的频率成分产生的，而图2中的第二组频率成分是由幅值A₂为0.8的频率成分产生的。

(3)当k＝3，A_i取1、0.8和0.6时，构造3组信号，每组信号含有3个有效频率成分，第一组信号的频率f_i为10Hz、20Hz和30Hz；第二组信号的频率f_i为40Hz、50Hz和60Hz；第三组信号的频率f_i为70Hz、80Hz和90Hz。而三组信号对应的相位分别取10，20，30，构造的信号和得到的特征值分布如图3所示；

从图3可以看出，当k＝3，A_i取1、0.8和0.6时，所构造的三组信号中，每组信号都含有三个有效频率成分。信号中的频率f_i(i＝1,2,3)幅值A_i大小对应相同，但频率f_i和相位不同。由三组信号的协方差矩阵C的特征值分布图中同样可看出，每组信号x_i(t)产生三组非零的特征值，每组特征值包含两个特征值λ_i与λ_i+1，一前一后紧密排列；而且每个信号所产的三组特征值大小与另外两个信号所产生的特征值的大小相同。

对比图3与图2、图1，可以发现，此三图中第一组特征值的大小基本相同。再由图1可知，第一组非零的特征值λ₁与λ₂是由幅值A₁取1的频率成分产生的；再对比图2与图3发现，第二组特征值的大小也基本相同，所以第二组非零的特征值λ₃与λ₄是由幅值A₂取0.8的频率成分产生的；而图3中的第三组频率成分λ₅与λ₆是由幅值A₃取0.6的频率成分产生的。

若继续增加信号的有效频率成分，可得到同样结果。因此，从图1、图2和图3对比得出以下结论：信号x(t)在Hankel矩阵的构造模式下，l＝min(m,n)，m为行数，n为列数，k为信号x(t)中的有效频率个数，在满足采样定理的条件下，如果l＞2k，则有：

(1)信号中每个频率成分最多只产生两个非零的特征值，一前一后紧密排列；

(2)信号有效特征值的个数只与信号的频率个数相关，而与幅值A_i、频率f_i和相位的大小均无关。

(3)信号有效特征值在协方差矩阵C的特征值分布图中的特征值排列顺序由信号频率的幅值A_i决定，某频率的幅值越大，它所产生的两个特征值就越大，排列越靠前。

根据以上频率与特征值的关系，提出一种基于PCA分离特征频率的方法，步骤如下：

(1)对于一个确定的信号x(t)，首先通过FFT变换滤出原始信号的直流分量，然后利用消除直流分量的信号构造Hankel矩阵X；

(2)求矩阵X的协方差矩阵C，并求其特征值λ_i，按从大到小排列为λ₁,λ₂…λ_m，构造协方差矩阵特征值分布图，并求得对应特征向量为α₁,α₂…α_m；

(3)根据特征值的λ_i分布，选择与某一具体频率对应的两个特征值及相应的特征向量进行重构，若某一频率在原始信号幅值谱中的幅值大小排序为k，选择协方差矩阵特征值分布图中的第2k-1与第2k个特征值对应的特征向量进行信号重构，得到一个新的矩阵；

(4)将重构后的矩阵叠加原矩阵的均值获得矩阵

(5)从矩阵中采用平均法恢复出信号该信号就是需要提取的特征频率分量。

实施例四：

为进一步证明信号中的频率成分与特征值分布之间的关系，即实施例三中的结论。构造信号x₄(t)，如下式：

式中，s(n)为标准差1.4的高斯白噪声，对此信号以采样频率1024Hz采集4096个数据点，结果如图4(a)所示，其幅值谱如图4(b)。

现在利用本发明提出的方法来分离这三个频率，原信号的特征值分布如图4(c)，可见除含有三组比较大特征值外，由于噪声s(n)的存在，x₄(t)的大部分特征值都是非零的特征值，噪声信号幅值相对较小，因此其所产生的特征值分布在有效特征值之后。

将图4(c)中的第一组特征值λ₁与λ₂进行重构，得到信号的幅值谱如图4(d)，可见频率为30Hz，幅值为0.9655，这与原始信号中的sin(60πt+20)频率成分一致。利用图4(c)中的第二组特征值λ₃与λ₄进行重构得到信号的幅值谱如图4(e)，可见频率为50Hz，幅值为0.8052，这与原始信号中的0.8sin(100πt+20)频率成分一致。第三组特征值λ₅与λ₆重构得到信号的幅值谱如图4(f)，可见频率为40Hz，幅值为0.587，这与原始信号中0.6sin(80πt+20)频率成分一致。

上述的重构结果进一步说明原始信号频率的幅值大小决定协方差矩阵C的特征值分布图中的排列顺序：原始信号中某一频率的幅值越大，其特征值在分布图中的排列顺序也就越靠前，且特征值成对出现，一前一后紧密排列。由于噪声的幅值较小，从图4(c)中可看出噪声产生的特征值分布在有效特征值之后。同时，从图4(d)、(e)和(f)可以看出重构后的信号幅值谱十分纯净，这也是一般带通滤波器难以实现的。

将采用三组特征值重构的信号相加，得到信号如图4(g)，可见结果与原始无噪声信号几乎完全重合，没有相位偏移。由式(8)可知，PCA算法从原始信号中提取的分量信号和原始信号是一种加减法的关系，这也是为什么不会产生相位偏移的原因。PCA算法不同于一般的滤波器，滤波器多是通过卷积运算实现的，例如小波滤波器和FIR滤波器等，如果滤波器不具有线性相位特性，则其分离出来的信号会出现相位失真。

实施例五：

信号频率幅值与特征值分布关系的证明：

设某一信号为将信号构造一个m×n的Hankel矩阵X，且矩阵X可表示为X＝[x₁，x₂…x_m]^T，其中x_i＝(x_i1，x_i2…x_in)由主成分分析的定义可知，矩阵X的协方差矩阵C为：

C＝E[(x-E(x))·(x-E(x))^T] 式(13)

式中，C为矩阵X的协方差矩阵；E(x)为向量x的均值；T为转置标识符。

根据公式(6)可知，构造如下公式：

Cα_i(Cα_i)^T＝λ_iα_i(λ_iα_i)^T 式(14)

式中，为恢复后的矩阵；λ_i为协方差矩阵C的第i个特征向量对应的特征值；α_i为协方差矩阵C的第i个特征向量；T为转置标识符。

根据α_iα_i ^T＝1可得：

|C|²＝λ_i ² 式(15)

式中，C为协方差矩阵；λ_i为协方差矩阵C的第i个特征向量对应的特征值。

将信号x(t)构造Hankel矩阵后带入式(13)；

可知协方差矩阵C的能量与a²成正比，a²越大则矩阵C的能量越大，根据式(15)可知λ_i也就越大。由实施例三可知，因为信号x(t)只含有一个频率成分，所以在协方差矩阵C的特征值分布图中仅有两个有效特征值，而且这两个特征值的大小与频率ω的幅值a成比例，进一步地说明频率成分的幅值越大，协方差矩阵C的特征值分布图中对应的特征值也越大。根据这一结论，只要确定原信号幅值谱中某一频率的幅值排序，便可以根据特征值分布重构出对应的频率成分，从而实现单个或多个特征频率的提取。例如在幅值谱中，某一频率的幅值大小排序为k，那么选择协方差矩阵特征值分布图中的第2k-1与第2k个特征值对应的特征向量进行信号重构，便可提取出此频率成分。

根据式(8)所示的加法关系，利用PCA算法还可实现“点阻”滤波，在原始信号中简单地减去PCA算法提取到的某一频率成分，便可在原始信号中消去这一频率成分，实现“点阻”滤波，实施例七将通过实例进一步说明。

实施例六：

轴心轨迹的提纯：在旋转机械设备中，通常不同故障类型对应着不同形状的轴心轨迹，可根据轴心轨迹形状判断出系统对应的故障类型，比如正常情况下，轴心轨迹通常是长短轴相差不大的椭圆；椭圆形轴心轨迹通常对应着转动部件不平衡或主轴轴线不直引起的摆度过大的故障；规则或不规则的花瓣形轴心轨迹往往对应着动静件碰磨的故障；内“8”字形的轴心轨迹对应着油膜涡动的故障；香蕉形或外“8”字形的轴心轨迹对应的是不对中故障。

虽然轴心轨迹的形状是判断设备故障的重要依据，但大型旋转机械运行时工况相对复杂，其主轴位移信号受电磁干扰或高频噪声干扰严重，从未提纯的轴心轨迹中难以得到任何有效信息。对于如何提纯轴心轨迹，目前使用较多的方法有小波变换、小波包、谐波小波以及EEMD等方法。此处根据实施例三得出的结论，提出将PCA算法应用到轴心轨迹的提纯上。

实施例七：

通过大型转子振动测试试验台检测实验数据，为实时监测转子的工作情况，在试验台转子A/B两侧的轴颈处，每侧以相互垂直的方式安装两个Kaman KD2306-1S电涡流位移传感器，如图5和图6所示，采用LMS Test.Lab进行数据采集。

轴心轨迹特征的提取：

本实施例仅对A端面处的D11与D12号电涡流位移传感器采集的数据进行分析。试验台转速为1080r/min，采样频率为2048Hz，采样点为4096点，采用FFT变换滤掉直流分量后的时域信号和幅值谱如图7所示。

从图7(b)中可以看出原始信号不仅存在噪声的干扰，而且还存在工频50Hz及其倍频的干扰。信号的能量主要集中在前五倍频上，表1为A端面D11与D12前5倍频的幅值大小及其在原始信号幅值谱中排序统计。

表1A端面D11与D12前5倍频的幅值大小及排序

由图7(a)中的主轴位移信号D11与D12合成的轴心轨迹如图8所示，其中原始信号D11作为x轴，D12信号作为y轴。可见根本无法辨识出轴心轨迹的形状，进而无法得知设备运行状态。现在采用本发明进行轴心轨迹的提纯。

采用前二倍频进行轴心轨迹的提纯：

从图7(b)和表1中可以看出，信号的基频和二倍频的幅值较大，首先提取信号的基频与二倍频进行轴心轨迹的合成。D11与D12信号的协方差矩阵特征值分布如图9所示，采用PCA分离特征频率的方法，结合表1中D11与D12幅值排序，分别对图9中的第一组与第二组特征值(第1和第2个、第3和第4个特征值)对应的特征向量进行重构，得到基频和二倍频的频域图如图10所示。由图10可见，采用PCA算法滤波后的效果十分明显，不仅滤掉了噪声的干扰，而且滤除了工频50Hz及其谐波的干扰。同时还可看出滤波后的幅值与表1中原始信号的基频和二倍频的幅值十分接近，说明频谱泄漏可略为不计。

利用PCA提取到的基频和二倍频信号合成的轴心轨迹如图11所示，可以看出采用PCA滤波后合成的轴心轨迹为明显的内陷香蕉状，且凹陷程度较大，说明近电机端(A端)存在不对中故障，可见PCA算法在轴心轨迹的提纯上有着十分理想的效果。

采用前五倍频进行轴心轨迹的提纯：

由于转子的轴心轨迹一般是多个故障复合的结果，所以继续采用本发明的方法对原始信号前五倍频进行提取，并合成轴心轨迹。

根据表1中倍频大小在原始信号幅值谱中的排序进行信号重构，例如D11中的3倍频排序为6，因此采用第11和第12个特征值进行信号重构，5倍频排序为11，采用第21和第22个特征值进行信号重构。采用PCA算法提取D11与D12的前5倍频的频谱图，如图12所示。从图12可以看出滤波后的信号十分理想，不仅消除了噪声的干扰，而且还成功地去除了工频50Hz及其谐波的干扰，尽管工频50Hz与3倍频信号(54Hz)十分接近。对比图12与图7以及表1可看出，滤波后各倍频的幅值与表1中原始信号的幅值十分的接近。

利用PCA算法得到的前5倍频时域信号如图13所示，可见时域图十分纯净且幅值波动较小。

采用图13中的时域信号合成轴心轨迹如图14，可以看出滤波后的轴心轨迹为明显的不规则花瓣形，说明近电机端(A端)存在动静件碰磨的故障。

结合图11和图14说明近电机端(A端)既存在不对中的故障，同时也存在着动静件碰磨的故障。这表明大型旋转平台的转子一般存在多个故障的复合。

从上面的实例可见，本发明的方法对单个特征频率的提取上效果十分显著，而一般的带通滤波器是很难做到的，滤波器总是会存在一定的带宽，并且还有过渡带的影响，除非两个频率相隔较远，否则滤波器在分离指定的特征频率时总是会受到临近频率的影响。而采用本发明提出的PCA信号分离方法，即使两个频率相差1Hz也能够不受影响地将它们提取出来。

PCA算法与小波算法在轴心轨迹提纯上的对比：

(1)PCA算法与小波包对比

首先，将本发明与小波包算法处理的结果进行对比，由于小波分析具有很强的波形识别和剔除噪声功能，所以小波包变换在轴心轨迹的提纯上，是常用的方法之一。

但小波包变换是采用隔点采样，即隔二抽一采样。这虽然达到无冗余存储且可以重构信号的目的，但随着分解层数的增加，各层、各频段序列的采样频率减半，随着数据点数减少，信号细节部分的失真问题也会越明显。同时还应注意在频段的选取上应使指定的特征频率尽量不要落在频段边界上，以免出现波形畸变。例如本例中若将D11的信号分解到第8层，可将3倍频54Hz与工频干扰50Hz信号剥离，但又会面临二倍频36Hz，分解到第8频段(即小波包节点[8,5])32-36Hz上，则2倍频36Hz恰好落在此频段的边界上，不利于分析。

利用小波包变换提取D11中的倍频信号，原始信号的采样点数4096，采样频率为2048Hz，采用Daubechies5小波进行6层小波包分解，提取前4倍频信号如图15所示，D11的前5倍频特征如表1。

采用小波包变换提取原始信号的基频及其倍频信号时，从图15中发现滤波后的倍频信号周围存在着较多噪声频率的干扰，更有甚者，小波带通滤波器将特征频率附近的噪声干扰进行了放大，并且对比图15与表1中D11信号的倍频幅值，发现小波滤波后的倍频信号幅值减小很多，存在能量泄露的现象。

同时，由图15(c)可知三倍频信号54Hz与工频50Hz较为接近，若采用Daubechies5小波进行6层小波分解50Hz与54Hz均会落到48-64Hz的频段内，此时每段的有效点数为64个；若采用8层小波分解，可将50Hz与54Hz信号剥离，但又会面临二倍频36Hz在频段边界上，而且有效点数仅为16个，分离后的信号会严重失真。而采用PCA算法不仅能够提取出特征频率而且还能够避免带通滤波器频带内的噪声干扰，显然在轴心轨迹的提纯上，PCA算法要明显优于小波包算法。

(2)PCA算法与谐波小波对比

谐波小波算法是Newland提出的一种改进小波算法，是旋转试验台常用的信号处理方法之一，能够对整个频带进行无限细分，规避了二进小波与二进小波包的隔二抽一采样不足，在旋转试验台轴心轨迹的上，可以很方便的提取基频及其倍频信号。

但谐波小波算法同样也是一种带通滤波器，在提取特征频率的同时不可避免的提取出特征频率附近的噪声频率，而且在频带的选取上应注意尽量避免分析频率落在频段边界上，以免出现波形畸变。同样，若将D11的信号分解到第8层，可将3倍频54Hz与工频干扰50Hz信号剥离，但同时又会面临二倍频36Hz恰好落在频带32-36Hz的边界上的问题，对于这一问题谐波小波难以解决。

利用谐波小波算法提取D11中的基频与二倍频信号，原始信号的采样点数4096，采样频率为2048Hz，采用谐波小波进行7层小波分解，提取前2倍频信号如图16所示。

将图15中(a)和(b)图与图16中的(a)和(b)对比并参照表1可发现，谐波小波在提取特征频率时要比Daubechies小波包明显很多，且不存在能量泄露的问题。

由图16可看出谐波小波在提取特征频率的同时，特征频滤周围还存在许多噪声频率，这是谐波小波带通滤波器不可以避免的，若采用增加分解层数来减弱特征频率周围的噪声干扰，往往会造成分析频率落在变频带上，更加不利于分析。同时从图16中还可见，采用谐波小波提取特征频率时，并没有放大特征频率周围的噪声信号，这点明显优于Daubechies小波。

利用谐波小波算法提取到的基频和二倍频信号合成的轴心轨迹如图17所示，比较清楚地看出采用谐波小波滤波后合成的轴心轨迹为内陷的香蕉状，这与采用PCA算法得到结果一致，但图11得到的结果比图17合成的轴向轨迹要明显的多。这说明发明在轴心轨迹的提纯更加适用。

PCA算法实现单个频率滤除：

以上探讨了PCA算法在轴心轨迹方面的应用，探讨本发明在单个频率滤除方面的应用。以D11的实验数据为例，由于D11中含有50Hz的工频干扰，采用PCA算法滤除50Hz的工频干扰。由图7(b)可知，原始信号中的50Hz的幅值大小排序为3，根据本发明结论采用特征值分布图中的第5和第6个特征值对应的特征向量进行重构，PCA从原始信号中提取由特征值重构的特征信号时，由式(8)可知，提纯出的信号和原始信号是一种加减法的关系，故从D11的原始数据中减去重构后的信号得到滤除50Hz后的信号，其幅值谱如图18所示。

从图18与表1中可以看出，本算法不仅滤除了原始信号50Hz的工频干扰，而且对原始信号中其它频率成分也未产生影响，这表明本方法在实现单个频率的去除上有很高的应用价值。

信号协方差矩阵中有效特征值的选择是PCA算法研究中的一个重点，本发明利用了有效特征值与信号频率和幅值之间的关系，并将PCA算法应用到轴心轨迹的提纯上，得出以下结论：

(1)在Hankel矩阵方式下，则信号中每个有效的频率成分最多只产生两个非零的特征值，一前一后的紧密排列；信号有效特征值的个数与信号的频率个数相关，与f_i、A_i和的取值无关；信号有效特征值在特征值分布图中的排列顺序仅与信号频率对应的幅值A_i相关，与f_i和的取值无关，对应频率的幅值越大所产生的特征值也越大，排列也就越靠前。

(2)基于每个频率信号产生两个有效的特征值和频率对应的幅值大小决定特征值的排列顺序这种特性，提出了一种采用PCA提取指定的特征频谱(单个频率)方法，或消除指定特征频率的方法，指出即使两个频率相差1Hz也能够不受影响地提取单个频率。本发明利用这种特性将PCA算法应用到轴心轨迹的提纯上，并进行了小波包与谐波小波包提纯效果的对比，指出了小波包与谐波小波包变换存在的不足，显示出PCA算法在轴心轨迹的提纯上效果显著。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (7)

1.一种基于主成分分析的特征频率提取方法，其特征在于包括以下步骤：

首先，两个位移传感器以相互垂直的方式采用一定的采样频率采集试验台的转子位移数据；

第二步，将采集到的信号x(t)，利用快速傅里叶变换滤出原始信号的直流分量，然后将消除直流分量后的信号构造Hankel矩阵X，具体如下；

式中，X为m×n阶的矩阵，x_ij为X向量中的第i行，第j的元素，其中1≤i≤m，1≤j≤n，m为矩阵的行数，n为矩阵的列数；

同样地，将X表示成列向量的形式，即X＝[x₁，x₂…x_i…x_m]^T，x_i表示矩阵X中第i行，包含n个元素的行向量，T表示向量的转置；

当信号x(t)的数据长度为L时，L为偶数时，m＝L/2，n＝L/2+1；L为奇数时，m＝(L+1)/2，n＝(L+1)/2；

第三步，求矩阵X的协方差矩阵C，如下式，并求其特征值λ_i，按从大到小排列为λ₁,λ₂…λ_m，同时，构造协方差矩阵特征值分布图，并求出对应特征值的特征向量为α₁,α₂…α_m；

<mrow> <mi>C</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mrow></mrow> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>m</mi> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>m</mi> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cov</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>m</mi> </msub> <mo>,</mo> <msub> <mi>x</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

式中：cov(x_i,x_j)＝E[(x_i-E(x_i))(x_j-E(x_j))]；E(x_i)为求向量x_i的均值，cov()求协方差运算；C为协方差矩阵；x_i表示矩阵X中第i行，包含n个元素的行向量，1≤i≤m；

根据特征值λ_i在协方差矩阵特征值分布图中的排序，选择与某一具体频率对应的两个特征值及相应的特征向量进行重构，若某一频率在原始信号幅值谱中的幅值大小排序为k，则选择协方差矩阵特征值分布图中的第2k-1与第2k个特征值对应的特征向量α_i进行信号重构，得到一个新的矩阵如下所示；

<mrow> <msup> <mover> <mi>X</mi> <mo>^</mo> </mover> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mi>k</mi> </mrow> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>&amp;alpha;</mi> <mi>i</mi> </msub> <msubsup> <mi>&amp;alpha;</mi> <mi>i</mi> <mi>T</mi> </msubsup> <mi>X</mi> </mrow>

然后，将重构后的矩阵叠加原始信号的均值并获得矩阵

一般，从矩阵中恢复出信号有两种方法；

方法1：采用平均法恢复信号

<mrow> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>k</mi> </munderover> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>-</mo> <mi>i</mi> </mrow> </msub> <mo>/</mo> <mi>k</mi> <mo>,</mo> <mn>1</mn> <mo>&amp;le;</mo> <mi>k</mi> <mo>&lt;</mo> <mi>m</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>-</mo> <mi>i</mi> </mrow> </msub> <mo>/</mo> <mi>m</mi> <mo>,</mo> <mi>m</mi> <mo>&amp;le;</mo> <mi>k</mi> <mo>&lt;</mo> <mi>n</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mi>k</mi> <mo>-</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>-</mo> <mi>i</mi> </mrow> </msub> <mo>/</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>-</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>n</mi> <mo>&lt;</mo> <mi>k</mi> <mo>&amp;le;</mo> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

式中，为近似重构矩阵中第i行，第j列的元素；m为原始矩阵X的行数，n为列数；i，j，m，n，k均为正整数；

方法2：采用伪逆矩阵的方式恢复信号

<mrow> <mover> <mi>x</mi> <mo>^</mo> </mover> <mo>=</mo> <mover> <mi>X</mi> <mo>^</mo> </mover> <msup> <mi>M</mi> <mo>+</mo> </msup> </mrow>

式中，为恢复后的矩阵；为恢复的一维信号；M⁺称为M的伪逆；M为m个阶数为n的单位矩阵构成；其中，M⁺＝(M^TM)^-1M^T∈R^mn×l，l＝m+n-1，m，n，l为正整数；

信号就是从原始信号x(t)中提取的特征频率成分。

2.根据权利要求1所述的基于主成分分析的特征频率提取方法，其特征在于：选取工频干扰50Hz在特征值分布图中对应的两个特征值及对应的特征向量进行重构，然后从原始信号中减去重构后的特征信号，便滤除了50Hz工频干扰的信号。

3.根据权利要求1所述的基于主成分分析的特征频率提取方法，其特征在于：完成特征频率的提取后，进行轴心轨迹的合成，将第一个加速度传感器采集到的信号进行特征频率的提取，将处理后的信号作为x轴；然后，将处理后的第二个加速度传感器采集的信号作为y轴，进行二维轴心轨迹的合成。

4.根据权利要求3所述的基于主成分分析的特征频率提取方法，其特征在于：提取后的信号为基频、二倍频、三倍频或1/2基频的组合。

5.根据权利要求3所述的基于主成分分析的特征频率提取方法，其特征在于：合成三维轴心轨迹，即将t作为z轴，合成位移量与时间的三维轴心轨迹。

6.根据权利要求1所述的基于主成分分析的特征频率提取方法，其特征在于：滤除原始信号中的直流分量的方法采用平均值法或高通滤波器法。

7.根据权利要求1所述的基于主成分分析的特征频率提取方法，其特征在于：

设有m个随机向量x₁，x₂…x_m，每个向量x_i含有n个样本，即x_i＝(x_i1，x_i2…x_in)，则可构造一个m×n阶的矩阵X：

式2中，X为m×n阶的矩阵；x_ij为X向量中的第i行，第j的元素，其中1≤i≤m，1≤j≤n；m为矩阵的行数，n为矩阵的列数；

X可表示为X＝[x₁，x₂…x_m]^T，对矩阵X进行主成分分析，可得l个新的变量指标y_i，i＝1,2…l(l≤m)表达式如下：

式3中，y_i为l个新的变量指标，y_i∈R^1×n，其中1≤i≤l≤n；x_i为矩阵X的行向量；X为m×n阶的矩阵；α_i为X的协方差矩阵第i个特征值由大到小对应的特征向量，α_ij为特征向量α_i内的第j个元素，1≤j≤n；T为转置标识符；i，j为正整数；

其中y_i∈R^1×n，由主成分分析的定义可知α_i＝(α_i1,α_i2,…,α_im)^T为X的协方差矩阵第i个特征值由大到小对应的特征向量，且α_i满足：

式4中，α_ij为特征向量α_i内的第j个元素，1≤i≤n，1≤j≤n，i，j为正整数；

式2所示矩阵X的协方差矩阵为：

式5中：cov(x_i,x_j)＝E[(x_i-E(x_i))(x_j-E(x_j))]，E(x_i)为求向量的均值，cov()求协方差运算；C为协方差矩阵；x_i为第i个变量的n个样本组成的行向量，1≤i≤m；

由主成分分析原理可知，协方差矩阵C的特征方程如下：

Cα_i＝λ_iα_i 式6

式6中，C为协方差矩阵；α_i为协方差矩阵C的第i个征向量；λ_i是矩阵C的特征值；

式6中，λ_i是矩阵C的特征值，α_i是与λ_i对应的特征向量，把特征值按由大到小的顺序排列λ₁＞λ₂＞…＞λ_l。根据公式1和公式3便可以实现数据降维，由原来的m个变量转换成新的l个变量；若要进行信号处理，则根据公式3进行信号的重构，两边左乘α_i，并进行求和，得：

式7中，I为单位矩阵；y_i为降维后的矩阵；m为X向量的行数；α_i为协方差矩阵C的第i个特征向量，T为转置标识符；

若根据累积贡献率L_k选择前k个主成分进行重构，可得到一个近似的矩阵：

式8中，为恢复后的矩阵；y_i为降维后的矩阵；k为有效特征值的个数；α_i为协方差矩阵C的第i个特征向量；T为转置标识符；

与原矩阵X相比，近似重构的矩阵保留了原始矩阵的大部分信息，却消除了原始矩阵的冗余信息；

获得重构矩阵后按矩阵构成方式恢复信号即可，通过矩阵的逆变换从重构矩阵得到恢复后的信号即：

式9中，为恢复后的矩阵；为恢复的一维信号；M⁺称为M的伪逆；M为m个阶数为n的单位矩阵构成。其中，M⁺＝(M^TM)^-1M^T∈R^mn×l，l＝m+n-1，m，n，l为正整数。

式中，M⁺＝(M^TM)^-1M^T∈R^mn×l，M⁺称为M的伪逆。通过求伪逆矩阵恢复出信号实际就是对矩阵各副对角线上的元素求平均，这与采用平均法恢复信号的结果一致，由于矩阵M为稀疏矩阵，尤其是当m与n的取值较大时，计算伪逆矩阵M⁺时所需要的内存和时间会成倍的增加，因而本文采用平均法从矩阵中恢复信号表达式如下：

式10中，为近似重构矩阵中第i行，第j列的元素；m为原始矩阵X的行数，n为列数；i，j，m，n，k均为正整数；根据公式十便可以恢复出信号