2. 齿轮箱故障诊断的时域特征提取与分类

任务说明

本案例使用齿轮箱振动加速度信号，通过滑动窗口截取信号片段，提取均值、方差、均方根、峰值因子、峭度、波形因子、脉冲因子和裕度因子共8个时域统计特征，构建特征向量后使用随机森林分类器识别正常、齿面磨损、断齿和轴承故障4种状态。数据集由本案例代码自动生成，包含4类模拟振动信号各200个样本，其中80%用于训练，20%用于测试。

时域特征定义

设信号片段为 \(x_1, x_2, \ldots, x_N\)，各特征计算方式如下：

均值：\(\bar{x} = \frac{1}{N}\sum_{i=1}^{N}x_i\)
方差：\(\sigma^2 = \frac{1}{N-1}\sum_{i=1}^{N}(x_i-\bar{x})^2\)
均方根：\(X_{rms} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}x_i^2}\)
峰值因子：\(C_f = X_{peak} / X_{rms}\)，其中 \(X_{peak} = \max|x_i|\)
峭度：\(K = \frac{1}{N}\sum_{i=1}^{N}(x_i-\bar{x})^4 / \sigma^4\)
波形因子：\(S_f = X_{rms} / ( \frac{1}{N}\sum_{i=1}^{N}|x_i| )\)
脉冲因子：\(I_f = X_{peak} / ( \frac{1}{N}\sum_{i=1}^{N}|x_i| )\)
裕度因子：\(CL_f = X_{peak} / ( \frac{1}{N}\sum_{i=1}^{N}\sqrt{|x_i|} )^2\)

代码实现

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, classification_report

# ---------- 1. 模拟齿轮箱振动信号 ----------
def simulate_gear_signal(fault_type, n_samples=200, fs=8192, duration=1.0):
    """生成4类齿轮箱振动信号，每类n_samples个样本"""
    t = np.linspace(0, duration, int(fs*duration), endpoint=False)
    shaft_freq = 30.0  # 轴频 30 Hz
    gear_mesh_freq = 30.0 * 23  # 啮合频率 690 Hz (23齿)
    signals = []
    labels = []
    for _ in range(n_samples):
        # 基础啮合振动
        x = 1.0 * np.sin(2*np.pi*gear_mesh_freq*t + 0.5*np.sin(2*np.pi*shaft_freq*t))
        # 随机噪声
        x += 0.15 * np.random.randn(len(t))
        if fault_type == 'normal':
            label = 0
        elif fault_type == 'wear':
            # 齿面磨损：啮合频率二次谐波幅值增大
            x += 0.6 * np.sin(2*np.pi*2*gear_mesh_freq*t)
            label = 1
        elif fault_type == 'broken_tooth':
            # 断齿：每转一周产生一次冲击
            impact = 2.5 * np.exp(-2000 * (t % (1/shaft_freq)))
            x += impact
            label = 2
        elif fault_type == 'bearing_fault':
            # 轴承故障：高频共振调制
            bearing_freq = 105.0  # 轴承故障特征频率
            x += 1.2 * np.sin(2*np.pi*3000*t) * (0.5 + 0.5*np.sin(2*np.pi*bearing_freq*t))
            label = 3
        else:
            raise ValueError('Unknown fault type')
        signals.append(x)
        labels.append(label)
    return np.array(signals), np.array(labels)

# ---------- 2. 时域特征提取 ----------
def extract_time_features(signal):
    """从单个信号片段提取8个时域统计特征"""
    N = len(signal)
    mean_val = np.mean(signal)
    var_val = np.var(signal, ddof=1)
    rms_val = np.sqrt(np.mean(signal**2))
    peak_val = np.max(np.abs(signal))
    # 峭度
    if var_val < 1e-12:
        kurtosis_val = 0.0
    else:
        kurtosis_val = np.mean((signal - mean_val)**4) / (var_val**2)
    # 整流均值
    abs_mean = np.mean(np.abs(signal))
    if abs_mean < 1e-12:
        waveform_factor = 0.0
        impulse_factor = 0.0
        margin_factor = 0.0
    else:
        waveform_factor = rms_val / abs_mean
        impulse_factor = peak_val / abs_mean
        # 裕度因子
        sqrt_abs_mean = np.mean(np.sqrt(np.abs(signal)))**2
        margin_factor = peak_val / sqrt_abs_mean if sqrt_abs_mean > 1e-12 else 0.0
    return np.array([mean_val, var_val, rms_val, peak_val,
                     kurtosis_val, waveform_factor, impulse_factor, margin_factor])

# ---------- 3. 构建数据集 ----------
fault_types = ['normal', 'wear', 'broken_tooth', 'bearing_fault']
all_features = []
all_labels = []
for ft in fault_types:
    sigs, labs = simulate_gear_signal(ft)
    for sig, lab in zip(sigs, labs):
        feat = extract_time_features(sig)
        all_features.append(feat)
        all_labels.append(lab)

X = np.array(all_features)
y = np.array(all_labels)
print(f'特征矩阵形状: {X.shape}, 标签分布: {np.bincount(y)}')

# ---------- 4. 训练与评估 ----------
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y)

clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)

# ---------- 5. 输出评估结果 ----------
acc = accuracy_score(y_test, y_pred)
print(f'测试集准确率: {acc:.4f}')
print('\n分类报告:')
print(classification_report(y_test, y_pred,
      target_names=fault_types, zero_division=0))

# ---------- 6. 特征重要性排序 ----------
feature_names = ['均值', '方差', '均方根', '峰值', '峭度',
                 '波形因子', '脉冲因子', '裕度因子']
importance = clf.feature_importances_
print('\n特征重要性:')
for name, imp in sorted(zip(feature_names, importance),
                        key=lambda x: x[1], reverse=True):
    print(f'  {name}: {imp:.4f}')

落地与运行方式

本案例依赖 numpy 和 scikit-learn，安装命令：

pip install numpy scikit-learn

将上述代码保存为 gear_fault_diagnosis.py，在终端执行：

python gear_fault_diagnosis.py

验证方法与判定标准

程序运行后应依次输出特征矩阵形状（800行×8列）、4类标签各200个的分布信息、测试集准确率、分类报告以及特征重要性排序。判定标准：

特征矩阵形状为 (800, 8)，标签分布为每类200个；
测试集准确率不低于0.90（由于模拟信号类间差异明显，随机森林通常能达到0.95以上）；
分类报告中每个类别的精确率、召回率和F1值均不低于0.85；
特征重要性中峭度和峰值因子应排在前列，因为断齿冲击和轴承故障调制会显著改变信号峰值分布。

若准确率低于0.90，可检查模拟信号参数是否被修改，或增加随机森林的n_estimators值后重新运行。

信号特征与分类逻辑说明

本案例的模拟信号设计使4类故障在时域特征上具有可区分性：正常状态信号平稳，峭度接近3；齿面磨损引入二次谐波使波形因子和均方根增大；断齿产生周期性冲击，显著提高峰值因子和峭度；轴承故障的高频调制使方差和脉冲因子上升。随机森林基于这些特征差异构建决策边界，实现故障类型识别。

📚 课程目录