基于最優直方圖的巖土參數概率模型的最大熵法研究

孔令奇

(成都工業學院材料與環境工程學院，四川 成都 610031)

在巖土工程的風險分析和可靠度分析中，巖土參數分布類型的不同將直接影響可靠度指標的計算結果，所以巖土參數概率特性的研究始終是一項基礎性的重要工作[1-4]。研究人員做了大量的工作，從早期的基于參數估計的正態和對數正態分布的提出[5-7]，到發現beta 分布擬合可以克服正態、對數正態分布無界分布的不足，將巖土參數基于參數估計的概率分布模型最終鎖定在beta 分布[8]。但隨著研究的深入，研究人員發現參數估計法擬合巖土參數概率分布時存在較大誤差[9]，故而提出非參數估計擬合方法，如正交多項式[10-11]、正態信息擴散[12-13]和最大熵法[14-15]等。

參考文獻[14]指出最大熵法避免了一些非參數估計法過分依賴樣本數據的缺點，最大熵法估計巖土參數概率分布精度和穩健性更高。參考文獻[15]采用信息熵推斷巖石力學參數概率分布模型，指出最大熵法根據樣本信息和統計方法推斷巖石力學參數的概率分布，具有更充分的數學和物理意義。

同時，通過對既有研究工作的總結分析發現，既有研究工作中一直忽略了一個基礎性問題，目前各種巖土參數概率分布擬合方法，不論是采用參數估計還是非參數估計法，都將直方圖作為評價概率分布擬合優劣的參考基準。目前頻率直方圖個數和區間的劃分主要是研究人員憑主觀決定。

為此，本文在既有研究基礎上完成如下工作：采用最大熵法擬合巖土參數測試樣本，根據巖土參數測試樣本的特點，確定最大熵法最佳的矩的階數；在最大熵法矩的階數確定的情況下，研究直方圖分組不同對巖土參數概率密度擬合曲線形狀和誤差的影響；借助判定系數，提出最優直方圖分組方法；通過工程實際測試樣本驗證本文方法的有效性和合理性，完善最大熵法在巖土參數概率密度函數擬合上的應用。

1 最大熵法

1.1 基本原理

最大熵原理的概率分布估計是一個優化問題[16-19]，如公式（1）。

式中：S 為信息的均值，mi為第i 階原點矩，N為原點矩的階數。f（z）是隨機變量z 的概率密度函數。若是用N+1 個參數形式表示的，則可用某種數值優化方法通過該模型求解出待定的參數。

通過式（1），改變f（z）可使熵達到最大值。S[f（z）]求最大熵分布的概率密度函數，可用經典微分方法求解式（1）。

設L 為拉格朗日函數，其拉格朗日乘子為λ0，λ1，…λN，于是：

令導數dL/df（z）等于零，則得：

歸并積分號下的各項，得到：

由于積分式等于零，因此其中一種情況是被積函數等于零，由此可得：

式（6）即為最大熵分布概率密度函數的解析表達式，也是式（1）的一個全城最優解。

1.2 乘子λi 求解

確定各個乘子λi（i=0，1，…，N），即可用最大熵分布表示變量z 的隨機特性。

由此可得λ0的計算公式：

為求乘子λ1，λx，…，λN，可將式（8）對λi微分，得

或

上式等號的右邊就是第i 階原點矩的負值，即

另外，若將式（9）對λi微分，又可得：

由式（12）和（13）可得：

上式就是求λ1，λ2，…，λN的方程。對式（14）作如下處理：

最后求

就可得到問題的解，其中Ri為殘差。

求解式(16)時，需要輸入一個初始點λi0(i=0，1，…，N)，有一個好的初始點，對算法滿意地收斂很重要。將式（1）的約束條件寫成如下數值形式：

式中：aj是數值積分乘子。求解式(17)線性方程組可得N+1 個未知的f(zj)，將它代入式(6)，則有f(zy)=exp（λαλ1zj+…λNzjN）；j=1，2，…，N+1 個方程，由此解得的λj(j=0，1，2，…，N)值即作為初值。

1.3 最大熵法最佳矩的階數

為了避免出現計算時發生溢出的可能性，可將z 值域變換到[0,1]之間來計算，為此，后文對巖土參數測試樣本進行了極差歸一化處理[21]，即：

一般，當概率密度函數曲線形狀比較復雜時，采用最大熵法擬合時矩的階數要取大一點，但是必須注意，如果樣本容量比較小，高階矩的統計值會因誤差較大而失去意義。

2 直方圖分組不同對擬合效果的影響

2.1 既有直方圖分組方法

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

其中，n 為樣本個數。

2.2 評價指標

K-S 檢驗法的檢驗量是整個取值范圍內的最大偏差值，在一定的顯著性水平下，可評估某一分布擬合的有效性，但不能提供某一分布擬合是否良好的絕對信息。為了評價擬合方法的優劣，引入一個無量綱的判定系數，它將擬合估計值與實際值對比分析，是直觀評價擬合效果的指標[21]。判定系數定義如下

式中：yi為實測值，為實測值的均值，Yi為擬合估計值。由表達式可知，若實測值yi與擬合估計值Yi之間誤差越小，則判定系數越接近1，擬合效果也就越好。

2.3 直方圖不同分組對擬合曲線影響

下面以兩組巖土參數測試樣本為例，說明直方圖分組不同對擬合曲線形狀和誤差的影響。選擇參考文獻[2]中提供的黏聚力和內摩擦角測試樣本，共81 組，文中取直方圖分組均為10，按經驗公式(18)計算直方圖分組為9。

采用最大熵法擬合兩組巖土參數測試樣本的概率分布，矩的階數取為5～6，直方圖分組不同時，擬合效果（判定系數）對比見表1，限于篇幅在此僅給出黏聚力樣本不同直方圖分組下擬合曲線對比如圖1-2 所示。

對于兩組巖土參數測試樣本，由表1、圖1-2 可見：對內摩擦角樣本（矩的階數為5～6）和黏聚力樣本（矩的階數為6），在直方圖人為分組取10 和按經驗公式取9 時，均未獲得最佳的擬合效果，特別是內摩擦角測試樣本在直方圖人為分組取10 時的判定系數較直方圖分組為8 的判定系數相差較大。

圖1 矩的階數為6 時黏聚力樣本最大熵法擬合曲線對比

表1 擬合效果對比

綜上所述，無論是通過擬合曲線的直觀對比，還是通過檢驗值的定量檢驗，最大熵法擬合優劣取決于直方圖的分組，如何確定巖土參數測試樣本的最優直方圖分組至關重要。

圖2 矩的階數為5 時黏聚力樣本最大熵法擬合曲線對比

3 基于最優直方圖的最大熵法擬合實例分析

通過上述分析可見，直方圖分組不同，判定系數大小不同，在確定的最大熵法矩的階數前提下，通過尋找最大的判定系數來確定最佳直方圖分組是簡單可行的方法。

仍以上述兩組巖土參數測試樣本為例，通過最優直方圖分組數的確定方法，確定兩組測試樣本的最優直方圖分組數。矩的階數取5～6 時，在最佳直方圖分組下的判定系數見表2。在最優直方圖分組下，矩的階數取為6 時，丙組巖土參數測試樣本的擬合曲線如圖3 所示。

圖3 最優直方圖分組下的擬合曲線

由表1、表2、圖1-3 可見，矩的階數為5～6 時，最優直方圖分組下的判定系數均大于表1 內的直方圖分組下的判定系數，最優直方圖分組下最大熵法擬合誤差最小，擬合效果最佳。

表2 最優直方圖分組下的擬合判定系數

為了說明確定最優直方圖分組的必要性及論文方法的普適性，選取參考文獻[13]提供的1 組摩擦因數小樣本，樣本個數為31，選取參考文獻[11]提供的1 組液限小樣本，樣本個數為26。圖4(a)為摩擦因數測試樣本在矩的階數取為5，直方圖分組為7 時（文獻中的直方圖分組）的最大熵法擬合曲線，判定系數為0.8866。圖4(b)為最優直方圖分組為6時的最大熵法擬合曲線，判定系數為0.8975。圖5(a)為液限測試樣本在矩的階數為5，直方圖分組為6時的最大熵法擬合曲線，判定系數為0.9068，圖5（b）為最優直方圖分組為5 時的最大熵法擬合曲線，判定系數為0.9437。

圖4 摩擦因數最大熵法擬合曲線

圖5 液限最大熵法擬合曲線

針對兩組大樣本和兩組小樣本，無論是通過擬合曲線的直觀對比，還是通過檢驗值的定量檢驗，均證明采用最優直方圖分組后的最大熵法擬合精度均大于既有文獻中直方圖分組下的擬合精度。

5 結論

首次提出直方圖分組的確定是最大熵法擬合巖土參數效果評價的前提和基礎，在已有研究基礎上，完善最大熵法對巖土參數概率密度函數擬合。既有的直方圖確定方法人為因素較大，不能客觀地反應樣本概率分布。

實測樣本的統計結果表明，采用最優直方圖分組下最大熵法擬合巖土參數概率分布的判定系數最大，擬合精度最高，結果驗證了基于最優直方圖分組的最大熵擬合方法的有效性和適應性。