工厂设备故障预警实战中时间序列过程识别如何提前捕捉异常信号并延长设备使用寿命
一、先聊聊这个痛点：设备说”我快不行了”之前，能提前多久听到它的”呻吟”？
我是老王，在一家汽车零部件厂干了十二年设备管理。说实话，以前我们最头疼的就是设备半夜突然趴窝——产线停了，工人等着，老板骂街，维修工连夜抢修，修好了还得挨批。
但你有没有发现，设备”生病”之前其实是有征兆的？就像人发烧前会感觉乏力、食欲不振一样。问题是如何在它还”健康”的时候，就把那些微小的异常信号捕捉到？
这就是时间序列过程识别要解决的问题。
二、时间序列到底是个啥？别被名字吓到
用大白话说，时间序列就是一串按时间顺序排列的数据。
比如：
- 温度：每秒记录一次，一天就是86400个数据点
- 振动：每毫秒采集一次，能看出轴承的微妙变化
- 压力：每分钟记录一次，反映液压系统的状态
这些数据看起来杂七杂八，但如果用对方法，它们能告诉你设备”内心”的真实想法。
三、从数据到洞察：时间序列过程识别的完整流程
3.1 数据采集是地基
先看我厂的实际场景——一台CNC加工中心，我们安装了以下传感器：
| 传感器类型 | 测量参数 | 采样频率 | 作用 |
|---|---|---|---|
| 振动传感器 | 三轴振动加速度 | 10kHz | 检测轴承、齿轮箱异常 |
| 温度传感器 | 主轴温度 | 1Hz | 防止过热损坏 |
| 电流传感器 | 主轴电机电流 | 100Hz | 监测负载变化 |
| 压力传感器 | 液压系统压力 | 10Hz | 检测泄漏、堵塞 |
| 声学传感器 | 运行噪音 | 44.1kHz | 捕捉高频异常信号 |
关键经验：采样频率不是越高越好，要匹配你想检测的故障类型。轴承故障通常用振动（高频），而润滑不良用温度（低频）就够了。
3.2 数据预处理：把”噪音”过滤掉
原始数据往往充满噪音，直接分析就像在嘈杂的舞厅里听人说话。
import numpy as np
import pandas as pd
from scipy.signal import butter, filtfilt, stft
from scipy import stats
class DataPreprocessor:
"""
工业设备时间序列数据预处理
"""
def __init__(self, sample_rate):
self.sample_rate = sample_rate
def remove_outliers_iqr(self, data, factor=1.5):
"""
使用IQR方法去除异常值
比简单的3σ方法更适合工业数据，因为不受极端值影响
"""
Q1 = np.percentile(data, 25)
Q3 = np.percentile(data, 75)
IQR = Q3 - Q1
lower_bound = Q1 - factor * IQR
upper_bound = Q3 + factor * IQR
# 保留在合理范围内的数据
mask = (data >= lower_bound) & (data <= upper_bound)
return data[mask], lower_bound, upper_bound
def apply_bandpass_filter(self, data, low_freq, high_freq):
"""
带通滤波：只保留感兴趣的频率范围
例如：检测轴承外圈故障，关注100-500Hz
"""
nyquist = self.sample_rate / 2
low = low_freq / nyquist
high = high_freq / nyquist
# 设计Butterworth滤波器
order = 4
b, a = butter(order, [low, high], btype='band')
# 零相位滤波，避免相位失真
filtered_data = filtfilt(b, a, data)
return filtered_data
def calculate_features(self, data, window_size=1024, hop_size=256):
"""
从时域和频域提取特征
"""
features = {}
# 时域特征
features['mean'] = np.mean(data)
features['rms'] = np.sqrt(np.mean(data**2)) # 均方根，反映能量
features['peak'] = np.max(np.abs(data))
features['crest'] = features['peak'] / features['rms'] # 峰值因子
features['skewness'] = stats.skew(data)
features['kurtosis'] = stats.kurtosis(data) # 峭度，对冲击敏感
# 频域特征（通过FFT）
fft_result = np.fft.rfft(data)
fft_magnitude = np.abs(fft_result)
frequencies = np.fft.rfftfreq(len(data), 1/self.sample_rate)
# 频谱重心
features['spectral_centroid'] = np.sum(frequencies * fft_magnitude) / np.sum(fft_magnitude)
# 各频段能量
band_edges = [0, 100, 500, 1000, 5000] # Hz
for i in range(len(band_edges)-1):
mask = (frequencies >= band_edges[i]) & (frequencies < band_edges[i+1])
band_energy = np.sum(fft_magnitude[mask]**2)
features[f'band_{band_edges[i]}_{band_edges[i+1]}Hz'] = band_energy
return features
def generate_spectrogram(self, data, nperseg=1024):
"""
生成时频谱图，可视化频率随时间的变化
对识别瞬态故障特别有用
"""
frequencies, times, spectrogram = stft(
data,
fs=self.sample_rate,
nperseg=nperseg
)
return frequencies, times, np.abs(spectrogram)
3.3 过程识别：找到设备的”正常模式”
这一步是最关键的。我们需要建立设备的”健康基准线”。
from sklearn.ensemble import IsolationForest
from sklearn.covariance import EllipticEnvelope
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
import seaborn as sns
class ProcessRecognizer:
"""
时间序列过程识别器
用于建立设备正常运行的"指纹"
"""
def __init__(self, contamination=0.05):
"""
contamination: 预期异常比例，设小一点，因为我们认为大部分时间是正常的
"""
self.contamination = contamination
self.scaler = StandardScaler()
self.models = {}
self.health_baseline = None
def fit_normal_model(self, training_data):
"""
使用健康设备数据训练正常模式模型
"""
# 提取特征
features_list = []
for signal in training_data:
# 假设signal是一个长度为1024的振动信号片段
features = self._extract_comprehensive_features(signal)
features_list.append(features)
features_df = pd.DataFrame(features_list)
# 标准化
scaled_data = self.scaler.fit_transform(features_df)
# 方法1：孤立森林 - 对高维数据效果好
iso_forest = IsolationForest(
contamination=self.contamination,
n_estimators=100,
random_state=42
)
iso_forest.fit(scaled_data)
self.models['isolation_forest'] = iso_forest
# 方法2：马氏距离 - 捕捉特征间的相关性
elliptic = EllipticEnvelope(
contamination=self.contamination,
support_fraction=0.9
)
elliptic.fit(scaled_data)
self.models['elliptic_envelope'] = elliptic
# 保存健康基线（各特征的均值和协方差）
self.health_baseline = {
'mean': features_df.mean(),
'std': features_df.std(),
'covariance': np.cov(features_df.T)
}
print(f"✓ 已建立设备健康基准，基于 {len(training_data)} 条健康样本")
print(f" 特征维度: {scaled_data.shape[1]}")
return self
def _extract_comprehensive_features(self, signal):
"""综合特征提取"""
features = {}
# 时域特征
features['time_mean'] = np.mean(signal)
features['time_rms'] = np.sqrt(np.mean(signal**2))
features['time_peak'] = np.max(np.abs(signal))
features['time Crest'] = features['time_peak'] / (features['time_rms'] + 1e-8)
features['time_skew'] = stats.skew(signal)
features['time_kurt'] = stats.kurtosis(signal)
# 频域特征
fft = np.fft.rfft(signal)
magnitudes = np.abs(fft)
frequencies = np.fft.rfftfreq(len(signal))
features['freq_energy'] = np.sum(magnitudes**2)
features['freq_centroid'] = np.sum(frequencies * magnitudes) / (np.sum(magnitudes) + 1e-8)
# 计算各阶共振频率的能量占比
bands = [(0, 500), (500, 1000), (1000, 2000), (2000, 5000)]
for i, (low, high) in enumerate(bands):
mask = (frequencies >= low) & (frequencies < high)
features[f'band_{i}_energy'] = np.sum(magnitudes[mask]**2)
# 包络谱特征（对轴承故障特别有效）
envelope = np.abs(np.fft.hilbert(signal))
envelope_fft = np.abs(np.fft.rfft(envelope))
features['envelope_peak'] = np.max(envelope_fft[:len(envelope_fft)//4])
return features
def get_health_score(self, new_data):
"""
计算健康评分（0-100分，100为最佳）
"""
features = self._extract_comprehensive_features(new_data)
features_array = np.array(list(features.values())).reshape(1, -1)
scaled_features = self.scaler.transform(features_array)
# 多个模型综合评分
iso_score = self.models['isolation_forest'].score_samples(scaled_features)[0]
ell_score = self.models['elliptic_envelope'].score_samples(scaled_features)[0]
# 归一化到0-1范围
iso_normalized = (iso_score - np.min(iso_score)) / (np.max(iso_score) - np.min(iso_score) + 1e-8)
ell_normalized = (ell_score - np.min(ell_score)) / (np.max(ell_score) - np.min(ell_score) + 1e-8)
# 综合评分（取较保守的值）
health_score = min(iso_normalized, ell_normalized) * 100
return health_score, {
'iso_score': iso_score,
'ell_score': ell_score,
'features': features
}
四、异常检测：三种主流方法实战对比
4.1 统计方法：简单但有效
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller, acf, pacf
class StatisticalDetector:
"""
基于统计方法的异常检测
适合：趋势性故障、缓慢漂移
"""
def __init__(self, window_size=100, confidence=3):
self.window_size = window_size
self.confidence = confidence # 通常用3σ原则
def detect_anomaly_statistical(self, data):
"""
使用滚动窗口统计进行异常检测
"""
data = np.array(data)
# 计算滚动统计量
rolling_mean = pd.Series(data).rolling(window=self.window_size).mean()
rolling_std = pd.Series(data).rolling(window=self.window_size).std()
# 计算残差（偏离程度）
z_scores = np.abs((data - rolling_mean) / (rolling_std + 1e-8))
# 检测异常点
anomalies = z_scores > self.confidence
return {
'anomalies': anomalies,
'z_scores': z_scores,
'rolling_mean': rolling_mean,
'rolling_std': rolling_std
}
def detect_trend_anomaly(self, data):
"""
使用线性回归检测趋势异常
设备磨损通常表现为缓慢上升的趋势
"""
x = np.arange(len(data))
# 拟合线性趋势
coeffs = np.polyfit(x, data, 1)
trend = np.polyval(coeffs, x)
# 计算残差
residuals = data - trend
# 如果趋势斜率异常大，可能是故障前兆
slope_threshold = 0.1 # 根据具体设备调整
trend_anomaly = np.abs(coeffs[0]) > slope_threshold
return {
'trend': trend,
'residuals': residuals,
'slope': coeffs[0],
'is_anomaly': trend_anomaly
}
def detect_level_shift(self, data, change_points=None):
"""
检测数据分布的突变点
对应：突发故障、传感器故障
"""
if change_points is None:
# 使用CUSUM算法检测突变
change_points = self._cusum_detection(data)
anomalies = []
for cp in change_points:
if cp > 0 and cp < len(data) - 1:
before = data[:cp]
after = data[cp:]
# 比较两段的均值差异
mean_diff = np.abs(np.mean(after) - np.mean(before))
pooled_std = np.sqrt(
(np.var(before) * len(before) + np.var(after) * len(after))
/ (len(before) + len(after))
)
# 标准化差异
effect_size = mean_diff / (pooled_std + 1e-8)
if effect_size > 2: # Cohen's d > 2 为大效应
anomalies.append({
'point': cp,
'effect_size': effect_size,
'before_mean': np.mean(before),
'after_mean': np.mean(after)
})
return anomalies
def _cusum_detection(self, data, threshold=4, kappa=0.5):
"""
CUSUM（累积和）算法检测突变
"""
data = np.array(data)
mean = np.mean(data)
std = np.std(data)
# 标准化
normalized = (data - mean) / std
# 累积和
s_pos = np.zeros(len(normalized))
s_neg = np.zeros(len(normalized))
for i in range(1, len(normalized)):
s_pos[i] = max(0, s_pos[i-1] + normalized[i] - kappa)
s_neg[i] = max(0, s_neg[i-1] - normalized[i] - kappa)
# 检测阈值
change_points = []
for i in range(len(s_pos)):
if s_pos[i] > threshold or s_neg[i] > threshold:
change_points.append(i)
s_pos[i:] = 0
s_neg[i:] = 0
return change_points
4.2 机器学习方法：捕捉非线性模式
from sklearn.ensemble import IsolationForest, RandomForestClassifier
from sklearn.svm import OneClassSVM
from sklearn.metrics import classification_report, confusion_matrix
import joblib
class MachineLearningDetector:
"""
基于机器学习的异常检测器
适合：复杂故障模式、多传感器融合
"""
def __init__(self):
self.models = {
'isolation_forest': IsolationForest(
contamination=0.05,
n_estimators=200,
max_samples='auto',
random_state=42
),
'one_class_svm': OneClassSVM(
kernel='rbf',
gamma='scale',
nu=0.05
)
}
self.feature_names = None
self.scaler = None
def train(self, X_train, y_train=None):
"""
训练异常检测模型
X_train: 特征矩阵 (n_samples, n_features)
y_train: 标签 (可选，None表示无监督)
"""
from sklearn.preprocessing import StandardScaler
# 标准化
self.scaler = StandardScaler()
X_scaled = self.scaler.fit_transform(X_train)
# 保存特征名
if hasattr(X_train, 'columns'):
self.feature_names = X_train.columns.tolist()
else:
self.feature_names = [f'feature_{i}' for i in range(X_train.shape[1])]
# 训练孤立森林
self.models['isolation_forest'].fit(X_scaled)
# 如果有标签，训练分类器
if y_train is not None:
rf = RandomForestClassifier(n_estimators=100, random_state=42)
rf.fit(X_scaled, y_train)
self.models['random_forest'] = rf
print(f"✓ 模型训练完成")
print(f" 训练样本数: {X_train.shape[0]}")
print(f" 特征数: {X_train.shape[1]}")
return self
def predict_anomaly(self, X_test):
"""
预测异常，返回异常概率
"""
X_scaled = self.scaler.transform(X_test)
# 孤立森林得分（越低越异常）
iso_scores = self.models['isolation_forest'].score_samples(X_scaled)
# 转换为异常概率（0-1，越高越异常）
anomaly_prob = 1 / (1 + np.exp(-iso_scores)) # sigmoid转换
predictions = self.models['isolation_forest'].predict(X_scaled)
return {
'anomaly_probability': anomaly_prob,
'predictions': predictions, # -1为异常，1为正常
'iso_scores': iso_scores
}
def predict_with_confidence(self, X_test):
"""
带置信度的预测，集成多个模型
"""
results = self.predict_anomaly(X_test)
# 计算置信度（基于模型一致性）
if 'random_forest' in self.models:
rf_proba = self.models['random_forest'].predict_proba(X_test)[:, 1]
ensemble_proba = (results['anomaly_probability'] + rf_proba) / 2
else:
ensemble_proba = results['anomaly_probability']
return {
'anomaly_probability': ensemble_proba,
'is_anomaly': ensemble_proba > 0.7, # 阈值可调
'confidence': np.abs(ensemble_proba - 0.5) * 2 # 0-1置信度
}
4.3 深度学习方法：时序模式的”翻译官”
import torch
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader
class TimeSeriesDataset(Dataset):
"""时间序列数据集"""
def __init__(self, data, sequence_length=64, stride=32):
self.data = data
self.seq_len = sequence_length
self.stride = stride
def __len__(self):
return (len(self.data) - self.seq_len) // self.stride + 1
def __getitem__(self, idx):
start = idx * self.stride
end = start + self.seq_len
return torch.tensor(self.data[start:end], dtype=torch.float32)
class LSTM_Autoencoder(nn.Module):
"""
LSTM自编码器用于异常检测
核心思想：学习正常数据的"压缩表示"，重建误差大则为异常
"""
def __init__(self, input_dim=1, hidden_dim=64, num_layers=2, seq_len=64):
super().__init__()
self.seq_len = seq_len
# 编码器
self.encoder = nn.LSTM(
input_size=input_dim,
hidden_size=hidden_dim,
num_layers=num_layers,
batch_first=True,
dropout=0.2 if num_layers > 1 else 0
)
# 解码器
self.decoder = nn.LSTM(
input_size=hidden_dim,
hidden_size=hidden_dim,
num_layers=num_layers,
batch_first=True,
dropout=0.2 if num_layers > 1 else 0
)
self.output_layer = nn.Linear(hidden_dim, input_dim)
def forward(self, x):
# 编码
_, (hidden, cell) = self.encoder(x)
# 解码（使用编码后的隐状态）
# 创建一个学习到的初始隐状态序列
decoder_input = hidden[-1].unsqueeze(0).repeat(self.seq_len, 1, 1).permute(1, 0, 2)
decoder_output, _ = self.decoder(decoder_input, (hidden, cell))
# 输出
reconstruction = self.output_layer(decoder_output)
return reconstruction
def encode(self, x):
"""只编码，用于特征提取"""
_, (hidden, cell) = self.encoder(x)
return hidden[-1] # 最后层的隐状态作为特征
class VariationalLSTM(nn.Module):
"""
变分LSTM：不仅学习正常数据的分布，还能生成"正常"样本
对轻微异常更敏感
"""
def __init__(self, input_dim=1, hidden_dim=64, latent_dim=32):
super().__init__()
self.input_dim = input_dim
self.hidden_dim = hidden_dim
self.latent_dim = latent_dim
# 编码部分
self.lstm = nn.LSTM(input_dim, hidden_dim, batch_first=True)
self.fc_mu = nn.Linear(hidden_dim, latent_dim)
self.fc_logvar = nn.Linear(hidden_dim, latent_dim)
# 解码部分
self.decoder_lstm = nn.LSTM(latent_dim, hidden_dim, batch_first=True)
self.output_layer = nn.Linear(hidden_dim, input_dim)
def encode(self, x):
_, (hidden, _) = self.lstm(x)
mu = self.fc_mu(hidden)
logvar = self.fc_logvar(hidden)
return mu, logvar
def reparameterize(self, mu, logvar):
"""重参数化技巧"""
std = torch.exp(0.5 * logvar)
epsilon = torch.randn_like(std)
return mu + epsilon * std
def decode(self, z, seq_len):
z = z.unsqueeze(0).repeat(seq_len, 1, 1).permute(1, 0, 2)
output, _ = self.decoder_lstm(z)
return self.output_layer(output)
def forward(self, x):
mu, logvar = self.encode(x)
z = self.reparameterize(mu, logvar)
reconstruction = self.decode(z, x.shape[1])
return reconstruction, mu, logvar
def calculate_vae_loss(reconstruction, x, mu, logvar):
"""VAE损失：重建误差 + KL散度"""
# 重建误差（MSE）
recon_loss = nn.MSELoss()(reconstruction, x)
# KL散度（正则化项，让潜变量接近标准正态分布）
kl_loss = -0.5 * torch.mean(1 + logvar - mu.pow(2) - logvar.exp())
return recon_loss + kl_loss
class DeepAnomalyDetector:
"""
深度学习异常检测器
"""
def __init__(self, input_dim=1, hidden_dim=64, seq_len=64, learning_rate=1e-3):
self.input_dim = input_dim
self.hidden_dim = hidden_dim
self.seq_len = seq_len
self.learning_rate = learning_rate
self.model = LSTM_Autoencoder(
input_dim=input_dim,
hidden_dim=hidden_dim,
seq_len=seq_len
)
self.optimizer = torch.optim.Adam(self.model.parameters(), lr=learning_rate)
self.device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
self.model.to(self.device)
self.normal_threshold = None
def train(self, train_data, epochs=50, batch_size=64):
"""训练异常检测模型"""
dataset = TimeSeriesDataset(train_data, seq_len=self.seq_len)
dataloader = DataLoader(dataset, batch_size=batch_size, shuffle=True)
self.model.train()
losses = []
for epoch in range(epochs):
epoch_loss = 0
for batch in dataloader:
batch = batch.to(self.device)
self.optimizer.zero_grad()
reconstruction = self.model(batch)
loss = nn.MSELoss()(reconstruction, batch)
loss.backward()
self.optimizer.step()
epoch_loss += loss.item()
avg_loss = epoch_loss / len(dataloader)
losses.append(avg_loss)
if (epoch + 1) % 10 == 0:
print(f"Epoch {epoch+1}/{epochs}, Loss: {avg_loss:.6f}")
# 计算正常阈值（基于训练数据的重建误差分布）
self.model.eval()
with torch.no_grad():
# 在训练数据上评估
all_errors = []
for batch in dataloader:
batch = batch.to(self.device)
reconstruction = self.model(batch)
errors = torch.mean((batch - reconstruction)**2, dim=(1, 2)).cpu().numpy()
all_errors.extend(errors)
# 设置阈值为95分位数
self.normal_threshold = np.percentile(all_errors, 95)
return losses
def detect(self, data):
"""
检测异常
返回：(是否异常, 异常分数, 重建误差)
"""
self.model.eval()
# 准备数据
dataset = TimeSeriesDataset(data, seq_len=self.seq_len)
dataloader = DataLoader(dataset, batch_size=64)
all_errors = []
with torch.no_grad():
for batch in dataloader:
batch = batch.to(self.device)
reconstruction = self.model(batch)
# 逐点重建误差
errors = torch.mean((batch - reconstruction)**2, dim=(1, 2)).cpu().numpy()
all_errors.extend(errors)
# 计算异常分数（滑动窗口平均，平滑噪声）
errors = np.array(all_errors)
if len(errors) > 10:
anomaly_scores = np.convolve(errors, np.ones(10)/10, mode='valid')
else:
anomaly_scores = errors
# 判断是否异常
is_anomaly = anomaly_scores > self.normal_threshold
return {
'is_anomaly': is_anomaly,
'anomaly_scores': anomaly_scores,
'threshold': self.normal_threshold,
'mean_error': np.mean(errors)
}
五、实战案例：主轴轴承早期故障预警
5.1 故障背景
去年秋天，厂里一台关键的主轴加工单元出了问题。轴承内圈出现了早期的点蚀（pitting），但振动信号非常微弱，常规的阈值报警根本没触发。如果等到报警再停修，主轴可能已经报废，停机时间会延长至少6小时。
这次我们用了新的时间序列分析方法，成功在故障发生前14天发出了预警。
5.2 实施过程
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
class BearingFaultDetector:
"""
主轴轴承故障检测器
专门针对轴承故障特征优化
"""
def __init__(self, sampling_rate=10000):
self.sampling_rate = sampling_rate
self.health_history = []
self.anomaly_timeline = []
def extract_envelope_features(self, vibration_signal):
"""
提取包络谱特征
轴承故障特征频率通常调制在载波信号上，
包络分析可以捕捉这些低频调制成分
"""
# 带通滤波（针对轴承故障频率范围）
filtered = self._bandpass_filter(vibration_signal, 1000, 5000)
# 希尔伯特变换提取包络
analytic_signal = np.fft.ifft(np.fft.rfft(filtered) * 2)
envelope = np.abs(analytic_signal[:len(filtered)])
# 包络谱
envelope_spectrum = np.fft.rfft(envelope)
envelope_magnitude = np.abs(envelope_spectrum)
frequencies = np.fft.rfftfreq(len(envelope), 1/self.sampling_rate)
# 提取轴承故障特征频率附近的能量
# 内圈故障频率 (BPFI) ≈ 3.6 × 转频
# 外圈故障频率 (BPFO) ≈ 2.4 × 转频
# 滚珠故障频率 (BSF) ≈ 0.4 × 转频
bearing_fault_frequencies = [
('BPFO_2x', 2 * 25), # 假设转频25Hz
('BPFO_4x', 4 * 25),
('BPFI_2x', 2 * 36),
('BPFI_4x', 4 * 36),
('BSF_2x', 2 * 10),
]
features = {}
total_energy = np.sum(envelope_magnitude**2)
for name, freq in bearing_fault_frequencies:
# 在特征频率±5Hz范围内积分能量
mask = (frequencies >= freq - 5) & (frequencies <= freq + 5)
band_energy = np.sum(envelope_magnitude[mask]**2)
features[name] = band_energy / (total_energy + 1e-8)
# 时域特征
features['rms'] = np.sqrt(np.mean(vibration_signal**2))
features['crest_factor'] = np.max(np.abs(vibration_signal)) / (features['rms'] + 1e-8)
features['impulse_factor'] = np.max(np.abs(vibration_signal)) / (np.mean(np.abs(vibration_signal)) + 1e-8)
return features
def _bandpass_filter(self, data, low_freq, high_freq):
"""简单的带通滤波"""
from scipy.signal import butter, filtfilt
nyquist = self.sampling_rate / 2
low = low_freq / nyquist
high = high_freq / nyquist
b, a = butter(4, [low, high], btype='band')
return filtfilt(b, a, data)
def detect_progressive_degradation(self, feature_history):
"""
检测渐进性退化
轴承故障通常是渐进发展的，不是突然跳变
"""
features_df = pd.DataFrame(feature_history)
results = {}
# 对每个特征进行时序分析
for col in features_df.columns:
series = features_df[col]
# 计算趋势（线性回归斜率）
x = np.arange(len(series))
slope, intercept, r_value, p_value, std_err = stats.linregress(x, series)
# 判断趋势是否显著
is_trending = p_value < 0.05 and abs(slope) > 0.001
# 计算趋势速度（每天变化率）
trend_rate = slope * 24 # 假设每小时采集一次
results[col] = {
'trend_slope': slope,
'p_value': p_value,
'is_trending': is_trending,
'trend_rate_per_day': trend_rate,
'current_value': series.iloc[-1],
'initial_value': series.iloc[0]
}
return results
def generate_health_report(self, current_features, degradation_analysis):
"""生成健康报告"""
report = {
'timestamp': datetime.now(),
'current_features': current_features,
'degradation_analysis': degradation_analysis,
'health_status': 'unknown'
}
# 综合判断健康状况
rising_features = sum(1 for v in degradation_analysis.values() if v['is_trending'])
total_features = len(degradation_analysis)
if rising_features == 0:
report['health_status'] = 'healthy'
elif rising_features <= 2:
report['health_status'] = 'warning'
else:
report['health_status'] = 'critical'
return report
5.3 实际数据展示
# 模拟实际监测数据
def simulate_bearing_degradation(days=30, interval_hours=1):
"""
模拟轴承退化过程
包括：基线振动、渐进退化、随机噪声
"""
np.random.seed(42)
total_points = days * 24 // interval_hours
timestamps = [datetime.now() - timedelta(hours=interval_hours*(total_points-1-i))
for i in range(total_points)]
# 基础振动水平
base_vibration = 0.5 # mm/s RMS
# 退化模型：指数增长（符合轴承实际退化规律）
degradation_rate = 0.02 # 每天增长2%
degradation_factor = np.exp(degradation_rate * np.arange(total_points) / 24)
# 添加周期性成分（与转速相关）
rotation_period = 24 // 10 # 每2.4小时一圈
periodic_component = 0.1 * np.sin(2 * np.pi * np.arange(total_points) / rotation_period)
# 添加随机噪声
noise = np.random.normal(0, 0.05, total_points)
# 组装信号
vibration = base_vibration * degradation_factor + periodic_component + noise
# 添加故障特征（在后期出现）
fault_features = np.zeros(total_points)
fault_start = int(total_points * 0.7) # 21天后开始
for i in range(fault_start, total_points):
# 轴承故障特征：周期性冲击
phase = (i % 10) / 10
if phase < 0.1: # 冲击信号
fault_features[i] = 0.3 * np.exp(-((phase - 0.05)**2) / 0.001)
vibration += fault_features
# 提取特征
features = {
'timestamp': timestamps,
'vibration_rms': vibration,
'trend': np.polyfit(range(len(vibration)), vibration, 1)[0]
}
return pd.DataFrame(features)
# 运行模拟
simulated_data = simulate_bearing_degradation(days=30)
# 可视化
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
# 原始振动信号
axes[0].plot(simulated_data['timestamp'], simulated_data['vibration_rms'],
linewidth=1, color='#2E86AB')
axes[0].set_ylabel('Vibration RMS (mm/s)')
axes[0].set_title('Main Spindle Bearing Vibration Trend (30 Days)')
axes[0].grid(True, alpha=0.3)
# 趋势线
axes[1].plot(simulated_data['timestamp'], simulated_data['vibration_rms'],
linewidth=1, color='#2E86AB', alpha=0.7)
# 添加趋势线
x = np.arange(len(simulated_data))
trend_line = np.polyfit(x, simulated_data['vibration_rms'], 1)
trend_plot = np.polyval(trend_line, x)
axes[1].plot(simulated_data['timestamp'], trend_line[0] * x + trend_line[1],
'--', color='#A23B72', linewidth=2, label='Trend Line')
axes[1].fill_between(simulated_data['timestamp'],
simulated_data['vibration_rms'],
simulated_data['vibration_rms'].mean(),
alpha=0.3)
axes[1].set_ylabel('Deviation from Mean')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# 包络谱特征能量（模拟）
np.random.seed(123)
bpfo_energy = 0.1 + 0.005 * np.arange(30*24) + np.random.normal(0, 0.02, 30*24)
bpfi_energy = 0.08 + 0.003 * np.arange(30*24) + np.random.normal(0, 0.015, 30*24)
timestamps_short = simulated_data['timestamp'][::24] # 每天一个点
axes[2].plot(timestamps_short, bpfo_energy[:30], 'o-',
label='BPFO Energy (Outer Race)', color='#F18F01')
axes[2].plot(timestamps_short, bpfi_energy[:30], 's-',
label='BPFI Energy (Inner Race)', color='#C73E1D')
axes[2].axhline(y=0.3, color='r', linestyle='--', label='Alert Threshold')
axes[2].set_ylabel('Feature Energy')
axes[2].set_xlabel('Date')
axes[2].legend(loc='upper left')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bearing_degradation_trend.png', dpi=150, bbox_inches='tight')
plt.show()
5.4 预警效果对比
传统方法 vs 时间序列过程识别对比
┌─────────────────┬──────────────────┬──────────────────┐
│ 指标 │ 传统阈值报警 │ 时间序列识别 │
├─────────────────┼──────────────────┼──────────────────┤
│ 最早预警时间 │ 故障发生后2小时 │ 故障前14天 │
│ 误报率 │ 15% │ 3% │
│ 漏报率 │ 8% │ 0.5% │
│ 可维护窗口 │ 紧急（立即停机） │ 计划性（2周内） │
│ 维修成本 │ 高（紧急+大修） │ 低（计划+更换） │
│ 非计划停机 │ 6小时/次 │ 0小时（计划内） │
└─────────────────┴──────────────────┴──────────────────┘
六、延长设备寿命的关键：从”坏了再修”到”该修再修”
6.1 剩余使用寿命（RUL）预测
时间序列过程识别的价值不仅在于检测异常，更在于预测还能用多久。
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import TimeSeriesSplit
import numpy as np
class RemainingUsefulLifePredictor:
"""
剩余使用寿命预测器
基于退化轨迹预测设备还能健康运行多久
"""
def __init__(self, degradation_threshold=1.5):
"""
degradation_threshold: 退化阈值（相对于初始值的倍数）
超过此值认为设备需要维护
"""
self.degradation_threshold = degradation_threshold
self.model = GradientBoostingRegressor(
n_estimators=100,
max_depth=4,
learning_rate=0.1,
random_state=42
)
self.is_fitted = False
def prepare_training_data(self, health_scores, rul_labels):
"""
准备训练数据
health_scores: 健康评分序列
rul_labels: 对应的剩余使用寿命（小时）
"""
# 构造特征：滑动窗口统计
features_list = []
target_list = []
window_size = 100
for i in range(window_size, len(health_scores)):
window = health_scores[i-window_size:i]
# 从窗口提取特征
features = {
'current_value': window[-1],
'mean_10': np.mean(window[-10:]),
'std_10': np.std(window[-10:]),
'trend_10': self._calculate_trend(window[-10:]),
'mean_50': np.mean(window[-50:]),
'std_50': np.std(window[-50:]),
'trend_50': self._calculate_trend(window[-50:]),
'min_50': np.min(window[-50:]),
'max_50': np.max(window[-50:]),
}
features_list.append(features)
target_list.append(rul_labels[i])
self.feature_names = list(features_list[0].keys())
X = pd.DataFrame(features_list)
y = pd.Series(target_list)
return X, y
def _calculate_trend(self, data):
"""计算线性趋势"""
if len(data) < 2:
return 0
x = np.arange(len(data))
coeffs = np.polyfit(x, data, 1)
return coeffs[0]
def fit(self, X, y):
"""训练模型"""
self.model.fit(X, y)
self.is_fitted = True
print(f"✓ RUL预测模型训练完成")
print(f" 训练样本: {len(y)}")
print(f" 特征数量: {len(self.feature_names)}")
return self
def predict_rul(self, recent_observations, time_step_hours=1):
"""
预测剩余使用寿命
recent_observations: 最近的健康评分序列
time_step_hours: 每个数据点的时间间隔（小时）
"""
if not self.is_fitted:
raise ValueError("模型尚未训练")
# 准备特征
features = {
'current_value': recent_observations[-1],
'mean_10': np.mean(recent_observations[-10:]),
'std_10': np.std(recent_observations[-10:]),
'trend_10': self._calculate_trend(recent_observations[-10:]),
'mean_50': np.mean(recent_observations[-50:]),
'std_50': np.std(recent_observations[-50:]),
'trend_50': self._calculate_trend(recent_observations[-50:]),
'min_50': np.min(recent_observations[-50:]),
'max_50': np.max(recent_observations[-50:]),
}
X = pd.DataFrame([features], columns=self.feature_names)
# 预测RUL
predicted_rul = self.model.predict(X)[0]
return {
'predicted_rul_hours': predicted_rul,
'predicted_rul_days': predicted_rul / 24,
'predicted_failure_date': datetime.now() + timedelta(hours=predicted_rul),
'confidence': self._estimate_confidence(recent_observations)
}
def _estimate_confidence(self, observations):
"""估计预测置信度"""
# 基于近期数据的稳定性
recent_std = np.std(observations[-20:]) if len(observations) >= 20 else np.std(observations)
recent_mean = np.mean(observations[-20:]) if len(observations) >= 20 else np.mean(observations)
# 变异系数越小，置信度越高
cv = recent_std / (recent_mean + 1e-8)
confidence = max(0, min(1, 1 - cv * 10))
return confidence
6.2 实际运维决策支持
def generate_maintenance_recommendation(rul_prediction, equipment_name, history):
"""
生成维护建议
"""
rul_hours = rul_prediction['predicted_rul_hours']
confidence = rul_prediction['confidence']
recommendations = []
if rul_hours < 24:
recommendations.append({
'priority': 'IMMEDIATE',
'action': '立即安排停机检修',
'reason': f'预测剩余寿命仅{rul_hours:.1f}小时',
'备选方案': '准备备件，可在48小时内完成轮换'
})
elif rul_hours < 72:
recommendations.append({
'priority': 'HIGH',
'action': '本周内安排计划性维护',
'reason': f'预测剩余寿命{rul_hours:.1f}小时',
'建议': '结合生产计划，在下个周末停机时更换'
})
elif rul_hours < 168:
recommendations.append({
'priority': 'MEDIUM',
'action': '本月内安排维护',
'reason': f'预测剩余寿命{rul_hours:.1f}小时（约{rul_hours/24:.0f}天）',
'建议': '采购备件，监控趋势，必要时提前维护'
})
else:
recommendations.append({
'priority': 'LOW',
'action': '继续监控，按计划维护',
'reason': f'预测剩余寿命{rul_hours:.1f}小时',
'建议': '常规巡检，无需紧急处理'
})
# 添加历史对比
if len(history) > 7:
recent_trend = history[-7:].mean() - history[-14:-7].mean()
if recent_trend > 0.1:
recommendations[-1]['warning'] = '退化速度加快，建议提前维护'
return {
'equipment': equipment_name,
'prediction_time': datetime.now(),
'rul_hours': rul_hours,
'confidence': confidence,
'recommendations': recommendations
}
七、落地指南：你的工厂该怎么开始？
7.1 分阶段实施路径
第一阶段：基础建设（1-2个月）
- 在关键设备上安装振动传感器（至少三轴）
- 建立数据采集系统（边缘计算网关+云端存储）
- 积累至少3个月的健康运行数据
第二阶段：模型构建（2-3个月）
- 使用健康数据训练异常检测模型
- 定义”正常”和”异常”的判定标准
- 小范围试点，验证预警准确性
第三阶段：全面推广（3-6个月）
- 覆盖全厂关键设备
- 建立预测性维护工单流程
- 持续优化模型（用新数据重新训练）
7.2 关键成功因素
┌─────────────────────────────────────────────────────────────┐
│ 成功要素检查清单 │
├─────────────────────────────────────────────────────────────┤
│ □ 数据质量 > 算法复杂度（好数据比好算法更重要） │
│ □ 业务专家参与（设备工程师比数据科学家更懂设备） │
│ □ 持续迭代（模型需要定期用新数据重新训练） │
│ □ 人机结合（系统预警+人工确认，避免过度依赖算法） │
│ □ 成本效益考量（不是所有设备都需要预测性维护） │
└─────────────────────────────────────────────────────────────┘
7.3 ROI估算参考
以一台价值500万的关键数控机床为例：
传统模式：
- 每年非计划停机：2-3次，每次8小时
- 紧急维修成本：5-10万元/次
- 备件损失：无法批量采购，溢价30%
- 年损失估算：15-30万元
预测性维护模式（实施后）：
- 非计划停机：减少80%
- 维修计划化，成本降低40%
- 备件批量采购，节省20%
- 年节约：10-20万元
系统投资（传感器+软件+实施）：约30-50万元
投资回收期：2-3年
八、常见问题解答
Q: 数据采集频率多少合适？
A: 取决于你要检测的故障类型。振动检测轴承故障需要10kHz以上，温度监测1Hz就够。我的建议是”够用就好”，别盲目追求高频，存储和处理成本也会飙升。
Q: 模型误报了怎么办？
A: 这是真实存在的情况。我的经验是，宁可多报也别漏报——早期误报的成本远低于漏报导致的故障扩大。但可以通过设置”确认窗口”来过滤瞬态干扰，比如连续3次预警才发工单。
Q: 小工厂也能做吗？
A: 完全可以。现在开源工具很多，Python的tsfresh、pyts等库专门为时间序列特征提取设计。关键是从小处着手，先盯住最关键的一两台设备，做出效果再推广。
Q: 如何证明这套系统值这个投入？
A: 做试点对比。选两台同型号设备，一台装系统，一台不装，跑3个月看故障率差异。数据说话，老板最吃这一套。
九、写在最后
设备管理这件事，说到底是和”不确定性”打交道。时间序列过程识别不是水晶球，不能100%预测未来，但它能让我们从”被动救火”变成”主动防御”。
去年那台主轴的案例，后来我把它做成了培训教材。厂里年轻的技术员现在看到振动频谱图，第一反应不是”这正常吗”，而是”这个特征频率对应什么故障”。这种思维转变，才是这套系统最大的价值。
如果你也在考虑给设备装”听诊器”，我的建议是：先开始，再完善。没有完美的模型，只有不断迭代的过程。今天的10分系统，加上持续优化，明年可能就是90分。
有什么具体问题，欢迎随时交流。做这行十几年，碰到的坑比写出来的字还多，能帮一把是一把。