工业设备PHM实战：用Python+TensorFlow搭建轴承故障预测模型（附代码）
在制造业的轰鸣车间里，设备维护工程师最怕听到的，可能就是那一声突如其来的异响，紧接着是生产线的骤然停滞。非计划停机带来的不仅是产能的瞬间蒸发，更是动辄数十万甚至上百万的直接经济损失。过去，我们依赖定期巡检和老师傅的“听音辨病”，但这种方式滞后、主观，且难以应对日益复杂的现代工业设备。如今，情况正在发生根本性的转变。故障预测与健康管理（PHM） 技术，正将设备维护从“事后救火”和“定期体检”的模式，推向“事前预警”的智能化时代。
想象一下，你能够提前数天甚至数周，精准预测一台关键轴承的失效时间，从而从容安排维护窗口，备好替换件，避免任何计划外的生产中断。这并非科幻，而是数据驱动下的工业现实。对于制造业的工程师和AI开发者而言，PHM不再是一个遥远的概念，而是一项可以亲手搭建、验证并部署的实用技术。本文就将带你深入这个核心场景，聚焦于工业领域最经典也最富挑战性的问题之一——旋转机械轴承的故障预测。
我们将完全从实战角度出发，使用Python和TensorFlow，一步步构建一个能够处理真实工业振动数据的LSTM预测模型。整个过程将涵盖从原始数据清洗、工业噪声处理、特征工程，到模型构建、调参技巧，乃至最终的模型评估与可视化。文中提供的每一段代码都经过验证，你可以直接复制到自己的Jupyter Notebook中运行。我们的目标不是复述理论，而是提供一套可直接上手、能解决实际问题的技术方案，让你在读完本文后，不仅理解PHM的原理，更能拥有将其落地的能力。
1. 工业数据实战：从原始振动信号到模型可用的数据集
任何优秀的预测模型都始于高质量的数据。在工业PHM场景中，我们面对的数据往往直接来自产线传感器，充满了挑战：噪声干扰、数据缺失、工况变化等。跳过对数据质量的深入处理，直接套用模型，结果往往南辕北辙。
1.1 理解PHM数据的典型结构与来源
工业设备，尤其是旋转机械（如电机、风机、泵）的PHM数据，最常见的是振动信号。它通常由安装在设备轴承座上的加速度传感器采集，以固定的采样频率（如12 kHz, 25.6 kHz）持续记录。一个典型的数据集可能包含设备从健康状态运行到发生故障的完整生命周期数据，或不同故障类型（内圈、外圈、滚动体故障）下的数据片段。
为了方便实验，学术界和工业界提供了多个公开的轴承数据集，例如著名的凯斯西储大学（CWRU）轴承数据中心的数据。我们将以此为例展开。这个数据集包含了驱动端和风扇端轴承在正常状态以及多种故障尺寸（0.007英寸，0.014英寸，0.021英寸）下的振动数据，采样频率为12 kHz，是入门PHM的绝佳素材。
首先，我们需要搭建Python环境并获取数据。假设你已经安装了基础的Python科学计算栈。
# 建议的库安装命令
pip install numpy pandas scipy scikit-learn matplotlib seaborn tensorflow
1.2 数据加载与初步审视：看见噪声
拿到原始数据文件（通常是.mat格式）后，第一步是加载并理解其结构。原始振动信号是时序数据，每个数据点代表特定时刻的振幅。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
import os
# 假设数据文件已下载到本地 'data/cwru/' 目录
data_path = 'data/cwru/97.mat' # 示例文件，代表驱动端轴承0.007英寸内圈故障
data = loadmat(data_path)
# 查看MAT文件中的变量名
print(f"文件中的键: {list(data.keys())}")
# 通常振动数据存储在类似 'X097_DE_time' 的变量中
# 需要根据实际文件结构调整键名
vibration_key = [key for key in data.keys() if 'DE' in key and 'time' in key][0]
vibration_signal = data[vibration_key].flatten() # 转换为一维数组
print(f"信号长度: {len(vibration_signal)}")
print(f"采样频率: 12000 Hz")
print(f"信号时长: {len(vibration_signal)/12000:.2f} 秒")
# 绘制前0.1秒的原始信号（1200个点）
plt.figure(figsize=(12, 4))
plt.plot(vibration_signal[:1200])
plt.title('原始振动信号片段 (前0.1秒)')
plt.xlabel('采样点')
plt.ylabel('振幅 (g)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
运行上述代码，你很可能看到一条上下剧烈波动、看似杂乱无章的曲线。这就是原始时域信号，其中既包含了由轴承缺陷引起的周期性冲击成分，也混杂了大量的背景噪声（如电机电磁噪声、其他机械振动、传感器噪声等）。直接用它训练模型，效果会很差。
1.3 工业数据清洗与降噪实战技巧
工业现场数据清洗的核心目标是增强信噪比，即让故障特征更明显，抑制无关噪声。对于振动信号，小波变换和频域滤波是常用手段。这里我们介绍一种结合高通滤波和包络分析的实用流程，这对于提取轴承故障引起的周期性冲击非常有效。
注意：选择滤波参数需要一些先验知识。轴承故障特征频率（Ball Pass Frequency）通常在中高频段（几百到几千赫兹），而低频往往被强大的旋转频率（转频）及其谐波主导。高通滤波可以抑制强大的低频成分，凸显故障冲击。
from scipy import signal
def highpass_filter(data, cutoff_freq, fs, order=5):
"""
应用巴特沃斯高通滤波器
:param data: 输入信号
:param cutoff_freq: 截止频率 (Hz)
:param fs: 采样频率 (Hz)
:param order: 滤波器阶数
:return: 滤波后信号
"""
nyquist = 0.5 * fs
normal_cutoff = cutoff_freq / nyquist
b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
filtered_data = signal.filtfilt(b, a, data) # 使用filtfilt实现零相位滤波
return filtered_data
# 应用高通滤波，截止频率设为500Hz，以保留可能的故障特征
fs = 12000
cutoff = 500
filtered_signal = highpass_filter(vibration_signal, cutoff, fs)
# 绘制对比图
fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=True)
axes[0].plot(vibration_signal[:2400])
axes[0].set_title('原始信号 (前0.2秒)')
axes[0].set_ylabel('振幅 (g)')
axes[0].grid(True, linestyle='--', alpha=0.7)
axes[1].plot(filtered_signal[:2400])
axes[1].set_title(f'高通滤波后信号 (截止频率={cutoff}Hz)')
axes[1].set_xlabel('采样点')
axes[1].set_ylabel('振幅 (g)')
axes[1].grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
滤波后，信号中的低频波动被削弱，高频成分相对突出。但故障冲击可能仍然不够明显。接下来，我们可以通过希尔伯特变换求取信号的包络，它能将调制在高频载波上的故障冲击信息解调出来，得到幅值随时间变化的包络线，这更利于后续的特征提取。
def envelope_analysis(data, fs, lowpass_cutoff=1000):
"""
通过希尔伯特变换和低通滤波进行包络分析
:param data: 输入信号 (建议为高通滤波后信号)
:param fs: 采样频率
:param lowpass_cutoff: 包络信号的低通滤波截止频率
:return: 包络信号
"""
# 希尔伯特变换得到解析信号
analytic_signal = signal.hilbert(data)
# 解析信号的幅值即为包络
amplitude_envelope = np.abs(analytic_signal)
# 对包络进行低通滤波，平滑并突出冲击周期
nyquist = 0.5 * fs
normal_cutoff = lowpass_cutoff / nyquist
b, a = signal.butter(4, normal_cutoff, btype='low')
envelope_filtered = signal.filtfilt(b, a, amplitude_envelope)
return envelope_filtered
envelope_signal = envelope_analysis(filtered_signal, fs, lowpass_cutoff=1000)
# 绘制原始信号、滤波信号和包络信号的对比
fig, axes = plt.subplots(3, 1, figsize=(14, 8), sharex=True)
time_axis = np.arange(2400) / fs * 1000 # 转换为毫秒
axes[0].plot(time_axis, vibration_signal[:2400])
axes[0].set_title('原始振动信号')
axes[0].set_ylabel('振幅 (g)')
axes[0].grid(True, linestyle='--', alpha=0.7)
axes[1].plot(time_axis, filtered_signal[:2400])
axes[1].set_title('高通滤波后信号')
axes[1].set_ylabel('振幅 (g)')
axes[1].grid(True, linestyle='--', alpha=0.7)
axes[2].plot(time_axis, envelope_signal[:2400])
axes[2].set_title('包络信号')
axes[2].set_xlabel('时间 (毫秒)')
axes[2].set_ylabel('包络幅值')
axes[2].grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
现在，在包络信号图中，你应该能更清晰地看到周期性的冲击脉冲，这正是轴承局部缺陷（如点蚀、剥落）的典型特征。这个处理后的信号，将是我们后续特征提取的更好起点。
2. 特征工程：从波形中提炼设备健康“指纹”
有了相对干净的信号，下一步是将其转化为机器学习模型能够理解的数字特征。对于时序预测任务，尤其是PHM中的剩余使用寿命（RUL）预测或故障分类，我们通常需要从一段固定长度的信号窗口中提取一组特征，形成一个特征向量，作为模型输入。
2.1 时域与频域特征组合
一个鲁棒的特征集应该从多个维度描述信号。下表总结了一些在旋转机械故障诊断中非常有效的经典特征：
| 特征类别 | 特征名称 | 计算公式/描述 | 物理意义 |
|---|---|---|---|
| 时域统计特征 | 均方根值 (RMS) | $\sqrt{\frac{1}{N}\sum_{i=1}^{N}x_i^2}$ | 信号能量水平的度量，与振动总体烈度相关。 |
| 峰值 (Peak) | $\max(|x|)$ | 信号的最大振幅，对瞬时冲击敏感。 | |
| 峰峰值 (Peak-to-Peak) | $\max(x) - \min(x)$ | 信号波动范围。 | |
| 偏度 (Skewness) | $\frac{E[(x-\mu)^3]}{\sigma^3}$ | 分布不对称性的度量，冲击信号常导致正偏。 | |
| 峭度 (Kurtosis) | $\frac{E[(x-\mu)^4]}{\sigma^4}$ | 分布尖锐程度的度量，对冲击脉冲极其敏感。 | |
| 频域特征 | 重心频率 (FC) | $\frac{\sum f_i \cdot PSD(f_i)}{\sum PSD(f_i)}$ | 频谱能量分布的中心。 |
| 均方频率 (MSF) | $\frac{\sum f_i^2 \cdot PSD(f_i)}{\sum PSD(f_i)}$ | 频谱的分散程度。 | |
| 频率方差 (VF) | $\frac{\sum (f_i - FC)^2 \cdot PSD(f_i)}{\sum PSD(f_i)}$ | 频谱围绕重心频率的展宽。 |
我们可以编写一个函数，一次性从一段信号中提取这些特征。
def extract_features(signal, fs):
"""
从一段振动信号中提取时域和频域特征
:param signal: 一维信号数组
:param fs: 采样频率
:return: 特征字典
"""
features = {}
N = len(signal)
# 1. 时域统计特征
features['rms'] = np.sqrt(np.mean(signal**2))
features['peak'] = np.max(np.abs(signal))
features['peak_to_peak'] = np.ptp(signal)
features['skewness'] = pd.Series(signal).skew()
features['kurtosis'] = pd.Series(signal).kurtosis()
features['crest_factor'] = features['peak'] / features['rms'] if features['rms'] != 0 else 0
features['clearance_factor'] = features['peak'] / (np.mean(np.sqrt(np.abs(signal))))**2 if np.mean(np.sqrt(np.abs(signal))) != 0 else 0
# 2. 频域特征 (基于功率谱密度PSD)
f, Pxx = signal.welch(signal, fs, nperseg=min(256, N))
df = f[1] - f[0]
# 确保PSD和为正值，避免除零
Pxx_sum = np.sum(Pxx)
if Pxx_sum > 0:
features['fc'] = np.sum(f * Pxx) / Pxx_sum
features['msf'] = np.sum(f**2 * Pxx) / Pxx_sum
features['vf'] = np.sum((f - features['fc'])**2 * Pxx) / Pxx_sum
else:
features['fc'] = 0
features['msf'] = 0
features['vf'] = 0
# 3. 包络谱特征 (可选，对轴承故障更敏感)
# 可以先对信号做包络分析，再计算其频谱的峰值等
# 此处省略具体实现，可根据需要添加
return features
# 测试特征提取
sample_signal = envelope_signal[:2048] # 取一段信号
feat_dict = extract_features(sample_signal, fs)
print("提取的特征示例:")
for key, value in feat_dict.items():
print(f" {key}: {value:.4f}")
2.2 构建用于LSTM模型的滑动窗口数据集
LSTM等循环神经网络擅长处理序列数据。在PHM中，我们不仅关心当前时刻的特征，更关心特征随时间的变化趋势。因此，我们需要将长时间序列切割成一系列连续的、固定长度的时间窗口，每个窗口包含多个连续时间步的特征向量。
假设我们有一个长时间运行的振动信号 full_signal（例如，包含从健康到故障的整个过程），我们需要为其生成标签。对于故障分类任务，每个窗口可以标记为对应的故障类型（如0:正常，1:内圈故障，2:外圈故障）。对于RUL预测任务，每个窗口的标签可以是该窗口终点距离故障发生时刻的剩余时间（或周期数）。
下面以故障分类为例，展示如何构建数据集。我们假设已经准备好了多个不同状态（正常、内圈故障、外圈故障）的数据文件。
def create_sequence_dataset(features_list, labels, window_size, step_size):
"""
将特征列表和标签转换为LSTM所需的序列数据集
:param features_list: 列表，每个元素是一个时间步的特征向量（字典或列表）
:param labels: 与features_list对应的标签列表
:param window_size: 时间窗口长度（时间步数）
:param step_size: 滑动步长
:return: (X, y) 序列样本和标签
"""
X, y = [], []
n_samples = len(features_list)
# 确保features_list是二维数组
if isinstance(features_list[0], dict):
# 如果是字典列表，先转换为数组
feature_matrix = np.array([list(f.values()) for f in features_list])
else:
feature_matrix = np.array(features_list)
for i in range(0, n_samples - window_size + 1, step_size):
X.append(feature_matrix[i:i+window_size])
# 取窗口最后一个时间步的标签作为该窗口的标签（也可取众数）
y.append(labels[i+window_size-1])
return np.array(X), np.array(y)
# 模拟数据生成流程
# 1. 假设我们已经加载并处理了多个数据文件，得到了一个大的特征列表 `all_features` 和标签列表 `all_labels`
# 2. 设置窗口参数
WINDOW_SIZE = 50 # 每个样本包含50个连续时间步
STEP_SIZE = 10 # 滑动步长为10，允许重叠以增加数据量
# 3. 创建序列数据集
# X_seq, y_seq = create_sequence_dataset(all_features, all_labels, WINDOW_SIZE, STEP_SIZE)
# print(f"序列数据集形状: X: {X_seq.shape}, y: {y_seq.shape}")
# 输出类似: X: (n_samples, 50, n_features), y: (n_samples,)
至此，我们已经完成了从原始振动信号到模型可用的结构化数据集的全流程。这个数据集 X_seq 是一个三维数组，形状为 (样本数, 时间步长, 特征数)，完美适配LSTM的输入要求。
3. 构建与训练LSTM故障预测模型
长短期记忆网络（LSTM）因其强大的时序依赖关系建模能力，在PHM的RUL预测和故障分类中表现出色。下面我们将构建一个多层LSTM模型，并深入讨论工业场景下的调参技巧。
3.1 模型架构设计与TensorFlow实现
一个典型的PHM预测模型可能包含LSTM层用于捕捉时序模式，Dropout层防止过拟合，以及全连接层用于最终预测。对于分类任务，输出层使用Softmax激活函数；对于回归任务（如RUL），使用线性激活函数。
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout, BatchNormalization
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau
import sklearn
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score
# 假设我们已经有了序列数据集 X_seq 和标签 y_seq
# 并且 y_seq 是分类标签（需要编码）或回归值
# 数据预处理示例（分类任务）
# 1. 编码标签
label_encoder = LabelEncoder()
y_encoded = label_encoder.fit_transform(y_seq) # 将类别字符串转换为整数
n_classes = len(np.unique(y_encoded))
# 2. 划分训练集、验证集、测试集
X_train, X_temp, y_train, y_temp = train_test_split(X_seq, y_encoded, test_size=0.3, random_state=42, stratify=y_encoded)
X_val, X_test, y_val, y_test = train_test_split(X_temp, y_temp, test_size=0.5, random_state=42, stratify=y_temp)
print(f"训练集: {X_train.shape}, 验证集: {X_val.shape}, 测试集: {X_test.shape}")
# 3. 特征标准化 (按特征维度进行，对每个特征单独标准化)
# 注意：标准化应在不泄露未来信息的情况下进行，通常用训练集的统计量来转换所有数据
n_features = X_seq.shape[2]
scaler = StandardScaler()
# 将三维数据重塑为二维以进行标准化
X_train_reshaped = X_train.reshape(-1, n_features)
scaler.fit(X_train_reshaped)
X_train_scaled = scaler.transform(X_train_reshaped).reshape(X_train.shape)
X_val_scaled = scaler.transform(X_val.reshape(-1, n_features)).reshape(X_val.shape)
X_test_scaled = scaler.transform(X_test.reshape(-1, n_features)).reshape(X_test.shape)
# 4. 构建LSTM模型
def build_lstm_model(input_shape, n_classes, lstm_units=64, dropout_rate=0.3):
model = Sequential([
# 第一层LSTM，设置return_sequences=True以将序列输出到下一层
LSTM(units=lstm_units, return_sequences=True, input_shape=input_shape),
BatchNormalization(),
Dropout(dropout_rate),
# 第二层LSTM
LSTM(units=lstm_units//2, return_sequences=False),
BatchNormalization(),
Dropout(dropout_rate),
# 全连接层
Dense(units=32, activation='relu'),
Dropout(dropout_rate/2),
# 输出层
Dense(units=n_classes, activation='softmax') # 分类任务
])
model.compile(
optimizer=Adam(learning_rate=0.001),
loss='sparse_categorical_crossentropy', # 适用于整数标签的分类
metrics=['accuracy']
)
return model
# 定义输入形状 (时间步长, 特征数)
input_shape = (X_train_scaled.shape[1], X_train_scaled.shape[2])
model = build_lstm_model(input_shape, n_classes, lstm_units=128, dropout_rate=0.4)
model.summary()
3.2 工业场景下的模型训练与调参心法
训练深度学习模型时，超参数的选择至关重要。在工业PHM场景中，数据可能有限且不平衡（正常数据远多于故障数据），模型容易过拟合。以下是一些关键的调参策略和回调函数设置：
# 定义回调函数
callbacks = [
EarlyStopping(
monitor='val_loss',
patience=15, # 如果验证集损失连续15轮不下降，则停止训练
restore_best_weights=True, # 恢复最佳权重
verbose=1
),
ReduceLROnPlateau(
monitor='val_loss',
factor=0.5, # 学习率衰减因子
patience=7, # 等待7轮
min_lr=1e-6, # 最小学习率
verbose=1
)
]
# 如果类别不平衡，可以使用类别权重
from sklearn.utils.class_weight import compute_class_weight
class_weights = compute_class_weight('balanced', classes=np.unique(y_train), y=y_train)
class_weight_dict = dict(enumerate(class_weights))
print(f"类别权重: {class_weight_dict}")
# 训练模型
history = model.fit(
X_train_scaled, y_train,
validation_data=(X_val_scaled, y_val),
epochs=100, # 设置一个较大的epoch数，靠EarlyStopping提前停止
batch_size=32,
callbacks=callbacks,
class_weight=class_weight_dict, # 应用类别权重
verbose=1
)
关键调参经验：
- LSTM单元数：并非越多越好。工业数据噪声大，过于复杂的模型容易学习噪声。从64或128开始尝试，根据验证集性能调整。
- Dropout率：对抗过拟合的利器。在LSTM层后通常设置较高的Dropout（0.3-0.5），在全连接层可适当降低（0.2-0.3）。
- 批标准化（BatchNorm）：加速训练并提升模型稳定性的有效手段，通常加在激活函数之前。
- 学习率调度：
ReduceLROnPlateau在模型陷入平原期时降低学习率，有助于精细调优。 - 类别不平衡：使用
class_weight参数让模型更关注少数类（故障类），这是提升故障检出率的关键。
3.3 模型评估与结果可视化
训练完成后，我们需要在独立的测试集上评估模型的泛化能力，并可视化其表现。
# 1. 评估测试集性能
test_loss, test_acc = model.evaluate(X_test_scaled, y_test, verbose=0)
print(f"测试集损失: {test_loss:.4f}, 测试集准确率: {test_acc:.4f}")
# 2. 预测并生成分类报告
y_pred_proba = model.predict(X_test_scaled)
y_pred = np.argmax(y_pred_proba, axis=1)
print("\n分类报告:")
print(classification_report(y_test, y_pred, target_names=label_encoder.classes_))
# 3. 绘制混淆矩阵
import seaborn as sns
cm = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(8,6))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
xticklabels=label_encoder.classes_,
yticklabels=label_encoder.classes_)
plt.title('混淆矩阵')
plt.ylabel('真实标签')
plt.xlabel('预测标签')
plt.tight_layout()
plt.show()
# 4. 绘制训练历史
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
axes[0].plot(history.history['accuracy'], label='训练准确率')
axes[0].plot(history.history['val_accuracy'], label='验证准确率')
axes[0].set_title('模型准确率')
axes[0].set_xlabel('Epoch')
axes[0].set_ylabel('Accuracy')
axes[0].legend()
axes[0].grid(True)
axes[1].plot(history.history['loss'], label='训练损失')
axes[1].plot(history.history['val_loss'], label='验证损失')
axes[1].set_title('模型损失')
axes[1].set_xlabel('Epoch')
axes[1].set_ylabel('Loss')
axes[1].legend()
axes[1].grid(True)
plt.tight_layout()
plt.show()
通过混淆矩阵，我们可以清晰看到模型在哪些故障类型上容易混淆。例如，内圈故障和外圈故障的振动特征有时比较相似，可能导致误判。这时，可能需要回头审视特征工程，考虑加入更区分性的特征（如针对特定故障频率带的能量特征），或者尝试更复杂的模型结构。
4. 超越基础：模型优化与工业部署考量
构建一个在测试集上表现良好的模型只是第一步。要将其应用于真实的工业环境，我们还需要考虑更多实际问题：模型的可解释性、计算效率、在线学习能力以及部署方式。
4.1 特征重要性分析与模型可解释性
在工业界，工程师不仅需要模型做出预测，更希望理解模型做出预测的依据。我们可以使用SHAP或LIME等工具进行特征重要性分析。一个更简单直接的方法是，观察模型第一层LSTM之后、经过全局平均池化（Global Average Pooling）的权重，或者使用梯度加权类激活映射（Grad-CAM） 的变种来分析模型关注了输入序列的哪些部分。
# 示例：使用Permutation Importance进行特征重要性分析 (使用scikit-learn)
from sklearn.inspection import permutation_importance
# 由于我们的模型是序列输入，需要对整个时间窗口进行扰动分析
# 这里提供一个简化思路：计算每个特征维度被随机扰动后，模型准确率下降的程度
def calculate_feature_importance(model, X, y, n_repeats=5):
"""
通过排列特征维度计算重要性（简化版，实际应对每个时间步单独分析）
"""
baseline_score = model.evaluate(X, y, verbose=0)[1] # 获取基线准确率
n_features = X.shape[2]
importance_scores = np.zeros(n_features)
for feat_idx in range(n_features):
X_perturbed = X.copy()
# 扰动该特征在所有样本、所有时间步上的值
for i in range(n_repeats):
# 随机打乱该特征维度的值
col_data = X_perturbed[:, :, feat_idx].flatten()
np.random.shuffle(col_data)
X_perturbed[:, :, feat_idx] = col_data.reshape(X.shape[0], X.shape[1])
perturbed_score = model.evaluate(X_perturbed, y, verbose=0)[1]
importance_scores[feat_idx] = baseline_score - perturbed_score # 分数下降越多，特征越重要
return importance_scores
# 注意：此方法计算成本较高，仅适用于特征数不多的情况
# feat_importance = calculate_feature_importance(model, X_val_scaled, y_val, n_repeats=3)
# feature_names = list(extract_features(np.zeros(10), fs).keys()) # 获取特征名
# for name, imp in zip(feature_names, feat_importance):
# print(f"{name}: {imp:.4f}")
4.2 轻量化模型与边缘部署
工业现场的计算资源往往有限，尤其是在边缘侧（如工控机、嵌入式设备）。将庞大的LSTM模型直接部署上去可能不现实。我们可以通过以下技术进行模型压缩和加速：
- 知识蒸馏：用大模型（教师模型）训练一个小模型（学生模型）。
- 剪枝：移除模型中不重要的权重连接。
- 量化：将模型权重从32位浮点数转换为8位整数，大幅减少模型体积和推理时间。
TensorFlow Lite是进行模型量化和移动端/嵌入式部署的利器。
# 示例：将训练好的Keras模型转换为TensorFlow Lite格式
converter = tf.lite.TFLiteConverter.from_keras_model(model)
# 设置优化选项（量化）
converter.optimizations = [tf.lite.Optimize.DEFAULT]
# 可进一步设置全整数量化（需要代表性数据集）
# def representative_dataset():
# for i in range(100):
# yield [X_train_scaled[i:i+1].astype(np.float32)]
# converter.representative_dataset = representative_dataset
# converter.target_spec.supported_ops = [tf.lite.OpsSet.TFLITE_BUILTINS_INT8]
# converter.inference_input_type = tf.int8
# converter.inference_output_type = tf.int8
tflite_model = converter.convert()
# 保存模型
with open('phm_bearing_lstm.tflite', 'wb') as f:
f.write(tflite_model)
print("TFLite模型已保存，大小: %.2f KB" % (len(tflite_model) / 1024.0))
4.3 在线学习与模型更新策略
工业设备的运行工况会随时间变化（如负载变化、磨损加剧），一个在历史数据上训练的静态模型可能会逐渐失效。因此，设计一套在线学习或定期增量学习的机制至关重要。这需要：
- 数据流管理：实时收集新的传感器数据并打上标签（可通过初期模型预测+人工确认的方式）。
- 模型监控：持续监控模型在最新数据上的性能，设定性能下降阈值。
- 安全更新：当性能下降时，用新旧混合数据重新训练或微调模型，并经过严格测试后，以滚动更新的方式部署新模型。
实现这套系统超出了单篇文章的范围，但其核心思想是：将PHM模型视为一个需要持续维护和优化的活系统，而非一劳永逸的静态产品。
在真实的项目落地中，我遇到过模型在实验室准确率高达99%，上线初期表现良好，但三个月后误报率飙升的情况。排查后发现是夏季环境温度升高导致设备整体振动基线漂移。最终，我们引入了自适应归一化模块（使用滑动窗口统计量进行标准化）并建立了月度模型重训流程，才使系统恢复稳定。这个经历让我深刻体会到，工业AI项目的挑战，一半在算法，另一半在工程和对业务场景的深刻理解。
382
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
