張道兵 ，楊小禮，朱川曲，王衛(wèi)軍

(1. 中南大學 土木工程學院，湖南 長沙，410075；2. 煤礦安全開采技術湖南省重點實驗室，湖南 湘潭，411201)

由于隧道襯砌結(jié)構的參數(shù)具有隨機性，其分布具有多樣性，且極限狀態(tài)功能函數(shù)具有高度非線性的特征，所以，在進行其可靠度計算中，要從參數(shù)估計、極限狀態(tài)功能函數(shù)和結(jié)構可靠度計算方法這3個方面進行考慮，任何一方面的較大誤差都會給計算結(jié)果帶來較大偏差，從而影響計算結(jié)果的準確性。在參數(shù)估計方面，最常用的確定變量概率分布的方法是先假設隨機變量的概率分布，然后，通過檢驗來確定拒絕或不拒絕假設，或者按照主觀經(jīng)驗來認定數(shù)據(jù)是否服從某種分布。這些做法給處理結(jié)果添加了主觀因素，使得結(jié)果與客觀事實相偏離。而運用最大熵原理可通過數(shù)據(jù)樣本估計出最接近真實情況、不添加任何主觀因素的概率密度分布，保證參數(shù)估計的準確性。在可靠度計算方法方面，目前主要有蒙特卡羅數(shù)值模擬法、隨機有限元法、響應面法[1]、矩法(包括一次二階矩法[2]、一次三階矩法[3]、二次二階矩法[4]、二次四階矩法等)、最優(yōu)化法共5類。針對隧道襯砌結(jié)構，其極限狀態(tài)功能函數(shù)是高度非線性的，若采用矩法計算結(jié)構的可靠度，需要求極限功能函數(shù)的導數(shù)，不僅計算量很大，而且當功能函數(shù)高度非線性時，矩法精度低，存在一定的誤差；若采用隨機有限元、響應面法和蒙特卡羅數(shù)值模擬法計算其結(jié)構可靠度，需編制專門的程序，計算較復雜，要求掌握編程語言的算法和編程技巧，難度比較大。即使在Matlab環(huán)境下實現(xiàn)蒙特卡羅法直接抽樣計算可靠度[5-6]，也需要掌握 Matlab的語法和編程技巧。到目前為止，運用最優(yōu)化原理求解可靠度指標的方法很多，如梯度下降法、懲罰函數(shù)法、序列二次規(guī)劃法、Lagrange乘子法、H-L法及其改進H-L法等，并對各種優(yōu)化進行了驗證。但這些最優(yōu)化法對算法和編程要求高，且效率低，程序通用性較差，對一般工程技術人員來說，具有相當?shù)碾y度。也有一些研究者直接利用 Matlab優(yōu)化工具箱實現(xiàn)最優(yōu)化方法計算結(jié)構可靠度指標，但由Matlab優(yōu)化工具箱中的所有函數(shù) (除 linprog之外)得到的僅是局部優(yōu)化解，要求所求模型必須為凸規(guī)劃模型[7]，計算結(jié)果才是準確的，而隧道襯砌結(jié)構可靠性模型不符合凸規(guī)劃條件，為此，本文作者根據(jù)可靠度指標的幾何意義，建立隧道襯砌結(jié)構可靠度指標計算的優(yōu)化模型，利用Microsoft Excel中的規(guī)劃求解器進行可靠度指標的優(yōu)化求解，彌補了上述各計算方法的不足，快速、有效地計算出隧道襯砌結(jié)構的可靠度指標，并采用蒙特卡羅100萬次直接抽樣法對結(jié)果進行了驗證，證實其結(jié)果精確度高，適用性好。在極限狀態(tài)功能函數(shù)方面，采用現(xiàn)行隧道設計規(guī)范的襯砌截面抗拉和抗壓檢算式，建立襯砌截面抗拉極限狀態(tài)方程和抗壓極限狀態(tài)方程。本文在該極限狀態(tài)功能函數(shù)的基礎上，采用最大熵原理計算得出襯砌結(jié)構參數(shù)的統(tǒng)計特征，并利用Microsoft Excel的規(guī)劃求解器直接由最優(yōu)化方法計算其結(jié)構可靠度，且采用蒙特卡羅法對結(jié)果進行驗證。

1 基于最大熵原理的概率密度函數(shù)的確定與擬合檢驗

1.1 概率密度函數(shù)的確定

最大熵原理的定義為：在所有滿足給定約束條件下的許多概率密度函數(shù)中，信息熵最大的密度函數(shù)就是最佳密度函數(shù)[8]。熵最大意味著人為假定最小，從熵作為不確定性的度量來看， 此時解包含的主觀成分最少，因而是最客觀的[9]。其表達式可寫為：

其滿足約束條件為：

式中：f(x)為概率密度函數(shù)；R為積分區(qū)間；m為所用矩的階數(shù)；mi為樣本的第i階原點矩。

采用經(jīng)典微積分方法，并引入拉格朗日算子，可求得最大熵密度函數(shù)解析式為：

式中：λ0，λ0，…，λm為拉格朗日乘子。

根據(jù)式(4)，只要確定參數(shù)λ就可以完全確定f(x)，為了便于數(shù)值求解，經(jīng)數(shù)學變換可得到：

式中：Ri為殘差，可以用數(shù)值計算的方法把它減小到接近于0。用非線性規(guī)劃求這些殘差平方和的最小值，就可以得到問題的解：

當R＜ε或所有的|Ri|＜ε時(其中，ε為規(guī)定的允許誤差)，則認為上式收斂。

1.2 擬合檢驗

根據(jù)最大熵原理求得的最大熵密度函數(shù)必須通過擬合良好性檢驗，才能從理論上確定是否能夠反映隨機變量真實的概率特性?？蛇x用精度高的K-S檢驗法(或A-D檢驗法)，把隨機變量的實測數(shù)據(jù)擬合到一個常用的較為合適的分布上，然后，通過對經(jīng)驗分布函數(shù)與假設理論分布函數(shù)的檢驗來確定拒絕或是接受假設。其檢驗步驟為[10]：

(1) 設由樣本矩求得的最大熵密度函數(shù)式f(x)積分得到總體的概率分布函數(shù)為F(x)，假設：F(x)=F0(x)；

(2) 設經(jīng)驗概率分布函數(shù)為Fn(x)，計算統(tǒng)計量

(3) 給定顯著性水平α，通常取α=0.10或α=0.05；

(4) 根據(jù)α和樣本容量n在K-S 檢驗表上查得臨界值Dn,α；

(5) 統(tǒng)計推斷：當Dn≤Dn,α時，則接受原假設分布；當Dn＞Dn,α時，否定原假設，選用另一種概率分布。

2 最優(yōu)化方法計算可靠度指標的數(shù)學模型與映射變換

2.1 數(shù)學模型

結(jié)構可靠度的最優(yōu)化方法是指從結(jié)構可靠度指標的幾何含義出發(fā)，采用可靠度指標的優(yōu)化數(shù)學模型，直接求解目標函數(shù)值(可靠度指標)的結(jié)構可靠度計算方法。設X1，X2，…，Xn是結(jié)構中n個任意分布的獨立隨機變量,由這些隨機變量表示的結(jié)構極限狀態(tài)方程為：

采用R-F(拉科維茨-菲斯萊法)法將非正態(tài)變量當量正態(tài)化,得到等效正態(tài)分布的均值及可靠度指標

根據(jù)結(jié)構可靠度指標的幾何意義，優(yōu)化模型可表示為：

約束條件為：

由此將求解可靠度指標β的問題轉(zhuǎn)化為有約束條件下的極小值問題。即在式(12)的約束條件下所確定的可行域內(nèi)搜索式(11)的最優(yōu)解，β的最小值即為可靠度指標。

2.2 分布的映射變換

可靠度指標β是指標準化正態(tài)坐標系中原點到極限狀態(tài)曲面的最短距離，因此，需要將非標準正態(tài)隨機變量Xi映射為標準正態(tài)隨機變量Yi。對于隧道襯砌結(jié)構參數(shù)中的正態(tài)分布和對數(shù)正態(tài)分布，可進行如下映射變換。

(1) 當Xi服從正態(tài)分布時，

(2) 當Xi服從對數(shù)正態(tài)分布時，lnXi服從正態(tài)分布，

2.3 Microsoft Excel規(guī)劃求解器實現(xiàn)最優(yōu)化方法計算可靠度指標的步驟

規(guī)劃求解是Excel工作表軟件的一個加載工具，用于求解復雜的方程及各類線性或非線性有約束優(yōu)化問題。其計算可靠度指標的步驟為：新建1個Excel 表格，首先在A單元格中依次輸入變量名稱，在B單元格中依次輸入對應變量符號，在C1單元格中輸入目標函數(shù)(在Excel規(guī)劃求解工具中，包含公式的目標單元格代表的是目標函數(shù))β的計算公式，在C2，C3，…，Ci單元格中輸入決策變量X1，X2，…，Xn(在Excel規(guī)劃求解工具中，可變單元格代表的是決策變量)的初值，一般可設初值為 0，如果變量為分母或取對數(shù)可設初值為1；然后，選擇 “工具” 欄中的 “規(guī)劃求解”選項，在“約束”欄中“添加”約束條件(用單元格表示)，并在“選項”按鈕中根據(jù)需要設置迭代方法、精度、允許誤差、收斂度；最后，在對話框中選擇目標單元格的值為“最小值”，單擊對話框右上的“求解” 按鈕，彈出 “規(guī)劃求解結(jié)果” 對話框，得到可靠指標β及正態(tài)空間中的驗算點Y*，得到結(jié)果，計算結(jié)束。

3 隧道襯砌結(jié)構的極限狀態(tài)方程與優(yōu)化模型

根據(jù)現(xiàn)行隧道設計規(guī)范的襯砌截面抗拉和抗壓檢算式，建立襯砌截面抗拉極限狀態(tài)方程和抗壓極限狀態(tài)方程，抗拉極限狀態(tài)屬正常使用極限狀態(tài)，而抗壓極限狀態(tài)則屬承載能力極限狀態(tài)。

(1) 當偏心矩e0≤0.2t(t為襯砌厚度)時，截面由抗壓強度控制承載能力，相應的抗壓極限狀態(tài)方程為：

式中：N極限為襯砌混凝土所能承受的極限軸力(即抗力)；N為計算所得的截面軸力(即載荷效應)；kPR為抗力計算模式不定性；b為縱向?qū)挾?，? m；t為截面厚度；α為偏心影響系數(shù)；Ra為混凝土抗壓強度。

偏心影響系數(shù)α的概率分布為正態(tài)分布，均值和變異系數(shù)計算公式為：

其優(yōu)化模型為：

(2) 當偏心矩e0＞0.2t時，截面由抗拉強度控制承載能力，相應的抗拉極限狀態(tài)方程為

式中：M為計算所得截面彎矩；RI為混凝土極限抗拉強度；kPS為荷載效應計算模式不定性；其余變量含義同上。

其優(yōu)化模型為：

上述極限狀態(tài)方程中，荷載效應N和M的可由蒙特卡羅有限單元法求得，其統(tǒng)計特征可由最大熵原理求得，其他隨機變量的值可通過大量的現(xiàn)場調(diào)查、試驗獲得，其統(tǒng)計特征由最大熵原理求得。

4 工程計算

以四川1個鐵路工程隧道某段為例。該段圍巖較軟，裂隙發(fā)育，巖體不完整。按照單線鐵路隧道襯砌標準進行設計，根據(jù)蒙特卡羅有限單元法、現(xiàn)場調(diào)查、測量、實驗室試驗獲得原始數(shù)據(jù)。

現(xiàn)以隧道襯砌厚度為例進行統(tǒng)計特征估計。采用探地雷達測得50個襯砌厚度，其樣本平均值為0.600 1 m, 標準差為0.090 2，通過前面所述方法求解其最大熵密度函數(shù)為：

其中：殘差平方和R為8.28×10-6。由K-S法進行檢驗, 假設樣本服從正態(tài)分布，統(tǒng)計量Dn=0.064 08。給定顯著性水平α= 0.10，樣本容量n=50，由K-S檢驗表上查得臨界值Dn,α=0.169 59。Dn＜Dn,α，所以，接受原假設，并可判定襯砌厚度為正態(tài)分布。

同理，依次可求得其它參數(shù)的統(tǒng)計特征，如表1所示。

表1 基本隨機變量統(tǒng)計特征Table 1 Probability property for basic random variable

按承載能力抗壓極限狀態(tài)和抗拉極限狀態(tài)分別建立隧道襯砌結(jié)構極限狀態(tài)方程，將上述參數(shù)代入式(18)和(20)，按照Microsoft Excel規(guī)劃求解器實現(xiàn)可靠度指標計算步驟，且設置迭代方法為牛頓法，精度為0.000 001，允許誤差0.1%，收斂度為0.000 1，可得：抗壓可靠度指標βra=4.474 3，抗拉可靠度指標βrl=2.108 8；在普通雙核計算機上運行的時間為0.6 s。

采用蒙特卡羅法 100萬次直接抽樣法計算結(jié)果為：βra=4.473 9，βrl=2.108 3。前者計算結(jié)果與此結(jié)果非常接近。并與文獻[11]中建議的目標可靠指標4.200 0和2.000 0比較，該隧道襯砌結(jié)構的可靠度能滿足工程可靠性要求。其主要參數(shù)(襯砌厚度t、襯砌厚度的變異系數(shù)Vt、混凝土抗拉極限Rl、混凝土抗拉極限的變異系數(shù)彎矩M的變異系數(shù)VM、軸力N的變異系數(shù)VN、混凝土抗壓極限Ra)與可靠度的關系曲線如圖1～7 所示(其中以式(14)進行計算，VN和Ra以式(12)進行計算)。

(1) 從圖1、圖2可知：隧道襯砌厚度對結(jié)構可靠度指標影響非常大，是影響結(jié)構可靠度的主要參數(shù)。

(2) 從圖3、圖4、圖7可知：混凝土抗拉極限與混凝土抗壓極限及其變異系數(shù)對可靠度指標有一定的影響，但不顯著。

圖1 襯砌厚度與可靠度指標關系曲線Fig.1 Relationship between lining thickness and reliability index

圖2 襯砌厚度的變異系數(shù)與可靠度指標關系曲線Fig.2 Relationship between variation coefficient of lining thickness and reliability index

圖3 混凝土抗拉極限與可靠度指標關系曲線Fig.3 Relationship between tensile ultimate stress of concrete and reliability index

圖4 混凝土抗拉極限的變異系數(shù)與可靠度指標關系曲線Fig.4 Relationship between variation coefficient of tensile ultimate stress of concrete and reliability index

圖5 彎矩的變異系數(shù)與可靠度指標關系曲線Fig.5 Relationship between variation coefficient of bending moment and reliability index

圖6 軸力的變異系數(shù)與可靠度指標關系曲線Fig.6 Relationship between variation coefficient axial force and reliability index

圖7 混凝土抗壓極限與可靠度指標關系曲線Fig.7 Relationship between pressure ultimate stress of concrete and reliability index

(3) 從圖5、圖6可知：彎矩和軸力的變異系數(shù)對可靠度指標的影響較大，在設計過程中，系統(tǒng)考慮彎矩和軸力的離散性和隨機性是必要的。

5 結(jié)論

(1) 為保證隧道襯砌結(jié)構可靠度計算精度，必須從參數(shù)估計、極限狀態(tài)功能函數(shù)和結(jié)構可靠度計算方法3個方面著手，任何一個方面出現(xiàn)較大誤差都會給結(jié)果造成較大偏差。

(2) 采用基于最大熵原理的方法計算隧道襯砌結(jié)構參數(shù)的概率密度函數(shù)，并通過擬合檢驗確定分布的正確性，可以滿足工程可靠度計算對參數(shù)的要求。

(3) Microsoft Excel規(guī)劃求解器實現(xiàn)最優(yōu)化方法計算可靠度指標不需要編程，只需要簡單的Excel操作基礎，即可進行可靠度計算，操作非常簡便，可滿足一般工程技術人員的要求。采用該方法計算隧道襯砌結(jié)構可靠度無需將狀態(tài)函數(shù)線性化，不受基本變量維數(shù)限制，其收斂速度快，計算效率高，且與蒙特卡羅100萬次直接抽樣法計算結(jié)果比較，具有很高的精度。

(4) 提高隧道支護結(jié)構的可靠度可從以下幾個方面進行：增加襯砌厚度；減小襯砌厚度的變異系數(shù)；增加混凝土抗拉、抗壓極限；減小混凝土抗拉、抗壓極限的變異系數(shù)。該方法還能廣泛適用于交通工程、水利工程、防護工程、巷道工程等隧道及地下工程領域的可靠度計算和分析。

[1]Bucher C G, Bourgund U. A fast and efficient response surface approach for structural reliability problems[J]. Structural Safety,1990, 7(1): 57-66.

[2]Karamchandani A, Cornell A C. Sensitivity estimation within first and second order reliability methods[J]. Structural Safety,1992, 11(2): 95-107.

[3]Tichy M. First-order three-moment reliability method[J].Structural Safety, 1994, 16: 189-200.

[4]Hohenbichle M, Rackwitz R. Improvement of second-order reliability estimates by importance sampling[J]. Journal of Engineering Mechanics, ASCE, 1994, 120(1): 2195-2203.

[5]馮曉波, 楊樺. 用MATLAB實現(xiàn)蒙特卡羅法計算結(jié)構可靠度[J]. 中國農(nóng)村水利水電, 2002(8): 50-51.FENG Xiao-bo, YANG Hua. The structural reliability calculation based on Monte-carlo and Matlab[J]. China Rural Water and Hydropower, 2002(8):50-51.

[6]朱川曲, 張道兵, 朱海燕.基于蒙特卡羅法的煤巷錨桿支護結(jié)構可靠性分析[J]. 中國安全科學學報, 2008, 18(4): 146-150.ZHU Chuan-qu, ZHANG Dao-bing, ZHU Hai-yan. Reliability analysis on bolt support structure of coal roadway based on Monte-carlo[J]. China Safety Science Journal, 2008, 18(4):146-150.

[7]張道兵, 朱川曲, 劉澤, 等. 基于 MATLAB優(yōu)化工具箱的煤巷頂板錨桿支護結(jié)構可靠性分析[J]. 中國安全科學學報,2007, 17(9): 131-134.ZHANG Dao-bing, ZHU Chuan-qu, LIU Ze, et al. Reliability analysis on bolt support structure of coal roadway roof based on Matlab’s optimization toolbox[J]China Safety Science Journal,2007, 17(9): 131-134.

[8]程亮, 童玲. 最大熵原理在測量數(shù)據(jù)處理中的應用[J]. 電子測量與儀器學報, 2009, 23(1): 47-51.CHENG Liang, TONG Ling. Measurement data processing based on maximum entropy method[J]. Journal of Electronic Measurement and Instrument, 2009, 23(1): 47-51.

[9]楊杰, 胡德秀, 吳中如. 基于最大熵原理的貝葉斯不確定性反分析方法[J]. 浙江大學學報, 2006, 40(5): 810-815.YANG Jie, HU De-xiu, WU Zhong-ru. Bayesian uncertainty inverse analysis method based on pome[J]. Journal of Zhejiang University, 2006, 40(5): 810-815.

[10]叢培江, 張燕. 最大熵原理在壩體混凝土斷裂韌度反演中應用[J]. 武漢理工大學學報, 2008, 30(1): 83-86.CONG Pei-jiang, ZHANG Yan. Application on maximum entropy principle in inversion analysis of fracture toughness of dam concrete[J]. Journal of Wuhan University of Technology,2008, 30(1): 83-86.

[11]譚忠盛. 隧道支護結(jié)構體系可靠度的理論研究及其工程應用[D]. 成都: 西南交通大學土木工程學院, 1998: 115-128.TAN Zhong-sheng. Theoretical research and engineering application on reliability of tunnel structure[D]. Chengdu:Southwest Jiaotong University. School of Civil Engineering,1998: 115-128.