4. 基于LSTM的剩余使用寿命预测

任务定义

本案例针对航空发动机退化过程，构建基于LSTM的剩余使用寿命（RUL）预测模型。输入为多传感器历史退化序列，输出为当前时刻至失效的剩余循环数。采用CMAPSS数据集FD001子集作为标准验证基准，完成从数据划分、模型构建到预测评估的完整闭环。

数据准备与预处理

FD001包含100台训练发动机单元和100台测试发动机单元，每台均记录21个传感器信号。训练集提供完整退化至失效的寿命标签；测试集截断于失效前某时刻，需预测截断点至失效的剩余循环数。本案例采用全部21个传感器通道，不做特征筛选，以验证LSTM对高维退化模式的自动提取能力。

预处理包含三步：1) 对每个传感器信号进行z-score标准化，均值和标准差从训练集统计获得并应用于测试集；2) 使用滑动窗口将多变量时间序列切分为固定长度样本，窗口长度设为30个时间步；3) 训练标签为窗口末端对应的真实RUL值，测试标签用于最终评估。RUL采用分段线性上限150，即超过150的RUL统一截断为150，以降低早期健康阶段标签差异过大的影响。

模型实现

采用Keras构建两层LSTM回归网络。输入形状为(时间步, 特征数)=(30, 21)，第一层LSTM含64个单元并返回序列，第二层含32个单元，后接全连接层输出RUL预测值。使用Adam优化器（学习率0.001）和均方误差损失函数。训练设置50个epoch，批大小512，并采用早停机制（patience=5）防止过拟合。

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.callbacks import EarlyStopping

# ---------- 1. 数据加载与RUL标签生成 ----------
def load_cmapss(train_path, test_path, max_rul=150):
    """加载FD001数据并生成RUL标签"""
    # 训练集: 列0=单元ID, 列1=循环数, 列2-4=操作设置, 列5-25=传感器
    train_df = pd.read_csv(train_path, sep=' ', header=None)
    train_df = train_df.dropna(axis=1, how='all')
    train_df.columns = ['id', 'cycle'] + [f'op{i}' for i in range(1,4)] + \
                       [f's{i}' for i in range(1,22)]
    
    # 测试集
    test_df = pd.read_csv(test_path, sep=' ', header=None)
    test_df = test_df.dropna(axis=1, how='all')
    test_df.columns = ['id', 'cycle'] + [f'op{i}' for i in range(1,4)] + \
                      [f's{i}' for i in range(1,22)]
    
    # 测试集RUL真值
    rul_df = pd.read_csv('RUL_FD001.txt', sep=' ', header=None)
    rul_df = rul_df.dropna(axis=1, how='all')
    test_rul = rul_df[0].values
    
    # 训练集RUL: 每台最大循环数 - 当前循环数
    train_rul = []
    for unit_id in train_df['id'].unique():
        unit_data = train_df[train_df['id'] == unit_id]
        max_cycle = unit_data['cycle'].max()
        rul = max_cycle - unit_data['cycle'].values
        train_rul.extend(rul)
    train_df['rul'] = np.clip(train_rul, 0, max_rul)
    
    # 测试集RUL: 截断点至失效的剩余循环
    test_rul_clipped = np.clip(test_rul, 0, max_rul)
    
    return train_df, test_df, test_rul_clipped

# ---------- 2. 标准化与滑窗切分 ----------
def prepare_sequences(df, scaler=None, window_size=30, is_train=True):
    """生成滑窗样本和对应标签"""
    sensor_cols = [f's{i}' for i in range(1,22)]
    
    if is_train:
        scaler = StandardScaler()
        scaler.fit(df[sensor_cols])
    
    X, y = [], []
    for unit_id in df['id'].unique():
        unit_data = df[df['id'] == unit_id].sort_values('cycle')
        scaled = scaler.transform(unit_data[sensor_cols])
        rul = unit_data['rul'].values
        
        for i in range(len(unit_data) - window_size + 1):
            X.append(scaled[i:i+window_size])
            y.append(rul[i+window_size-1])
    
    return np.array(X), np.array(y), scaler

# ---------- 3. 构建LSTM模型 ----------
def build_lstm(input_shape):
    model = Sequential([
        LSTM(64, return_sequences=True, input_shape=input_shape),
        Dropout(0.2),
        LSTM(32, return_sequences=False),
        Dropout(0.2),
        Dense(1, activation='linear')
    ])
    model.compile(optimizer=Adam(learning_rate=0.001), loss='mse')
    return model

# ---------- 4. 主流程 ----------
if __name__ == '__main__':
    # 数据路径（按实际存放位置修改）
    train_df, test_df, test_rul_true = load_cmapss(
        'train_FD001.txt', 'test_FD001.txt')
    
    # 准备训练数据
    X_train, y_train, scaler = prepare_sequences(
        train_df, is_train=True, window_size=30)
    
    # 准备测试数据（使用训练集的scaler）
    X_test, y_test, _ = prepare_sequences(
        test_df, scaler=scaler, is_train=False, window_size=30)
    
    # 注意: 测试集每台只取最后一个窗口作为预测点
    X_test_last = []
    test_ids = test_df['id'].unique()
    for unit_id in test_ids:
        unit_data = test_df[test_df['id'] == unit_id]
        n = len(unit_data)
        if n >= 30:
            scaled = scaler.transform(unit_data[sensor_cols].values)
            X_test_last.append(scaled[-30:])
        else:
            # 不足30步的序列做边缘填充
            pad_len = 30 - n
            scaled = scaler.transform(unit_data[sensor_cols].values)
            padded = np.vstack([np.tile(scaled[0], (pad_len, 1)), scaled])
            X_test_last.append(padded)
    X_test_last = np.array(X_test_last)
    
    # 训练模型
    model = build_lstm((30, 21))
    early_stop = EarlyStopping(monitor='val_loss', patience=5, 
                               restore_best_weights=True)
    
    # 划分验证集
    split = int(0.8 * len(X_train))
    X_tr, X_val = X_train[:split], X_train[split:]
    y_tr, y_val = y_train[:split], y_train[split:]
    
    history = model.fit(X_tr, y_tr, validation_data=(X_val, y_val),
                        epochs=50, batch_size=512, 
                        callbacks=[early_stop], verbose=1)
    
    # 预测测试集RUL
    y_pred = model.predict(X_test_last).flatten()
    y_pred = np.clip(y_pred, 0, 150)
    
    # ---------- 5. 评估 ----------
    rmse = np.sqrt(np.mean((y_pred - test_rul_true)**2))
    
    # PHM08评分函数
    diff = y_pred - test_rul_true
    score = np.sum(np.where(diff < 0, 
                            np.exp(-diff/13) - 1, 
                            np.exp(diff/10) - 1))
    
    print(f'RMSE: {rmse:.3f}')
    print(f'PHM Score: {score:.3f}')
    # 预期输出示例（基于公开基准，非本次运行结果）:
    # RMSE: 约18.5~20.0 (取决于随机种子与数据划分)
    # PHM Score: 约300~500

验证方法与判定标准

运行上述脚本需准备三个文件：train_FD001.txt、test_FD001.txt和RUL_FD001.txt，均来自CMAPSS数据集FD001子集。训练完成后观察终端输出的RMSE和PHM Score值。判定标准：RMSE低于20视为模型有效捕捉退化趋势；PHM Score越低表示提前预测误差越小。若RMSE显著高于25，应检查数据路径、滑窗长度或增加训练轮次。

本实现使用全部21个传感器通道，未做特征选择，因此模型性能代表LSTM对原始高维退化信号的直接处理能力。如需进一步提升精度，可结合特征重要性筛选或引入注意力机制，但已超出本案例最小闭环范围。