别再只盯着EMD了！用Python手把手实现LMD（局部均值分解）处理轴承故障信号

别再只盯着EMD了！用Python手把手实现LMD处理轴承故障信号

轴承故障诊断领域，信号分解方法一直是研究热点。经验模态分解（EMD）作为经典算法，虽然应用广泛，但在处理非平稳信号时存在端点效应、模态混叠等问题。局部均值分解（LMD）作为一种新兴的自适应信号处理方法，通过乘积函数（PF）的分解方式，能更准确地提取信号特征。本文将用Python完整实现LMD算法，并展示其在轴承故障诊断中的实际应用效果。

1. 为什么LMD更适合轴承故障诊断

轴承振动信号具有明显的非平稳特性，传统傅里叶变换难以有效分析。EMD方法虽然能自适应分解信号，但在实际应用中存在几个关键问题：

• 端点效应：信号两端会出现发散现象
• 虚假分量：产生无物理意义的模态分量
• 包络问题：过包络或欠包络导致分解不准确

LMD通过局部均值和包络估计的迭代计算，将信号分解为若干乘积函数（PF）之和。每个PF分量都是包络信号与纯调频信号的乘积，具有更明确的物理意义。

关键优势对比：

2. LMD算法核心原理与实现步骤

LMD的核心思想是通过迭代过程将信号分解为多个PF分量。让我们深入理解其数学原理和实现步骤。

2.1 算法数学基础

给定信号x(t)，LMD分解过程可表示为：

x(t) = ΣPFₚ(t) + uₖ(t)

其中PFₚ(t)是第p个乘积函数，uₖ(t)是残余分量。

每个PF分量又可表示为：

PF(t) = a(t) × s(t)

a(t)是包络函数，s(t)是纯调频信号。

2.2 Python实现关键步骤

以下是LMD分解的核心步骤实现：

import numpy as np from scipy.signal import argrelextrema from scipy.interpolate import interp1d def lmd(signal, max_sifting=10, tol=0.001): """ 局部均值分解(LMD)实现 :param signal: 输入信号 :param max_sifting: 最大筛分次数 :param tol: 收敛容差 :return: PF分量列表 """ residual = signal.copy() pfs = [] while True: # 步骤1：提取局部极值 maxima_indices = argrelextrema(residual, np.greater)[0] minima_indices = argrelextrema(residual, np.less)[0] if len(maxima_indices) < 2 or len(minima_indices) < 2: break # 步骤2：计算局部均值和包络估计 mean_env = calculate_mean_envelope(residual, maxima_indices, minima_indices) # 步骤3-6：迭代筛分过程 pf, residual = sifting_process(residual, mean_env, max_sifting, tol) pfs.append(pf) return pfs, residual

完整实现还需要以下几个辅助函数：

def calculate_mean_envelope(signal, max_indices, min_indices): """计算局部均值和包络估计""" # 合并极值点并排序 extrema = np.sort(np.concatenate([max_indices, min_indices])) # 计算相邻极值点间的平均值和幅值 means = (signal[extrema[:-1]] + signal[extrema[1:]]) / 2 amps = np.abs(signal[extrema[:-1]] - signal[extrema[1:]]) / 2 # 创建插值函数 t = np.arange(len(signal)) mean_func = interp1d(extrema[:-1], means, kind='linear', fill_value='extrapolate') amp_func = interp1d(extrema[:-1], amps, kind='linear', fill_value='extrapolate') return mean_func(t), amp_func(t)

3. 完整LMD实现与轴承信号分析

现在我们将上述代码整合成完整的LMD实现，并应用于轴承故障信号分析。

3.1 完整LMD类实现

class LMDAnalyzer: def __init__(self, max_sifting=10, tol=0.001, max_pfs=8): self.max_sifting = max_sifting self.tol = tol self.max_pfs = max_pfs def decompose(self, signal): """完整LMD分解流程""" residual = signal.copy() pfs = [] for _ in range(self.max_pfs): pf, residual = self._extract_pf(residual) if pf is None: break pfs.append(pf) return pfs, residual def _extract_pf(self, signal): """提取单个PF分量""" h = signal.copy() a = np.ones_like(signal) for _ in range(self.max_sifting): # 获取极值点 max_idx = argrelextrema(h, np.greater)[0] min_idx = argrelextrema(h, np.less)[0] if len(max_idx) < 2 or len(min_idx) < 2: return None, signal # 计算局部均值和包络 mean, env = self._calculate_local_mean_env(h, max_idx, min_idx) # 更新h和a h = h - mean a = a * env # 检查是否纯调频 if np.all(np.abs(h) <= 1 + self.tol) and \ np.all(np.abs(env - 1) < self.tol): break pf = a * h residual = signal - pf return pf, residual def _calculate_local_mean_env(self, signal, max_idx, min_idx): """计算局部均值和包络""" extrema = np.sort(np.concatenate([max_idx, min_idx])) t = np.arange(len(signal)) # 计算局部均值 means = (signal[extrema[:-1]] + signal[extrema[1:]]) / 2 mean_func = interp1d(extrema[:-1], means, kind='linear', fill_value='extrapolate') # 计算包络估计 amps = np.abs(signal[extrema[:-1]] - signal[extrema[1:]]) / 2 amp_func = interp1d(extrema[:-1], amps, kind='linear', fill_value='extrapolate') return mean_func(t), amp_func(t)

3.2 轴承故障信号分析实例

让我们用实际轴承数据测试我们的LMD实现：

import matplotlib.pyplot as plt from scipy.io import loadmat # 加载轴承故障数据 data = loadmat('bearing_fault.mat') vibration = data['vibration'].flatten() fs = 12000 # 采样频率12kHz # LMD分解 lmd_analyzer = LMDAnalyzer(max_sifting=15, tol=0.005) pfs, residual = lmd_analyzer.decompose(vibration) # 可视化结果 plt.figure(figsize=(12, 8)) plt.subplot(len(pfs)+1, 1, 1) plt.plot(vibration) plt.title('原始振动信号') for i, pf in enumerate(pfs): plt.subplot(len(pfs)+1, 1, i+2) plt.plot(pf) plt.title(f'PF {i+1}') plt.tight_layout() plt.show()

典型分析结果特征：

1. 外圈故障：PF1分量中会出现明显的周期性冲击
2. 内圈故障：高频PF分量中可见调制现象
3. 滚动体故障：多个PF分量均可能出现特征频率

4. LMD与EMD实际效果对比

为了直观展示LMD的优势，我们使用同一段轴承故障信号对比两种方法的分解结果。

4.1 分解效果对比

from PyEMD import EMD # EMD分解 emd = EMD() imfs = emd.emd(vibration) # 绘制对比图 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10)) # EMD结果 for i, imf in enumerate(imfs[:5]): ax1.plot(imf + i*0.5, label=f'IMF {i+1}') ax1.legend() ax1.set_title('EMD分解结果') # LMD结果 for i, pf in enumerate(pfs[:5]): ax2.plot(pf + i*0.5, label=f'PF {i+1}') ax2.legend() ax2.set_title('LMD分解结果') plt.tight_layout() plt.show()

4.2 关键指标对比分析

我们通过几个量化指标来客观比较两种方法：

端点效应评估：

def evaluate_edge_effect(signal, components): """评估端点效应""" edge_ratio = 0.1 # 评估两端10%的区域 n = len(signal) edge_indices = int(n * edge_ratio) edge_energies = [] center_energies = [] for comp in components: edge_energy = np.sum(comp[:edge_indices]**2) + np.sum(comp[-edge_indices:]**2) center_energy = np.sum(comp[edge_indices:-edge_indices]**2) edge_energies.append(edge_energy) center_energies.append(center_energy) return np.mean(np.array(edge_energies) / np.array(center_energies))

计算效率对比：

import time # EMD计算时间 start = time.time() imfs = emd.emd(vibration) emd_time = time.time() - start # LMD计算时间 start = time.time() pfs, _ = lmd_analyzer.decompose(vibration) lmd_time = time.time() - start print(f"EMD计算时间: {emd_time:.3f}秒") print(f"LMD计算时间: {lmd_time:.3f}秒")

实际测试结果：

5. 工程应用中的优化技巧

在实际工程应用中，我们发现以下几个技巧可以显著提升LMD的分析效果：

5.1 参数调优指南

1. 最大筛分次数：

   • 一般设置在10-20次之间
   • 过小可能导致分解不充分
   • 过大会增加计算量，可能引入噪声

2. 收敛容差：

   • 典型值0.001-0.01
   • 信号幅值较大时可适当放宽
   • 对微弱特征检测应设置更严格

3. 极值点处理：

   • 对噪声较大的信号可先进行平滑处理
   • 极值点过少时可适当延长信号

5.2 常见问题解决方案

问题1：分解得到的PF分量过多或过少

解决方案：调整max_pfs参数，或检查信号是否包含太多噪声

问题2：端点效应仍然明显

解决方案：对信号进行镜像延拓后再分解

问题3：计算时间过长

优化代码：使用numba加速关键循环

from numba import jit @jit(nopython=True) def fast_local_mean(signal, extrema): # 使用numba加速的局部均值计算 means = np.zeros(len(extrema)-1) for i in range(len(extrema)-1): means[i] = (signal[extrema[i]] + signal[extrema[i+1]]) / 2 return means

5.3 轴承故障诊断完整流程

基于LMD的轴承故障诊断典型流程：

1. 数据采集：获取振动信号，确保足够采样率
2. 信号预处理：
   • 去直流分量
   • 必要时带通滤波
3. LMD分解：得到PF分量
4. 特征提取：
   • 包络谱分析
   • 能量熵计算
5. 故障识别：根据特征频率判断故障类型

包络谱分析示例代码：

from scipy.signal import hilbert def envelope_spectrum(pf, fs): """计算PF分量的包络谱""" analytic_signal = hilbert(pf) env = np.abs(analytic_signal) env_fft = np.abs(np.fft.rfft(env)) freqs = np.fft.rfftfreq(len(env), 1/fs) return freqs, env_fft # 分析PF1的包络谱 freqs, spec = envelope_spectrum(pfs[0], fs) plt.plot(freqs, spec) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude') plt.title('Envelope Spectrum of PF1') plt.show()

在实际项目中，LMD算法展现出了比传统EMD更稳定的性能。特别是在处理高速轴承的振动信号时，PF分量能更清晰地揭示故障特征频率，减少了人工干预的需要。一个实用的建议是：对于不同型号的轴承，可以先采集正常状态下的信号作为基准，这样故障诊断时对比效果会更明显。