PythonAI算法
基于物理机理引导和自编码器融合的机械早期故障诊断
本文介绍了一种结合物理机理与数据驱动的机械早期故障诊断方法。利用 Wiener 去噪和自定义 bTSTFT 变换提取振动信号的时频特征，构建双通道输入。采用卷积自编码器在健康数据上无监督学习正常模式，通过重构误差（MSE）检测异常。引入 CUSUM 累积控制图和滚动方差评分提升鲁棒性。实验表明，该方法能在超帧 460~470 触发报警，比传统指标提前约 60~80 帧，适用于边缘硬件部署。
本文介绍了一种结合物理机理与数据驱动的机械早期故障诊断方法。利用 Wiener 去噪和自定义 bTSTFT 变换提取振动信号的时频特征，构建双通道输入。采用卷积自编码器在健康数据上无监督学习正常模式，通过重构误差（MSE）检测异常。引入 CUSUM 累积控制图和滚动方差评分提升鲁棒性。实验表明，该方法能在超帧 460~470 触发报警，比传统指标提前约 60~80 帧，适用于边缘硬件部署。
首先，利用 Wiener 去噪的振动信号，通过自定义的 bTSTFT 变换提取每个超帧的 256×256 幅值和相位矩阵，构建双通道时频表示。然后，采用卷积自编码器在健康数据上无监督学习正常振动模式的重构；在测试阶段，计算每个超帧的重构误差（MSE）作为异常指标。为提升检测鲁棒性，对 MSE 序列应用 CUSUM 累积和控制图，以及基于滚动方差的早期评分，两者均能有效捕捉退化初期信号的细微变化。实验结果表明，该方法在超帧 460470 附近即触发持续报警，比传统峭度、RMS 等指标提前约 6080 帧（约 1~1.5 分钟）。算法将物理先验（BPFO 谐波抑制、Wiener 去噪）与数据驱动（bTSTFT+ 自编码器）有机结合，在保持计算可行性的前提下（每超帧约 1 秒，可部署于边缘硬件），实现了极高的检测灵敏度和提前量，为旋转机械的预测性维护提供了可靠的技术方案。
从磁盘加载预计算的 bTSTFT 幅值和相位矩阵，每个超帧（由 5 个连续 1 秒帧构成）生成 256×256 的二维时频表示。 对幅值和相位分别进行 Z-score 标准化，使不同量级的特征统一尺度，避免网络训练偏向幅值或相位。 将标准化后的幅值和相位堆叠为双通道输入，形状为 (980, 256, 256, 2)。
根据轴承全寿命周期，将前 300 个超帧（对应原帧 50300）标记为健康状态，用作训练集。
将中间 300800 超帧（对应原帧 300~800）作为测试集，包含退化初期至失效阶段。
从训练集中随机抽取 20% 作为验证集，用于早停和模型选择。
编码器：由两层卷积 + 最大池化组成，逐步压缩空间维度，提取高层次抽象特征（通道数 32→64→128）。 解码器：通过上采样和卷积恢复原始尺寸，输出层采用线性激活以保留幅值和相位的真实动态范围。 损失函数为均方误差（MSE），优化器为 Adam，学习率 0.001。
仅使用健康数据训练自编码器，使其学习重构正常振动模式的 bTSTFT 表示。 设置早停（patience=8）、学习率衰减（factor=0.5）和模型检查点，防止过拟合并保存最佳模型。
对测试集每个超帧，用训练好的自编码器重构，计算像素级 MSE 作为异常分数。 同时计算训练集上的 MSE，确定健康状态的百分位数带（如 P15~P85）作为正常波动范围参考。
以健康 MSE 的均值为目标值，容许偏差 k=0.5×健康标准差，决策阈值 h=5×健康标准差。 对测试集 MSE 序列逐点计算正向累积和，一旦累积和超过 h，判定为持续异常，记录首次报警帧号。 CUSUM 能有效平滑单点噪声，检测到持续的偏离趋势。
对测试集 MSE 序列计算滚动方差（窗口 10），突出变异增强现象（退化初期信号波动加剧）。 计算健康段滚动方差的均值和标准差，设定阈值=均值 +4×标准差。 检测首次超过阈值的超帧，作为另一种互补的早期报警信号。
绘制 MSE 曲线、CUSUM 曲线和方差评分曲线，标注健康带和报警线。
对比传统时域指标（峭度、RMS）的报警时间（约 530540 帧），验证本方法提前至 460470 帧。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.layers import Input, Conv2D, MaxPooling2D, UpSampling2D
from tensorflow.keras.models import Model
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau, ModelCheckpoint
from scipy.fft import fft, fftfreq
# ==================== 1. 加载预计算的 bTSTFT 张量 ====================
# bTSTFT 是自定义宽带时间连续小波变换，包含幅值和相位信息
# 输入：980 个超帧（每个超帧由 5 个连续 1 秒帧构成，步长 1 帧），每个超帧生成 256×256 的幅值和相位矩阵
data_array_btstft_mag = np.load('data_array_bTSTFT_f5_bearing1_focus_Wiev10_v30_mag.npy') # 幅值 (980,256,256)
data_array_btstft_phas = np.load('data_array_bTSTFT_f5_bearing1_focus_Wiev10_v30_phas.npy') # 相位 (980,256,256)
# 对幅值和相位分别进行标准化（Z-score），使数据分布稳定
mag_flat = data_array_btstft_mag.flatten().reshape(-1, 1)
phas_flat = data_array_btstft_phas.flatten().reshape(-1, 1)
scaler_mag = StandardScaler()
scaler_phas = StandardScaler()
mag_norm_flat = scaler_mag.fit_transform(mag_flat)
phas_norm_flat = scaler_phas.fit_transform(phas_flat)
mag_norm = mag_norm_flat.reshape(data_array_btstft_mag.shape)
phas_norm = phas_norm_flat.reshape(data_array_btstft_phas.shape)
# 将标准化后的幅值和相位堆叠为双通道输入 (980,256,256,2)
input_data = np.stack([mag_norm, phas_norm], axis=-1)
# ==================== 2. 划分训练集和测试集 ====================
# 训练集：健康超帧（原帧索引 50~300），测试集：待检测超帧（原帧索引 300~800）
healthy_start, healthy_end = 50, 300
test_start, test_end = 300, 800
X_train = input_data[healthy_start:healthy_end] # 训练数据 (250,256,256,2)
X_test = input_data[test_start:test_end] # 测试数据 (500,256,256,2)
# 从训练集中再划分验证集（20%），用于早停和调参
from sklearn.model_selection import train_test_split
X_train, X_val, _, _ = train_test_split(X_train, X_train, test_size=0.2, random_state=10)
# ==================== 3. 构建卷积自编码器 ====================
input_layer = Input(shape=(256, 256, 2))
# 编码器：逐步下采样，提取深层特征
x = Conv2D(32, (5,5), activation='relu')(input_layer)
x = MaxPooling2D(2,2)(x)
x = Conv2D(64, (3,3), activation='relu')(x)
x = MaxPooling2D(2,2)(x)
encoded = Conv2D(128, (3,3), activation='relu')(x)
# 解码器：上采样恢复原始尺寸，输出线性激活以保留动态范围
x = Conv2D(64, (3,3), activation='relu')(encoded)
x = UpSampling2D(2)(x)
x = Conv2D(32, (5,5), activation='relu')(x)
x = UpSampling2D(2)(x)
decoded = Conv2D(2, (3,3), activation='linear')(x)
autoencoder = Model(input_layer, decoded)
autoencoder.compile(optimizer=Adam(learning_rate=0.001), loss='mse')
# ==================== 4. 训练自编码器（仅用健康数据） ====================
early_stop = EarlyStopping(monitor='val_loss', patience=8, restore_best_weights=True)
reduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=6, min_lr=1e-6)
checkpoint = ModelCheckpoint('best_autoencoder.keras', monitor='val_loss', save_best_only=True)
history = autoencoder.fit(
X_train, X_train, epochs=40, batch_size=16,
validation_data=(X_val, X_val), shuffle=False,
callbacks=[early_stop, reduce_lr, checkpoint], verbose=1
)
# ==================== 5. 计算重构误差（MSE）并检测异常 ====================
# 对测试集进行预测并计算每个超帧的 MSE
recon_test = autoencoder.predict(X_test)
mse_per_surframe = np.mean((X_test - recon_test)**2, axis=(1,2,3)) # (500,)
# 计算训练集上的 MSE，用于确定健康状态的百分位数带
recon_train = autoencoder.predict(X_train)
train_mse = np.mean((X_train - recon_train)**2, axis=(1,2,3))
p15_threshold = np.percentile(train_mse, 15) # 健康带下限
p85_threshold = np.percentile(train_mse, 85) # 健康带上限
# 绘制 MSE 曲线，观察早期偏离
indices = np.arange(test_start, test_start + len(mse_per_surframe))
plt.figure()
plt.plot(indices, mse_per_surframe, label='MSE Reconstruction')
plt.axhspan(p15_threshold, p85_threshold, color='green', alpha=0.15, label='Healthy band (P15-P85)')
plt.axvline(470, color='red', linestyle='--', label='Dip ~470')
plt.title('Reconstruction MSE on Test Set')
plt.legend()
plt.show()
# ==================== 6. CUSUM 累积和控制图 ====================
# 参考健康 MSE 的均值和标准差，设置 CUSUM 参数
mse_ref = np.mean(train_mse) # 健康均值
k = 0.5 * np.std(train_mse) # 容许偏差（斜率）
h = 5 * np.std(train_mse) # 决策阈值
cusum_pos = 0
cusum = [0]
for m in mse_per_surframe:
cusum_pos = max(0, cusum_pos + (m - mse_ref) - k)
cusum.append(cusum_pos)
cusum = cusum[:len(mse_per_surframe)] # 对齐长度
# 找出首次超过阈值 h 的超帧索引
crossing = np.where(np.array(cusum) > h)[0]
if len(crossing) > 0:
first_alarm = crossing[0] + test_start
print(f"CUSUM 首次报警在超帧 {first_alarm}")
# ==================== 7. 基于 MSE 方差的早期评分 ====================
# 对 MSE 序列计算滚动方差（窗口 10），突出变异增强
rolling_var = pd.Series(mse_per_surframe).rolling(window=10, min_periods=1).var()
# 计算健康段方差均值和阈值（均值 +4 倍标准差）
healthy_var = rolling_var[:healthy_end - test_start] # 对应健康部分
thresh_var = np.mean(healthy_var) + 4 * np.std(healthy_var)
# 检测首次超过阈值的帧
alarm_var = np.where(rolling_var > thresh_var)[0]
if len(alarm_var) > 0:
first_alarm_var = alarm_var[0] + test_start
print(f"方差评分首次报警在超帧 {first_alarm_var}")
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
生成新的随机RSA私钥和公钥pem证书。 在线工具，RSA密钥对生成器在线工具，online
基于 Mermaid.js 实时预览流程图、时序图等图表，支持源码编辑与即时渲染。 在线工具，Mermaid 预览与可视化编辑在线工具，online
解析常见 curl 参数并生成 fetch、axios、PHP curl 或 Python requests 示例代码。 在线工具，curl 转代码在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online