4. 特征工程（上）：时域特征提取与频域特征提取

各位同学，欢迎来到特征工程环节。说实话，很多做寿命预测的朋友，模型选得再好，数据预处理做得再漂亮，如果特征没提好，那基本就是白忙活。我见过太多人，上来就扔一堆原始数据给模型，结果模型学了半天，学到的全是噪声。

特征工程，说白了就是帮模型“划重点”。今天我们先聊时域和频域这两大块。这是最基础、最常用的特征提取手段，也是我入行时最先啃透的部分。

一、时域特征提取：从波形里“抠”信息

时域特征，就是直接从原始信号的时间序列里算出来的统计量。你想想看，一个轴承从健康到故障，它的振动波形会怎么变？振幅会变大，波形会变得“毛糙”，偶尔还会出现很大的冲击。这些变化，都能用一些数值指标捕捉到。

1. 均值与方差

均值反映信号的直流分量或整体趋势。对于旋转机械，如果均值发生漂移，往往意味着传感器出了问题，或者设备发生了结构性变化。

方差反映信号的波动程度。设备磨损越严重，振动幅度通常越大，方差自然就上去了。

import numpy as np

def extract_time_features(signal):
    """
    提取时域特征
    signal: 一维数组，振动信号
    """
    mean_val = np.mean(signal)          # 均值
    var_val = np.var(signal)            # 方差
    rms_val = np.sqrt(np.mean(signal**2))  # 均方根值
    
    return mean_val, var_val, rms_val

# 举个栗子
vibration_data = np.random.randn(1000) * 0.5 + 10  # 模拟振动信号
mean, var, rms = extract_time_features(vibration_data)
print(f"均值: {mean:.3f}, 方差: {var:.3f}, 均方根: {rms:.3f}")

2. 峰值与峭度

峰值就是信号的最大绝对值。这个指标对早期故障特别敏感——设备刚开始出现裂纹时，偶尔会蹦出一个很大的冲击，峰值一下子就上去了。

峭度（Kurtosis）是个很有意思的指标。它衡量信号分布的“尖锐程度”。健康设备的振动信号接近正态分布，峭度值在3左右。一旦出现故障，信号中会出现很多“尖刺”，峭度值就会飙升。

def kurtosis(signal):
    """
    计算峭度
    """
    n = len(signal)
    mean = np.mean(signal)
    fourth_moment = np.sum((signal - mean)**4) / n
    variance = np.var(signal)
    
    # 峭度 = 四阶中心矩 / 方差^2
    kurt = fourth_moment / (variance**2)
    return kurt

# 模拟一个带冲击的信号
healthy_signal = np.random.randn(1000)  # 健康信号，峭度≈3
faulty_signal = np.copy(healthy_signal)
faulty_signal[::50] += 5  # 每50个点加一个冲击

print(f"健康信号峭度: {kurtosis(healthy_signal):.2f}")
print(f"故障信号峭度: {kurtosis(faulty_signal):.2f}")

3. 其他常用时域特征

除了上面这些，还有几个指标我经常用：

特征名称	计算公式	物理含义
波形因子	RMS / 整流平均值	反映波形形状，故障时通常增大
峰值因子	峰值 / RMS	检测冲击性故障，早期故障敏感
脉冲因子	峰值 / 整流平均值	类似峰值因子，但更稳定
裕度因子	峰值 / 方根幅值	对磨损类故障敏感

二、频域特征提取：换个角度看信号

时域特征虽然直观，但有个问题——它把不同频率成分混在一起了。你想想看，一个轴承故障，它的特征频率可能是几百赫兹；而齿轮啮合频率可能是几千赫兹。时域上根本分不清。这时候，就得请出频域分析这个神器了。

1. FFT（快速傅里叶变换）

FFT是频域分析的基石。它把时域信号分解成不同频率的正弦波分量。说白了，就是告诉你“信号里都有哪些频率的成分，每个成分有多强”。

import numpy as np
import matplotlib.pyplot as plt

def compute_fft(signal, sampling_rate):
    """
    计算FFT并返回频率和幅值
    """
    n = len(signal)
    fft_vals = np.fft.fft(signal)
    fft_abs = np.abs(fft_vals[:n//2])  # 取正频率部分
    freqs = np.fft.fftfreq(n, 1/sampling_rate)[:n//2]
    
    return freqs, fft_abs

# 模拟一个包含50Hz和120Hz成分的信号
fs = 1000  # 采样率1000Hz
t = np.arange(0, 1, 1/fs)
signal = 2 * np.sin(2*np.pi*50*t) + 1.5 * np.sin(2*np.pi*120*t)

freqs, magnitude = compute_fft(signal, fs)

# 找到幅值最大的前3个频率
top_indices = np.argsort(magnitude)[-3:][::-1]
print("主要频率成分:")
for idx in top_indices:
    print(f"  {freqs[idx]:.1f} Hz, 幅值: {magnitude[idx]:.2f}")

2. 功率谱密度（PSD）

PSD比FFT更进一步。它描述的是信号功率在频率上的分布。对于随机振动信号，PSD比FFT的幅值谱更稳定、更有统计意义。

常用的PSD估计方法有：

• 周期图法：直接对信号做FFT然后取平方。简单，但方差大。
• Welch法：把信号分段，每段加窗，然后平均。方差小，是我最推荐的方法。
• 自相关法：先算自相关函数，再做FFT。适合短数据。

from scipy import signal as scipy_signal

def compute_psd_welch(signal, fs, nperseg=256):
    """
    使用Welch法计算功率谱密度
    """
    freqs, psd = scipy_signal.welch(
        signal, 
        fs=fs, 
        nperseg=nperseg,
        window='hann'
    )
    return freqs, psd

# 计算PSD
freqs_psd, psd_vals = compute_psd_welch(signal, fs)

# 提取频域特征
# 1. 频谱峰值（最大功率对应的频率）
peak_freq = freqs_psd[np.argmax(psd_vals)]

# 2. 重心频率（功率加权平均频率）
centroid = np.sum(freqs_psd * psd_vals) / np.sum(psd_vals)

# 3. 频带能量（比如0-200Hz的能量）
band_energy = np.trapz(psd_vals[(freqs_psd >= 0) & (freqs_psd <= 200)], 
                       freqs_psd[(freqs_psd >= 0) & (freqs_psd <= 200)])

print(f"峰值频率: {peak_freq:.1f} Hz")
print(f"重心频率: {centroid:.1f} Hz")
print(f"0-200Hz频带能量: {band_energy:.3f}")

3. 频域特征的实际应用

在实际项目中，我通常会提取以下频域特征：

1. 频谱峰值及其对应频率：直接反映主要振动源的频率。比如轴承外圈故障频率。
2. 重心频率：反映能量集中的频段。设备退化时，重心频率往往会向高频移动。
3. 频率方差：反映频谱的分散程度。故障越严重，频谱越“散”。
4. 特定频带能量：比如提取1倍频、2倍频、3倍频附近的能量。这对齿轮箱故障诊断特别有用。

三、时域与频域特征的融合策略

说了这么多，你可能要问：到底该用哪些特征？我的建议是：

• 初期探索：把所有时域和频域特征都算出来，做个相关性分析，剔除高度相关的特征。
• 特征筛选：用随机森林或XGBoost的特征重要性，选出Top 10-20个特征。
• 最终组合：保留5-8个时域特征 + 5-8个频域特征，构成最终的特征集。

记住，特征不是越多越好。我见过有人提取了100多个特征，结果模型过拟合得一塌糊涂。少而精，才是王道。

---