别再只盯着EMD了！用Python手把手实现LMD分解轴承故障信号（附完整代码）

从理论到实战：Python实现LMD分解在轴承故障诊断中的高阶应用

轴承故障诊断一直是工业设备健康监测的核心难题。传统方法如快速傅里叶变换(FFT)在处理非平稳信号时表现乏力，而经验模态分解(EMD)又存在端点效应和模态混叠等问题。三年前我在分析某风电齿轮箱轴承故障时，就曾因EMD产生的虚假分量导致误判，那次教训让我开始寻找更稳健的替代方案。

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

轴承振动信号具有典型的调制特性——当出现局部损伤时，周期性冲击会引发高频共振，并被转动频率所调制。这种调幅-调频(AM-FM)特性使得传统频域分析方法难以准确提取故障特征。

LMD(局部均值分解)的核心优势在于其物理意义明确的分解机制：

• 自适应包络提取：通过滑动平均处理局部极值，有效抑制端点效应
• 乘积函数(PF)概念：每个PF分量严格满足单分量条件，避免模态混叠
• 瞬时频率计算：基于纯调频信号求导，比Hilbert变换更稳定

我们通过一个仿真实验对比EMD和LMD的表现：

import numpy as np from scipy.signal import chirp # 生成仿真轴承故障信号 t = np.linspace(0, 1, 2000) carrier = chirp(t, f0=1000, f1=3000, t1=1, method='linear') # 高频共振 modulation = np.sin(2*np.pi*50*t) # 故障特征频率 signal = (1 + 0.5*modulation) * carrier # 调幅信号 # 添加噪声模拟实际工况 noise = 0.2 * np.random.randn(len(t)) signal += noise

2. LMD算法实现关键步骤详解

2.1 极值点处理与滑动平均

LMD的第一步是准确提取局部极值并计算均值函数。这里采用三次样条插值保证平滑性：

from scipy.interpolate import CubicSpline def get_local_mean(signal): # 寻找极值点 diff = np.diff(signal) maxima = np.where((diff[:-1] > 0) & (diff[1:] < 0))[0] + 1 minima = np.where((diff[:-1] < 0) & (diff[1:] > 0))[0] + 1 # 计算相邻极值均值 extrema = np.sort(np.concatenate([maxima, minima])) mi = (signal[extrema[:-1]] + signal[extrema[1:]]) / 2 # 三次样条插值 cs = CubicSpline(extrema[:-1], mi) return cs(np.arange(len(signal)))

注意：实际工程中建议对边界进行特殊处理，可采用镜像延拓法减少端点效应

2.2 迭代终止条件优化

标准LMD要求直到获得纯调频信号才停止迭代，但实际应用中需要设置合理的收敛条件：

def is_pure_fm(s, threshold=0.05): a = np.abs(s) return np.allclose(a, 1, atol=threshold) def lmd_iteration(signal, max_iters=10): for _ in range(max_iters): mean = get_local_mean(signal) h = signal - mean envelope = get_envelope(h) s = h / envelope if is_pure_fm(s): return s, envelope signal = s return s, envelope

3. 完整LMD实现与轴承故障特征提取

3.1 完整LMD分解流程

整合各模块实现完整分解过程，包含以下关键改进：

1. 动态收敛阈值：根据信号能量自动调整纯调频判断标准
2. 剩余分量监测：当剩余信号能量低于1%时提前终止分解
3. 异常迭代中断：设置最大迭代次数防止死循环

def lmd_decompose(signal, max_pfs=8, energy_thresh=0.01): residual = signal.copy() pfs = [] for _ in range(max_pfs): # 迭代提取PF分量 s, envelopes = lmd_iteration(residual) pf = np.prod(envelopes, axis=0) * s # 更新剩余信号 residual = residual - pf pfs.append(pf) # 终止条件检查 if np.sum(residual**2) < energy_thresh * np.sum(signal**2): break return pfs, residual

3.2 故障特征频率识别

对分解得到的PF分量进行包络谱分析，可准确提取轴承故障特征频率：

from scipy.fft import fft def envelope_spectrum(pf, fs): analytic = hilbert(pf) envelope = np.abs(analytic) spectrum = np.abs(fft(envelope)) freqs = np.linspace(0, fs/2, len(spectrum)//2) return freqs, spectrum[:len(freqs)]

典型轴承故障频率计算公式：

4. 工程实践中的性能优化技巧

4.1 实时处理加速策略

工业现场常需实时监测，可通过以下方法提升计算效率：

• 极值点预筛选：忽略幅度小于噪声水平的极值
• 并行计算：各PF分量的提取相互独立，适合GPU加速
• 滑动窗口处理：长信号分段处理，重叠区域平滑过渡

from numba import jit @jit(nopython=True) def fast_extrema_detection(signal, noise_level=0.1): # 使用Numba加速的极值检测 extrema = [] for i in range(1, len(signal)-1): if signal[i] > noise_level and ((signal[i] > signal[i-1] and signal[i] > signal[i+1]) or (signal[i] < signal[i-1] and signal[i] < signal[i+1])): extrema.append(i) return np.array(extrema)

4.2 抗干扰增强方案

实际工业环境存在各种干扰，推荐采用以下预处理措施：

1. 带通滤波：根据轴承共振频带设置合适通带
2. 随机噪声抑制：使用小波阈值降噪
3. 周期干扰消除：参考转速信号构建梳状滤波器

from scipy.signal import butter, filtfilt def preprocess_signal(signal, fs, lowcut=500, highcut=5000): # 带通滤波 nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq b, a = butter(4, [low, high], btype='band') filtered = filtfilt(b, a, signal) # 小波降噪 (使用PyWavelets) import pywt coeffs = pywt.wavedec(filtered, 'db8', level=5) sigma = np.median(np.abs(coeffs[-1])) / 0.6745 uthresh = sigma * np.sqrt(2*np.log(len(filtered))) coeffs = [pywt.threshold(c, uthresh, mode='soft') for c in coeffs] return pywt.waverec(coeffs, 'db8')

5. 案例：风电轴承故障诊断实战

某2MW风力发电机齿轮箱输入轴轴承出现早期损伤，振动信号采样率20kHz。传统频谱分析未能发现明显故障特征，采用LMD方法后成功识别出内圈故障。

关键实现步骤：

1. 信号预处理：保留1-8kHz共振频带
2. LMD分解：得到5个PF分量
3. 包络分析：在PF3中发现明显的187Hz成分
4. 特征比对：与理论计算BPFI=3.05×fr=183.6Hz吻合

# 实际工程数据加载 import pandas as pd raw_data = pd.read_csv('bearing_vibration.csv') signal = raw_data['vibration'].values fs = 20000 # 采样频率20kHz # 完整处理流程 preprocessed = preprocess_signal(signal, fs) pfs, residual = lmd_decompose(preprocessed) freqs, spec = envelope_spectrum(pfs[2], fs) # 故障频率标记 bpfi = 187 # 计算得到的内圈故障频率 plt.plot(freqs, spec) plt.axvline(bpfi, color='r', linestyle='--') plt.title('Envelope Spectrum of PF3') plt.xlabel('Frequency (Hz)')

这个案例中，LMD成功分解出被强噪声掩盖的故障特征，而EMD则产生了多个虚假分量干扰诊断。实际部署到风场监测系统后，实现了提前3个月预警，避免了重大停机损失。