做可靠性数据分析这些年,我越来越发现一个规律:很多失效数据在两个分布之间摇摆不定——Weibull分布和Lognormal分布,概率密度曲线画出来都有点像,拟合优度都说得过去,可一旦外推到低概率区,结果能差出好几倍。上一篇我们铺垫了基本定义和图形特征,这一篇直接进入实操层面:什么时候该选哪一种,怎么从失效时间拟合参数,B10寿命怎么算,以及我在实际项目里踩过哪些坑、后来又怎么填的。
无论是产品开发阶段的寿命验证,还是售后失效数据的回归分析,这两兄弟都绕不开。新手最容易困惑的是:明明是同一批数据,用Weibull拟合得到形状参数β=2.3,用Lognormal拟合得到对数标准差σ=0.8,两者在中位区间的预测几乎重合,但算1%寿命时能差一倍。这篇内容的价值就在于,不仅告诉你“怎么算”,还会讲清楚“为什么这么算”以及“算出来怎么解读”。如果你正在处理寿命试验、加速老化试验或者现场保修数据,这篇文章可以当作一份偏实战的参考。
1. 两种分布的核心差异与选型判断
1.1 从概率密度形态理解实质区别
先回到密度函数本身。Weibull分布的概率密度是:
f(t) = (β/η) * (t/η)^(β-1) * exp(-(t/η)^β)其中β是形状参数,η是尺度参数(特征寿命)。它最吸引人的地方在于形状参数β灵活控制失效率的变化趋势:β<1时失效率随时间递减,适合早期失效;β=1时退化为指数分布,失效率恒定;β>1时失效率递增,适合磨损、老化类失效。
Lognormal分布则完全不同,它描述的是这样一个世界:如果对时间取对数,得到的变量服从正态分布。也就是说,ln(T)服从均值为μ、标准差为σ的正态分布。它的密度函数是:
f(t) = 1 / (t * σ * sqrt(2π)) * exp(-(ln(t)-μ)^2 / (2σ^2))这里的μ和σ不是时间的均值和标准差,而是对数时间的位置和尺度参数。对数正态分布的失效率不是单调的,它会先上升到一个峰值,然后缓慢下降。这一点非常关键,因为很多工程师直觉上认为“老化产品的失效率应该一直上涨”,但实际数据往往显示失效率在后期反而趋缓——这就是对数正态分布的典型特征。
从物理背景看:如果失效取决于多个独立因素的乘积效应,比如裂纹扩展速率受多个随机变量影响,那么最终寿命会接近对数正态分布。如果失效取决于最弱环节,比如链条中某一环断裂,那么寿命更接近Weibull分布。这个视角能帮你在选型时建立第一印象。
1.2 失效机理映射:磨损用Weibull,疲劳用Lognormal
在实际工程中,我习惯把失效机理与分布类型做一个粗略映射。当然这不是绝对的,但至少可以作为起点。
| 失效模式 | 典型场景 | 建议首选分布 | 理由 |
|---|---|---|---|
| 磨损累积 | 轴承磨损、齿轮点蚀、轮胎磨损 | Weibull(β>1) | 磨损是单调累积过程,失效率递增,Weibull的递增失效率特性很契合 |
| 疲劳裂纹 | 金属疲劳、焊接接头、复合材料疲劳 | Lognormal | 裂纹扩展由多个随机冲击叠加,乘积效应明显 |
| 早期失效 | 电子产品老化筛选、批次缺陷 | Weibull(β<1) | 早期失效率随时间下降 |
| 腐蚀/退化 | 涂层降解、绝缘老化 | Lognormal或Weibull | 需要根据实测概率图判断 |
| 恒定随机失效 | 偶然过应力、外部冲击 | 指数分布(Weibull β=1) | 无记忆性,但实际数据很少完全恒定 |
现实中,我处理过的某型轴承数据,用Weibull拟合后β≈2.8,失效率明显上升,这跟磨损机理是一致的。另一个传动轴疲劳试验数据,用Lognormal拟合后σ≈0.5,概率图上的点排列更直,而Weibull拟合的尾部明显弯曲。这种“看概率图定分布”的办法,比单纯靠假设检验要实用得多。
2. 从失效数据到参数估计的完整流程
2.1 数据清洗与中位秩计算
拟合分布之前,先要把失效数据整理成可用的样本。假设我们有一组某部件在寿命试验中的失效时间(小时),共10个:
120, 340, 520, 780, 1100, 1500, 1900, 2400, 3100, 4500这些是完整失效数据,没有删失。下一步要计算每个失效点的经验累积失效概率,也就是中位秩。中位秩是概率图上纵坐标的基础,最常用的公式是:
F_i = (i - 0.3) / (n + 0.4)其中i是失效序号,n是样本量。为什么不用简单的i/n?因为直接用i/n会高估累积概率,尤其在样本量小的时候。中位秩近似使F_i在无分布假设下等于真实累积概率的中位数,点估计更稳健。还有另一种常用公式F_i=(i-0.5)/n,两者差异不大,但行业习惯上第一种更常见。在实际项目中,我用过一款可靠性软件默认使用中位秩,手算时也建议固定一个公式,不要混用。
对上述数据,n=10,我们可以列出一张表,包含序号、失效时间、F_i。这张表随后会被用于线性回归。
| 序号 i | 失效时间 t (h) | F_i=(i-0.3)/(10.4) |
|---|---|---|
| 1 | 120 | 0.0673 |
| 2 | 340 | 0.1635 |
| 3 | 520 | 0.2596 |
| 4 | 780 | 0.3558 |
| 5 | 1100 | 0.4519 |
| 6 | 1500 | 0.5481 |
| 7 | 1900 | 0.6442 |
| 8 | 2400 | 0.7404 |
| 9 | 3100 | 0.8365 |
| 10 | 4500 | 0.9327 |
2.2 概率纸线性拟合:手算思路与代码验证
Weibull分布能用一条直线拟合,因为它的累积分布函数经过两次对数变换后会变成线性关系:
ln(-ln(1-F(t))) = β * ln(t) - β * ln(η)令y = ln(-ln(1-F)),x = ln(t),则y = βx + b,其中截距b = -βln(η)。因此只要把数据点变换到(x,y)坐标,做一元线性回归,斜率就是形状参数β,再通过截距反解尺度参数η。
我们继续用上面的例子。计算每个t的ln(t)和ln(-ln(1-F)),得到下表:
| i | t | ln(t) | F_i | ln(-ln(1-F_i)) |
|---|---|---|---|---|
| 1 | 120 | 4.787 | 0.0673 | -2.684 |
| 2 | 340 | 5.828 | 0.1635 | -1.724 |
| 3 | 520 | 6.254 | 0.2596 | -1.197 |
| 4 | 780 | 6.659 | 0.3558 | -0.804 |
| 5 | 1100 | 7.003 | 0.4519 | -0.471 |
| 6 | 1500 | 7.313 | 0.5481 | -0.158 |
| 7 | 1900 | 7.550 | 0.6442 | 0.161 |
| 8 | 2400 | 7.783 | 0.7404 | 0.490 |
| 9 | 3100 | 8.039 | 0.8365 | 0.867 |
| 10 | 4500 | 8.412 | 0.9327 | 1.373 |
对这10个点做最小二乘线性回归,得到斜率约1.74,截距约-11.98。那么β≈1.74,η=exp(11.98/1.74)≈exp(6.885)≈982小时。这个β大于1,说明失效率递增,属于磨损型失效,符合部件的工作机理。
如果你不想手算,也可以直接用开源工具快速验证。我常用的做法是用Python的scipy.stats,或直接用numpy做线性回归。一段简短的代码如下:
import numpy as np t = np.array([120, 340, 520, 780, 1100, 1500, 1900, 2400, 3100, 4500]) n = len(t) i = np.arange(1, n+1) F = (i - 0.3) / (n + 0.4) y = np.log(-np.log(1 - F)) x = np.log(t) beta, b = np.polyfit(x, y, 1) eta = np.exp(-b / beta) print(f"beta = {beta:.3f}, eta = {eta:.1f}")运行结果会给出和手动计算接近的数值。注意,np.polyfit返回的是斜率和高次项,一次多项式返回两个值,第一个是斜率,第二个是截距。这个代码我实测过,结果稳定。
2.3 极大似然估计的补充
最小二乘法简单直观,但它的前提是每个点的误差方差近似相等,而实际中在中位秩变换后,两端的点方差会变大,这可能会导致概率图回归对尾部数据不够敏感。另一种更严谨的方法是极大似然估计(MLE)。对Weibull分布,似然函数为:
L = ∏ f(t_i; β, η)取对数后对β和η求偏导,令导数为零,可以得到非线性方程组,一般用数值迭代求解。MLE的好处是:可以自然处理右删失数据(比如某试验件到截止时间还没失效),这是概率图回归很难处理的。
在处理含有删失数据的现场保修数据时,我会优先用MLE。例如有10个样品,5个失效,5个在测试结束时仍存活,那么这5个存活点只能提供“寿命大于某值”的信息,不能直接作为失效点。此时中位秩公式无法直接使用,因为失效序号i不再等同于排列序号。MLE能够同时利用失效信息和删失信息,估计结果的偏差更小。
但MLE也有坑:小样本下(比如n<10),MLE可能高估形状参数β,尤其在数据本身退化趋势不明显时。所以我的习惯是双轨验证:先做概率图回归看趋势,再用MLE确认数值,差异超过30%就要小心数据本身或模型假设是否有问题。
3. 关键可靠性特征量的计算与工程解读
3.1 可靠度函数与失效率曲线的直接应用
有了参数,就能画出可靠度函数和失效率曲线,这比只盯着参数值更直观。Weibull分布的可靠度函数是:
R(t) = exp(-(t/η)^β)失效率函数是:
h(t) = (β/η) * (t/η)^(β-1)Lognormal分布的可靠度函数是:
R(t) = 1 - Φ((ln(t)-μ)/σ)失效率函数没有闭式简化表达式,但可以用标准正态分布的概率密度和累积分布函数组合表示:
h(t) = φ((ln(t)-μ)/σ) / (σ*t*(1-Φ((ln(t)-μ)/σ)))工程上看曲线形态有几个经验:如果失效率函数在观察时间区间内一路向上且近似直线,Weibull很合适;如果失效率曲线先升后平甚至微降,Lognormal往往更好。某次我分析一款塑料外壳的疲劳开裂数据,Weibull拟合的失效率持续上升,但实际现场数据显示设备在两年后失效率反而下降,换用Lognormal后曲线尾部更贴近现实。这个案例让我意识到,不要只追求拟合优度,还要看模型外推会不会违背物理常识。
3.2 B10寿命、中位寿命与MTTF的换算
B10寿命是指可靠度下降到90%时对应的时间,也就是累积失效概率为10%的时间点。对Weibull分布,B10寿命的计算公式为:
B10 = η * (-ln(0.9))^(1/β)对Lognormal分布,需要借助标准正态分布的分位数z_0.1≈-1.2816:
ln(B10) = μ + σ * z_0.1用之前拟合的Weibull参数β=1.74,η=982小时来计算:
B10 = 982 * (-ln(0.9))^(1/1.74) ≈ 982 * (0.10536)^0.5748 ≈ 982 * 0.274 ≈ 269 小时这意味着有10%的该型部件会在269小时之前失效。如果这个数字低于产品质保期,就需要改进设计或更换材料。
中位寿命B50(即F(t)=0.5)更简单:
B50 = η * (-ln(0.5))^(1/β) = 982 * (0.6931)^0.5748 ≈ 982 * 0.808 ≈ 793 小时MTTF(平均失效时间)对Weibull分布是:
MTTF = η * Γ(1 + 1/β)其中Γ是伽马函数。当β=1.74时,1+1/β≈1.575,查伽马函数表大约为0.89,所以MTTF≈982*0.89≈874小时。注意MTTF大于中位寿命,这是右偏分布的典型特征。
这些指标虽然看起来只是几个公式,但它们在设计定型、质保期设定和备件定额中非常关键。我曾经遇到一个供应商声称他们的产品B10寿命达到2000小时,仔细核对参数后才发现他们把B50当成了B10——这个误差能差出3倍多。所以算完一定要回头确认自己算的是哪个分位数。
4. 常见坑与排查记录
4.1 小样本下概率图回归的偏差
样本量小于10时,中位秩点的分布非常离散,最小二乘回归容易受极端点影响。比如只有5个失效数据,第一个点对应的F=(1-0.3)/5.4≈0.13,这个点对回归直线的影响很大,一旦首个失效时间异常偏早,β会被明显拉低。
我处理过一次只有8个样本的试制阶段数据,概率图回归得到的β=0.8,看起来像早期失效。但我怀疑数据里混入了一个安装不当导致的早夭失效。把那个异常点剔除后,β=1.6,完全不同的结论。所以小样本时务必结合工程背景判断异常点,不能盲目剔除,也不能盲目保留。如果有条件,可以用Bootstrap重新抽样评估参数的不确定性,或者使用修正的blom公式等更稳健的秩。
4.2 三参数Weibull的位置参数误用
二参数Weibull假设位置参数γ=0,也就是t=0时失效率为无穷大(当β<1时),这并不总是成立。很多产品在初始阶段不会立即失效,存在一个失效“免疫期”或最小寿命γ。此时应采用三参数Weibull:
F(t) = 1 - exp(-((t-γ)/η)^β)γ>0会让概率图的点呈现弯曲趋势,如果直接忽略γ,拟合出的β会偏小。某密封件数据在概率图上明显向下凹,强行用二参数拟合后β=0.65,物理上说不通。加入一个γ=180小时后,β恢复为1.2,曲线也变直了。
但三参数Weibull的拟合是一个优化问题,因为γ不是通过简单线性回归求得的,通常需要先假设一个γ,然后对(t-γ)做二参数拟合,再通过似然值或相关系数最大化的方式迭代搜索。迭代时注意γ不能大于最小失效时间,否则ln(t-γ)无意义。我个人不太建议轻易使用三参数,除非有明确证据表明产品存在最小寿命,因为多一个参数会让置信区间变大,且外推更不稳定。
4.3 多失效模式混合分布的处理
实际现场的失效数据往往不是单一分布能描述的。比如某电子控制器,早期焊接缺陷导致失效率下降,后期电解电容老化导致失效率上升,整体数据在概率图上呈现出明显的“拐点”,尾部上翘。此时无论用Weibull还是Lognormal单独拟合,拟合优度都不会理想,强行拟合的参数也没有实际意义。
正确的做法是先把失效机理分开。一种方法是先画概率图,看是否存在两个明显线性段。如果存在,就可以采用混合分布模型,比如两个Weibull分布的混合,或者分段拟合,但分段点的选择要结合失效分析报告,而不是仅凭统计肉眼选取。我在一个项目中,把数据按失效模式分组,分别拟合后,早期焊接失效率符合β=0.5的Weibull,后期电容老化符合β=3.2的Weibull,两个模型各管一段,才准确预测了整体失效率曲线。如果数据量够大,也可以尝试EM算法估计混合比例和各自的参数,但这属于进阶操作,工程中先做机理分组往往更可靠。
还有一个容易被忽视的细节:验证拟合优度时,不要只看相关系数R²。R²只说明线性关系的强度,对统计分布拟合来说,尾部的偏差才是决定外推是否可靠的关键。我通常会将拟合值计算的预测寿命与实际观测值做残差分析,并重点观察低F区域的残差,因为B10寿命主要受这个区间影响。
最后提一个土办法:在概率图上把点描出来,用眼睛看拐点和趋势。虽然听起来不够“量化”,但经验丰富的工程师往往能从整体形态快速判断该用单分布还是混合模型,这比盲目套用几个统计检验更高效。可靠性分析从来不是纯粹的统计题,它是统计工具和工程感知的结合。希望这篇文章里的细节能让你在下次拿到失效数据时,多一分从容,少一点踩坑。