CN111854930A - 一种基于先验预估的振动信号工频干扰压制方法 - Google Patents

Abstract

一种基于先验预估的振动信号工频干扰压制方法，涉及频谱信号分析领域，解决现有方法存在信噪频谱重叠时损害有用信号，工频频率不稳定时压噪效果不理想，以及不能满足使用需要等问题，首先对振动信号做FFT变换得到频谱图，并确定振动信号中工频干扰的频率及工频干扰频率幅值大小的排列，然后对振动信号构造的Hankel矩阵进行主成分分析，提取与工频干扰频率幅值大小的排列后对应的主成分进行重构后恢复矩阵得到先验工频干扰，再通过FFT确定的频率及先验工频干扰求得的振幅和初相位构造出预估工频干扰，最后用振动信号减去预估工频干扰即可得到去噪后振动信号。本发明方法对信号细节保护能力强，可有效提高振动信号的信噪比。

Description

一种基于先验预估的振动信号工频干扰压制方法

技术领域

本发明涉及频谱信号分析领域，具体涉及一种基于先验预估的振动信号工频干扰压制方法，该方法用于轴承、转子等旋转机械设备振动信号的工频干扰压制。

背景技术

振动测试技术是在轴承、转子等旋转机械设备上布置位移、速度或加速度等传感器获取旋转机械设备的振动信号，以此来监测旋转机械设备的工作状态。但由于工业电等因素的影响，使传感器采集到的振动信号中存在工频干扰，造成有用信号被工频干扰信号淹没，进而增加了旋转机械设备工作状态判断和故障特征提取的难度。同时严重影响了旋转机械设备的状态判断和故障诊断。

目前压制工频干扰信号的方法主要有频域法和时域法两大类。频域法主要有：

现有文献1《基于奇异向量频谱的工频干扰消除方法》，华南理工大学学报(自然科学版)，2017,45(5):16-23，该方法运用奇异值向量的频谱去除工频干扰，虽然该方法在时域上可有效滤除大量噪声，但通过对压噪后的信号进行频谱分析发现，该方法在信噪频谱混叠时易损害有用信号。

现有文献2《基于谐波小波和SVD的可倾瓦滑动轴承工频干扰抑制》，华南理工大学学报(自然科学版),2016,44(10):1-7，该方法运用谐波小波和奇异值分解联合的方法压制工频干扰，但该方法在工频频率不稳定时压噪效果不理想。

现有文献3《基于独立分量分析的振动信号中工频干扰消除》，昆明理工大学学报(理工版)，2009,34(3):24-28，该方法通过独立分量分析法压制工频干扰，但该方法受限于观测信号源数目及信号的方向性。

现有文献4《滚动轴承磨损区域静电监测技术及寿命预测方法研究》，南京航空航天大学,2013，该方法中运用陷波法压制工频干扰，但该方法不适应于信噪频谱混叠的情况。

现有文献5《基于谱插值和奇异值差分谱的滚动轴承静电监测信号去噪方法》，航空动力学报,2014,28(9):1996-2002，该方法运用谱插值法压制工频干扰，但该方法同样在信噪频谱混叠区域易滤除部分有用信号。

现有文献6《余弦函数自适应法识别与消除单频干扰》，石油地球物理勘探,2011,46(1):64-69，通过正余弦函数逼近法消除工频干扰，但该方法受限于工频频率的稳定性。

时域法主要有：

现有文献7《基于PCA最优阶数的并行震源工频噪声压制》，仪器仪表学报,2018,39(1):15-23，通过分析主成分阶数与工频干扰频率之间的内在关系去除工频噪声，但该方法同样存在信噪频率混叠时损害有用信号的现象。

现有专利文件1，公开号为CN106908840A，名称为《基于主成分分析的地震资料工频干扰自动识别与压制方法》通过主成分分析方法压制工频干扰，但该方法只在信噪频谱重叠不严重时压噪效果理想。

现有专利文件2，公开号为CN107219555A，名称为《基于主成分分析的并行震源地震勘探资料强工频噪声压制方法》同样受限于信噪频谱混叠的情况。

现有专利文件3，公开号为CN106895985A，名称为《高速线材轧机的故障信号降噪重构特征识别方法》通过傅里叶变换搜索表征工频干扰奇异值的方法压制工频干扰，但该方法只适合于强能量工频干扰的环境。

综上，上述方法都不适应于振动信号工频干扰的压制。

发明内容

本发明为解决现有对工频干扰信号的压制方法存在信噪频谱重叠时损害有用信号，工频频率不稳定时压噪效果不理想，以及采用干扰奇异值的方法压制工频干扰不能满足使用需要等问题，提供一种基于先验预估的振动信号工频干扰压制方法。

一种基于先验预估的振动信号工频干扰压制方法，该方法由以下步骤实现：

步骤一、对振动信号进行采集并进行主成分分析，获得所述振动信号的主成分矩阵；并进行主成分重构，获得主成分重构矩阵Y：

步骤二、将步骤一获得的主成分重构矩阵Y的第一行和最后一列首尾相连，获得先验工频干扰x′(i)；

步骤三、设定步骤二获得的先验工频干扰x′(i)中元素的最大值为A_max，最小值为A_min，则预估工频干扰信号y(i)的振幅A为：

A＝(A_max+A_min)/2

步骤四、构造预估工频干扰信号y₁(i)的正弦函数，构造表达式为：

y₁(i)＝A·sin(2·π·f·i·Δt)

式中Δt为采样间隔，i为振动信号x(i)中采样点序号，f为振动信号x(i)中工频干扰的频率；

步骤五、以先验工频干扰x′(i)为参考信号，计算先验工频干扰x′(i)和正弦函数y₁(i)的互相关函数为：

式中，N为采样点数，

为互相关函数，j为偏移量，j＝0,1,…,L-1，

为x′(i)的平均值，

为y₁(i)的平均值；

步骤六、将步骤五中的互相关函数

中最大值对应的偏移量作为校正偏移量j^max；所述预估工频干扰信号y₁(i)的表达式则为：

y(i)＝A·sin(2·π·f·(i+j^max)·Δt))

步骤七、求取振动信号x(i)与步骤六获得的预估工频干扰信号y₁(i)的差值，获得压制工频干扰后的信号。

本发明的有益效果：本发明所述的振动信号工频干扰压制方法，首先对振动信号做FFT变换得到频谱图，根据频谱图确定振动信号中工频干扰的频率及工频干扰频率幅值大小的排列顺序，然后对振动信号构造的Hankel矩阵进行主成分分析，提取与工频干扰频率幅值大小的排列顺序对应的主成分进行重构后恢复矩阵得到先验工频干扰，再通过FFT确定的频率及先验工频干扰求得的振幅和初相位构造出预估工频干扰，最后用振动信号减去预估工频干扰即可得到去噪后振动信号。从而实现了振动信号工频干扰的压制。

本发明所述的振动信号工频干扰压制方法，能够实现在振动信号中压制工频干扰的目的，该算法可有效滤除工频干扰，同时在信噪频谱混叠时进行压噪不会损害有用信号，可有效提高振动信号信噪比。

附图说明

图1为本发明所述的一种基于先验预估的振动信号工频干扰压制方法的流程图；

图2为本发明所述的一种基于先验预估的振动信号工频干扰压制方法的振动信号频谱图；

图3为压噪前振动信号时域图；

图4为压噪前后振动信号时域对比图。

具体实施方式

结合图1至图4说明本实施方式，基于先验预估的振动信号工频干扰压制方法，在本实施方式中，振动信号包含变频信号和工频干扰,其中变频信号：sin(2·π·i·Δt+42·(i·Δt)²)，工频干扰：0.8·sin(2·π·51·i·Δt+1.34),采样间隔Δt＝0.001,i＝1,2,…,N，N＝1201，振动信号：x(i)＝sin(2·π·i·Δt+42·(i·Δt)²)+0_.8·sin(2·π·51·i·Δt+1_.34)。

步骤一、对振动信号进行采集并进行主成分分析，获得所述振动信号的主成分矩阵；并进行主成分重构，获得主成分重构矩阵Y：具体过程为：

A、设振动信号为x(i)，i＝1,2,...,N，i为采样点序号，N为采样点数，对x(i)进行快速傅里叶变换(FFT)得到频谱图，根据频谱图确定x(i)中工频干扰的频率f及工频干扰频率幅值大小的次序k，k＝1,2,...,N/2,本例中N＝1201,f＝51,k＝1；

B、根据相空间理论对x(i)，构造Hankel矩阵如下：

记矩阵H的行数为m，列数为n，若N为偶数，则m＝(N/2)+1，n＝N/2，若N为奇数，则m＝n＝(N+1)/2，本例中N为奇数，m＝n＝601；

C、计算H的协方差矩阵Γ，如公式：

其中H^T为H的转置矩阵，“·”表示矩阵乘法；

D、用奇异值分解法，计算协方差矩阵Γ的特征值矩阵Λ和特征向量矩阵R则存在公式：

Γ＝R·Λ·R^T (3)

其中Λ为特征值由大到小排列的对角矩阵，Λ＝diag[λ₁,λ₂,…,n]，λ₁,λ₂，λ₃，...，λ_n为特征值，R为各个特征值对应的特征向量矩阵，R^T为R的转置矩阵且满R^T·R＝R·R^T＝E，其中E为单位矩阵。

E、H经线性映射，得到主成分矩阵：

Φ＝R^T·H (4)

F、保留Φ的第2k-1，2k行主成分，其余各行置零后得到矩阵为Φ′，并进行主成分重构，则得到重构矩阵：

Y＝R·Φ′ (5)

本例中保留Φ的1、2行主成分。

步骤二、Y的第一行和最后一列首尾相连得到先验工频干扰x′(i)；

步骤三、x′(i)中大于零的元素中元素最大值为A_max，最小值为A_min，则预估工频干扰y(i)的振幅A为：

A＝(A_max+A_min)/2 (6)

步骤四、构造预估工频干扰y(i)正弦函数，构造表达式为：

y₁(i)＝A·sin(2·π·f·i·Δt) (7)

其中Δt为采样间隔，“·”表示相乘；本例中Δt＝0.001,i＝1,2,…,1201。

步骤五、以x′(i)为参考信号，计算x′(i)和y₁(i)的互相关函数为：

其中

为互相关函数，j为偏移量，j＝0,1,…,N-1，

为x′(i)的平均值，

为y₁(i)的平均值；本例中，

步骤六、将

中最大值对应的偏移量作为校正偏移量j^max，本例中j^max＝3；m、预估工频干扰y₁(i)的表达式为：

y(i)＝A·sin(2·π·f·(i+j^max)·Δt)) (9)

本例中y(i)＝0.8369·sin(2·π·51·((i+3)·0.001))

步骤七、压制工频干扰后的信号为：

x₁(i)＝x(i)-y(i) (10)

结合图2和图4说明本实施方式，图2为振动信号频谱图，基中工频干扰频率为51Hz，工频干扰频率幅值排第1。图3为压噪前振动信号时域图；图4为压噪前后振动信号时域对比图。可以看出，图3中的振动信号中包含大量工频干扰，而经图4的对比可以看出，压噪后振动信号变得平滑，清晰。

Claims (3)

1.一种基于先验预估的振动信号工频干扰压制方法，其特征是：该方法由以下步骤实现：

步骤一、对振动信号进行采集并进行主成分分析，获得所述振动信号的主成分矩阵；并进行主成分重构，获得主成分重构矩阵Y：

步骤二、将步骤一获得的主成分重构矩阵Y的第一行和最后一列首尾相连，获得先验工频干扰x′(i)；

步骤三、设定步骤二获得的先验工频干扰x′(i)中元素的最大值为A_max，最小值为A_min，则预估工频干扰信号y(i)的振幅A为：

A＝(A_max+A_min)/2

步骤四、构造预估工频干扰信号y₁(i)的正弦函数，构造表达式为：

y₁(i)＝A·sin(2·π·f·i·Δt)

式中Δt为采样间隔，i为振动信号x(i)中采样点序号，f为振动信号x(i)中工频干扰的频率；

步骤五、以先验工频干扰x′(i)为参考信号，计算先验工频干扰x′(i)和正弦函数y₁(i)的互相关函数为：

式中，N为采样点数，

为互相关函数，j为偏移量，j＝0,1,…,L-1，

为x′(i)的平均值，

为y₁(i)的平均值；

步骤六、将步骤五中的互相关函数

中最大值对应的偏移量作为校正偏移量j^max；所述预估工频干扰信号y₁(i)的表达式则为：

y(i)＝A·sin(2·π·f·(i+j^max)·Δt))

步骤七、求取振动信号x(i)与步骤六获得的预估工频干扰信号y₁(i)的差值，获得压制工频干扰后的信号。

2.根据权利要求1所述的一种基于先验预估的振动信号工频干扰压制方法，其特征在于：步骤一中获得振动信号的主成分矩阵的过程为：

步骤一一、设定振动信号为x(i)，i＝1,2,...,N，对x(i)进行FFT变换得到频谱图，根据频谱图确定x(i)中工频干扰的频率f及工频干扰频率幅值大小的次序k，k＝1,2,...,N/2；

步骤一二、根据相空间理论对x(i)，构造Hankel矩阵H如下：

所述H的行数为m，列数为n，若N为偶数，则m＝(N/2)+1，n＝N/2，若N为奇数，则m＝n＝(N+1)/2；

步骤一三、计算矩阵H的协方差矩阵Γ，用公式表示为：

式中，H^T为H的转置矩阵，·为矩阵乘法；

步骤一四、采用奇异值分解法，计算步骤一三中的协方差矩阵Γ的特征值矩阵Λ和特征向量矩阵R，则有公式：

Γ＝R·Λ·R^T

式中，Λ为特征值由大到小排列的对角矩阵，Λ＝diag[λ₁,λ₂,…,λ_n]，λ₁,λ₂,λ₃,…,λ_n为特征值，R为各个特征值对应的特征向量矩阵，R^T为R的转置矩阵且满R^T·R＝R·R^T＝E，所述E为单位矩阵；

步骤一五、H经线性映射，获得主成分矩阵：

Φ＝R^T·H。

3.根据权利要求1所述的一种基于先验预估的振动信号工频干扰压制方法，其特征在于：步骤一中，保留主成分矩阵Φ的第2k-1，2k行主成分，其余各行置零后的矩阵为Φ′，并进行主成分重构，则获得重构矩阵Y：

Y＝R·Φ′。

Images