用Python做CMM测量数据可视化与SPC统计分析（实战）
三坐标测量机（CMM, Coordinate Measuring Machine）批量测量后产生大量数据，传统做法是导出Excel看偏差彩图，但要做统计过程控制（SPC）、趋势分析、过程能力指数Cpk计算，Excel力不从心。本文用Python完整实现CMM数据导入、清洗、可视化、SPC分析全流程，代码可直接复用。
一、CMM测量数据特点
CMM测量数据有以下特征：
- 多变量：单件零件可能有几十到几百个特征
- 批量：批量生产时单批可能几百件
- 结构化：每条记录含零件ID、特征名、名义值、实测值、偏差、上下公差
- 时序性：按测量顺序排列，可做时间序列分析
根据ISO 22514-2: 2017标准，CMM数据可用于计算过程能力指数（Cp、Cpk、Pp、Ppk）。
二、技术栈
Python 3.10+
pandas 2.0+ # 数据处理
numpy 1.24+ # 数值计算
matplotlib 3.7+ # 可视化
scipy 1.10+ # 统计分布
openpyxl 3.1+ # Excel读写
安装：
pip install pandas numpy matplotlib scipy openpyxl
三、数据导入与清洗
3.1 数据格式约定
PC-DMIS导出的Excel数据通常格式：
3.2 导入代码
import pandas as pd
import numpy as np
def load_cmm_data(file_path):
"""加载CMM测量数据"""
df = pd.read_excel(file_path)
df['Deviation'] = df['Measured'] - df['Nominal']
df['Tol_Upper'] = df['Nominal'] + df['Tol_Plus']
df['Tol_Lower'] = df['Nominal'] + df['Tol_Minus']
df['Pass'] = (df['Measured'] >= df['Tol_Lower']) & \
(df['Measured'] <= df['Tol_Upper'])
return df
df = load_cmm_data("measurement_data.xlsx")
print(f"总记录数: {len(df)}")
print(f"合格率: {df['Pass'].mean():.2%}")
3.3 异常值检测
用3σ准则剔除异常测量数据：
def detect_outliers(df, feature_name, threshold=3):
"""3σ法则检测异常值"""
subset = df[df['Feature'] == feature_name]
mean = subset['Deviation'].mean()
std = subset['Deviation'].std()
z_score = (subset['Deviation'] - mean) / std
outliers = subset[abs(z_score) > threshold]
return outliers
outliers = detect_outliers(df, 'Hole_D1')
print(f"异常值数量: {len(outliers)}")
四、过程能力指数Cpk计算
4.1 公式定义
根据ISO 22514-2: 2017：
• Cp（过程能力指数）= (USL - LSL) / (6σ)
• Cpk（过程能力指数偏移）= min((USL - μ), (μ - LSL)) / (3σ)
其中：USL=上规格限，LSL=下规格限，μ=测量均值，σ=标准差
Cpk判读标准：
Cpk ≥ 1.67：过程能力优秀
1.33 ≤ Cpk < 1.67：过程能力充分
1.00 ≤ Cpk < 1.33：过程能力勉强
Cpk < 1.00：过程能力不足
4.2 Python实现
def calc_cpk(df, feature_name):
"""计算指定特征的Cpk"""
subset = df[df['Feature'] == feature_name]
usl = subset['Tol_Upper'].iloc[0]
lsl = subset['Tol_Lower'].iloc[0]
mean = subset['Measured'].mean()
std = subset['Measured'].std(ddof=1)
cp = (usl - lsl) / (6 * std)
cpk = min((usl - mean), (mean - lsl)) / (3 * std)
return {
'feature': feature_name,
'usl': usl,
'lsl': lsl,
'mean': mean,
'std': std,
'cp': cp,
'cpk': cpk,
'sample_size': len(subset)
}
result = calc_cpk(df, 'Hole_D1')
print(f"Cpk = {result['cpk']:.3f}, Cp = {result['cp']:.3f}")
五、可视化分析
5.1 偏差直方图+正态分布拟合
import matplotlib.pyplot as plt
import scipy.stats as stats
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
def plot_histogram_with_normal(df, feature_name):
"""绘制直方图+正态分布拟合曲线"""
subset = df[df['Feature'] == feature_name]
deviations = subset['Deviation'].values
fig, ax = plt.subplots(figsize=(10, 6))
n, bins, patches = ax.hist(deviations, bins=20,
density=True, alpha=0.7,
color='steelblue',
edgecolor='black')
# 正态分布拟合
mu, sigma = stats.norm.fit(deviations)
x = np.linspace(bins[0], bins[-1], 100)
ax.plot(x, stats.norm.pdf(x, mu, sigma),
'r-', linewidth=2, label=f'正态拟合 μ={mu:.4f}, σ={sigma:.4f}')
# 公差限
tol_plus = subset['Tol_Plus'].iloc[0]
tol_minus = subset['Tol_Minus'].iloc[0]
ax.axvline(x=tol_plus, color='red', linestyle='--',
label=f'上公差 +{tol_plus}')
ax.axvline(x=tol_minus, color='blue', linestyle='--',
label=f'下公差 {tol_minus}')
ax.set_xlabel('偏差 (mm)')
ax.set_ylabel('概率密度')
ax.set_title(f'{feature_name} 偏差分布')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{feature_name}_histogram.png', dpi=150)
plt.show()
plot_histogram_with_normal(df, 'Hole_D1')
5.2 SPC控制图（Xbar-R图）
def plot_xbar_r_chart(df, feature_name, subgroup_size=5):
"""绘制Xbar-R控制图"""
subset = df[df['Feature'] == feature_name]
measurements = subset['Measured'].values
# 分组
n_groups = len(measurements) // subgroup_size
groups = measurements[:n_groups * subgroup_size].reshape(
n_groups, subgroup_size)
group_means = groups.mean(axis=1)
group_ranges = groups.max(axis=1) - groups.min(axis=1)
# 控制限系数（n=5时A2=0.577, D3=0, D4=2.114）
A2, D3, D4 = 0.577, 0, 2.114
xbar_bar = group_means.mean()
r_bar = group_ranges.mean()
ucl_x = xbar_bar + A2 * r_bar
lcl_x = xbar_bar - A2 * r_bar
ucl_r = D4 * r_bar
lcl_r = D3 * r_bar
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
# Xbar图
ax1.plot(group_means, 'bo-', markersize=6)
ax1.axhline(y=xbar_bar, color='g', linestyle='-',
label=f'CL={xbar_bar:.4f}')
ax1.axhline(y=ucl_x, color='r', linestyle='--',
label=f'UCL={ucl_x:.4f}')
ax1.axhline(y=lcl_x, color='r', linestyle='--',
label=f'LCL={lcl_x:.4f}')
ax1.set_ylabel('子组均值')
ax1.set_title(f'{feature_name} Xbar-R 控制图')
ax1.legend()
ax1.grid(True, alpha=0.3)
# R图
ax2.plot(group_ranges, 'bo-', markersize=6)
ax2.axhline(y=r_bar, color='g', linestyle='-',
label=f'CL={r_bar:.4f}')
ax2.axhline(y=ucl_r, color='r', linestyle='--',
label=f'UCL={ucl_r:.4f}')
ax2.axhline(y=lcl_r, color='r', linestyle='--',
label=f'LCL={lcl_r:.4f}')
ax2.set_xlabel('子组编号')
ax2.set_ylabel('子组极差')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{feature_name}_xbar_r.png', dpi=150)
plt.show()
plot_xbar_r_chart(df, 'Hole_D1')
5.3 趋势分析图
def plot_trend(df, feature_name):
"""绘制测量趋势图"""
subset = df[df['Feature'] == feature_name].reset_index(drop=True)
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(subset.index + 1, subset['Measured'],
'bo-', markersize=5, alpha=0.7)
ax.axhline(y=subset['Tol_Upper'].iloc[0], color='r',
linestyle='--', label='上公差')
ax.axhline(y=subset['Tol_Lower'].iloc[0], color='r',
linestyle='--', label='下公差')
ax.axhline(y=subset['Nominal'].iloc[0], color='g',
linestyle='-', alpha=0.5, label='名义值')
# 滚动平均
rolling_mean = subset['Measured'].rolling(window=10).mean()
ax.plot(subset.index + 1, rolling_mean,
'r-', linewidth=2, label='10件滚动平均')
ax.set_xlabel('测量序号')
ax.set_ylabel('实测值 (mm)')
ax.set_title(f'{feature_name} 测量趋势')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{feature_name}_trend.png', dpi=150)
plt.show()
plot_trend(df, 'Hole_D1')
六、批量分析所有特征
def batch_analysis(df, output_excel):
"""批量分析所有特征，输出汇总报告"""
features = df['Feature'].unique()
results = []
for feat in features:
result = calc_cpk(df, feat)
results.append(result)
result_df = pd.DataFrame(results)
result_df['Status'] = result_df['cpk'].apply(
lambda x: '优秀' if x >= 1.67 else
('充分' if x >= 1.33 else
('勉强' if x >= 1.00 else '不足')))
result_df.to_excel(output_excel, index=False)
return result_df
summary = batch_analysis(df, 'cpk_summary.xlsx')
print(summary[['feature', 'cpk', 'cp', 'Status']])
输出示例：
feature cpk cp Status
0 Hole_D1 1.5234 1.6124 充分
1 Hole_D2 1.0987 1.2156 勉强
2 Hole_D3 0.8765 0.9234 不足
3 Plane_A 2.1234 2.2456 优秀
七、关键统计概念速查
指标 公式 说明
均值μ Σx/n 数据集中趋势
标准差σ sqrt(Σ(x-μ)²/(n-1)) 数据离散程度
Cp (USL-LSL)/(6σ) 过程能力潜在指数
Cpk min((USL-μ),(μ-LSL))/(3σ) 过程能力实际指数
Pp (USL-LSL)/(6s) 过程性能潜在指数
Ppk min((USL-x̄),(x̄-LSL))/(3s) 过程性能实际指数
💡 Cp和Cpk的区别：Cp用样本标准差（短期），Ppk用总体标准差（长期）。
八、扩展应用
- 自动报警：当Cpk连续3天低于1.33时触发邮件报警。
- MES集成：通过REST API将Cpk数据上传到MES系统。
- 机器学习：用Isolation Forest检测测量异常模式。
- Web可视化：用Streamlit搭建交互式SPC看板。
Streamlit快速搭建示例：
import streamlit as st
st.title("CMM测量数据SPC看板")
uploaded = st.file_uploader("上传Excel", type=['xlsx'])
if uploaded:
df = load_cmm_data(uploaded)
feature = st.selectbox("选择特征", df['Feature'].unique())
if st.button("分析"):
result = calc_cpk(df, feature)
st.write(f"Cpk: {result['cpk']:.3f}")
plot_histogram_with_normal(df, feature)
st.pyplot()
九、注意事项
样本量：Cpk计算样本量≥25件才有统计意义，建议≥50件。
正态性检验：Cpk假设数据服从正态分布，需先用Shapiro-Wilk检验。
测量系统分析：在做Cpk前先做MSA，确保测量系统GR&R≤10%。
数据时效：Cpk基于历史数据，过程调整后需重新计算。
Shapiro-Wilk正态性检验代码：
from scipy.stats import shapiro
def check_normality(data, alpha=0.05):
stat, p = shapiro(data)
if p > alpha:
print(f"数据服从正态分布 (p={p:.4f})")
else:
print(f"数据不服从正态分布 (p={p:.4f})，Cpk结果需谨慎使用")
十、参考资料
- ISO 22514-2: 2017 Statistical methods — Process performance
- ISO 10360-2: 2009 CMM验收标准
- AIAG SPC参考手册第二版
- Python官方文档：pandas、matplotlib、scipy
总结
本文完整介绍了CMM测量数据从导入、清洗、Cpk计算到可视化的Python实现方案。核心要点：
-
- pandas处理结构化数据高效
-
- scipy.stats做正态分布拟合
-
- matplotlib绘制SPC控制图
-
- Streamlit快速搭建交互看板
浙公网安备 33010602011771号