正態(tài)信息擴散法推斷巖土參數概率分布的窗寬優(yōu)化和區(qū)間截尾問題研究

孔令奇，李翠娟

(1.成都工業(yè)學院 材料與環(huán)境工程學院，四川 成都 610031；2.西南交通大學 土木工程學院，四川 成都 610031)

巖土參數概率分布的研究一直是一項基礎性研究工作，具有極其重要的價值和意義[1?3]。國內外研究人員已做了大量工作[4?7]，提出了多種巖土參數概率分布擬合方法，主要包括參數估計[8?9]和非參數估計2種[10?12]。基于非參數估計的主要有最大熵法、正交多項式等方法[1?3,10?12]。宮鳳強等[13]系統(tǒng)考察了經典分布擬合、最大熵法、一般多項式法、正交多項式法和正態(tài)信息擴散法在推斷巖土參數概率分布的優(yōu)劣，最終確定正態(tài)信息擴散法為相對最優(yōu)的推斷方法。宮鳳強等[10]提出正態(tài)信息擴散法推斷小樣本巖土力學參數概率密度函數的有效性。朱喚珍等[1]進一步證明正態(tài)信息擴散法在大樣本巖土參數概率分布推斷中的有效性，但該方法在應用的過程中存在窗寬選擇困難和截尾誤差的問題[14]。朱喚珍等[1]給出了最優(yōu)窗寬的計算方法，但忽略了窗寬選擇與頻率直方圖的匹配關系，且對于巖土參數測試樣本來說，由于受經濟技術條件限制，巖土體參數的試驗數據比較有限[15]，難于滿足最優(yōu)窗寬的取值條件。宮鳳強等[16]給出了正態(tài)信息擴散法截尾分布的處理方法，正態(tài)信息擴散法采用Gaussian作為核函數時，正態(tài)分布屬于無界分布，擬合曲線的首端和尾端向兩邊延展，給出的c3取值方法是針對樣本區(qū)間的截尾方法，但其擴大了正態(tài)信息擴散分布的擬合邊界。本文在既有研究基礎上提出正態(tài)信息擴散法窗寬確定方法，在不同直方圖分組基礎上，獲得測試樣本“匹配窗寬”；基于正態(tài)信息擴散分布首尾延展特性，結合巖土參數實測數據的分布特征，提出“收窄”正態(tài)信息擴散分布取值區(qū)間的方法；通過工程實際測試樣本驗證本文方法的有效性和合理性，完善正態(tài)信息擴散法在巖土參數概率密度函數擬合上的應用。

1 正態(tài)信息擴散法窗寬確定

1.1 正態(tài)信息擴散法的概率密度函數

為了避免區(qū)間劃分對概率密度函數曲線形狀的影響，正態(tài)信息擴散法考慮以任意一點x為中心的鄰域，統(tǒng)計落入這個區(qū)間的樣本點個數，區(qū)域的位置隨著x的變動而變動，任意一點x處的概率密度函數值由區(qū)域內的樣本點數決定。

設X1,X2,…,Xn是取自總體X的樣本，則樣本的概率密度函數為[1,17-19]

其中，K()為核函數；h為窗寬，即每個區(qū)間的長度；n為樣本個數。

1.2 既有窗寬確定方法

窗寬是正態(tài)信息擴散法推斷概率密度函數中一個重要的指標。文獻[1]和[17]提出通過求解積分均方誤差MISE(mean integrated squared error)函數的最小值獲得最優(yōu)窗寬，即

其中：f(x)為樣本總體的實際概率密度。

當K(x)滿足如下條件時：

1)K(x)在區(qū)間[?1,1]取值，且要求對稱；

2)∫K(x)dx=1；

3)∫xK(x)dx=0；

4)∫x2K(x)dx=σ2k>0，σk為標準差。

則

但獲得式(4)還有一個最關鍵的前提條件，即滿足h→0，nh→+∞時[17]，式(4)才近似相等。對于巖土參數測試樣本來說，由于受經濟技術條件限制，巖土體參數的試驗數據比較有限[15]，難于滿足此條件。所以對巖土參數測試樣本來說，采用正態(tài)信息擴散法擬合其概率密度函數時，通過式(4)計算出的h值并非最佳窗寬。

在此僅以文獻[20]提供的內摩擦角φ的測試樣本為例(樣本81組，具體樣本值見文獻[20])，詳細說明窗寬取值對擬合精度的影響，更多測試樣本的分析結果見1.3.2節(jié)。當直方圖分組數取為10時，采用不同窗寬擬合內摩擦角φ測試樣本的概率密度函數，擬合曲線如圖1所示。

圖1 內摩擦角樣本不同窗寬下的擬合曲線Fig.1 Fitting curves of cohesion samples under different window widths

K-S檢驗法的檢驗量是整個取值范圍內的最大偏差值，在一定的顯著性水平下，可評估某一分布擬合的有效性，但不能提供某一分布擬合是良好的絕對信息[21]。均方差可從全局角度評估正態(tài)信息擴散法對巖土參數概率密度函數的擬合質量，故本文采用均方差(Mean Squared Error，MSE)指標檢驗擬合曲線的優(yōu)劣[22]

式中：fi為理論頻率；S*(xi)為頻率直方圖中每個子區(qū)間的相對頻率；m為區(qū)間個數。

內摩擦角樣本在不同窗寬下采用正態(tài)信息擴散法推斷其概率密度函數時的均方差指標統(tǒng)計值如表1所示。

表1 內摩擦角樣本的檢驗結果Table 1 Test results of internal friction angle samples

由圖1和表1明顯可見，當采用式(4)計算的窗寬擬合內摩擦角樣本的概率密度函數時，擬合誤差最大，均方差指標為0.084 3，隨著窗寬取值的減小，擬合誤差在h/3時達到最小，均方差指標為0.035 7，繼續(xù)減小窗寬，擬合誤差增大。所以，對于內摩擦角測試樣本，其最佳的匹配窗寬是h/3，h并非最佳窗寬。

1.3 匹配窗寬確定

1.3.1 直方圖分組

理論和實測數據均已顯示選擇窗寬h做為正態(tài)信息擴散法的最優(yōu)窗寬是不恰當的，對既有研究工作總結發(fā)現(xiàn)，目前對巖土參數概率密度函數曲線的擬合不論是基于參數估計還是非參數估計法，均是根據測試樣本獲得直方圖，而直方圖繪制依賴于區(qū)間的劃分，區(qū)間劃分的不同直接影響曲線的形狀和擬合誤差的大小。在文獻[1]中，作者給出了最佳窗寬的計算公式，但忽略了窗寬選擇與直方圖的匹配關系，最優(yōu)窗寬是相對的，所以本文稱其為“匹配窗寬”。

仍以內摩擦角的測試樣本為例，當直方圖分組數取為5和20時，其正態(tài)信息擴散分布的擬合曲線如圖2～3所示。直方圖分組為5時檢驗結果統(tǒng)計值如表2所示。直方圖分組為20時檢驗結果統(tǒng)計值如表3所示。

表2 分組為5時檢驗結果Table 2 Test result when grouping is 5

表3 分組為20時檢驗結果Table 3 Test result when grouping is 20

圖2 分組為5時擬合曲線對比Fig.2 Comparison of fitting curves when grouping is 5

圖3 分組為20時擬合曲線對比Fig.3 Comparison of fitting curves when grouping is 20

由圖1～3明顯可見，當直方圖分組數取為5時，最佳匹配窗寬為h/2；當直方圖分組數取為10時，最佳匹配窗寬為h/3；當直方圖分組數取為20時，最佳匹配窗寬為h/5。測試樣本直方圖分組數越大，應采用較小的匹配窗寬；直方圖分組數越小，應采用較大的匹配窗寬。所以，采用正態(tài)信息擴散法推斷巖土參數概率分布的前提是確定測試樣本直方圖的分組數。

文獻[23]指出，確定直方圖分組數的原則是分組的結果能正確反映數據的分布規(guī)律，組數應根據數據多少來確定。組數過少，會掩蓋數據的分布規(guī)律；組數過多，使數據過于零星分散，也無法顯示樣本分布狀況。

在概率論與數理統(tǒng)計中，對于正態(tài)分布總體的隨機變量，其直方圖子區(qū)間的劃分與樣本數量有最佳關系，取分組數m=1.87(n-1)2/5[21]。但大量的研究工作已表明，巖土參數測試樣本離散性嚴重，呈偏態(tài)分布。所以，對于工程中非正態(tài)分布的隨機變量總體，文獻[24]指出，當樣本個數大于50時，可將直方圖繪制時的分組數m取為

其中，n為樣本個數。

對于巖土參數測試樣本來說，由于受經濟技術條件限制，巖土體參數的試驗數據比較有限[15]，樣本個數通常不足50個。文獻[23]給出當樣本總數在50～100時，直方圖分組數取6～10。文獻[25]給出當樣本數在100以內時，直方圖分組數一般分5～12組。

綜上所述，本文按如下分組：

當樣本個數大于50時，按式(6)進行直方圖分組；當樣本個數小于50時，直方圖分組數取為5或6(依具體樣本選擇)。

1.3.2 匹配窗寬確定方法

由上文的分析可知，直方圖分組數的不同，反映了樣本不同程度的波動性和離散性，與之匹配的窗寬值也應是不同的。

根據式(4)計算的窗寬h對波動性較小的巖土參數測試樣本適用，但對具有一定波動性的巖土參數測試樣本，顯然較小的窗寬才合適。

為此，本文提出巖土參數測試樣本窗寬選擇方法為：取h為初始值，依次改變窗寬的大小為h/2，h/3，h/4…，誤差評價指標均方差最小為目標，確定最佳的匹配窗寬h′。

以文獻[20]粉質黏土壓縮指數(55組樣本)、黏聚力c(81組樣本)2組大樣本和文獻[26?28]提供的5組小樣本為例，驗證本文方法的有效性和準確性。5組小樣本巖土參數分別為：內摩擦角正切值(24組樣本)[26]、液限(26組樣本)[26]、黏聚力(21組樣本)[26]、350號混凝土斷裂韌度(35組樣本)[27]和標準風壓(25組樣本)參數[28]。擬合曲線如圖4～10所示，MSE誤差指標如表4所示。

表4 MSE誤差評價指標Table 4 MSE error evaluation index

圖4 粉質黏土壓縮指數樣本擬合曲線比較Fig.4 Comparison of fitting curves of silty clay compression index samples

圖5 黏聚力c樣本擬合曲線比較Fig.5 Comparison of fitting curves of cohesion c samples

圖6 內摩擦角正切值樣本擬合曲線比較Fig.6 Comparison of sample fitting curves of tangent value of internal friction angle

圖7 液限樣本擬合曲線比較Fig.7 Comparison of fitting curves of liquid limit samples

圖8 黏聚力樣本擬合曲線比較Fig.8 Comparison of fitting curves of cohesive force samples

圖9 350號混凝土斷裂韌度樣本擬合曲線比較Fig.9 Comparison of fitting curves of No350 concrete fracture toughness samples

圖10 標準風壓樣本擬合曲線比較Fig.10 Comparison of fitting curves of standard wind pressure samples

由圖4～10和表4明顯可見：無論是通過擬合曲線的直觀對比，還是通過均方差指標的定量檢驗，都說明本文提出的匹配窗寬下的信息擴散分布所得到的概率密度函數優(yōu)于h窗寬下的概率密度函數。說明了本文方法在正態(tài)信息擴散法窗寬選擇中的正確性和適用性。正態(tài)信息擴散法推斷巖土參數測試樣本概率分布的匹配窗寬與測試樣本的個數和樣本的離散程度相關，每組測試樣本的匹配窗寬不同，應單獨計算。

正如文獻[1]所述，正態(tài)信息擴散法能較為理想的推斷巖土參數測試樣本的概率分布，但存在一個不容忽視的問題，由圖1～10也可見，那就是擬合曲線存在上、下截尾誤差。

2 擬合邊界確定

不論對哪一種巖土參數概率密度函數推斷方法，通常尾部樣本點分散，取值概率很小，擬合過程中，如果上界選擇過大，概率密度擬合曲線的尾部就會出現(xiàn)波動；如果上界選擇過小，擬合曲線的尾部會上翹，因為在這種界限附近，樣本點尚未明顯減少，使得這個點以外的點上的取值概率得以累加。如果樣本下界選得過大，擬合曲線在首端的取值概率較大，會使曲線變得平坦。

2.1 既有擬合邊界確定方法

文獻[16]提出了幾種正態(tài)信息擴散法區(qū)間取值情況：

1)最大、最小值作為區(qū)間邊界。如圖1～10可見，如果選擇最大、最小值作為正態(tài)信息擴散分布法的擬合邊界，顯然出現(xiàn)截尾誤差，累積概率分布小于1。

2)文中對比分析了3σ，4σ，c33和c3型區(qū)間取值方法的優(yōu)劣，最終確定c3型為最合理的截尾區(qū)間。

c3型區(qū)間取值方法為：以[μ-3σ,μ+3σ]為基礎，參考偏度系數c進行調整，當c<0，左端邊界取μ-(3-c)σ，減小下限值；當c>0，右端邊界取μ+(3+c)σ，增大上限值。

正態(tài)信息擴散分布法推斷巖土參數概率分布時首尾延展，超出樣本邊界，出現(xiàn)截尾誤差，繼續(xù)擴大樣本取值邊界，會使截尾誤差增大。以文獻[20]的黏聚力c測試樣本(81組)為例，樣本最大、最小值為[7,43]，測試樣本的偏度系數為0.367大于0，分布右偏，采用c3型截尾區(qū)間為[?3.98,49.3]，顯然正態(tài)信息擴散分布的擬合區(qū)間范圍超出樣本的最大、最小值區(qū)間，且下限為負，這與巖土參數的實際分布不符。

2.2 正態(tài)信息擴散法擬合邊界取值

在一般概率密度擬合方法中，通常取樣本隨機變量z∈[zmin,zmax]，下界值取稍小于zmin，上界值取稍大于。但該方法不能直接應用于正態(tài)信息擴散估計法。比如在確定最大熵法及最佳平方逼近法積分邊界時，上、下界限均是在樣本最大、最小值基礎上增大上界或減小下界，這樣取值是將擬合曲線向首尾兩端擴展，包絡樣本邊界。但正態(tài)信息擴散估計法不同，采用Gaussian作為核函數時，正態(tài)分布屬于無界分布，擬合曲線的首端和尾端向兩邊無限延展，遠超出樣本邊界。所以，正態(tài)信息擴散估計法中首尾邊界的確定方法應與最大熵法等相反，收窄首尾邊界，在樣本最大、最小值基礎上，減小上界、增大下界。

為了確定正態(tài)信息擴散法的擬合邊界，關鍵是如何減小上界、增大下界。目前沒有經驗值、沒有統(tǒng)一的標準，更無文獻做過系統(tǒng)的研究。

正如文獻[16]所述，首尾區(qū)間的選擇應以實測數據的分布范圍為依據，擬合曲線包含樣本區(qū)間范圍，不能取得太寬；確定分布區(qū)間取值范圍的標準是使該區(qū)間內累積概率值接近1。最簡單、直接的方法是：確定樣本正態(tài)信息擴散分布匹配窗寬的基礎上，試探減小上界值、增大下界值，直到正態(tài)信息擴散分布擬合區(qū)間下界值稍大于zmin，上界值稍小于zmax。

為了節(jié)約試探取值時間，本文提出借助樣本隊列曲線法。以文獻[20]給出的黏聚力c樣本(81組)為例，繪制樣本的隊列曲線如圖11所示。

圖11 黏聚力c樣本隊列曲線Fig.11 Cohesion c sample queue curve

由圖11可見，黏聚力c首、尾部樣本點分散，取值概率很小，正因如此，才造成圖12所示的概率密度擬合曲線尾部出現(xiàn)波動，擬合曲線首尾與樣本實際概率分布不吻合。為了實現(xiàn)擬合曲線在樣本最大、最小值處的概率值為0(或趨于0)，應減小上界，增大下界。根據圖11，對于黏聚力c測試樣本，樣本原取值邊界為[7,43]，將該測試樣本從小到大排列，排列之后可見樣本主要取值在第5～79個樣本之間，舍去樣本最大值，上界取第80個的樣本值38，下界舍去前4個離散的小樣本，取第5個樣本值9為下界。改進樣本邊界(即減小樣本上界，增大樣本下界)后的擬合曲線與原樣本邊界(即取樣本最大、最小值作為區(qū)間邊界)擬合曲線對比如圖12所示。

圖12 黏聚力c樣本擬合曲線對比Fig.12 Comparison of cohesion c sample fitting boundary

原樣本邊界下的正態(tài)信息擴散分布的擬合曲線邊界為[1.47,48.53]，遠超出樣本邊界；改進樣本邊界后正態(tài)信息擴散分布的擬合曲線邊界為[3.47,43.53]，與樣本值邊界更接近。原樣本邊界下正態(tài)信息擴散分布的累積概率值為0.997 1，K-S檢驗值為0.048 3；改進樣本邊界后累積概率值為0.999 9，接近1，K-S檢驗值為0.050 9，小于臨界值0.1511。

上文已說明，文獻[16]提出的c3型區(qū)間確定方法下，擬合區(qū)間與實際樣本區(qū)間不匹配，本文不再與之對比。為了進一步證明本文方法的有效性和普適性，對上文的內摩擦角和5組小樣本為例，繪制樣本擬合曲線對比圖，如圖13～18所示。區(qū)間取值和誤差檢驗值如表5所示。

由圖13～18及表5可見：收窄樣本區(qū)間后，正態(tài)信息擴散分布不僅包絡樣本區(qū)間，且與樣本的實際分布區(qū)間更接近；收窄樣本區(qū)間后在滿足擬合誤差的前提下，截尾誤差變小，樣本累積概率值更逼近1。增大的誤差來自擬合曲線的首、尾，實測樣本在首、尾樣本點分散，在繪制頻率直方圖時應在首尾處不等份分組，使樣本首尾處的概率值逐漸減小趨于0。

表5 樣本邊界、累積概率及誤差檢驗統(tǒng)計值Table 5 Sample boundary,cumulative probability and error test statistics

圖13 內摩擦角樣本擬合曲線對比Fig.13 Comparison of boundary fitting of internal friction angle samples

圖14 內摩擦角正切值樣本擬合曲線對比Fig.14 Comparison of fitting boundary of tangent value samples of internal friction angle

圖15 液限樣本擬合曲線對比Fig.15 Comparison of fitting boundary of liquid limit sample

圖16 黏聚力樣本擬合曲線對比Fig.16 Comparison of cohesive sample fitting boundary

圖17 350號混凝土斷裂韌度樣本擬合曲線對比Fig.17 Comparison of fitting boundary of No.350 concrete fracture toughness sample

圖18 標準風壓樣本擬合曲線對比Fig.18 Comparison of fitting boundary of standard wind pressure samples

3 結論

1)合適的樣本直方圖分組是正態(tài)信息擴散法推斷巖土參數概率密度函數的前提，測試樣本直方圖分組數越大，應采用較小的匹配窗寬；直方圖分組數越小，應采用較大的匹配窗寬。在直方圖分組數確定的情況下，以均方差最小為目標，逐步減小窗寬h的“匹配窗寬”確定方法簡單有效。

2)提出的減小上界、增大下界的正態(tài)信息擴散分布樣本邊界確定方法，可有效減小正態(tài)信息擴散法推斷巖土參數概率分布時出現(xiàn)的截尾誤差問題。依據樣本隊列曲線，觀測首尾樣本點分布，可直觀判定樣本上界減小量和下界增大量。改進樣本邊界后的正態(tài)信息擴散分布區(qū)間內的累積概率更接近1。

正態(tài)信息擴散法的窗寬選擇和截尾誤差是研究人員無法回避的問題，本文的研究工作使窗寬選擇和減小截尾誤差問題又向前推進了一步，下一步將研究通過對大量巖土參數測試樣本統(tǒng)計分析，探尋正態(tài)信息擴散法推斷巖土參數概率分布中截尾誤差取值的理論公式。