物理机理嵌入和自适应学习的机械早期故障诊断（Python）

首先基于轴承故障频率BPFO及谐波的理论知识设计多频带阻滤波器从健康帧中提取纯净的背景噪声轮廓随后采用 Wiener 滤波对原始振动信号进行去噪在保留故障特征的同时大幅降低噪声水平。去噪后的信号逐帧计算标准差STD经标准化后输入 LSTM 网络利用其时间序列建模能力学习健康状态下 STD 的演变规律。在测试阶段LSTM 预测的 STD 与真实值进行比较当预测值超过健康段统计阈值均值1.5倍标准差时触发早期预警。实验结果显示该方法能在帧 458 附近检测到异常比传统时域指标提前约 1.5 分钟显著提升了故障检测的及时性。算法充分结合了物理先验BPFO 谐波抑制与数据驱动LSTM 时序学习的优势具有计算高效、可解释性强、适合边缘部署的特点为旋转机械的预测性维护提供了可靠的技术方案。算法步骤物理机理噪声估计根据轴承参数滚子数、节径、接触角等计算理论故障频率 BPFO ≈ 236 Hz。选取若干健康帧如帧 100–105对每个健康帧应用多频带阻 Butterworth 滤波器精确抑制 BPFO 及其整数倍谐波直至 10 kHz尤其对 900–1100 Hz 共振区域采用更宽阻带以覆盖实测强峰。对滤除谐波后的信号进行高通滤波100 Hz剔除极低频机械干扰。将多个健康帧的噪声估计取平均获得稳定的背景噪声轮廓。Wiener 滤波去噪将噪声轮廓平铺至与原始 16 分钟信号等长作为噪声参考信号。使用 Welch 法估计原始信号和噪声的功率谱密度PSD。计算频域 Wiener 增益G PSD_signal / (PSD_signal PSD_noise)并设置增益下限如 0.1防止过度抑制。将增益插值到原始信号的 FFT 频率点对频谱加权后反变换得到去噪信号。将去噪信号重塑为 [984 帧 × 20480 样本] 的数组供后续特征提取和建模。时域特征提取对每一帧去噪信号计算标准差STD形成 984 点的 STD 序列。对 STD 序列进行标准化Z-score使其均值为 0、方差为 1便于 LSTM 训练。LSTM 预测模型构建与训练取健康段帧 50–250的标准化 STD 作为训练数据。使用滑动窗口窗口大小 10 帧构造输入-输出对用前 10 帧预测第 11 帧的 STD 值。搭建单层 LSTM50 单元 全连接输出层的回归模型以均方误差为损失函数Adam 优化器训练 100 轮并划分 20% 验证集。异常检测与早期预警对测试段帧 300–800应用训练好的 LSTM 模型逐窗口预测 STD 值。计算健康段 STD 的均值与标准差设定预警阈值 均值 1.5 × 标准差。比较模型预测值与实际值当预测值超过阈值时认为出现异常记录首次超阈值的帧号作为早期故障预警点。实验表明该方法在帧 458 附近即触发预警比传统指标峭度、RMS的显著变化约帧 540提前约 1.5 分钟。import numpy as np from scipy import signal from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Input from sklearn.preprocessing import StandardScaler # 从健康帧例如 frame 100,101,102,103,105提取噪声并平均 healthy_frames [100, 101, 102, 103, 105] noise_profiles [] for idx in healthy_frames: frame data_array[idx, 0, :] # 取通道 1轴承 1 frame_stop multi_bandstop_butterworth(frame) # 抑制 BPFO 谐波 # 进一步高通滤波100 Hz去除极低频机械干扰 frame_highpass highpass_butterworth(frame_stop, fs20480.0, cutoff100.0) noise_profiles.append(frame_highpass) # 对多个健康帧的噪声估计取平均得到稳定的背景噪声轮廓 noise_avg np.mean(noise_profiles, axis0) # 2. Wiener 滤波去噪 def wiener_filter(signal_in, noise_in, fs20480.0): 基于功率谱密度的 Wiener 滤波器。 原理在频域计算增益 G P_signal / (P_signal P_noise)对信号频谱加权 从而在抑制噪声的同时保留信号特征。 参数 signal_in : 待去噪的原始信号长序列如整个 16 分钟数据 noise_in : 噪声估计信号应与 signal_in 等长或通过平铺实现 fs : 采样率 返回 filtered_signal : 去噪后的信号 # 计算信号和噪声的功率谱密度Welch 法长段估计 f_signal, psd_signal signal.welch(signal_in, fsfs, nperseglen(signal_in)//4) f_noise, psd_noise signal.welch(noise_in, fsfs, nperseglen(noise_in)//4) # 将噪声 PSD 插值到信号的频率轴上 psd_noise_interp np.interp(f_signal, f_noise, psd_noise, leftpsd_noise[0], rightpsd_noise[-1]) # Wiener 增益加小常数防止除零 gain psd_signal / (psd_signal psd_noise_interp 1e-6) gain np.maximum(gain, 0.1) # 设置增益下限避免过度抑制 # 将信号变换到频域 fft_signal np.fft.rfft(signal_in) freqs_fft np.fft.rfftfreq(len(signal_in), 1/fs) # 将增益插值到 FFT 频率点 gain_fft np.interp(freqs_fft, f_signal, gain) # 应用增益并反变换 fft_filtered fft_signal * gain_fft filtered_signal np.fft.irfft(fft_filtered) return filtered_signal # 将噪声轮廓平铺至与原始信号等长984 帧 × 20480 样本 noise_duplicated np.tile(noise_avg, 984) raw_16min data_array[:, 0, :].flatten() # 原始 16 分钟信号 wiener_16min wiener_filter(raw_16min, noise_duplicated) # 重塑为帧数组 [984, 20480] 供后续分析 data_array_Wiener wiener_16min.reshape(984, 20480) np.save(data_array_Wiener.npy, data_array_Wiener) # 保存去噪结果供 Part II 使用 # 3. LSTM 预测 STD 实现早期预警 # 计算每帧的标准差STD作为时域特征 n_files data_array_Wiener.shape[0] std_sequence np.zeros(n_files) for i in range(n_files): std_sequence[i] np.std(data_array_Wiener[i, :]) # 标准化 STD 序列零均值单位方差利于 LSTM 训练 scaler StandardScaler() std_norm scaler.fit_transform(std_sequence.reshape(-1, 1)).flatten() # 选取健康帧50–250训练 LSTM 预测器 healthy_start, healthy_end 50, 250 window 10 # 用过去 10 帧预测下一帧 n_samples healthy_end - healthy_start - window 1 X_healthy np.array([std_norm[i:iwindow] for i in range(healthy_start, healthy_start n_samples)]) y_healthy std_norm[healthy_start window : healthy_end 1] # 注意对齐 X_healthy X_healthy.reshape(X_healthy.shape[0], window, 1) # 构建 LSTM 模型单层 50 单元输出层为 Dense(1) model Sequential() model.add(Input(shape(window, 1))) model.add(LSTM(50, return_sequencesFalse)) model.add(Dense(1)) model.compile(optimizeradam, lossmse) model.fit(X_healthy, y_healthy, epochs100, validation_split0.2, verbose0) # 在中段300–800 帧进行预测检测异常 mid_start, mid_end 300, 800 mid_slice std_norm[mid_start:mid_end] n_mid len(mid_slice) - window 1 mid_seq np.zeros((n_mid, window, 1)) for j in range(n_mid): mid_seq[j] mid_slice[j:jwindow].reshape(window, 1) pred model.predict(mid_seq).flatten() # 计算健康段 STD 的均值和标准差设定预警阈值均值 1.5 倍标准差 healthy_std std_norm[healthy_start:healthy_end] threshold np.mean(healthy_std) 1.5 * np.std(healthy_std) # 找出预测值超过阈值的帧索引 alarm_indices np.where(pred threshold)[0] if len(alarm_indices) 0: first_alarm_frame alarm_indices[0] mid_start window print(f首次预警发生在第 {first_alarm_frame} 帧) else: print(无预警)担任《Mechanical System and Signal Processing》《中国电机工程学报》《宇航学报》《控制与决策》等期刊审稿专家擅长领域信号滤波/降噪机器学习/深度学习时间序列预分析/预测设备故障诊断/缺陷检测/异常检测参考文章物理机理嵌入和自适应学习的机械早期故障诊断Python - 哥廷根数学学派的文章https://zhuanlan.zhihu.com/p/2006771585596551304