最大熵法密度函數的求解及在樁板墻穩定分析中的應用

劉春剛

(中鐵第一勘察設計院集團有限公司，西安 710043)

1 最大熵法概率密度函數的確定和擬合

1.1 概率密度函數的確定

在概率論中，熵的定義［1］如式(1)所示

式中，S為熵;f(x)為概率密度函數。最大熵法的理論基礎是最大信息熵原理:在給定所有滿足約束條件的概率密度函數中，信息熵最大的概率密度函數為最佳的概率密度函數。設x為一個連續型隨機變量，其概率密度函數f(x)滿足以下條件

式中，Mi為x的i階原點矩，可通過統計樣本計算確定。這樣最大熵法就轉化為在式(2)、式(3)約束下求式(1)最大值問題。構造拉格朗日函數，如式(4)、式(5)所示

式中，λ0，λ1，λ2，…，λn為待定系數，由式(2)、式(3)構成的非線性方程組求解。

式(7)給出了n階矩法概率密度函數表達式，韋征等［3］對n值的選取問題進行的研究，以一個單層網殼結構可靠度分析為驗證算例，比較了最大熵理論三階矩算法、四階矩算法及五階矩算法的計算結果。通過對比發現三階矩算法得到的概率密度不能正確地反應結構位移響應變化規律，與四階矩算法相比，五階矩算法尾部誤差較大，已經不能有效地提高概率密度函數、超越概率和可靠指標的精度，并且求解待定系數的方程組過于復雜，采用數值解法求解方程組時對初始值要求嚴格，不容易得到能使方程收斂的初始值。最終韋征等建議n值取4，即采用四階矩算法。

1.2 擬合檢驗

根據最大熵原理求得函數只有通過擬合檢驗才能確定其性能的好壞，即能否最大程度反映隨機變量本身概率特性。由于總體分布未知，在檢驗過程中自然就不會牽扯到總體分布的參數，對于這種總體參數未知的檢驗，被稱為“非參數檢驗”，常用的非參數檢驗方法有K-S檢驗法、變量值隨機性檢驗等，以K-S檢驗為例對最大熵法求得概率密度函數進行檢驗，其具體過程如下:

(1)原假設根據最大熵法得到的最大概率密度函數積分求得概率分布函數F0(X);

(2)利用樣本數據計算個樣本數據點的累計概率得到檢驗累計概率分布函數S0(X);

我國高鐵走出國家的形勢雖然良好，有著很大的潛力，但一些挑戰是不可避免的，激烈的競爭、國際形勢、未知的不利因素都是阻礙高鐵發展的重要因素。

(3)計算F0(X)和S0(X)在相應的變量值點x上的差D(x)，得到差最大值Dmax;

(4)對于給定的檢驗水平α和樣本容量n，查Kolmogorov分布的分位數表，得到臨界值Dn;

(5)判斷:若Dmax≤Dn時，則接受原假設，否則否定原假設。

2 概率密度函數的求解

最大熵法擬合概率密度函數的核心和難點在于以下非線性方程組的求解

式中，Mi統計樣本 x的 i階原點矩;λi為待定系數。

方程組的求解需要解決三方面的問題:①各方程積分域的選擇;②方程組簡化;③方程組求解方法。

2.1 積分域選取

方程組式(8)中各方程都含有一個不定積分項，直接求解非常困難，需要先把各方程的不定積分計算出來，去掉方程中的積分符號，簡化方程組。由于被積函數比較復雜，需要人為的選擇一個積分區域，用被積函數在該區域上的定積分代替原不定積分，然后采用數值方法進行積分運算。趙國藩［2］建議積分區域選為［μ-10σ，μ+10σ］并用梯形積分公式進行積分運算。通過算例對比發現積分區域的選取是否合適對擬合的好壞有重要影響。以均值為25，標準差為11.2極值Ⅰ型分布為例，當積分區域分別選擇［0，400］和［0，80］時擬合的概率密度函數如圖(1)、圖(2)所示。

由圖1、圖2對比發現，積分區域取［0，400］時得不到正確的擬合函數，擬合的效果很差，最大誤差為0.76。而積分區域取［0，80］時，最大誤差量級為10-2，二者相差76倍。通過算例驗證發現，積分區域應控制在［μ-10σ，μ+10σ］以內為宜，同時積分區間端點處對應的概率密度函數值不宜超過10-4。

圖1 積分域［0，400］

圖2 積分域［0，80］

2.2 方程組簡化

觀察第一個方程，作以下變換

即

即

式(12)代入式(13)得

此時，最大熵法系數方程組就由式(8)化簡為式(14)，即將原五元的方程組簡化為四元的方程組，提高了方程組收斂速度和穩定性。由方程組(14)得到λ1、λ2、λ3、λ44 個系數后，按照式(10)求解 λ0，從而得到所有的待定系數。

2.3 方程組求解方法

非線性方程組(8)比較復雜，求解起來比較麻煩，已有許多文獻對方程組的求解方法進行了探索研究。趙國藩等［2］利用循環中點求積Newton算法對非線性方程組牛頓算法進行改進，改進后的新算法具有收斂快、穩定、對初始值無限制，成功率高的優點。此外，還有一些文獻［4-5］應用遺傳算法、粒子群算法等優化算法求解，都取得了一定的成果。但是這些計算方法對數學知識要求較高，并且需要編寫較為復雜的程序才能實現。

Matlab具有強大的數值計算能力，并且整合了許多優化算法、非線性方程組數字解法程序作為可以直接調用的內部函數，求解方程組比較方便。在求解實例時發現直接調用Matlab優化工具箱中的遺傳算法程序去求解方程組(14)，不需要設置初始值，但是得到的解精度并不高，調用ga函數求解方程組一般上所能達到的精度為0.1～0.01，很難得到精度更高的解。在Matlab中還可以用Fsolve函數求解非線性方程組，該方法可以得到精度為10-7以及更好的解，但是使用Fsolve函數對初始解的要求非常嚴格，初始解選擇不恰當將導致程序不收斂。比較上述兩種算法的優缺點后，將兩種方法有機的結合起來，先使用遺傳算法得到一個精度為0.1—0.01的解，然后將這個解作為方程組的初始解調用Fsolve函數得到精度更好的解。這種方法前后兩步都是直接調用Matlab內置函數，實現起來非常簡單，并且求解成功率很高。本文分別用遺傳算法和遺傳算法與Fsolve相結合的方法對四種分布進行了擬合，四種分布的參數見表1。

表1 分布參數

分別使用上述兩種計算方法對上述四種分布進行計算，計算結果如圖3～圖10所示，其中圖3～圖6為遺傳算法計算結果;圖7～圖10為遺傳算法和Fslove函數相結合的計算結果。

對以上模擬結果的精度進行對比分析，結果如表2所示。

圖3 分布1遺傳算法擬合結果

圖4 分布2遺傳算法擬合結果

圖5 分布3遺傳算法擬合結果

圖6 分布4遺傳算法擬合結果

圖7 分布1改進算法擬合

圖8 分布2改進算法擬合

圖9 分布3改進算法擬合

圖10 分布4改進算法擬合

表2 各分部計算精度對比

由表2可以看出，單純用遺傳算法計算的結果精度不太高，概率密度函數擬合的誤差在0.3～0.1之間，將這樣的擬合概率密度函數代替真實密度函數進行可靠度計算會帶來較大的誤差。而采用遺傳算法和Fsolve函數相結合的算法得到的擬合函數最大擬合誤差均在0.03以下，精度較高。其中正態分布(分布1)擬合精度由原來的0.25提高至10-5，提高了2 400倍;分布3的最大擬合誤差為0.02，代數與遺傳算法相比精度提高了9倍。因此用遺傳算法得到擬合函數的待定系數后，再以此結果為初始解用Fsolve進行一次求解，可以有效地提高計算精度。解決方程組(8)求解困難的問題，并且該方法直接調用Matlab內置函數，不需要編寫實現遺傳算法和非線性方程組數值解法的程序，實際操作起來簡便易行。

3 樁板墻的極限狀態方程與優化模型

3.1 極限狀態方程

以路塹式樁板墻為研究對象，建立其可靠度分析方程，極限狀態方程的一般形式為

式中，R為抗力;S為荷載。

樁板墻破壞主要是樁發生破壞，其破壞形式主要有受彎破壞、受剪破壞和錨固段內側土體受壓破壞三種形式，因此根據這三種形式求其極限方程。

按照式(15)極限狀態方程的形式把受彎和受剪的極限方程分別表示為

式中，Ωpw、Ωpj分別為受彎和受剪計算模式的隨機變量;Mmax、Qmax分別為樁板墻受到的最大彎矩和最大剪力;Mu、Qu分別為樁板墻極限彎矩和極限剪力。

以“m”法為例

式中，A3、A4、B3、B4、C3、C4、D3、D4為關于樁埋深的函數;x0為滑面處的位移;M0、Q0分別為樁錨固段頂點的彎矩和剪力。對于式(20)和式(21)分別求一階偏導為零，即可得到Mmax和Qmax表達式［6］。

3.2 隨機變量的確定

綜上所述，式(16)、式(17)為失效模式下的偏微分方程。對于抗力項隨機變量，在式(16)中主要考慮縱向鋼筋受拉強度fyu、縱向鋼筋截面面積As、混凝土軸心抗壓強度fcu、樁的截面寬度b和有效高度h0以及分項系數Ωpw;在式(17)中主要考慮箍筋受拉強度fyvu、混凝土抗拉強度ftu、樁的截面寬度b和有效高度h0以及分項系數Ωpj。對于荷載項由于計算復雜性，暫時直接按照分項系數法直接計算。

4 算例驗證

以路塹式樁板墻為例，樁長23 m，樁間距6 m，樁采用2 m×3 m的T形樁，設計懸臂端長8 m，錨固段長15 m，懸臂端主要位于膨脹土和粉質黏土，錨固段主要位于圓粒土、強風化片巖和弱風化片巖中，局部位于粉質黏土中，其中設計土體容重取γ=18 kN/m3，綜合φ=28°，內力計算按照“m”法，其中膨脹土取值3 000 kPa/m2，圓粒土和強風化片巖取值12 000 kPa/m2，弱風化片巖120 000 kPa/m，土壓力分布按照三角形計算。樁底支撐采用自由端。

對于荷載的計算，利用鐵一院研發的路基助手計算，其荷載分項系數按照K=1.35考慮，對于材料和樁的結構相關參數如表3所示。

表3 計算資料

對于彎矩利用最大熵原理求解得到，γ0=-0.919，γ1=0.020，γ2=-0.500，γ3=-3.03，γ4=0，失效概率為3.34e-8，可靠指標 β=5.399;同樣對于剪力 γ0=-0.743，γ1=0.031，γ2=-0.459，γ3=-2.975，γ4=0.04，失效概率為7.95e-7，可靠指標β=4.864。

5 結論

(1)通過選擇方程的積分區域和對方程組的簡化，然后利用matlab中Fsolve函數和ga工具箱相結合的方式，實現了對基于最大熵函數概率密度函數方程組簡單和精確求解。

(2)通過對路塹式樁板墻計算，然后利用熵函數分別得到受彎、受剪和錨固段樁側土體可靠指標均值，而且一般只需要較少樣本即可得到相對精確的概率密度函數和可靠指標。

(3)樁板墻的破壞屬于延性破壞，按照文獻［7］的規定，可靠度標準不低于3.2，所以計算出來的彎矩可靠指標5.399和剪力可靠指標4.864均滿足要求。

［1］ E.T.JaYNES.Information Theory and Stastical Mechanics［J］.The Physical Review，1957，106(104):620-630.

［2］ 趙國番.工程結構可靠性理論與應用［M］.大連理工大學出版社，1996:137-142.

［3］ 韋征，葉繼紅，沈士釗.最大熵法可靠度理論在工程中的應用［J］.振動與沖擊，2007，26(6):146-151.

［4］ 謝世坤，段芳，李強征，等.非線性方程組求解的三種Newton法比較［J］.井岡山學院學報(自然科學)，2006，27(8):8-11.

［5］ 陳長憶，葉永春.基于粒子群算法的非線性方程組求解［J］.計算機應用與軟件，2006，23(5):137-139.

［6］ 新型支擋結構設計與工程實例［M］.北京:人民交通出版社，2010:275-296.

［7］ 中華人民共和國建設部.GB50068—2001 建筑結構可靠度統一標準［S］.北京:中國建筑工業出版社，2001.