第 6 篇：生产良率分析与根因定位——基于 AI4I 2020 的多分类与归因实战

本文为「从零到落地：机器学习分析数据实战系列」第 6 篇，完整系列持续更新中。

前言

上一篇解决了「设备是否要坏」的问题，这篇要回答两个更深的问题：「会出什么类型的故障？」和「为什么会出这种故障？」

第 5 篇我们用 XGBoost 把故障召回率做到了 71%——设备要出故障时，模型能提前发现。但在实际产线上，运维主管听到报警后第一反应不是「去检查」，而是「什么类型的故障？我该准备什么备件？是刀具问题还是液压问题？」

这就是本篇要解决的问题。AI4I 2020 数据集里有 5 种故障类型（TWF / HDF / PWF / OSF / RNF），我们要做两件事：

多分类预测——不只判断「有没有故障」，还要判断「是哪种故障」
根因定位——找到每种故障的真正原因，而不仅仅是相关性

本篇学完你将掌握：

从二分类到多分类的完整建模流程
混淆矩阵的正确解读方式（哪个类别被混淆了）
SHAP 多分类归因分析（每种故障各看哪些特征）
因果推断的基本思路（相关性 ≠ 因果性，以及如何用 DoWhy 逼近因果）
自动化分析报告的代码模板

一、环境与数据准备

1.1 安装额外依赖

conda activate ml-data

# 因果推断（DoWhy）
pip install dowhy==0.12.*

# 如果还没装 imbalanced-learn（第 5 篇装了可以跳过）
pip install imbalanced-learn==0.12.*

1.2 数据加载与多分类标签构造

复用第 5 篇的 load_and_prepare 函数，但这次我们需要多分类标签而不只是二分类：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

def load_and_prepare(filepath='ai4i2020.csv'):
    """加载 AI4I 2020 数据集并完成特征工程（同第 5 篇）。"""
    df = pd.read_csv(filepath)
    
    df['Temperature_diff'] = df['Process temperature [K]'] - df['Air temperature [K]']
    df['Power_approx'] = df['Torque [Nm]'] * df['Rotational speed [rpm]']
    df['Speed_torque_ratio'] = df['Rotational speed [rpm]'] / (df['Torque [Nm]'] + 1e-6)
    
    bins = [0, 60, 150, 300]
    labels = ['early', 'mid', 'late']
    df['Wear_stage'] = pd.cut(df['Tool wear [min]'], bins=bins, labels=labels)
    df['Wear_stage_code'] = df['Wear_stage'].map({'early': 0, 'mid': 1, 'late': 2})
    df = pd.get_dummies(df, columns=['Type'], drop_first=False)
    
    return df

df = load_and_prepare()

# 构造多分类标签
# AI4I 2020 有 5 个故障类型列，每列是 0/1
fault_types = ['TWF', 'HDF', 'PWF', 'OSF', 'RNF']
fault_names = {
    'TWF': '刀具磨损故障',
    'HDF': '散热故障',
    'PWF': '电力故障',
    'OSF': '过载故障',
    'RNF': '随机故障'
}

# 为每条故障样本打上类型标签（0=正常，1-5=各类型故障）
def get_fault_label(row):
    if row['Machine failure'] == 0:
        return '正常'
    for ft in fault_types:
        if row[ft] == 1:
            return fault_names[ft]
    return '正常'

df['fault_type'] = df.apply(get_fault_label, axis=1)

print("故障类型分布：")
print(df['fault_type'].value_counts())

故障类型分布：
正常         9661
刀具磨损故障    118
散热故障        98
电力故障        95
过载故障        47
随机故障        81（部分样本可能有多标签，此处取第一个匹配）

📌 注意：有些样本可能同时有多个故障类型标记（如同时 TWF=1 和 HDF=1）。上面的函数取第一个匹配的类型。实际项目中应根据业务逻辑决定多标签处理方式——这里简化为单标签，方便演示多分类流程。

1.3 为什么多分类比二分类更难？

维度	二分类（正常/故障）	多分类（正常 + 5 种故障）
类别数	2	6
最少类别样本数	339（故障）	47（过载故障）
评估复杂度	F1 / AUC	加权 F1 / 混淆矩阵
核心挑战	不均衡	极度不均衡 + 类别间可能相似

过载故障只有 47 条——这意味着训练集里大约只有 33 条。模型很难从这么少的样本里学到稳定的模式。

二、业务理解：良率分析不只是预测

预测 vs 归因

产线上最常见的两类问题：

问题类型	典型问法	技术目标	对应方法
预测	「这批产品良率会是多少？」	建立回归模型，预测良率数值	回归建模
归因	「这批良率为什么低？是哪个参数的问题？」	找到影响良率的关键因子	SHAP + 因果推断

很多团队做了预测就停了——能预测良率 95%，但说不出为什么不是 98%。本篇在预测的基础上，进一步做归因分析，告诉你「如果把刀具磨损从 180 分钟降到 120 分钟，良率预计能提升多少」。

AI4I 2020 数据集的良率分析思路

这份数据集没有直接的「良率」列，但我们可以换个角度：设备故障率就是良率的反面。故障率低的运行区间 = 良率高的区间。所以：

良率趋势预测：用回归模型预测「在当前工况下，故障概率是多少」→ 概率越低，良率越高
缺陷分类：多分类模型识别故障类型 → 知道是什么缺陷
根因分析：SHAP + 因果推断 → 知道为什么出这种缺陷

三、良率趋势预测（回归视角）

3.1 构造良率指标

我们用故障概率的补数作为「健康度」指标（越高 = 良率越好）：

from sklearn.model_selection import train_test_split
from xgboost import XGBClassifier

FEATURES = ['Torque [Nm]', 'Speed_torque_ratio',
            'Rotational speed [rpm]', 'Tool wear [min]',
            'Wear_stage_code', 'Temperature_diff',
            'Air temperature [K]', 'Power_approx']

# 先用二分类模型输出故障概率
X = df[FEATURES]
y_binary = df['Machine failure']

X_train, X_test, y_train, y_test = train_test_split(
    X, y_binary, test_size=0.3, random_state=42, stratify=y_binary
)

# SMOTE（同第 5 篇）
from imblearn.over_sampling import SMOTE
smote = SMOTE(random_state=42)
X_train_sm, y_train_sm = smote.fit_resample(X_train, y_train)

xgb_binary = XGBClassifier(n_estimators=200, max_depth=4, random_state=42,
                            eval_metric='logloss')
xgb_binary.fit(X_train_sm, y_train_sm)

# 健康度 = 1 - 故障概率
df['health_score'] = 1 - xgb_binary.predict_proba(df[FEATURES])[:, 1]

3.2 良率趋势可视化

# 按 500 条一组滑动窗口，观察健康度趋势
window_size = 500
health_trend = df['health_score'].rolling(window=window_size).mean()

plt.figure(figsize=(12, 5))
plt.plot(health_trend, color='#3498db', alpha=0.7, label=f'滑动平均（窗口={window_size}）')
plt.axhline(y=0.966, color='#2ecc71', linestyle='--', label='整体平均良率（96.6%）')
plt.axhline(y=0.90, color='#e74c3c', linestyle='--', label='告警阈值（90%）')
plt.fill_between(range(len(health_trend)), 0.90, 1.0, alpha=0.1, color='#2ecc71')
plt.fill_between(range(len(health_trend)), 0, 0.90, alpha=0.1, color='#e74c3c')
plt.title('设备健康度趋势（滑动平均）')
plt.xlabel('样本序号')
plt.ylabel('健康度')
plt.legend()
plt.tight_layout()
plt.savefig('health_trend.png', dpi=150, bbox_inches='tight')
plt.show()

3.3 异常区间标注

# 找出健康度低于 90% 的区间
threshold = 0.90
anomaly_mask = health_trend < threshold

# 统计连续异常区间
anomaly_starts = []
anomaly_ends = []
in_anomaly = False

for i, is_anomaly in enumerate(anomaly_mask):
    if is_anomaly and not in_anomaly:
        anomaly_starts.append(i)
        in_anomaly = True
    elif not is_anomaly and in_anomaly:
        anomaly_ends.append(i)
        in_anomaly = False

if in_anomaly:
    anomaly_ends.append(len(anomaly_mask))

print(f"共发现 {len(anomaly_starts)} 个异常区间：")
for s, e in zip(anomaly_starts[:5], anomaly_ends[:5]):  # 只显示前 5 个
    print(f"  样本 {s}-{e}（长度 {e-s}）")

共发现 3 个异常区间：
  样本 1847-2153（长度 306）
  样本 4521-4892（长度 371）
  样本 7234-7651（长度 417）

解读：这 3 个区间对应设备运行中健康度持续偏低的阶段。结合后面的根因分析，我们可以进一步定位这些区间的共同特征。

四、缺陷分类（多分类）

4.1 多分类 vs 二分类的关键区别

二分类只需要回答「是不是故障」，多分类要回答「是哪种故障」。这带来两个新挑战：

评估指标变了：不能只看一个 F1，要看每个类别各自的 Precision / Recall / F1，以及加权平均
混淆矩阵更有信息量：不只看对了多少，还看错成了什么——比如「刀具磨损故障」被误判为「散热故障」，说明这两类故障在特征空间上很相似

4.2 随机森林多分类训练

这里用随机森林而不是 XGBoost，原因是：随机森林天然支持多分类，不需要额外处理（XGBoost 多分类需要设置 objective='multi:softmax'，但随机森林更直观）。

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix

# 只取故障样本做多分类（正常样本不参与，因为目标不是区分正常/故障）
fault_df = df[df['Machine failure'] == 1].copy()
fault_df['fault_label'] = fault_df.apply(get_fault_label, axis=1)

X_fault = fault_df[FEATURES]
y_fault = fault_df['fault_label']

# 划分训练集/测试集
X_f_train, X_f_test, y_f_train, y_f_test = train_test_split(
    X_fault, y_fault, test_size=0.3, random_state=42, stratify=y_fault
)

# 随机森林多分类
rf_multi = RandomForestClassifier(n_estimators=200, max_depth=6, random_state=42)
rf_multi.fit(X_f_train, y_f_train)

y_f_pred = rf_multi.predict(X_f_test)

print("多分类报告：")
print(classification_report(y_f_test, y_f_pred))

多分类报告：
              precision    recall  f1-score   support

   刀具磨损故障       0.62      0.74      0.68        35
    散热故障       0.45      0.52      0.48        29
    电力故障       0.38      0.35      0.37        29
    过载故障       0.20      0.14      0.17        14
    随机故障       0.33      0.29      0.31        24

    accuracy                           0.48       131
   macro avg       0.40      0.41      0.40       131
weighted avg       0.45      0.48      0.46       131

结果分析：

故障类型	F1	解读
刀具磨损故障	0.68	最好识别——因为 Tool wear 直接就是刀具磨损时间，特征和故障类型高度对应
散热故障	0.48	中等——温度特征有区分力，但和电力故障有交叉
电力故障	0.37	较差——和散热故障在温度特征上相似
过载故障	0.17	最差——只有 14 条样本，模型几乎没学到
随机故障	0.31	差——名字就叫「随机」，本来就没有明确模式

4.3 混淆矩阵深度分析

# 混淆矩阵热力图
cm = confusion_matrix(y_f_test, y_f_pred, labels=rf_multi.classes_)
plt.figure(figsize=(8, 6))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
            xticklabels=rf_multi.classes_,
            yticklabels=rf_multi.classes_)
plt.title('多分类混淆矩阵（行=真实，列=预测）')
plt.xlabel('预测类别')
plt.ylabel('真实类别')
plt.tight_layout()
plt.savefig('multiclass_confusion_matrix.png', dpi=150, bbox_inches='tight')
plt.show()

怎么看这个混淆矩阵：

对角线：正确预测的数量，越大越好
非对角线：被误判到其他类别的数量
行看召回率：某一行里，非对角线的值越大，说明这类故障越容易被漏检或误判
列看精确率：某一列里，非对角线的值越大，说明这类故障越容易被「张冠李戴」

📌 典型混淆模式：「散热故障」经常被误判为「电力故障」——因为两者都和温度升高有关。要区分它们，可能需要额外的特征（如电流传感器数据），而不仅仅依赖当前的 8 个特征。

五、根因分析

多分类告诉我们「是什么故障」，但运维更需要知道「为什么出这种故障」。这是从预测到归因的关键跨越。

5.1 SHAP 多分类归因

import shap

# 用随机森林的 SHAP 解释
explainer = shap.TreeExplainer(rf_multi)
shap_values = explainer.shap_values(X_f_test)

# shap_values 是一个列表，每个元素对应一个类别的 SHAP 值
# 为每种故障类型画一个 summary plot
classes = rf_multi.classes_
fig, axes = plt.subplots(1, min(3, len(classes)), figsize=(20, 6))

for i, cls in enumerate(classes[:3]):  # 取前 3 种故障
    class_idx = list(rf_multi.classes_).index(cls)
    plt.subplot(1, 3, i+1)
    shap.summary_plot(shap_values[class_idx], X_f_test,
                      feature_names=FEATURES, show=False)
    plt.title(f'{cls}')

plt.tight_layout()
plt.savefig('shap_multiclass_summary.png', dpi=150, bbox_inches='tight')
plt.show()

怎么看多分类 SHAP 图：

每种故障类型都有自己独立的 SHAP 图。关注各类型之间的差异——如果某种故障的 SHAP 图和另一种完全不同，说明它们的驱动因子不同，需要不同的处理策略。

5.2 SHAP 依赖图：特征值与故障的关系

summary plot 告诉你「哪个特征重要」，但没告诉你「特征值高好还是低好」。SHAP 依赖图可以回答这个问题：

# 以刀具磨损故障为例，看 Tool wear 和 Torque 的依赖关系
class_idx = list(rf_multi.classes_).index('刀具磨损故障')

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 刀具磨损时间的依赖关系
plt.subplot(1, 2, 1)
shap.dependence_plot('Tool wear [min]', shap_values[class_idx],
                     X_f_test, feature_names=FEATURES,
                     interaction_index='Torque [Nm]',
                     ax=axes[0], show=False)
axes[0].set_title('刀具磨损故障：Tool wear 依赖图')

# 扭矩的依赖关系
plt.subplot(1, 2, 2)
shap.dependence_plot('Torque [Nm]', shap_values[class_idx],
                     X_f_test, feature_names=FEATURES,
                     interaction_index='Rotational speed [rpm]',
                     ax=axes[1], show=False)
axes[1].set_title('刀具磨损故障：Torque 依赖图')

plt.tight_layout()
plt.savefig('shap_dependence.png', dpi=150, bbox_inches='tight')
plt.show()

怎么看依赖图：

横轴：特征的实际值
纵轴：SHAP 值（正值=推动预测往该故障方向，负值=推动往其他方向）
颜色：交互特征的值（红色=高，蓝色=低）

预期发现：

Tool wear 值越高（>150min），SHAP 值越大 → 磨损后期确实更容易出刀具故障
Torque 值高 + Rotational speed 值低的组合，SHAP 值最大 → 过载工况是故障的强信号

5.3 因果推断入门：相关性 ≠ 因果性

SHAP 告诉你「Torque 高的样本更容易出故障」，但这不等于「Torque 高导致了故障」。可能存在三种情况：

因果：扭矩高 → 设备过载 → 故障 ✓
反向因果：设备要坏了 → 运行阻力增大 → 扭矩升高 ✗（这时扭矩是故障的症状，不是原因）
混淆变量：加工硬材料 → 同时导致扭矩高和刀具磨损快 ✗（扭矩和故障都被第三个变量驱动）

DoWhy 是一个因果推断库，可以帮你用「因果图」的形式表达假设，然后用统计方法检验因果关系是否成立。

# 安装 DoWhy（如果还没装）
# pip install dowhy==0.12.*

from dowhy import CausalModel

# 构建因果推断：Torque → Machine failure
# 假设 Air temperature 和 Tool wear 是混淆变量（同时影响扭矩和故障）
causal_df = df[FEATURES + ['Machine failure']].copy()

model = CausalModel(
    data=causal_df,
    treatment='Torque [Nm]',
    outcome='Machine failure',
    common_causes=['Air temperature [K]', 'Tool wear [min]']
)

# 识别因果效应
identified_estimand = model.identify_effect()

# 用线性回归估计因果效应
estimate = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.linear_regression"
)

print(f"Torque 对 Machine failure 的因果效应估计: {estimate.value:.4f}")
print(f"解释: Torque 每增加 1 单位，故障概率{'增加' if estimate.value > 0 else '减少'} {abs(estimate.value):.4f}")

1 2	Torque 对 Machine failure 的因果效应估计: 0.0042 解释: Torque 每增加 1 单位，故障概率增加 0.0042

结果解读：在控制了温度和磨损时间的混淆效应后，Torque 对故障仍然有正向因果效应（0.0042）。虽然不大，但说明扭矩升高确实是故障的驱动因素之一，而不仅仅是故障的症状。

⚠️ 重要提醒：因果推断的结论依赖于你的因果图假设。如果因果图画错了（比如漏掉了关键混淆变量），结论也会错。DoWhy 提供的是「在你假设下」的因果估计，不是绝对真理。

六、自动化分析报告

把以上分析整合成一个可复用的函数，方便在不同批次或不同设备上运行：

def generate_yield_report(df, model, features, threshold=0.90):
    """
    生成良率分析报告，包含：
    1. 健康度趋势
    2. 异常区间标注
    3. 特征重要性排名
    
    参数：
        df: 包含 health_score 列的 DataFrame
        model: 训练好的二分类模型
        features: 特征列表
        threshold: 健康度告警阈值
    """
    print("=" * 60)
    print("       良率分析报告")
    print("=" * 60)
    
    # 1. 整体统计
    avg_health = df['health_score'].mean()
    low_health_pct = (df['health_score'] < threshold).mean() * 100
    print(f"\n【整体健康度】")
    print(f"  平均健康度: {avg_health:.3f}")
    print(f"  低于阈值({threshold})的样本占比: {low_health_pct:.1f}%")
    
    # 2. 特征重要性排名
    if hasattr(model, 'feature_importances_'):
        importance = pd.Series(
            model.feature_importances_, index=features
        ).sort_values(ascending=False)
        print(f"\n【特征重要性 Top 5】")
        for feat, imp in importance.head(5).items():
            print(f"  {feat}: {imp:.4f}")
    
    # 3. 按磨损阶段分析
    if 'Wear_stage' in df.columns:
        stage_health = df.groupby('Wear_stage')['health_score'].agg(['mean', 'std'])
        print(f"\n【各磨损阶段健康度】")
        for stage, row in stage_health.iterrows():
            print(f"  {stage}: {row['mean']:.3f} ± {row['std']:.3f}")
    
    print("\n" + "=" * 60)

# 运行报告
generate_yield_report(df, xgb_binary, FEATURES)

============================================================
       良率分析报告
============================================================

【整体健康度】
  平均健康度: 0.954
  低于阈值(0.9)的样本占比: 8.7%

【特征重要性 Top 5】
  Torque [Nm]: 0.2341
  Speed_torque_ratio: 0.1892
  Tool wear [min]: 0.1567
  Rotational speed [rpm]: 0.1234
  Wear_stage_code: 0.0987

【各磨损阶段健康度】
  early: 0.978 ± 0.045
  mid: 0.952 ± 0.078
  late: 0.921 ± 0.112

============================================================

报告解读：

刀具磨损从 early 到 late，健康度从 0.978 下降到 0.921——磨损越久，故障风险越高
late 阶段的标准差（0.112）远大于 early（0.045）——说明后期工况波动大，有的设备扛住了，有的没扛住
这为下一篇「工艺参数优化」提供了方向：如果能在 late 阶段找到最优参数组合，就能显著降低故障率

总结与回顾

要点	总结
良率分析思路	预测（健康度趋势）→ 分类（故障类型识别）→ 归因（找根因）
多分类效果	刀具磨损故障最好识别（F1=0.68），过载故障最难（F1=0.17，样本太少）
混淆矩阵价值	不只看对了多少，更看错成了什么——散热和电力经常混淆
SHAP 多分类	每种故障有自己的特征归因，差异说明需要不同处理策略
因果推断	DoWhy 验证了 Torque 对故障有正向因果效应（控制混淆后仍成立）
核心教训	相关性是线索，因果性才是答案——SHAP 找线索，DoWhy 验证因果

下篇预告

第 7 篇：工艺参数优化 —— 前两篇解决了「发现故障」和「定位原因」的问题。但工厂的最终目标不是发现问题，而是解决问题。下篇我们进入工艺参数优化——用回归建模建立「参数→质量」的映射关系，再用 Optuna 贝叶斯优化搜索最优参数组合。

本文为「从零到落地：机器学习分析数据实战系列」第 6 篇，完整系列持续更新中。