宮鳳強,王天成,黃天朗
(1.中南大學資源與安全工程學院,湖南長沙,410083;2.東南大學土木工程學院,江蘇南京,211189)
在工程可靠度分析中,工程隨機參數的最優(yōu)概率密度或分布函數模型直接影響可靠度結果的準確性[1-2],因此,對最優(yōu)概率模型推斷方法的研究始終是一項基礎性工作[3-4]。在很多工程可靠度領域,經常遇到大量的極值型參數。所謂極值型參數,主要是考慮一系列最大(小)值參數,如建(構)筑物抗震設計中需要考慮的年最大地震荷載、海洋巖土工程中需要考慮的年最大風荷載等。當采用可靠度方法分析上述參數時,獲取它們的最佳概率分布函數是首要步驟。目前,很多學者通常采用極值型分布對該類參數進行擬合檢驗。例如,張延年等[5]統(tǒng)計了中國159 個代表性城市在1951—2008年的歷年最大風速,選取Gumbel(極值Ⅰ型)分布進行擬合分析;李凱平等[6]利用極值Ⅲ型分布擬合地震烈度;程思軍[7]利用現場試驗數據,對局部腐蝕鋼筋的最大腐蝕深度進行概率統(tǒng)計,揭示Gumbel 分布可作為局部腐蝕鋼筋的最大腐蝕深度概率模型。另外,一些研究者認為其他經典分布(正態(tài)分布、對數正態(tài)分布等)也可以用于擬合極值型參數的概率模型,如莫華美[8]研究最大積雪深度的概率分布時,認為對數正態(tài)分布代替極值Ⅰ型分布更具優(yōu)越性;段忠東等[9]認為極值風速的最優(yōu)分布為威布爾分布。從上述研究可知,極值型參數的最優(yōu)概率分布并不一定是極值型分布,但是均從經典分布范圍中選取。這種擬合方法存在以下2個不易解決的根本問題:1)區(qū)間不匹配。經典分布的定義區(qū)間一般為無限區(qū)間或半無限區(qū)間,但極值型工程參數的分布區(qū)間為有限數值區(qū)間,因此,利用經典分布擬合此類參數時在理論上無法滿足累積概率等于1的要求。2)經典分布無法反映工程參數實際分布的隨機波動性。經典分布均為單峰值型分布,但極值型工程參數的分布可能存在多個峰值,呈現一定波動性,因此,發(fā)展能夠滿足有限區(qū)間上累積概率等于1并能夠體現實際分布隨機波動性的概率分布推斷方法十分必要。在這方面,很多專家提出了很多種推斷方法。如蘇永華等[10]提出了利用一般多項式推斷巖土參數概率分布的方法;LI等[4]提出了正交多項式方法;張道兵等[11]運用最大熵原理估計了隧道襯砌結構參數的密度函數;李松輝[12]引入了車輛荷載效應的截尾分布模型。此外,HUANG等[13]在自然災害風險分析中引入了信息擴散理論,目前該方法在巖土工程中得到了廣泛應用。宮鳳強等[14]利用正態(tài)信息擴散方法推斷了小樣本巖土參數的概率分布,取得了較好的擬合效果。周道成等[15]用正態(tài)信息擴散方法確定了河冰抗壓強度概率分布,證明了正態(tài)信息擴散分布更加接近河冰抗壓強度的真實分布,并優(yōu)于經典分布的擬合方法;徐志軍等[16]基于正態(tài)信息擴散理論,擬合了路基沉降的概率分布,并通過工程實例驗證了此方法的正確性和有效性;黃達等[17-18]將正態(tài)信息擴散原理廣泛用于推斷巖體抗剪強度參數、巖體質量參數和粗粒土參數的單變量邊緣分布,取得了良好效果;宮鳳強等[19]利用正態(tài)信息擴散方法研究了Mohr-Coulomb 強度準則參數概率模型推斷方法,其概率分布更加接近參數的實際分布規(guī)律。根據上述研究,本文作者提出基于正態(tài)信息擴散原理的極值型工程參數概率分布的推斷方法,并考察該方法是否滿足上述2個根本條件。在研究過程中,以樣本個數分別為25和61的年最大標準風壓樣本和內摩擦角樣本為例,分別采用正態(tài)信息擴散分布和Gumbel 分布進行擬合,對比分析2種分布擬合所得結果的優(yōu)良性。在此基礎上,利用模擬樣本,考察樣本個數對擬合精度的影響,并將正態(tài)信息擴散分布應用于2組核管道最大腐蝕深度的預測。
極值理論研究的主要問題是極大(小)值的極限分布問題,是統(tǒng)計數學中的重要分支。設X1,X2,…,Xn為獨立同分布的隨機變量,分布函數為F(x)(稱為底分布),對自然數n,令Mn=max{X1,X2,…,Xn},mn=min{X1,X2,…,Xn},分別表示n個隨機變量的最大值與最小值,則[20]
若已知分布函數F(x),則可以根據式(1)和(2)精確求出最大值和最小值的分布函數。但在實際應用中,隨機變量的分布函數F(x)往往是未知的,因此,很難直接用于統(tǒng)計分析。為求得原始樣本的分布函數,必須考慮n→∞時分布函數的情況即極值漸進分布,簡稱極值型分布。
極值統(tǒng)計理論已經證明對于3種類型的連續(xù)型隨機變量原始分布,根據其分布尾部的不同形式,當n→∞時,其極限分布只有3種類型,包括極值Ⅰ型(Gumbel)分布、極值Ⅱ型(Frechet)分布、極值Ⅲ型(Weibull)分布。其中,Gumbel 分布的概率密度函數和概率分布函數表達式見式(3)和(4),其他2種分布的概率函數表達式見文獻[20]。
設隨機變量O的概率密度函數為f(x),定義μ(x)在(-∞,+ ∞)上的1個波雷爾可測函數,Δn為常數且大于0,則稱[13]
為總體密度f(x)的1 個擴散估計;μ(x)為擴散函數;Δn為總體密度f(x)的窗寬;n為隨機變量的樣本數;xi為樣本觀測值。根據信息擴散過程,μ(x)由下式確定:
式(6)滿足式(5)對擴散函數的要求,且與概率論中正態(tài)分布的密度函數形式一致,稱為正態(tài)信息擴散函數。將式(5)和式(6)聯合求解,則隨機變量O的概率密度函數f(x)的正態(tài)信息擴散估計為
Gumbel 分布函數的定義區(qū)間一般為(-∞,+ ∞),而樣本工程物理力學參數的取值通常大于0 且其范圍為有限區(qū)間。當用Gumbel 分布函數擬合這些數據時,必然存在區(qū)間不匹配問題。為了解決這個問題,基于“3σ”統(tǒng)計原則并考慮Gumbel 分布的不對稱性,確定采用考慮偏度c的“3σ”分布區(qū)間截尾方法,具體的區(qū)間截尾公式見表1。
表1 截尾區(qū)間確定方法Table 1 Truncated interval determination method
選擇高大釗[22]列舉的年最大標準風壓和張蕾等[23]列舉的內摩擦角這2 種工程參數為例進行研究,具體數據見表2。
為了后續(xù)擬合檢驗計算的需要,對實際工程樣本參數(如最值、平均值、標準差及偏度)進行統(tǒng)計,并計算相應的截尾區(qū)間左、右端點值,具體結果見表3。
為了對比正態(tài)信息擴散分布與Gumbel 分布的擬合效果,從概率分布曲線、K-S檢驗值和累積概率3個方面進行具體分析。
圖1(a)和圖1(b)所示分別為2 組實際工程樣本的累積概率分布曲線。從圖1(a)可見:NID累積概率分布曲線與Gumbel 累積概率分布曲線在前半段幾乎重合;在后半段,NID累積概率分布曲線完全沿著經驗分布折線的變化趨勢延伸,而Gumbel 累積概率分布曲線與經驗分布折線的偏差較大。圖2(a)和圖2(b)所示分別為2組實際工程樣本的概率密度函數曲線。從圖2(a)可見,年最大標準風壓樣本的直方圖反映了實際分布具有很大的波動性,而Gumbel 分布曲線呈單峰值型分布,顯然無法很好地刻畫實際分布的波動特性;NID概率密度分布曲線具有多峰值特點,更加貼近直方圖的分布趨勢。同理,對圖2(b)的分析結果類似。綜上分析可知NID概率密度分布比Gumbel分布更優(yōu)。
表4 所示為2 組實際工程樣本的K-S 檢驗值計算結果。根據K-S檢驗方法的特點,對2種工程參數的分布檢驗均在顯著性水平為0.05(置信水平為95%)時進行。從表4 可見:在同一置信水平下,2組樣本的K-S 檢驗法臨界值分別為0.264 0 和0.174 0;NID 分布檢驗值的計算結果分別為0.127 2和0.043 8,而Gumbel分布檢驗值的計算結果分別為0.132 4 和0.054 5。顯然,這2 組樣本中NID 分布和Gumbel 分布的檢驗結果均通過臨界值檢驗,但Gumbel 分布的檢驗值均比NID 分布的檢驗值高,這表明NID分布具有更優(yōu)的擬合效果。
表2 實際工程樣本數據Table 2 Sample data of actual engineering
表3 樣本數據統(tǒng)計信息Table 3 Statistical information of sample data
圖1 實際工程樣本的概率分布函數曲線比較Fig.1 Comparison of probability distribution function curves for actual engineering samples
圖2 實際工程樣本的概率密度函數曲線比較Fig.2 Comparison of probability density function curves for actual engineering samples
表4 實際工程樣本K-S檢驗值和累積概率計算結果Table 4 Calculation results of K-S test values and cumulative probability for actual engineering samples
從表4還可見:根據確定后的截尾區(qū)間計算整個區(qū)間上的累積概率,對于Gumbel 分布,2 個實例樣本所得到的累積概率分別為0.997 2和0.995 7,而NID 分布的累積概率均為1.000 0。上述結果表明在截尾區(qū)間下,Gumbel 分布無法滿足累積概率等于1的先決條件。
在實際工程中,所獲得的樣本很多都是小樣本,很難考察不同樣本個數對上述2種擬合方法在擬合精度上的影響。為了研究樣本個數對2種分布擬合精度的影響,基于Monte-Carlo 方法抽取不同個數的樣本,分別利用上述2種分布進行擬合,并根據K-S 檢驗值與累積概率2 種指標,評判2 種分布的擬合效果。
以年最大標準風壓樣本為例,并以α=0.001 5和β=17 602.505 8(即9 602.505 8+8 000)的Gumbel分布函數作為母函數,利用Monte-Carlo 方法生成8 組樣本,個數分別為 15,20,30,50,100,200,500 和1 000。表5 所示為模擬樣本數據的統(tǒng)計結果(由于生成的模擬數據較多,文中沒有給出具體的模擬樣本數據)。
以樣本個數分別為15,50,200 和1 000 的4組模擬樣本為例進行說明,繪制概率分布函數曲線,見圖3。從圖3 可以看出:不論樣本個數較少還是較多,相對于Gumbel 分布曲線,NID 分布曲線和經驗分布曲線都更加接近。圖3所示概率概率分布曲線還表明:不論樣本數據的實際分布完全符合經典分布還是有一定的波動性,正態(tài)信息擴散方法均可以很好地刻畫樣本的實際分布。
參考4.2節(jié)的檢驗對比過程,得出模擬數據的K-S 檢驗值和累積概率的計算結果,見表6。從表6 可見:對于每一組模擬樣本而言,NID 分布和Gumbel 分布均通過了置信水平為95%的臨界值檢驗,但Gumbel 分布的檢驗結果均大于NID 分布的檢驗結果,這說明在不同的樣本個數下,NID分布的擬合優(yōu)良性更佳。此外,在截尾區(qū)間下,Gumbel 分布的累積概率始終小于1.000 0,但NID分布的累積概率恒等于1.000 0,并不受樣本個數變化的影響。
根據表6 繪制模擬樣本個數遞增下Gumbel 分布與NID 分布的K-S 檢驗值變化曲線,如圖4 所示。從圖4 可見:在樣本個數遞增時,臨界值和NID 分布的檢驗值均逐漸遞減并趨于收斂,而Gumbel 分布的檢驗值在n=20 和n=500 的點附近存在先增后減的情況,即存在一定波動性。
同樣,根據表6中的數據繪制模擬樣本個數遞增下Gumbel 分布與NID 分布的累積概率變化曲線,如圖5 所示。由圖5 可知:在樣本個數遞增時,Gumbel 分布的累積概率波動范圍較大且毫無規(guī)律性,而NID 分布的累積概率始終為1.000 0,與樣本個數的變化無關。
為了進一步說明正態(tài)信息擴散分布推斷方法的可行性與有效性,以2個核管道腐蝕深度的預測問題作為工程實例進行對比分析[24-25]。
1)某被腐蝕核管道在一定時間內的平面觀測深度實測值分別為2.82,2.96,3.08,3.09,3.18,3.19,3.22,3.32,3.33,3.52,3.58,3.61,3.62,3.95 和4.12 mm。周國強等[26]認為其分布符合極值Ⅱ型分布并預測該核管道最大腐蝕深度為4.840 3 mm。用正態(tài)信息擴散分布擬合該組數據并計算相應的K-S 檢驗值,結果顯示:極值Ⅱ型分布的檢驗結果(即以核管道最大腐蝕深度為變量的檢測結果,量綱一參數)為0.115 0,NID分布的檢驗結果為0.106 4。這2 種分布的檢驗值均小于臨界值0.338 0,顯然,上述分布均通過了臨界值檢驗。根據張博庭[27]提出的有限比較法可知,極值Ⅱ型分布的有限比較結果大于NID 分布的有限比較結果,這說明NID分布更適合作為此核管道腐蝕深度的概率分布。以NID 分布作為最大腐蝕深度的預測模型,計算得到該核管道最大腐蝕深度為4.284 9 mm。由此可知:選擇擬合效果更好的NID分布作為預測模型,其最大腐蝕深度的預測結果會更加精確。
表5 年最大標準風壓的模擬數據統(tǒng)計信息Table 5 Statistical information for simulated data of annual maximum standard wind pressure
圖3 模擬樣本的概率分布函數曲線隨樣本個數增加的比較Fig.3 Comparison of probability distribution function curves for simulated samples with the increase of sample numbers
圖4 模擬樣本的K-S檢驗值隨樣本個數增加的比較Fig.4 Comparison of K-S test values for simulated samples with the increase of sample numbers
圖5 模擬樣本的累積概率隨樣本個數增加的比較Fig.5 Comparison of cumulative probability values for simulated samples with increase of sample numbers
為了更加清晰地看出NID 分布的擬合效果,繪制NID 分布和極值Ⅱ型分布的擬合曲線,如圖6所示。由圖6 可知:極值Ⅱ型分布的累積概率未達到1,而NID 分布的累積概率為1,這說明NID 分布函數能夠較好地擬合核管道腐蝕深度的數據點分布。圖7 所示為NID 分布與極值Ⅱ型分布的累積概率偏差的變化情況。從圖7可明顯看出:NID分布的最大偏差小于極值Ⅱ型分布的最大偏差,避免了局部偏差過大的情況,這表明NID 分布能夠更好地擬合真實分布;另外,核管道腐蝕深度模型的優(yōu)化能夠更好地提高最大腐蝕深度的預測精度。
圖6 實例1的NID分布和極值Ⅱ型分布擬合曲線Fig.6 Fitting curves of NID and extremum Ⅱ-type distributions for example 1
圖7 實例1的擬合誤差折線Fig.7 Fitting error poly-lines for example 1
2)某核電站設冷水系統(tǒng)除淤管道不銹鋼部分運行1 a后的超聲測量厚度,在同一條件下,其實測值分別為 2.9,3.7,3.9,3.9,4.6,4.7,4.8,5.0,5.4,5.5,5.7,5.9,6.0,6.3,6.3,7.0 和8.1 mm。周國強等[26]認為其分布符合極值Ⅲ型分布并預測該核管道最大腐蝕深度為9.3 mm。按照實例1 的計算過程,進行K-S 檢驗值計算,結果顯示:極值Ⅲ型分布的管道不銹鋼部分腐蝕深度檢驗結果為0.120 1,NID分布的檢驗結果為0.085 9。這2 種分布的檢驗值均小于臨界值0.318 0。同理,這2種分布均通過了K-S檢驗,但極值Ⅲ型分布的有限比較結果大于NID 分布的比較結果。顯然,NID分布比極值Ⅲ型分布更適合作為此核管道腐蝕深度的概率分布。利用NID 分布預測該核管道最大腐蝕深度為8.7 mm,表明由極值Ⅲ型分布預測的最大腐蝕深度比NID分布的大。與實例1進行比較可得出同樣的結論,即選擇擬合效果更好的NID分布作為預測模型,其最大腐蝕深度的預測結果會更加精確。
同理,繪制該樣本的NID 分布和極值Ⅲ型分布的擬合曲線圖,如圖8 所示(其中,極值變量指核管道最大腐蝕深度)。由圖8 可知:核管道腐蝕深度樣本數據點大致分布在NID分布曲線的兩側,表明NID分布有較好的擬合效果。圖9所示為該樣本的擬合誤差折線圖。從圖9可知:極值Ⅲ型分布的最大偏差大于NID 分布的最大偏差,說明NID分布更加接近真實分布。
圖8 實例2的NID分布和極值Ⅲ型分布擬合曲線Fig.8 Fitting curves of NID and extremum Ⅲ-type distributions for example 2
圖9 實例2的擬合誤差折線Fig.9 Fitting error poly-lines for example 2
1)利用正態(tài)信息擴散分布和極值型分布分別對年最大標準風壓和內摩擦角2組工程樣本進行擬合,正態(tài)信息擴散分布的K-S 檢驗值比極值Ⅰ型(Gumbel)分布的K-S檢驗值??;正態(tài)信息擴散分布的累積概率等于1,而Gumbel 分布的累積概率小于1。正態(tài)信息擴散分布的曲線能夠更好地反映真實樣本的直方圖或經驗分布折線的波動變化。
2)在不同樣本個數下,正態(tài)信息擴散分布與Gumbel 分布均通過了K-S 檢驗,但正態(tài)信息擴散分布的檢驗值均低于Gumbel 分布的檢驗值。隨著樣本個數增加,正態(tài)信息擴散分布的K-S檢驗結果具有更快的收斂速度和更好的收斂穩(wěn)定性;正態(tài)信息擴散分布的累積概率始終等于1,并不受樣本個數的影響,而Gumbel 分布的累積概率有較大的波動范圍并恒小于1。
3)正態(tài)信息擴散分布預測模型比極值Ⅱ型分布或極值Ⅲ型分布的模型精度更高;由正態(tài)信息擴散分布預測的最大腐蝕深度比極值Ⅱ型分布或極值Ⅲ型分布的預測結果更加精確。
4)正態(tài)信息擴散方法在推斷極值型工程參數的概率密度函數方面具有分布參數唯一、計算過程簡便、累積概率恒定、擬合檢驗值低、對樣本個數適應性強的特點。