徐家進(jìn)

(力中國(guó)際融資租賃有限公司，廣州 510620）

《疲勞應(yīng)用統(tǒng)計(jì)學(xué)》［1］是高鎮(zhèn)同先生在30多年前寫的一本書，將統(tǒng)計(jì)學(xué)正式引入了解決各種疲勞問(wèn)題的工程實(shí)踐之中，建立了疲勞統(tǒng)計(jì)學(xué)這一學(xué)科［2］。

疲勞統(tǒng)計(jì)學(xué)是疲勞可靠性的理論基礎(chǔ)，對(duì)于解決各類可靠性問(wèn)題，特別是航空工業(yè)領(lǐng)域，其意義無(wú)論怎樣強(qiáng)調(diào)都不為過(guò)。在疲勞統(tǒng)計(jì)學(xué)中，有一個(gè)非常獨(dú)特而重要的統(tǒng)計(jì)分布就是威布爾分布。事實(shí)上，威布爾分布是疲勞統(tǒng)計(jì)學(xué)的一個(gè)最主要的內(nèi)容之一。威布爾分布一個(gè)非常重要的特點(diǎn)就是具有最小壽命這個(gè)參數(shù)。但高鎮(zhèn)同先生指出，最小壽命這個(gè)提法是不夠全面的，因?yàn)榘凑胀紶柗植嫉亩x，應(yīng)用于疲勞壽命時(shí)可明確為100%的可靠度之下的安全壽命。此提法是科學(xué)的，以后在疲勞統(tǒng)計(jì)學(xué)中就將這個(gè)參數(shù)改稱為安全壽命。而這個(gè)參數(shù)對(duì)于像飛機(jī)這樣產(chǎn)品的定壽起了關(guān)鍵性的作用。問(wèn)題是威布爾分布的數(shù)學(xué)形式比較復(fù)雜，特別是在當(dāng)時(shí)計(jì)算機(jī)的性能、使用和普及遠(yuǎn)遠(yuǎn)不如現(xiàn)在的情況下，使用起來(lái)比較困難。盡管目前已有了不少關(guān)于威布爾分布應(yīng)用的軟件，但絕大部分是二參數(shù)威布爾分布［3］，恰恰將安全壽命這個(gè)參數(shù)歸零了。這樣的簡(jiǎn)化顯然是不合情理的。同時(shí)，高鎮(zhèn)同先生強(qiáng)調(diào)威布爾分布是一種全態(tài)分布，既包括正態(tài)分布，也包括偏態(tài)分布。但因?yàn)橥紶柗植嫉膹?fù)雜性影響了它的推廣和使用。因此，疲勞統(tǒng)計(jì)學(xué)智能化的當(dāng)務(wù)之急就是要解決威布爾分布3個(gè)參數(shù)的估計(jì)問(wèn)題。本文提出了與已有估計(jì)威布爾分布3個(gè)參數(shù)的2個(gè)方法不同的智能化解決方案

1 疲勞統(tǒng)計(jì)學(xué)的特點(diǎn)

疲勞統(tǒng)計(jì)學(xué)是將統(tǒng)計(jì)學(xué)原理應(yīng)用到結(jié)構(gòu)疲勞這個(gè)領(lǐng)域中去。三參數(shù)威布爾分布可以更好地描寫結(jié)構(gòu)疲勞壽命的特點(diǎn)，在疲勞統(tǒng)計(jì)學(xué)中具有非常重要的地位。但三參數(shù)威布爾分布并不為大家所熟知，這是因?yàn)槠鋽?shù)學(xué)形式比較復(fù)雜，而且通過(guò)樣本來(lái)確定其3個(gè)參數(shù)比較困難，應(yīng)用起來(lái)也就有一定的難度［4-6］。自然也與對(duì)威布爾分布的意義認(rèn)識(shí)不足有關(guān)，為此有必要先簡(jiǎn)要介紹一下威布爾分布的特點(diǎn)。

2 威布爾分布的特點(diǎn)

威布爾分布概率密度函數(shù)［1］為

式中：N為疲勞壽命；N0為安全壽命；b為形狀參數(shù)；Na為特征壽命；xa為特征參數(shù)；x0為位置參數(shù)；且0≤x0≤x＜∞。

為方便可設(shè)：

式中：λ為尺度或比例參數(shù)。

從統(tǒng)計(jì)分布角度看，由尺度參數(shù)λ代替特征參數(shù)更有一般性。這是因?yàn)閤0＝0時(shí)，λ＝xa。即λ可看成x0＝0時(shí)的特征參數(shù)。因此，將式(3）代替式(1b）作為一般的三參數(shù)威布爾分布表達(dá)式。三參數(shù)分別為形狀參數(shù)b、尺度參數(shù)λ和位置參數(shù)x0。

進(jìn)一步，若設(shè)x0＝0，式(3）則為兩參數(shù)威布爾分布概率密度函數(shù)。

再設(shè)λ＝1就成了“標(biāo)準(zhǔn)”(一個(gè)參數(shù)）的威布爾分布概率密度函數(shù)。

對(duì)于標(biāo)準(zhǔn)化的威布爾分布，b分別為0.5，2，3.5，5時(shí)，威布爾分布概率密度函數(shù)如圖1所示。

由威布爾分布概率密度函數(shù)及圖1可看出其有如下特點(diǎn)：

圖1 標(biāo)準(zhǔn)威布爾分布概率密度函數(shù)Fig.1 Standard Weibull distribution probability density function

1）b為威布爾分布的形狀參數(shù)是名至實(shí)歸的。0＜b＜1時(shí)類似1／x函數(shù)，而1＜b＜3是左偏分布，3＜b＜4則近似正態(tài)分布，b＞4則為右偏分布。此外，當(dāng)b＝1時(shí)為指數(shù)分布，b＝2時(shí)為瑞利分布。這也是高鎮(zhèn)同先生為什么稱威布爾分布為“全態(tài)分布”的主要原因。

2）威布爾密度函數(shù)有3個(gè)參數(shù)，比正態(tài)分布多1個(gè)，這是其優(yōu)點(diǎn)(可描寫安全壽命）但也是其缺點(diǎn)(數(shù)學(xué)形式比較復(fù)雜）的來(lái)源。正態(tài)分布能很好地描寫對(duì)稱狀態(tài)的數(shù)據(jù)，問(wèn)題是在現(xiàn)實(shí)生活中真正對(duì)稱的數(shù)據(jù)只是一種理想狀態(tài)，因此正態(tài)分布往往只是一階近似。而威布爾分布既能描寫對(duì)稱狀態(tài)(盡管是近似的），又能描寫左偏或右偏的數(shù)據(jù)，故其應(yīng)用范圍要比正態(tài)分布大得多。

上文對(duì)威布爾分布做了定性的描述，現(xiàn)進(jìn)一步對(duì)其做定量的研究。為方便，開始時(shí)僅對(duì)標(biāo)準(zhǔn)威布爾分布進(jìn)行比較詳細(xì)的研究，而二參數(shù)和三參數(shù)的威布爾分布均可得到類似的結(jié)果。

1）威布爾分布概率密度函數(shù)是滿足歸一化條件的。即要證明：

可設(shè)z＝xb，式(6）左邊可變?yōu)?/p>

對(duì)于兩參數(shù)和三參數(shù)的威布爾分布同樣可得到這個(gè)結(jié)論，只要分別設(shè)：

2）容易證明，威布爾分布函數(shù)為

事實(shí)上，和1）相同可設(shè)z＝xb：

這個(gè)分布函數(shù)恰恰就是破壞率，而可靠度P為

式(9a）給出了可靠度的函數(shù)形式，對(duì)于應(yīng)用有重大意義。即xp為與可靠度P相應(yīng)x的值。兩參數(shù)和三參數(shù)威布爾分布分別為

若此時(shí)取xp＝xa，xp－x0＝(x0＋λ）－x0＝λ，則

這表明任何威布爾分布當(dāng)其值為特征參數(shù)時(shí)可靠度一定為36.8%，這是特征參數(shù)xa的一個(gè)物理意義。而取xp＝x0時(shí)P＝100%，也是參數(shù)x0為100%可靠度的安全壽命的由來(lái)。順便指出，對(duì)于正態(tài)分布不可能有這個(gè)概念，因其左側(cè)與水平軸不相交。

3）當(dāng)b＞1時(shí)密度函數(shù)是凸函數(shù)，因此必有一個(gè)極(大）值。事實(shí)上：

令f′(x）＝0，即有b－1－bxb＝0。

再求f(x）的二階導(dǎo)數(shù)：

因f″(xmax）＝b exp(－xb）xb－3b(1－b）＜0，即證明了xmax為f(x）的極大值，且f(x）是凸函數(shù)。由式(11a）不難發(fā)現(xiàn)，當(dāng)b越大時(shí)，ln［(b－1）／b］→0，xmax→1，峰值集中在x＝1處。

而對(duì)于兩參數(shù)和三參數(shù)威布爾分布概率密度函數(shù)的極大值點(diǎn)位置分別為

4）威布爾分布的數(shù)學(xué)期望按照定義且設(shè)z＝xb可得

再注意到伽馬函數(shù)的定義：

類似地，對(duì)于兩參數(shù)和三參數(shù)威布爾分布的數(shù)學(xué)期望分別為

5）威布爾分布的方差。

同樣按照定義可有

由式(13）可得

類似地，對(duì)于兩參數(shù)和三參數(shù)威布爾分布的方差分別為

很明顯，三參數(shù)的威布爾分布方差的期望值與x0沒(méi)有直接關(guān)系。

6）關(guān)于威布爾分布和正態(tài)分布相似性的討論。前面已提到當(dāng)3＜b＜4時(shí)2種分布是相似的。正態(tài)分布意味著完全對(duì)稱性，即均值、中值、峰值三者合一，而威布爾分布則不具備這個(gè)特點(diǎn)。不過(guò)在一定條件下，可讓這3個(gè)值中的某2個(gè)值一致，如可令均值和中值相等，因此在這個(gè)意義上可認(rèn)為與正態(tài)分布相似。均值可由式(12a）給出，而中值則可由式(9a）中P＝50%時(shí)給出，0.5＝exp(－），x50為50%可靠度的x值，可由式(12a）代替。

式(16）為超越方程，用Excel來(lái)解較麻煩。利用Python能很快得到結(jié)果。但不管哪一個(gè)方法都是利用牛頓二分法。為此可將式(16）改寫為

利用Excel可得b＝3.44(精確到10－5）。

Python中運(yùn)行的結(jié)果為：k(對(duì)分次數(shù)）＝19，E(精度）＝1×10－8，b＝3.439 54。

Excel的優(yōu)點(diǎn)在于直觀，但在計(jì)算過(guò)程中需要人的介入，且效率較低。Python效率高，代碼一旦調(diào)試成功，全部工作都可自動(dòng)完成，精度也高。

順便指出，這個(gè)結(jié)論雖然是從標(biāo)準(zhǔn)威布爾分布推出的，但對(duì)于兩參數(shù)和三參數(shù)分布都是適用的，因?qū)ΨQ性僅僅與形狀參數(shù)b有關(guān)。另外按照文獻(xiàn)［7］，威布爾分布的偏態(tài)系數(shù)為

不難看出，這個(gè)偏態(tài)系數(shù)也僅僅與b有關(guān)。同時(shí)，偏態(tài)系數(shù)為零時(shí)也應(yīng)該是對(duì)稱性最好的時(shí)候。亦是和正態(tài)分布最接近的時(shí)候，這對(duì)于比較這2種分布的異同意義重大，為此有必要求出這個(gè)b值。即要解如下方程：

式(19）為一個(gè)關(guān)于b的超越方程，可利用相同Python代碼得：k(對(duì)分次數(shù)）＝26，E(精度）＝1×10－8，b＝3.602 35。這與前文的b＝3.44比較接近，但并不重合，這也表明威布爾分布不可能完全對(duì)稱。

3 疲勞統(tǒng)計(jì)學(xué)的智能化

20世紀(jì)80年代，計(jì)算機(jī)在中國(guó)剛剛開始起步和普及，做一些簡(jiǎn)單的統(tǒng)計(jì)工作沒(méi)有問(wèn)題，但在工程實(shí)踐中應(yīng)用較少。一方面是因?yàn)橛布恍?，另外一方面懂得編程的工程技術(shù)人員不多。因此，各種現(xiàn)成的圖表，如正態(tài)概率坐標(biāo)紙、威布爾分布概率坐標(biāo)紙等還是廣為使用。這個(gè)方法雖然看起來(lái)比較簡(jiǎn)單實(shí)用，但誤差較大，且必須將數(shù)據(jù)對(duì)數(shù)化，這從數(shù)學(xué)的角度看是一種空間變換，而這樣一變換很可能會(huì)將這2種分布變得看起來(lái)不可區(qū)分。更加重要的是，這種方法難以將擬合程度量化。只能夠看起來(lái)擬合不錯(cuò)，但實(shí)際上并非如此。所謂智能化，就是要將得到的數(shù)據(jù)直接通過(guò)計(jì)算機(jī)來(lái)得到人們要求的統(tǒng)計(jì)結(jié)果，即不僅是由計(jì)算機(jī)畫出直觀的圖表，而且還有給出判斷的定性與定量的結(jié)果?，F(xiàn)舉例說(shuō)明。

例1取一組試件疲勞壽命如表1所示［1］。

表1 一組疲勞壽命數(shù)據(jù)［1］Table 1 Aset of fatigue life data［1］

表1中：cycle表示循環(huán)次數(shù)，而均值Nav、標(biāo)準(zhǔn)差s及中值Nm皆由Excel得出。若這些數(shù)據(jù)是服從正態(tài)分布的，則相應(yīng)的正態(tài)分布參數(shù)的估計(jì)值可由式(20）給出：

式中：“＾”表示估計(jì)值。

于是，相應(yīng)的正態(tài)分布密度函數(shù)為

現(xiàn)若假定它們服從三參數(shù)威布爾分布，則求其3個(gè)參數(shù)，由式(9c）、式(12b）及式(15c）分別得到

式中：N0為估計(jì)安全壽命；Na為估計(jì)特征壽命。

從式(24）得

再由式(23）得

再將式(26）～式(28）代入式(22）可得

由以上參數(shù)和式(1a）得出相應(yīng)的威布爾分布概率密度函數(shù)如下：

利用Excel，將由該數(shù)據(jù)得到的正態(tài)分布式(22）和威布爾分布式(30）的可靠度在圖2中做一個(gè)比較。

但從圖2中很難回答哪一個(gè)分布的可靠度估計(jì)更好一些。為此考慮到可靠度估計(jì)值［8-9］(也可參考文獻(xiàn)［1］）：

圖2 正態(tài)和威布爾分布可靠度比較(例1）Fig.2 Comparison of normal and Weibull distribution reliability(Example 1）

式中：i為數(shù)據(jù)(觀測(cè)值）由小到大排列的序數(shù)；n為數(shù)據(jù)的個(gè)數(shù)。

不過(guò)仍然可從另外一個(gè)角度來(lái)看正態(tài)分布是否被滿足，即卡方檢驗(yàn)，但此法相當(dāng)麻煩，可參考文獻(xiàn)［1］。

例2利用文獻(xiàn)［1］表12-3中的數(shù)據(jù)繪制表2。

表2 100個(gè)試件在同一應(yīng)力條件下疲勞壽命數(shù)據(jù)Table 2 Fatigue life data of 100 specimens under the same stress condition

對(duì)于例2，仍然可和例1一樣畫出正態(tài)分布和威布爾分布的可靠度擬合圖，以及給出兩者對(duì)于理想可靠度的決定系數(shù)［6］。不過(guò)，因?yàn)閿?shù)據(jù)較多，采用Excel很困難，但用Python處理較容易。Python代碼運(yùn)行結(jié)果為：Nav＝5.315，s＝1.289，Nm＝5.07，Cs＝1.021，k(對(duì)分次數(shù)）＝20，E(精度）＝1×10－7，b＝1.526，N0＝3.39，Na＝5.53(102cycle）。

驗(yàn)算：Nav1＝5.315，s＝1.289，Nm1＝5.07，正態(tài)分布與理想可靠度的決定系數(shù)為0.979 14，威布爾分布與理想可靠度的決定系數(shù)為0.988 44。

在此要注意這里用的是決定系數(shù)［10］而不是相關(guān)系數(shù)來(lái)區(qū)分?jǐn)M合水平，主要是因?yàn)閳D3的橫坐標(biāo)變成了疲勞壽命，這樣其與可靠度的關(guān)系不再是線性關(guān)系。計(jì)算結(jié)果表明，威布爾分布的決定系數(shù)仍然比正態(tài)分布大，盡管不是大很多。最重要的是，這些數(shù)據(jù)是偏態(tài)的，用正態(tài)分布已經(jīng)不能很好描述。而威布爾分布恰恰彌補(bǔ)了這個(gè)不足，如圖3所示。不過(guò)，也必須指出用這個(gè)解析法來(lái)求威布爾分布的3個(gè)參數(shù)并非沒(méi)有瑕疵。例如，經(jīng)過(guò)計(jì)算，若這些數(shù)據(jù)滿足威布爾分布，那么安全壽命是3.39，但開始3個(gè)數(shù)據(jù)是小于這個(gè)值的，則計(jì)算過(guò)程沒(méi)有問(wèn)題，很可能是因?yàn)檫@個(gè)方法先將形狀參數(shù)b找出來(lái)后再計(jì)算N0和Na，就無(wú)法保證N0一定會(huì)小于輸入的數(shù)據(jù)的最小值。故需要對(duì)威布爾分布參數(shù)算法做進(jìn)一步的研究。

圖3 正態(tài)和威布爾分布可靠度比較(例2）Fig.3 Comparison of reliability between Normal and Weibull distribution(Example 2）

4 高鎮(zhèn)同法

第3節(jié)指出計(jì)算威布爾分布的解析算法還存在一些問(wèn)題，如計(jì)算出來(lái)的安全壽命比實(shí)際的數(shù)據(jù)還要大，說(shuō)明了這個(gè)算法是不自冾的。這里存在2種可能性：一是這些數(shù)據(jù)并不那么符合威布爾分布；二是上述計(jì)算威布爾參數(shù)的算法還存在問(wèn)題，即用樣本估計(jì)的均值和方差值來(lái)代替總體相應(yīng)的均值和方差存在較大誤差。為此，可從另外一個(gè)角度來(lái)計(jì)算這些參數(shù)，如用最小二乘法。為簡(jiǎn)單，從二參數(shù)的威布爾分布開始，由式(9c）可有

其估計(jì)值可利用式(31），式(32）在取二次自然對(duì)數(shù)后可得

于是，可通過(guò)最小二乘法來(lái)求出b、d及λ。但問(wèn)題并非那么簡(jiǎn)單，因?yàn)橛枚?shù)威布爾分布的前提是設(shè)位置參數(shù)或安全壽命x0＝0。在實(shí)際情況下，x0≠0，而且恰恰是因?yàn)檫@個(gè)x0≠0才顯示出威布爾分布的優(yōu)越性。事實(shí)上，系數(shù)b和λ及相應(yīng)的相關(guān)系數(shù)r都是x0的函數(shù)，可通過(guò)解析法來(lái)求出使得相關(guān)系數(shù)r的極大值，但這種方法推導(dǎo)較麻煩，容易出錯(cuò)。用Python在0≤x0＜xmin區(qū)間中(這里的xmin就取給定數(shù)據(jù)的最小值）可直接將r關(guān)于x0的圖畫出來(lái)，如圖4所示。然后用Python智能地將相關(guān)系數(shù)最大的r的x0找出來(lái)。理論上分為兩步，但在實(shí)際編的代碼中是一氣呵成的，即以x0為變量，找出r的極大值同時(shí)將b和λ確定，這就是高鎮(zhèn)同法。具體表述如下：

圖4 安全壽命與相關(guān)系數(shù)的關(guān)系(例3）Fig.4 Safety life versus correlation coefficient(Example 3）

1）輸入原始數(shù)據(jù)，若原始數(shù)據(jù)沒(méi)有排好序，則先排序。

2）利用Python中的scipy可直接通過(guò)以給定精度的間隔來(lái)遍歷x0的可能值區(qū)間［0，xmin］，以求出使得相關(guān)系數(shù)最大的x0，即x0max。

3）注意到，在scipy中計(jì)算相關(guān)系數(shù)時(shí)，事實(shí)上是先求出最小二乘法中直線方程的相應(yīng)系數(shù)b和d，推出λ＝exp(－d／b）。因此一旦定出了x0max，則威布爾分布相應(yīng)的參數(shù)b與λ也就同時(shí)得到了。

例3現(xiàn)利用例1的數(shù)據(jù)在Python代碼上進(jìn)行試算，可得到如下結(jié)果：N＝｛350，380，400，430，450，470，480，500，520，540，550，570，600，610，630，650，670，730，770，840｝，Nav＝557.0，s＝132.152，Nm＝545.0，Cs＝0.408，r＝0.999 218，bm＝2.040，λ＝320.98，N0＝276.60，Cs＝0.605，Nav1＝560.97，s＝146.04，Nm1＝544.79，解析法威布爾分布與理想可靠度的相關(guān)系數(shù)為0.999 20，高鎮(zhèn)同法威布爾分布與理想可靠度的相關(guān)系數(shù)為0.999 22。

例4高鎮(zhèn)同法與解析法最大不同是b值不同(前者b＝2.21，這里b＝2.040，相差8%），從而標(biāo)準(zhǔn)差s不同(132 vs 146，相差10%）。但相關(guān)系數(shù)r(0.999 20 vs 0.999 22）只有10－5，在這個(gè)意義上兩者并沒(méi)有什么大的差別。因解析法對(duì)例2的解是不自冾的，對(duì)Python的代碼稍加修改可得如下結(jié)果：Nav＝5.32，s＝1.29，Nm＝5.07，km＝2 780，r＝0.993 763，bm＝2.147，λ＝2.87，N0＝2.780，Nav1＝5.32；s1＝1.25，Nm1＝5.20，正態(tài)分布與理想可靠度的決定系數(shù)為0.979 11，威布爾分布與理想可靠度的決定系數(shù)為0.989 03。

圖5形象地給出了例4中的數(shù)據(jù)如何使用高鎮(zhèn)同法。先求出相關(guān)系數(shù)和安全壽命N0的關(guān)系曲線，再找出使得相關(guān)系數(shù)最大的安全壽命。

圖5 安全壽命與相關(guān)系數(shù)的關(guān)系(例4）Fig.5 Safety life versus correlation coefficient(Example 4）

由此可見(jiàn)，高鎮(zhèn)同法的結(jié)果與前面的解析法明顯不同，即不存在超過(guò)安全壽命N0的原始數(shù)據(jù)。幾個(gè)參數(shù)都較符合，即使是標(biāo)準(zhǔn)差也只有3%的誤差，且高鎮(zhèn)同法得到的相關(guān)系數(shù)明顯大于解析法。在此意義上，可認(rèn)為高鎮(zhèn)同法是優(yōu)于解析法的。

例5值得注意的是，高鎮(zhèn)同法與文獻(xiàn)［11］中提出的“確定威布爾分布三參數(shù)相關(guān)系數(shù)優(yōu)化法”還是有所不同的，主要體現(xiàn)在：文獻(xiàn)［11］中的數(shù)學(xué)推導(dǎo)較麻煩，相應(yīng)的代碼也會(huì)較復(fù)雜，仍然沒(méi)有充分發(fā)揮出利用計(jì)算機(jī)智能的優(yōu)勢(shì)。將文獻(xiàn)［11］中表3的數(shù)據(jù)放入同樣的Python代碼即得如下結(jié)果：N＝｛325，376.3，387.3，447.5，456.3，524.3，574.4，635.1｝，Nav＝465.77，s＝105.763，Nm＝451.9，Cs＝0.302，r＝0.992 976，bm＝1.747，λ＝250.00，N0＝251.84，Cs＝0.823，Nav1＝474.52，s1＝131.51，Nm1＝454.54。

圖6給出了例5如何使用高鎮(zhèn)同法(參照文獻(xiàn)［11］中的表3數(shù)據(jù)繪制圖6）。而得到的結(jié)果與文獻(xiàn)［11］中的結(jié)果幾乎一樣，但要注意在文獻(xiàn)［11］中得到的相關(guān)系數(shù)是負(fù)的，而這里是正的，主要是因?yàn)樵谖墨I(xiàn)［11］中設(shè)Y＝－ln［ln(1／p）］。

圖6 安全壽命與相關(guān)系數(shù)的關(guān)系(文獻(xiàn)［11］中的表3）Fig.6 Safety life versus correlation coefficient(Table 3 of Ref.［11］）

由此可見(jiàn)，本文方法相對(duì)簡(jiǎn)單，容易掌握，不會(huì)出現(xiàn)矛盾的結(jié)果。這表明高鎮(zhèn)同法具有相當(dāng)?shù)膬?yōu)越性，值得推廣。

5 高鎮(zhèn)同法的應(yīng)用

高鎮(zhèn)同法對(duì)于擬合三參數(shù)疲勞性能曲線也是可以應(yīng)用的。因?yàn)檫@與求威布爾分布的三參數(shù)在數(shù)學(xué)上是沒(méi)有區(qū)別的。即用相關(guān)系數(shù)最大化來(lái)定出位置的同時(shí)，再由最小二乘法求另2個(gè)參數(shù)，即可智能化地同時(shí)得到3個(gè)參數(shù)。

例6文獻(xiàn)［7］中給出了一個(gè)三參數(shù)冪函數(shù)表達(dá)式：

式中：S0、m和C為材料常數(shù)；Smax為最大應(yīng)力。

文獻(xiàn)［7］雖然也是用最小二乘法，利用相關(guān)系數(shù)最大化來(lái)求出S0，再去求另外2個(gè)參數(shù)，但較麻煩?，F(xiàn)介紹高鎮(zhèn)同法如何用應(yīng)在例6中。設(shè)S0為已知，在式(37）兩邊取10為底的對(duì)數(shù)：

再設(shè)：Y ＝lg(Smax－S0），X ＝lg N。

即可得

按照高鎮(zhèn)同法，可將S0根據(jù)需要而定的精度遍歷區(qū)間［0，Smax］，以求出使得相關(guān)系數(shù)最大的S0，同時(shí)得出相應(yīng)的b和d，再由式(39）定出m和C。結(jié)果與文獻(xiàn)［7］幾乎相同。以文獻(xiàn)［7］第11章例4中的數(shù)據(jù)為例，繪制表3。

表3 例6的一組數(shù)據(jù)Table 3 Exam ple 6 aset of data

其Python代碼與高鎮(zhèn)同法類似(要注意此時(shí)相關(guān)系數(shù)是負(fù)的），可得如下結(jié)果：km＝27 089，r＝－0.999 14，S0＝270.89，m ＝2.146，C＝9.599 4×106。

由圖7可知，高鎮(zhèn)同法得到的曲線與實(shí)際值擬合得很好。圖中：S0＝270.89，m＝2.146 4，C＝9.599×106，r＝－0.999 14。其結(jié)果和文獻(xiàn)［7］的結(jié)果：r＝－0.999 1，S0＝270.89，m＝2.142 5，C＝9.444 5×106，除了C相差比較大一點(diǎn)(不超過(guò)1.7%），其他參數(shù)的相對(duì)誤差都不到10－3。但用高鎮(zhèn)同法計(jì)算較方便。

圖7 S-N擬合曲線Fig.7 S-Ncurve fiting

例7以文獻(xiàn)［7］中第11章例5中的數(shù)據(jù)為例，繪制表4。

表4 例7的一組數(shù)據(jù)Table 4 Example 7 aset of data

表4中數(shù)據(jù)服從如下經(jīng)驗(yàn)公式：

式中：a為試件出現(xiàn)裂縫的長(zhǎng)度；C、m、a0(也稱為初始裂縫長(zhǎng)度）都是與材料有關(guān)的常數(shù)。

很明顯，式(41）與式(37）是不一樣的，物理意義更加不同。取對(duì)數(shù)后：

式中：Y＝lg(a－a0），X＝lg N。

同樣用高鎮(zhèn)同法，以視需要而定的精度使a0遍歷區(qū)間［0，amin］，可求出使得相關(guān)系數(shù)r最大的a0。再根據(jù)式(43）定出C和m。將Python的代碼做出適當(dāng)?shù)男薷目傻萌缦陆Y(jié)果：r＝0.992 03，a0＝0.173，m＝0.343，C＝6.989 6×103。

圖8表明，由高鎮(zhèn)同法計(jì)算出來(lái)的結(jié)果與實(shí)際值擬合相當(dāng)好。圖中：a0＝0.173，m＝0.343 10，C＝6.989 6×103，r＝0.992 03。其結(jié)果與文獻(xiàn)［7］的結(jié)果：r＝0.992 0，a0＝0.176，m ＝0.336 5，C＝6 976相差非常小。就相關(guān)系數(shù)而言，其差別可以忽略不計(jì)，如果僅從圖形上看，兩者也難以得出有多大差別。而對(duì)于其他系數(shù)，相對(duì)誤差都不超過(guò)2%。即相關(guān)系數(shù)10－4的差異可能會(huì)引起擬合曲線各參數(shù)的差異擴(kuò)大上百倍。這個(gè)結(jié)論也符合三參數(shù)威布爾分布，與取對(duì)數(shù)后的線性化是有關(guān)的。

圖8 a-N擬合曲線Fig.8 a-N curve fitting

最后，要對(duì)高鎮(zhèn)同法的由來(lái)再做一個(gè)簡(jiǎn)要的說(shuō)明。①高鎮(zhèn)同法的基礎(chǔ)嚴(yán)格來(lái)講應(yīng)該是最優(yōu)擬合，最小二乘法只是其特例。在文獻(xiàn)［11］中已給出了一些很好的例子。②在文獻(xiàn)［11］的基礎(chǔ)上，不少學(xué)者［12-15］都做了不少改進(jìn)工作，取得了一定的成績(jī)，特別是回看文獻(xiàn)［15］離高鎮(zhèn)同法只有一步之遙，很可惜沒(méi)有再深入下去。③為了解決解析法出現(xiàn)的不自冾問(wèn)題，不通過(guò)直接求相關(guān)系數(shù)極大值，而是利用scipy這個(gè)軟件庫(kù)將所有可能的相關(guān)系數(shù)求出來(lái)之后通過(guò)排序?qū)⒆畲笾登蟪鰜?lái)，同時(shí)也將相應(yīng)的威布爾參數(shù)求了出來(lái)。且從數(shù)學(xué)角度看此法對(duì)于負(fù)的相關(guān)系數(shù)也是同樣適用的，很快就能推廣到三參數(shù)疲勞性能曲線的擬合中使用。

6 結(jié)論

1）威布爾分布是一種全態(tài)分布，比正態(tài)分布更具有一般性，可以相信在大數(shù)據(jù)時(shí)代將發(fā)揮出更大的作用。特別是應(yīng)用到疲勞統(tǒng)計(jì)學(xué)中具有100%可靠度的安全壽命，簡(jiǎn)稱為安全壽命這個(gè)參數(shù)，其涉及可靠性的實(shí)際應(yīng)用領(lǐng)域中具有重大意義。

2）將疲勞統(tǒng)計(jì)學(xué)智能化的一個(gè)重要切入點(diǎn)是將估計(jì)威布爾分布三參數(shù)智能化。威布爾分布的數(shù)學(xué)形式雖然比較復(fù)雜，但完全可用像高鎮(zhèn)同法(充分利用Python的特點(diǎn)而創(chuàng)造出來(lái)的一種新算法）智能化地克服圖解法和解析法的弱點(diǎn)，很方便地解決威布爾分布三參數(shù)的計(jì)算問(wèn)題，不僅如此，同時(shí)可解決廣義三參數(shù)疲勞性能曲線的擬合問(wèn)題。

3）疲勞統(tǒng)計(jì)學(xué)智能化不僅是實(shí)際應(yīng)用中必不可少的一個(gè)工具，在理論研究中也是功不可沒(méi)。例如，如何確定威布爾分布的無(wú)偏性和對(duì)稱性，若沒(méi)有計(jì)算機(jī)智能的幫助其計(jì)算復(fù)雜性也是令人卻步的。

4）本文提倡的疲勞統(tǒng)計(jì)學(xué)中的智能化主要是指：①不再需要各種各樣的有關(guān)統(tǒng)計(jì)函數(shù)值的表，都可通過(guò)計(jì)算機(jī)得到；②所有的數(shù)據(jù)處理和圖表都可由計(jì)算機(jī)自動(dòng)完成；③只要給出(疲勞壽命）樣本數(shù)據(jù)通過(guò)計(jì)算機(jī)就可判斷該數(shù)據(jù)的總體服從的最佳分布(正態(tài)分布或威布爾分布），并同時(shí)給出該分布的參數(shù)。特別是高鎮(zhèn)同法可使得威布爾分布的應(yīng)用更加方便，對(duì)于其普及具有重大意義。

當(dāng)然上面說(shuō)的智能化還是處于初級(jí)階段，這里有兩層意思：一是疲勞統(tǒng)計(jì)學(xué)的智能化水平本身還不夠高，二是使用的人還不夠廣泛，盡管Python已經(jīng)非常容易學(xué)習(xí)和使用。隨著計(jì)算機(jī)在疲勞統(tǒng)計(jì)學(xué)中的使用越益廣泛和深入，可以相信達(dá)到比較成熟的智能化階段并非遙不可及。

致謝感謝力中國(guó)際融資租賃有限公司萬(wàn)偉浩先生對(duì)于本文及有關(guān)研究工作的支持。