可靠性分析威布尔三参数估计方法比较分析Word格式.docx
《可靠性分析威布尔三参数估计方法比较分析Word格式.docx》由会员分享，可在线阅读，更多相关《可靠性分析威布尔三参数估计方法比较分析Word格式.docx（14页珍藏版）》请在冰豆网上搜索。
因此,国内外一直
有人在致力于相关研究。
早期的做法通常是先由作图法得到位置参数,然后线性回归分析得形状参数和尺度参数
[8,9]
然而小样本作图法误差大,相关
研究人员找到了新的参数估计方法。
傅惠民[10]
提
出了相关系数优化法。
张秀之
[11]
将概率权重矩法
应用到海风预测方面。
然各参数估计方法适用范围和精度不同,工程选用上就有一定困难,因此对各方法进行比较,确定其适用范围和优缺点很有必要。
本文在综合分析基础上,对目前可靠性研究方面比较常用的相关系数优化法[10,12]
、概率权重矩
法
[11,1315]
与本文提出的参数拟合割线优化法进行
了比较研究,得出相关结论对在不同样本数据条件下,选用何种参数估计方法进行拟合,有一定的指导作用。
1威布尔的三参数分布模型
[115]
设随机变量X服从三参数威布尔分布,记为（,,
f（x=
-1
exp-
x
0,x<
（1
X的累积分布函数,即寿命分布函数为
F（x=1-exp-
（2
X的逆累积分布函数为
x=+-l1-F1/（3式中为形状参数,>
0;
为尺度参数,>
为位置参数,0。
2威布尔分布参数估计方法
21相关系数优化法[10,12]
相关系数优化法是一种十分有效和通用的参数估计方法。
其基本思想是先以线性相关系数为目标函数确定位置参数,然后通过线性回归分析求得形状参数和尺度参数[10],具体如下。
由X的累积分布函数可知
F（x=1-exp-
（x（4
对（4式进行变换,取对数得
ln
（1-F=-
x-
-lln
1
（1-F
=
ln-lx-（5
令
=Y,lx-=X（6ln=A,=B（7
可得Y=A+BX（8由于变量X和Y之间成线性关系,所以,可以根据已知的一组试验数据（xi,F（xi,通过式（6,变换得到新的数据（Xi,Yi,再由线性回归分析确定出待定参数A、B和系数R（X,Y分别为
A=Y-BX,B=LXY/LXX,R=LXY/XXYY（9a
Y=1
n
!
n
i=1
Yi,X=
Xi,LXX=!
ni=1X2i-nX2,
LXY=!
n1XiYi-nXY,LYY!
n1Y2i-nY2（9b
由上式可见,由于Xi,LXX,LXY均是含位置参数的函数,A,B,R也是与位置参数有关的函数。
又因所要求的位置参数必须使线性相关系数R的绝对值取最大值,所以可以由下面的方法求得位置参数。
线性相关系数R（X,Y为关于X,Y的函数,即为位置参数的函数,R对求导有
R=
LXY
XXLYY
i=1
XiYi-nXY
（!
X2i-nX2（!
Y2i-nY2
d
=0∀
2
=0（10化简整理得方程
LXY
dLXY
-
2LXX
dLXX
=0（11解上面的微分方程就可得到位置参数的值。
显然,由于式（11无法直接求解,所以必须用数值分析的方法求解得位置参数的允许精度范围内的数值解。
求得,由式（9即可得A,B,R,从而由式（7可求得形状参数和尺度参数的值。
22参数拟合最优化割线法
在文献[10]中,作者先建立含有相关系数的函数,通过逐步搜索的方法求解相关系数,再进行线性回归分析,最后求得三参数。
与之相比,本文采用一种新方法进行参数拟合。
由于文献[10]中逐步搜索的方法有一定的误差,这里提出新方法,参数拟合割线优化法能有效地减少计算误差,提高计算的精度,过程如下。
设由试验得到一组疲劳寿命的样本为（T1,T2,T3,#,Tn,先求出其相应的顺序统计量为（t1∃,t2∃,#,tn∃,样本容量为n,显然,位置参数%t1∃。
对（4式进行变换,得到一个是线性关系的变量Y与X（8式,的最优解是两变量X、Y要有（Xi,Yi最接近一条直线,此时才是最佳拟合效果[10,16]。
由式（11可知,位置参数的值可由解微分方,
6118科学技术与工程10卷
F（=LXYdLXYd-2LXXd
LXXd
LXYd（SXY-XYd-2LXXd（SXX-X2
;
d（SXY-XYd=-n!
ni=1Yiti--n!
ni=1Yi
ti-d（SXX-X2
d=-2n!
ni=1Xiti--2n!
ni=1Xi
ti-
=F（=0
（12由式（12知,需对参数求解。
在此采用和文献[10]相比更优的割线法
[17]
来
迭代求解方程!
（=0。
先设的初始值为1和0;
第j-1步的迭代值为j-1;
第j步的迭代值为j;
则第j+1步的迭代值和迭代终止条件可得知为j+1=j
（j（j-j-1
j-!
（j-
j=0,1,2,#!
（j+1%∀。
（13
式（3中∀为精度,j=j+1时方程!
（
=0即得解,退出迭代。
综上所述,用参数拟合割线优化法进行Weibull分布的三个参数计算步骤归结如下:
（1令0=t1∃
-01;
（2选择c的取值（c可取大于1的任意实数,确定递推步长h=t1∃/c,并且令j=0;
（3计算j+1=j–h;
（4若!
（j+1!
（j>
0,令j=j+1,返回（3;
（5若!
（j<
0,终止递推,且令1=j,0=j+1;
（6将1和0的结果代入（13式的迭代过程,求得;
（7由（7式求得和。
23基于MATLAB的超过概率权重矩法
对文献[11]中的方法进行优化,抛弃文献[11]中的一些近似公式,提出精度更高的基于MATLAB的超过概率权重矩法,并编写出MATLAB计算程序。
由实验得样本数据直接用MATLAB迭代即可得出所需参数,具体方法如下。
x随机变量x的概率权重矩
[7,11,14,15]
为:
Mi,j,k=E{xi
[F（x]j
[1-F（x]k
}
&
[x（F]iFj[1-F]k
dF
（14
为了避免观测值的高次乘方造成较大的抽样误差,一般取i=1。
当j和k为整数,当j=0或k=0时,概率权重矩Mi,j,k可表示为Mi,0,k（k阶超过概率权重矩或Mi,j,0（j阶不及概率权重矩。
样本（X1,X2,
X3,#,X
转变为其样本顺序统计量（x1,x2,#,
xn,则其超过概率权重矩Ms
1,0,k可重写为
M
s
1,0,k
E
xi（1-Pik
P（15
其中!
Pi为xi的频率,Pi为xi的累积频率,因此我们可知Pi的求取公式
[7,11]
Pj
i=
（n-1（n-2#（n-j
Pi=
（16
而由于本文使用超过概率权重矩M1,0,k,所以公式需转换为（1-Pi为1-Pi
k
（n-ik（
n-1k
（n-k-1（n-k-2#（n-k-i+1
（n-1（n-2#（n-i+1
（17
这样威布尔分布的k阶超过概率权重矩式为
M1,0,k=
1+k+
#（1+
1（1+k1+
1（18
式中,#为gamma函数,为便于计算这里取k=0,1,3,得M1,0,1,M1,0,2,M1,0,3解出威布尔分布三参数为
4（M1,0,3M1,0,0-M2
1,0,1
4M1,0,3+M1,0,0-4M1,0,1=
M1,0,0-
#lM1,0,0-2M1,0,1M1,0,1-2M1,0,3
ln2
ln2ln
M1,0,0-2M1,0,12（M1,0,1-2M1,0,3
（19
Pi[7,11]
为
611925期郭必柱,等:
Pi=
i=1,2,3,#,n（20易得观测样本的超过概率权重矩为
Ms1,0,0=1n!
ni=1xi
Ms1,0,1=n!
xi1-i-0n
Ms1,0,3=1
xi1-
3
（21
将MATLAB用于概率权重矩来估计分布参数,提高了估计的精度,特别是即使没有厚实数学基础的人员也可方便地直接由样本数据得到系统的可靠性或零件的寿命等分析结果。
过程为先用样本数据（sample计算样本概率权重矩,然后估计出分布参数。
下面为其主要程序代码,读者输入样本数据即得参数。
sample=[?
]%输入样本;
S=sort（sample;
num=length（sample;
pwm（1=mean（S;
forj=1:
num-1
b=0;
fori=1:
num
ifnum-i>
=j
a=prod（1:
num-i/prod（1:
j/prod（1:
num-i-j;
else
a=0;
end
b=b+a*S（i;
b=b/num;
b=b*prod（1:
j*prod（1:
num-j-1/prod（1:
num-1;
pwm（j+1=b;
pwm
4*（pwm（4*pwm（1-pwm（2*pwm（2/（4*pwm（4+pwm（1-4*pwm（2
a=4*（pwm（4*pwm（1-pwm（2*pwm（2/（4*pwm（4+pwm（1-4*pwm（2
（pwm（1-a/gamma（log（（pwm（1-2*pwm（2/（pwm（2-2*pwm（4/log（2
log（2/log（（pwm（1-2*pwm（2/（2*（pwm（2-2*pwm（4
4*（-（（*pwm（1-4*pwm（2
b=（pwm（1-a/gamma（log（（pwm（1-2*pwm（2/
（pwm（2-2*pwm（4/log（2
c=log（2/log（（pwm（1-2*pwm（2/（2*（pwm（2-2*pwm（4
%a,b,c分别表示威布尔分布三参数,,。
由程序可知,计算中未采用文献[7]中的（21式和（22式,而是直接将（19式通过MATLAB程序得出结果。
3计算实例与结果分析
算例1
某材料缺口试样（Kr在应力水平Smax=2009MPa（应力比R=075作用下,测得的一组试件的疲劳寿命数据为:
325000,376300,387300,447500,456300,524300,574400,635100（单位:
周。
用相关系数优化法求得Weibull分布的三参数为=2516647,=2502278,=17495。
参数拟合割线优化法采用中位秩拟合法时=2666836,=2274186,=17826;
采用平均秩拟合法=2515826,=2502784,=17502。
基于MATLAB的超过概率权重矩法求得=2566192,=2357030,=18995。
算例2
某一疲劳试验材料在同一应力水平下,测得一组试件的疲劳寿命（单位:
千周,如表1所示。
表1疲劳寿命试验数据
序数
疲劳寿命
Ni/kC
6120科学技术与工程10卷
参数进行估计,从算例1和算例2的计算结果和相关研究的大量计算结果都表明这3种算法都有相当的精度,可在工程上应用。
算例2还采用Theil不等系数及相关系数对计算结果进行了检验,结果对照列于表2中。
表2三种参数估计方法结果对比
算法位置
参数
尺度
形状
Theil不等
系数
相关系数
相关系数
优化法
27647053211485204140057163099841
参数拟合割
线优化法
27635973212263204240020679099967
改进的概率
权重矩法
27907113138252214010022331099915
算例1参数拟合割线优化法同时给出了中位秩平均秩的结果。
但为了便于比较,表2中相关系数优化法、参数拟合割线优化法都是按照中位秩算法计算样本分布。
下面是分析结果。
（13种方法中相关系数优化法可适用于截尾试验数据;
由（10式表明,相关系数的大小只与位置参数的估计值有关。
由表2可看出,相关系数优化法使得试验数据的拟合已具有非常好的线性度,本文提出的参数拟合割线优化法与之相比有更好的拟合效果,相关系数更大,Theil不等系数也更小,采用割线优化法有更高的精度。
而且与相关系数优化法相比,后者采用割线法求解参数方程比前者通过逐步搜索的方法有更快的收敛速度。
（2概率权重矩法直接计算各参数,与前两者相比更方便简单,文献[7,11]中采用一些近似公式后,不需要迭代计算,但它不适合极少子样情况。
当样本数量极小,如小于10时有研究表明,为负值,显然不符合事实。
本文提出的改进的概率权重矩法,通过应用MATLAB进行迭代计算,提高了精度,结果接近真值。
如算例1,为小样本的情况,计算结果与前面两个算法基本一致,可知其可行性。
因此,概率权重矩法在疲劳寿命等可靠性分析中具有重要实用价值。
尤其对于小样本,概率权重的Weibull分布参数准确度会更高[14]。
（3不同容量的样本实例计算表明,随着样本容量的增加,各估计方法之间的差异越来越小,在样本容量较大时,各种估计方法均可使用;
在样本容量较小时,各种估计方法的差异较大,基于MATLAB的超过概率权重矩法在中等样本数时精度更高,而且更便于工程人员应用,只需要样本数就可快速得出相应参数和可靠性指标。
因此,使用者需根据实际情况选择合适的估计方法。
4结论
（1提出确定Weibull分布参数的新方法,参数拟合割线优化法,并用它进行工程实例计算和与其他方法比较,确定其优劣程度以及适用范围。
（2对文献[11]中概率权重矩法加以改进,提出了更高精度,而且又更适合于工程使用的基于MATLAB的概率权重矩法。
给出了计算样本概率重矩的MATLAB语言程序,直接由样本数据可运行得出三参数值。
程序过程为:
根据试验样本值,计算各阶样本超过概率加权矩值Ms1,0,k,而样本概率加权矩是总体概率加权矩的无偏估计（Ms1,0,k=M（k,然后根据（19式可得疲劳寿命Weibull分布的参数。
（3对威布尔分布参数估计的3种方法进行比较研究,通过相关系数、Theil不等系数比较了各种方法的差异,得出各自相应的适用范围,对正确选用三参数威布尔分布参数求解方法有一定指导意义。
参考文献
1WeibullWAstatisticaldistributionfunctionofwideapplicabilityJournalofAppliedMechanics,1951;
73（9:
293297
2HasumiT,AkimotoT,AizawaYTheWeibulllogWeibulldistributionforinteroccurrencetimesofearthquakesPhysicaA,2009;
388（4:
491498
3高镇同疲劳应用统计学北京:
国防工业出版社,1986
4TiryakioluM,HudakDOnestimatingWeibullmodulusbythelinearregressionmethodJournalofMaterialsScience,2007;
42（24:
6121
25期郭必柱,等:
6122科学技术与工程10卷概率权重矩法及其在Weibull分布参数估计中的应5WeibullWFatiguetestingandanalysisofresultsNewYorkMacmil.:
lanCompany,19616LittleREJebeEHStatisticaldesignoffatigueexpermentsLondon,i,EnglandAppliedSciencePublishers1975:
7GreenwoodJA,LandwehrJM,MatalasNC,etalProbabilityweightedmomentsdefinitionandrelationtoparametersofseveraldistribu:
tionsexpressableininverseforWaterResourcesResearch,1979;
m15（5:
1049–10548高木升社,19779徐灏概率疲劳沈阳:
东北大学出版社,1994可靠性基础数学广五所,译北京:
国防工业出版11张秀之用海洋预报,1994;
11（3:
56–6112史景钊,蒋国良用相关系数法估计威布尔分布的位置参数河南农业大学学报,1995;
29（2:
167–17113宋松柏P∋分布参数的概率权重矩法S函数计算200828（5:
1;
14邓建5上海:
同济大水文,地下结构可靠度理论与应用研究学,
升级会员 