可靠性工程统计：让失效变得"可预测"工程世界里有一个残酷的真相——所有东西都会坏。问题不是"会不会坏"，而是"什么时候坏

工程世界里有一个残酷的真相——所有东西都会坏。问题不是"会不会坏"，而是"什么时候坏、坏得多快、能不能提前知道"。可靠性工程统计，正是那门试图把"失效"这件事讲清楚、算明白的学问。它横跨概率论、统计推断与工程实践，是航空、汽车、半导体、电力等行业质量管理的核心工具。

一、失效率与可靠度函数：描述"死亡"的语言

要研究失效，先得有一套描述它的数学语言。设产品寿命为随机变量 $T$ ，其概率密度函数为 $f(t)$ ，累积分布函数为 $F(t) = P(T \leq t)$ 。

可靠度函数（Reliability Function）描述的是产品在时刻 $t$ 仍然存活的概率：

$R(t) = P(T > t) = 1 - F(t)$

失效率函数（Hazard Function，也叫风险函数）则更有工程直觉——它描述的是"已经活到时刻 $t$ 的产品，在接下来极短时间内失效的条件概率密度"：

$h(t) = \frac{f(t)}{R(t)}$

这两者之间有一个优雅的关系。由于 $f(t) = -R'(t)$ ，可以推出：

$R(t) = \exp\!\left(-\int_0^t h(u)\,du\right)$

工程师们还喜欢用累积失效率函数（Cumulative Hazard Function）：

$H(t) = \int_0^t h(u)\,du$

于是 $R(t) = e^{-H(t)}$ ，形式极为简洁。

浴盆曲线：失效率的"一生"

现实中大多数产品的失效率 $h(t)$ 呈现出著名的浴盆曲线（Bathtub Curve）形状：

早期失效期：失效率高且递减，对应制造缺陷、材料瑕疵
偶然失效期：失效率近似恒定，产品处于"最佳状态"
耗损失效期：失效率递增，老化、疲劳开始主导

这条曲线不只是示意图，它直接决定了产品的质保策略、维修计划和报废时机。

二、两大寿命分布：指数与威布尔

指数分布：最简单的"无记忆"模型

指数分布是可靠性分析的入门模型，其失效率为常数 $\lambda$ ：

$h(t) = \lambda$

对应的可靠度函数：

$R(t) = e^{-\lambda t}$

概率密度函数：

$f(t) = \lambda e^{-\lambda t}, \quad t \geq 0$

平均寿命（MTTF，Mean Time To Failure）为 $1/\lambda$ 。

指数分布有个迷人又危险的性质——无记忆性： $P(T > s+t \mid T > s) = P(T > t)$ 。产品"不记得"自己已经用了多久，每一刻都像刚出厂一样。这在电子元器件的偶然失效阶段是合理假设，但对于有磨损、疲劳的机械零件，就过于乐观了。

威布尔分布：工程师的瑞士军刀

威布尔（Weibull）分布是可靠性工程中用途最广的分布，没有之一。它有两个参数：形状参数 $\beta$ （shape）和尺度参数 $\eta$ （scale，也叫特征寿命）。

失效率函数：

$h(t) = \frac{\beta}{\eta}\left(\frac{t}{\eta}\right)^{\beta-1}$

可靠度函数：

$R(t) = \exp\!\left[-\left(\frac{t}{\eta}\right)^\beta\right]$

概率密度函数：

$f(t) = \frac{\beta}{\eta}\left(\frac{t}{\eta}\right)^{\beta-1}\exp\!\left[-\left(\frac{t}{\eta}\right)^\beta\right]$

威布尔分布的神奇之处在于 $\beta$ 的物理含义：

$\beta$ 值	失效率形态	对应阶段
$\beta < 1$	递减	早期失效（婴儿期）
$\beta = 1$	恒定（退化为指数分布）	偶然失效期
$\beta > 1$	递增	耗损失效期
$\beta \approx 3.5$	近似正态分布	典型磨损失效

一个分布，通过调整 $\beta$ ，就能模拟浴盆曲线的三个阶段——这正是它被称为"工程师的瑞士军刀"的原因。

平均寿命（MTTF）为：

$\text{MTTF} = \eta \cdot \Gamma\!\left(1 + \frac{1}{\beta}\right)$

其中 $\Gamma(\cdot)$ 是伽马函数。

三、加速寿命试验：用时间换信息

现代产品的设计寿命动辄十年、二十年，没有人等得起。加速寿命试验（Accelerated Life Testing，ALT）的思路是：在更高应力（温度、电压、振动等）下加速产品失效，再用物理-统计模型把结果"折算"回正常工作条件。

Arrhenius 模型：温度加速的经典

最常用的加速模型是 Arrhenius 模型，来自化学反应速率理论。它认为失效率与温度的关系为：

$\lambda(T) = A \cdot \exp\!\left(-\frac{E_a}{k_B T}\right)$

其中：

$E_a$ 是激活能（单位：eV），反映失效机理对温度的敏感程度
$k_B = 8.617 \times 10^{-5}\ \text{eV/K}$ 是玻尔兹曼常数
$T$ 是绝对温度（K）

加速因子（Acceleration Factor）定义为：

$AF = \frac{\lambda(T_{\text{stress}})}{\lambda(T_{\text{use}})} = \exp\!\left[\frac{E_a}{k_B}\left(\frac{1}{T_{\text{use}}} - \frac{1}{T_{\text{stress}}}\right)\right]$

实践中，工程师在多个应力水平下做试验，用最大似然估计（MLE）同时拟合威布尔参数和加速模型参数，再外推到使用条件。

四、系统可靠性：串联与并联

单个零件的可靠性固然重要，但真实系统是由成百上千个零件组成的网络。系统可靠性研究的是：零件的可靠性如何"组合"成系统的可靠性。

串联系统：最脆弱的那个决定一切

串联系统中，任意一个零件失效，系统就失效。若各零件独立，系统可靠度为：

$R_{\text{series}}(t) = \prod_{i=1}^{n} R_i(t)$

串联系统的可靠性由最薄弱环节决定。零件越多，系统可靠性越低——这就是为什么复杂系统的设计要尽量减少串联环节。

并联系统：冗余是工程师的护身符

并联系统中，所有零件都失效，系统才失效（主动冗余）。系统可靠度为：

$R_{\text{parallel}}(t) = 1 - \prod_{i=1}^{n}\left[1 - R_i(t)\right]$

并联引入了冗余，是提升系统可靠性最直接的手段。飞机的多发动机、航天器的备份计算机，都是这个原理的工程体现。

混合系统

现实中的系统往往是串并联的混合体。分析时通常先识别最小割集（Minimal Cut Sets）或最小路集（Minimal Path Sets），再逐步化简。

五、Python 代码：把概念画出来

理论讲完，该让代码说话了。下面的代码涵盖了上述所有核心概念的可视化。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from scipy.stats import weibull_min, expon
from scipy.special import gamma

# ── 全局绘图风格 ──────────────────────────────────────────────
plt.rcParams.update({
    'font.family': 'DejaVu Sans',
    'axes.titlesize': 13,
    'axes.labelsize': 11,
    'legend.fontsize': 9.5,
    'axes.spines.top': False,
    'axes.spines.right': False,
    'figure.dpi': 130,
})

t = np.linspace(0.001, 4, 500)

# ════════════════════════════════════════════════════════════════
# 图1：威布尔分布 — PDF / 可靠度 / 失效率
# ════════════════════════════════════════════════════════════════
fig1, axes1 = plt.subplots(1, 3, figsize=(15, 4.5))
fig1.suptitle("Weibull Distribution — PDF, Reliability & Hazard Rate", fontsize=14, fontweight='bold', y=1.02)

betas  = [0.5, 1.0, 1.5, 3.0, 5.0]
eta    = 1.0
colors = ['#e74c3c', '#e67e22', '#2ecc71', '#3498db', '#9b59b6']
labels = [f'β={b}' for b in betas]

for b, c, lbl in zip(betas, colors, labels):
    # scipy weibull_min: shape=β, scale=η, loc=0
    rv = weibull_min(c=b, scale=eta)
    pdf_vals = rv.pdf(t)
    rel_vals = rv.sf(t)          # survival function = R(t)
    haz_vals = pdf_vals / rel_vals

    axes1[0].plot(t, pdf_vals,  color=c, lw=2,   label=lbl)
    axes1[1].plot(t, rel_vals,  color=c, lw=2,   label=lbl)
    axes1[2].plot(t, np.clip(haz_vals, 0, 8), color=c, lw=2, label=lbl)

titles = ['Probability Density f(t)', 'Reliability R(t)', 'Hazard Rate h(t)']
ylabels = ['f(t)', 'R(t)', 'h(t)']
for ax, ttl, yl in zip(axes1, titles, ylabels):
    ax.set_title(ttl)
    ax.set_xlabel('Time t')
    ax.set_ylabel(yl)
    ax.legend()
    ax.set_xlim(0, 4)
    ax.set_ylim(bottom=0)

plt.tight_layout()
plt.savefig('weibull_distributions.png', bbox_inches='tight')
plt.show()

# ════════════════════════════════════════════════════════════════
# 图2：浴盆曲线（混合威布尔模拟）
# ════════════════════════════════════════════════════════════════
fig2, ax2 = plt.subplots(figsize=(9, 4.5))
fig2.suptitle("Bathtub Curve — Simulated via Mixed Weibull", fontsize=14, fontweight='bold')

t_bath = np.linspace(0.01, 10, 1000)

# 三段：早期(β<1) + 偶然(β=1) + 耗损(β>1)，加权叠加
def weibull_hazard(t, beta, eta):
    return (beta / eta) * (t / eta) ** (beta - 1)

h_infant = 0.6 * weibull_hazard(t_bath, beta=0.4, eta=0.8)
h_random = 0.15 * weibull_hazard(t_bath, beta=1.0, eta=5.0)
h_wearout = 0.05 * weibull_hazard(t_bath, beta=4.5, eta=8.0)
h_total = h_infant + h_random + h_wearout

ax2.fill_between(t_bath, h_total, alpha=0.15, color='steelblue')
ax2.plot(t_bath, h_total, color='steelblue', lw=2.5, label='Total h(t)')
ax2.plot(t_bath, h_infant,  '--', color='#e74c3c', lw=1.5, alpha=0.8, label='Infant mortality')
ax2.plot(t_bath, h_random,  '--', color='#2ecc71', lw=1.5, alpha=0.8, label='Random failures')
ax2.plot(t_bath, h_wearout, '--', color='#e67e22', lw=1.5, alpha=0.8, label='Wear-out')

ax2.annotate('Infant\nMortality', xy=(0.5, 0.55), fontsize=9, color='#e74c3c',
             ha='center', style='italic')
ax2.annotate('Useful Life\n(Constant Rate)', xy=(4.5, 0.08), fontsize=9, color='#2ecc71',
             ha='center', style='italic')
ax2.annotate('Wear-out', xy=(8.8, 0.35), fontsize=9, color='#e67e22',
             ha='center', style='italic')

ax2.set_xlabel('Time t')
ax2.set_ylabel('Hazard Rate h(t)')
ax2.set_ylim(0, 0.85)
ax2.legend(loc='upper right')
plt.tight_layout()
plt.savefig('bathtub_curve.png', bbox_inches='tight')
plt.show()

# ════════════════════════════════════════════════════════════════
# 图3：Arrhenius 加速寿命模型
# ════════════════════════════════════════════════════════════════
fig3, (ax3a, ax3b) = plt.subplots(1, 2, figsize=(12, 4.5))
fig3.suptitle("Arrhenius Accelerated Life Model", fontsize=14, fontweight='bold')

kb   = 8.617e-5   # eV/K
Ea   = 0.7        # 激活能，典型半导体值
T_use = 298       # 25°C 使用温度

# 左图：AF vs 应力温度
T_stress = np.linspace(330, 450, 200)  # 57°C ~ 177°C
AF = np.exp((Ea / kb) * (1/T_use - 1/T_stress))

ax3a.plot(T_stress - 273.15, AF, color='#c0392b', lw=2.5)
ax3a.fill_between(T_stress - 273.15, AF, alpha=0.1, color='#c0392b')
ax3a.set_xlabel('Stress Temperature (°C)')
ax3a.set_ylabel('Acceleration Factor (AF)')
ax3a.set_title(f'AF vs Stress Temp  (Ea={Ea} eV, T_use=25°C)')
ax3a.set_yscale('log')
ax3a.grid(True, alpha=0.3)

# 右图：不同 Ea 下的 AF
T_s_pt = np.array([60, 85, 105, 125, 150]) + 273.15
Ea_vals = [0.3, 0.5, 0.7, 1.0, 1.2]
colors_ea = ['#1abc9c','#3498db','#9b59b6','#e74c3c','#e67e22']

for ea, col in zip(Ea_vals, colors_ea):
    af_pts = np.exp((ea / kb) * (1/T_use - 1/T_s_pt))
    ax3b.plot(T_s_pt - 273.15, af_pts, 'o-', color=col, lw=1.8,
              markersize=6, label=f'Ea={ea} eV')

ax3b.set_xlabel('Stress Temperature (°C)')
ax3b.set_ylabel('Acceleration Factor (AF)')
ax3b.set_title('AF for Different Activation Energies')
ax3b.set_yscale('log')
ax3b.legend()
ax3b.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('arrhenius_model.png', bbox_inches='tight')
plt.show()

# ════════════════════════════════════════════════════════════════
# 图4：系统可靠性 — 串联 vs 并联
# ════════════════════════════════════════════════════════════════
fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle("System Reliability: Series vs Parallel", fontsize=14, fontweight='bold')

t_sys = np.linspace(0, 5, 500)
beta_comp, eta_comp = 2.0, 2.0
rv_comp = weibull_min(c=beta_comp, scale=eta_comp)
R_single = rv_comp.sf(t_sys)

# 串联：n 个相同零件
ax = axes4[0]
ax.set_title('Series System  (n identical components)')
for n, col in zip([1, 2, 3, 5, 10], colors):
    R_series = R_single ** n
    ax.plot(t_sys, R_series, color=col, lw=2, label=f'n={n}')
ax.axhline(0.9, color='gray', lw=1, ls=':', alpha=0.6)
ax.text(4.8, 0.91, 'R=0.9', fontsize=8, color='gray', ha='right')
ax.set_xlabel('Time t')
ax.set_ylabel('System Reliability R_sys(t)')
ax.legend()
ax.set_ylim(0, 1.02)

# 并联：n 个相同零件
ax = axes4[1]
ax.set_title('Parallel System  (n identical components)')
for n, col in zip([1, 2, 3, 5, 10], colors):
    R_parallel = 1 - (1 - R_single) ** n
    ax.plot(t_sys, R_parallel, color=col, lw=2, label=f'n={n}')
ax.axhline(0.9, color='gray', lw=1, ls=':', alpha=0.6)
ax.text(4.8, 0.91, 'R=0.9', fontsize=8, color='gray', ha='right')
ax.set_xlabel('Time t')
ax.set_ylabel('System Reliability R_sys(t)')
ax.legend()
ax.set_ylim(0, 1.02)

plt.tight_layout()
plt.savefig('system_reliability.png', bbox_inches='tight')
plt.show()

# ════════════════════════════════════════════════════════════════
# 图5：威布尔概率纸（Weibull Probability Plot）
# ════════════════════════════════════════════════════════════════
fig5, ax5 = plt.subplots(figsize=(8, 5.5))
fig5.suptitle("Weibull Probability Plot (Parameter Estimation)", fontsize=14, fontweight='bold')

np.random.seed(42)
# 模拟两组失效数据：β=1.8, η=3.0 和 β=3.5, η=2.0
datasets = [
    (1.8, 3.0, '#3498db', 'Dataset A  (β=1.8, η=3.0)'),
    (3.5, 2.0, '#e74c3c', 'Dataset B  (β=3.5, η=2.0)'),
]

for beta_true, eta_true, col, lbl in datasets:
    n_samples = 30
    data = eta_true * np.random.weibull(beta_true, n_samples)
    data_sorted = np.sort(data)

    # 中位秩（Median Rank）近似：(i - 0.3) / (n + 0.4)
    i_rank = np.arange(1, n_samples + 1)
    F_i = (i_rank - 0.3) / (n_samples + 0.4)

    # 线性化：ln(ln(1/(1-F))) vs ln(t)
    y_vals = np.log(np.log(1 / (1 - F_i)))
    x_vals = np.log(data_sorted)

    # 线性回归拟合
    coeffs = np.polyfit(x_vals, y_vals, 1)
    beta_est  = coeffs[0]
    eta_est   = np.exp(-coeffs[1] / coeffs[0])
    x_fit = np.linspace(x_vals.min() - 0.3, x_vals.max() + 0.3, 100)
    y_fit = np.polyval(coeffs, x_fit)

    ax5.scatter(x_vals, y_vals, color=col, s=30, alpha=0.7, zorder=5)
    ax5.plot(x_fit, y_fit, color=col, lw=2,
             label=f'{lbl}\n  → β̂={beta_est:.2f}, η̂={eta_est:.2f}')

ax5.set_xlabel('ln(t)')
ax5.set_ylabel('ln(ln(1 / (1−F(t))))')
ax5.legend(loc='upper left', fontsize=9)
ax5.grid(True, alpha=0.3)

# 添加概率刻度参考线
for p, lbl_p in [(0.1,'10%'), (0.5,'50%'), (0.9,'90%'), (0.99,'99%')]:
    y_ref = np.log(np.log(1/(1-p)))
    ax5.axhline(y_ref, color='gray', lw=0.8, ls='--', alpha=0.5)
    ax5.text(ax5.get_xlim()[0] if ax5.get_xlim()[0] > -5 else -0.5,
             y_ref + 0.05, lbl_p, fontsize=7.5, color='gray')

plt.tight_layout()
plt.savefig('weibull_probability_plot.png', bbox_inches='tight')
plt.show()

print("\n✅ 所有图表已生成：")
print("   1. weibull_distributions.png")
print("   2. bathtub_curve.png")
print("   3. arrhenius_model.png")
print("   4. system_reliability.png")
print("   5. weibull_probability_plot.png")

六、关键概念速查表

把上面所有核心公式整理成一张表，方便随时查阅：

概念	核心公式	关键参数含义
可靠度函数	$R(t) = 1 - F(t)$	存活概率
失效率函数	$h(t) = f(t)/R(t)$	条件失效强度
指数分布	$R(t) = e^{-\lambda t}$	$\lambda$ ：常数失效率
威布尔分布	$R(t) = \exp[-(t/\eta)^\beta]$	$\beta$ ：形状， $\eta$ ：特征寿命
Arrhenius 加速因子	$AF = \exp[(E_a/k_B)(1/T_u - 1/T_s)]$	$E_a$ ：激活能
串联系统	$R_{sys} = \prod R_i$	短板效应
并联系统	$R_{sys} = 1 - \prod(1-R_i)$	冗余增益

结语

可靠性工程统计的魅力，在于它把"失效"这件看似随机、令人沮丧的事，变成了可以建模、可以预测、可以主动干预的工程问题。威布尔分布的 $\beta$ 参数，像一个诊断仪，告诉你产品正处于生命周期的哪个阶段；加速寿命试验，让工程师在几周内看到产品二十年后的命运；串并联模型，则把系统级的可靠性设计变成了一道可以优化的数学题。

真正有趣的地方在于：这套理论的背后，是对"时间"本身的一种统计哲学——我们无法阻止熵增，但我们可以精确地描述它，并在它到来之前做好准备。