疲勞統計學智能化中的高鎮同法

徐家進

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

《疲勞應用統計學》［1］是高鎮同先生在30多年前寫的一本書，將統計學正式引入了解決各種疲勞問題的工程實踐之中，建立了疲勞統計學這一學科［2］。

疲勞統計學是疲勞可靠性的理論基礎，對于解決各類可靠性問題，特別是航空工業領域，其意義無論怎樣強調都不為過。在疲勞統計學中，有一個非常獨特而重要的統計分布就是威布爾分布。事實上，威布爾分布是疲勞統計學的一個最主要的內容之一。威布爾分布一個非常重要的特點就是具有最小壽命這個參數。但高鎮同先生指出，最小壽命這個提法是不夠全面的，因為按照威布爾分布的定義，應用于疲勞壽命時可明確為100%的可靠度之下的安全壽命。此提法是科學的，以后在疲勞統計學中就將這個參數改稱為安全壽命。而這個參數對于像飛機這樣產品的定壽起了關鍵性的作用。問題是威布爾分布的數學形式比較復雜，特別是在當時計算機的性能、使用和普及遠遠不如現在的情況下，使用起來比較困難。盡管目前已有了不少關于威布爾分布應用的軟件，但絕大部分是二參數威布爾分布［3］，恰恰將安全壽命這個參數歸零了。這樣的簡化顯然是不合情理的。同時，高鎮同先生強調威布爾分布是一種全態分布，既包括正態分布，也包括偏態分布。但因為威布爾分布的復雜性影響了它的推廣和使用。因此，疲勞統計學智能化的當務之急就是要解決威布爾分布3個參數的估計問題。本文提出了與已有估計威布爾分布3個參數的2個方法不同的智能化解決方案

1 疲勞統計學的特點

疲勞統計學是將統計學原理應用到結構疲勞這個領域中去。三參數威布爾分布可以更好地描寫結構疲勞壽命的特點，在疲勞統計學中具有非常重要的地位。但三參數威布爾分布并不為大家所熟知，這是因為其數學形式比較復雜，而且通過樣本來確定其3個參數比較困難，應用起來也就有一定的難度［4-6］。自然也與對威布爾分布的意義認識不足有關，為此有必要先簡要介紹一下威布爾分布的特點。

2 威布爾分布的特點

威布爾分布概率密度函數［1］為

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

為方便可設：

式中：λ為尺度或比例參數。

從統計分布角度看，由尺度參數λ代替特征參數更有一般性。這是因為x0＝0時，λ＝xa。即λ可看成x0＝0時的特征參數。因此，將式(3）代替式(1b）作為一般的三參數威布爾分布表達式。三參數分別為形狀參數b、尺度參數λ和位置參數x0。

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

再設λ＝1就成了“標準”(一個參數）的威布爾分布概率密度函數。

對于標準化的威布爾分布，b分別為0.5，2，3.5，5時，威布爾分布概率密度函數如圖1所示。

由威布爾分布概率密度函數及圖1可看出其有如下特點：

圖1 標準威布爾分布概率密度函數Fig.1 Standard Weibull distribution probability density function

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

2）威布爾密度函數有3個參數，比正態分布多1個，這是其優點(可描寫安全壽命）但也是其缺點(數學形式比較復雜）的來源。正態分布能很好地描寫對稱狀態的數據，問題是在現實生活中真正對稱的數據只是一種理想狀態，因此正態分布往往只是一階近似。而威布爾分布既能描寫對稱狀態(盡管是近似的），又能描寫左偏或右偏的數據，故其應用范圍要比正態分布大得多。

上文對威布爾分布做了定性的描述，現進一步對其做定量的研究。為方便，開始時僅對標準威布爾分布進行比較詳細的研究，而二參數和三參數的威布爾分布均可得到類似的結果。

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

可設z＝xb，式(6）左邊可變為

對于兩參數和三參數的威布爾分布同樣可得到這個結論，只要分別設：

2）容易證明，威布爾分布函數為

事實上，和1）相同可設z＝xb：

這個分布函數恰恰就是破壞率，而可靠度P為

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

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

這表明任何威布爾分布當其值為特征參數時可靠度一定為36.8%，這是特征參數xa的一個物理意義。而取xp＝x0時P＝100%，也是參數x0為100%可靠度的安全壽命的由來。順便指出，對于正態分布不可能有這個概念，因其左側與水平軸不相交。

3）當b＞1時密度函數是凸函數，因此必有一個極(大）值。事實上：

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

再求f(x）的二階導數：

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

而對于兩參數和三參數威布爾分布概率密度函數的極大值點位置分別為

4）威布爾分布的數學期望按照定義且設z＝xb可得

再注意到伽馬函數的定義：

類似地，對于兩參數和三參數威布爾分布的數學期望分別為

5）威布爾分布的方差。

同樣按照定義可有

由式(13）可得

類似地，對于兩參數和三參數威布爾分布的方差分別為

很明顯，三參數的威布爾分布方差的期望值與x0沒有直接關系。

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

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

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

Python中運行的結果為：k(對分次數）＝19，E(精度）＝1×10－8，b＝3.439 54。

Excel的優點在于直觀，但在計算過程中需要人的介入，且效率較低。Python效率高，代碼一旦調試成功，全部工作都可自動完成，精度也高。

順便指出，這個結論雖然是從標準威布爾分布推出的，但對于兩參數和三參數分布都是適用的，因對稱性僅僅與形狀參數b有關。另外按照文獻［7］，威布爾分布的偏態系數為

不難看出，這個偏態系數也僅僅與b有關。同時，偏態系數為零時也應該是對稱性最好的時候。亦是和正態分布最接近的時候，這對于比較這2種分布的異同意義重大，為此有必要求出這個b值。即要解如下方程：

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

3 疲勞統計學的智能化

20世紀80年代，計算機在中國剛剛開始起步和普及，做一些簡單的統計工作沒有問題，但在工程實踐中應用較少。一方面是因為硬件不行，另外一方面懂得編程的工程技術人員不多。因此，各種現成的圖表，如正態概率坐標紙、威布爾分布概率坐標紙等還是廣為使用。這個方法雖然看起來比較簡單實用，但誤差較大，且必須將數據對數化，這從數學的角度看是一種空間變換，而這樣一變換很可能會將這2種分布變得看起來不可區分。更加重要的是，這種方法難以將擬合程度量化。只能夠看起來擬合不錯，但實際上并非如此。所謂智能化，就是要將得到的數據直接通過計算機來得到人們要求的統計結果，即不僅是由計算機畫出直觀的圖表，而且還有給出判斷的定性與定量的結果?，F舉例說明。

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

表1 一組疲勞壽命數據［1］Table 1 Aset of fatigue life data［1］

表1中：cycle表示循環次數，而均值Nav、標準差s及中值Nm皆由Excel得出。若這些數據是服從正態分布的，則相應的正態分布參數的估計值可由式(20）給出：

式中：“＾”表示估計值。

于是，相應的正態分布密度函數為

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

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

從式(24）得

再由式(23）得

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

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

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

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

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

式中：i為數據(觀測值）由小到大排列的序數；n為數據的個數。

不過仍然可從另外一個角度來看正態分布是否被滿足，即卡方檢驗，但此法相當麻煩，可參考文獻［1］。

例2利用文獻［1］表12-3中的數據繪制表2。

表2 100個試件在同一應力條件下疲勞壽命數據Table 2 Fatigue life data of 100 specimens under the same stress condition

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

驗算：Nav1＝5.315，s＝1.289，Nm1＝5.07，正態分布與理想可靠度的決定系數為0.979 14，威布爾分布與理想可靠度的決定系數為0.988 44。

在此要注意這里用的是決定系數［10］而不是相關系數來區分擬合水平，主要是因為圖3的橫坐標變成了疲勞壽命，這樣其與可靠度的關系不再是線性關系。計算結果表明，威布爾分布的決定系數仍然比正態分布大，盡管不是大很多。最重要的是，這些數據是偏態的，用正態分布已經不能很好描述。而威布爾分布恰恰彌補了這個不足，如圖3所示。不過，也必須指出用這個解析法來求威布爾分布的3個參數并非沒有瑕疵。例如，經過計算，若這些數據滿足威布爾分布，那么安全壽命是3.39，但開始3個數據是小于這個值的，則計算過程沒有問題，很可能是因為這個方法先將形狀參數b找出來后再計算N0和Na，就無法保證N0一定會小于輸入的數據的最小值。故需要對威布爾分布參數算法做進一步的研究。

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

4 高鎮同法

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

其估計值可利用式(31），式(32）在取二次自然對數后可得

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

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

1）輸入原始數據，若原始數據沒有排好序，則先排序。

2）利用Python中的scipy可直接通過以給定精度的間隔來遍歷x0的可能值區間［0，xmin］，以求出使得相關系數最大的x0，即x0max。

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

例3現利用例1的數據在Python代碼上進行試算，可得到如下結果：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，解析法威布爾分布與理想可靠度的相關系數為0.999 20，高鎮同法威布爾分布與理想可靠度的相關系數為0.999 22。

例4高鎮同法與解析法最大不同是b值不同(前者b＝2.21，這里b＝2.040，相差8%），從而標準差s不同(132 vs 146，相差10%）。但相關系數r(0.999 20 vs 0.999 22）只有10－5，在這個意義上兩者并沒有什么大的差別。因解析法對例2的解是不自冾的，對Python的代碼稍加修改可得如下結果：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，正態分布與理想可靠度的決定系數為0.979 11，威布爾分布與理想可靠度的決定系數為0.989 03。

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

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

由此可見，高鎮同法的結果與前面的解析法明顯不同，即不存在超過安全壽命N0的原始數據。幾個參數都較符合，即使是標準差也只有3%的誤差，且高鎮同法得到的相關系數明顯大于解析法。在此意義上，可認為高鎮同法是優于解析法的。

例5值得注意的是，高鎮同法與文獻［11］中提出的“確定威布爾分布三參數相關系數優化法”還是有所不同的，主要體現在：文獻［11］中的數學推導較麻煩，相應的代碼也會較復雜，仍然沒有充分發揮出利用計算機智能的優勢。將文獻［11］中表3的數據放入同樣的Python代碼即得如下結果：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如何使用高鎮同法(參照文獻［11］中的表3數據繪制圖6）。而得到的結果與文獻［11］中的結果幾乎一樣，但要注意在文獻［11］中得到的相關系數是負的，而這里是正的，主要是因為在文獻［11］中設Y＝－ln［ln(1／p）］。

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

由此可見，本文方法相對簡單，容易掌握，不會出現矛盾的結果。這表明高鎮同法具有相當的優越性，值得推廣。

5 高鎮同法的應用

高鎮同法對于擬合三參數疲勞性能曲線也是可以應用的。因為這與求威布爾分布的三參數在數學上是沒有區別的。即用相關系數最大化來定出位置的同時，再由最小二乘法求另2個參數，即可智能化地同時得到3個參數。

例6文獻［7］中給出了一個三參數冪函數表達式：

式中：S0、m和C為材料常數；Smax為最大應力。

文獻［7］雖然也是用最小二乘法，利用相關系數最大化來求出S0，再去求另外2個參數，但較麻煩。現介紹高鎮同法如何用應在例6中。設S0為已知，在式(37）兩邊取10為底的對數：

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

即可得

按照高鎮同法，可將S0根據需要而定的精度遍歷區間［0，Smax］，以求出使得相關系數最大的S0，同時得出相應的b和d，再由式(39）定出m和C。結果與文獻［7］幾乎相同。以文獻［7］第11章例4中的數據為例，繪制表3。

表3 例6的一組數據Table 3 Exam ple 6 aset of data

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

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

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

例7以文獻［7］中第11章例5中的數據為例，繪制表4。

表4 例7的一組數據Table 4 Example 7 aset of data

表4中數據服從如下經驗公式：

式中：a為試件出現裂縫的長度；C、m、a0(也稱為初始裂縫長度）都是與材料有關的常數。

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

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

同樣用高鎮同法，以視需要而定的精度使a0遍歷區間［0，amin］，可求出使得相關系數r最大的a0。再根據式(43）定出C和m。將Python的代碼做出適當的修改可得如下結果：r＝0.992 03，a0＝0.173，m＝0.343，C＝6.989 6×103。

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

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

最后，要對高鎮同法的由來再做一個簡要的說明。①高鎮同法的基礎嚴格來講應該是最優擬合，最小二乘法只是其特例。在文獻［11］中已給出了一些很好的例子。②在文獻［11］的基礎上，不少學者［12-15］都做了不少改進工作，取得了一定的成績，特別是回看文獻［15］離高鎮同法只有一步之遙，很可惜沒有再深入下去。③為了解決解析法出現的不自冾問題，不通過直接求相關系數極大值，而是利用scipy這個軟件庫將所有可能的相關系數求出來之后通過排序將最大值求出來，同時也將相應的威布爾參數求了出來。且從數學角度看此法對于負的相關系數也是同樣適用的，很快就能推廣到三參數疲勞性能曲線的擬合中使用。

6 結論

1）威布爾分布是一種全態分布，比正態分布更具有一般性，可以相信在大數據時代將發揮出更大的作用。特別是應用到疲勞統計學中具有100%可靠度的安全壽命，簡稱為安全壽命這個參數，其涉及可靠性的實際應用領域中具有重大意義。

2）將疲勞統計學智能化的一個重要切入點是將估計威布爾分布三參數智能化。威布爾分布的數學形式雖然比較復雜，但完全可用像高鎮同法(充分利用Python的特點而創造出來的一種新算法）智能化地克服圖解法和解析法的弱點，很方便地解決威布爾分布三參數的計算問題，不僅如此，同時可解決廣義三參數疲勞性能曲線的擬合問題。

3）疲勞統計學智能化不僅是實際應用中必不可少的一個工具，在理論研究中也是功不可沒。例如，如何確定威布爾分布的無偏性和對稱性，若沒有計算機智能的幫助其計算復雜性也是令人卻步的。

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

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

致謝感謝力中國際融資租賃有限公司萬偉浩先生對于本文及有關研究工作的支持。