劉 丞，范文亮,2，余書君，李正良,2

(1.重慶大學(xué)土木工程學(xué)院，重慶 400045；2.山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué))，重慶 400045)

工程結(jié)構(gòu)可靠度分析的主要目標(biāo)是估計(jì)包含各類不確定性隨機(jī)變量的實(shí)際工程結(jié)構(gòu)的失效概率[1-2]。針對此類問題，已發(fā)展了一次可靠度方法(FORM)[3-6]、二次可靠度方法(SORM)[7-10]、代理模型方法[11-13]、矩方法[14-15]和蒙特卡羅模擬方法(MCS)[16-18]等一系列結(jié)構(gòu)可靠度分析方法。

迄今為止，一次可靠度方法仍是應(yīng)用最為廣泛的可靠度分析方法。其中，HASOFER 和LIND[3]首先提出了一類驗(yàn)算點(diǎn)迭代算法-驗(yàn)算點(diǎn)法，通過在標(biāo)準(zhǔn)正態(tài)空間中搜索從原點(diǎn)到極限狀態(tài)曲面的最小距離來獲得驗(yàn)算點(diǎn)坐標(biāo)，結(jié)構(gòu)的失效概率可通過該最小距離解析計(jì)算。RACKWITZ 和FLESSLER[5]將該算法進(jìn)一步擴(kuò)展到了非正態(tài)隨機(jī)變量的一般情況，即HLRF 方法。若引入Rosenblatt 變換[19]或Nataf 變換[20-21]，亦可將其方便地拓展至相關(guān)變量情形。上述一次可靠度方法的迭代過程簡單、易于實(shí)現(xiàn)，但需求解迭代過程中所有迭代點(diǎn)處的函數(shù)值和梯度值，且對于隱式功能函數(shù)的梯度值計(jì)算尚需引入差分法[22]，計(jì)算效率不很理想。引入代理模型是改善一次可靠度方法效率的有效途徑之一。目前，將代理模型中的Kriging 模型與一次可靠度方法相結(jié)合主要有兩種方式：其一，通過隨機(jī)選定的點(diǎn)集建立Kriging 模型，并以此模型近似計(jì)算迭代點(diǎn)處的函數(shù)值和梯度值，從而實(shí)現(xiàn)可靠度分析[23]；其二，通過隨機(jī)選定的初始點(diǎn)集建立初始Kriging 模型，以此模型計(jì)算梯度值并確定新的迭代點(diǎn)，并將迭代點(diǎn)的計(jì)算結(jié)果用于更新Kriging 模型，從而實(shí)現(xiàn)驗(yàn)算點(diǎn)的搜索[24]。仔細(xì)考察不難發(fā)現(xiàn)，前者可視為Kriging 模型和一次可靠度方法的簡單組合，由于僅由隨機(jī)點(diǎn)集進(jìn)行Kriging 模型的一次建模，計(jì)算精度依賴于該Kriging 模型對功能函數(shù)的近似精度，此外由于點(diǎn)集的隨機(jī)性必然導(dǎo)致計(jì)算結(jié)果的隨機(jī)性；后者則將Kriging 模型和一次可靠度的迭代過程實(shí)現(xiàn)了較好的融合，但初始點(diǎn)集的隨機(jī)性依然會(huì)導(dǎo)致計(jì)算效率的隨機(jī)性，且若初始Kriging 模型精度不高時(shí)可能會(huì)弱化一次可靠度方法的收斂性能，此外該方法未能充分利用Kriging 模型的統(tǒng)計(jì)特性。

圖1 示出了一次可靠度分析的典型迭代收斂過程(即后文中算例1 中工況1 的迭代收斂結(jié)果)，仔細(xì)考察不難發(fā)現(xiàn)，整個(gè)迭代過程可分為兩個(gè)階段：1)全局搜索階段，即前期迭代階段，此階段相鄰迭代點(diǎn)相距較遠(yuǎn)；2)局部搜索階段，即后期迭代階段，此階段迭代點(diǎn)在一個(gè)較小的局部區(qū)域內(nèi)波動(dòng)。若以全局搜索階段的迭代點(diǎn)為基礎(chǔ)建立Kriging 模型，且充分利用該模型的統(tǒng)計(jì)特性進(jìn)行Kriging 模型的主動(dòng)學(xué)習(xí)，進(jìn)而實(shí)現(xiàn)局部搜索階段迭代點(diǎn)的函數(shù)值及梯度值的自適應(yīng)計(jì)算，將有助于改善一次可靠度方法的計(jì)算效率。為此，基于上述思路，本文提出了一種基于主動(dòng)學(xué)習(xí)Kriging模型的改進(jìn)一次可靠度方法。

圖1 一次可靠度方法迭代過程Fig.1 Iterative process of FORM

1 通用一次可靠度方法

假設(shè)結(jié)構(gòu)可靠度問題的功能函數(shù)為：

式中：X={X1,X2,…,Xn}為基本隨機(jī)向量；n為隨機(jī)變量的數(shù)量。對Xi引入等概率變換，即：

式中：Ui為與Xi對應(yīng)的標(biāo)準(zhǔn)正態(tài)變量；Φ(·)為其累積分布函數(shù)；為Xi的累積分布函數(shù)Fi(·)的反函數(shù)。于是，式(1)可改寫為：

式中：U={U1,U2,…,Un}，且Ui和Uj的相關(guān)系數(shù)ρU,ij可由Xi和Xj的相關(guān)系數(shù)ρX,ij依下式確定，即[25]：

式中：μi和σi分別為Xi的均值和標(biāo)準(zhǔn)差；Hk(·)為Hermite 多項(xiàng)式；φ(·)為標(biāo)準(zhǔn)正態(tài)變量的概率密度函數(shù)；mi為由Xi的分布確定的常數(shù)。

式(3)所對應(yīng)的可靠指標(biāo)β 和驗(yàn)算點(diǎn)u＊為[20]：

顯然，u*=和β 需由式(6)～式(8)迭代確定。后文中將在迭代過程中由式(8)更新的點(diǎn)稱為迭代點(diǎn)。由式(6)和式(7)可見，在每一迭代步中均涉及h(U)在迭代點(diǎn)處的函數(shù)值和梯度值。

若h(U)為隱式功能函數(shù)，則u*處的梯度值可由差分法計(jì)算如下：

2 基于主動(dòng)學(xué)習(xí)Kriging 模型的改進(jìn)一次可靠度方法

考慮到一次可靠度方法在迭代搜索驗(yàn)算點(diǎn)時(shí)的收斂特性，可將迭代過程分為全局搜索階段和局部搜索階段兩部分。文中以第k次迭代點(diǎn)u(k)與第k-1 次迭代點(diǎn)u(k-1)的距離||u(k)-u(k-1)||達(dá)到某一較小的閾值c作為兩個(gè)階段的劃分準(zhǔn)則，若||u(k)-u(k-1)||≥c，則屬于全局搜索階段，否則屬于局部搜索階段。文中取c=0.1。

2.1 全局搜索階段

在全局搜索階段，以均值點(diǎn)u(1)為初始迭代點(diǎn)，按照第1 節(jié)的通用一次可靠度方法確定迭代點(diǎn)，直至達(dá)到||u(k)-u(k-1)||<0.1。且將此階段最終的k值記為q。

2.2 局部搜索階段

在局部搜索階段，將全局搜索階段的迭代點(diǎn){u(1),u(2),…,u(q-1)}及在迭代點(diǎn)處計(jì)算梯度值的差分點(diǎn)作為初始訓(xùn)練點(diǎn)集，基于此建立初始Kriging 模型，根據(jù)Kriging 模型在后續(xù)迭代點(diǎn)處的預(yù)測精確性自適應(yīng)地更新模型并借助更新后的Kriging 模型提供預(yù)測值和梯度值，將Kriging 模型的更新與通用一次可靠度方法的迭代融合在一起，直至收斂。

2.2.1 初始Kriging 模型

為簡便，將初始訓(xùn)練點(diǎn)集記為u={u(1),u(2),···,u(M)}，相應(yīng)函數(shù)值其中M=(n+ 1) (q- 1)，由此可得h(U)的Kriging 模型(U)為[26-28]：

式中：f(U)為U的多項(xiàng)式基函數(shù)，文中取f(U)=[1,U]；P為回歸系數(shù)向量；m(U)為方差為σ2的零均值高斯過程，其協(xié)方差滿足：

式中：Rθ(u(i),u(j))為含參數(shù)θ 的相關(guān)函數(shù)方程，文中采用高斯相關(guān)函數(shù)，即：

由式(10)可給出h(U)在任意點(diǎn)u處的預(yù)測值和標(biāo)準(zhǔn)差，即：

式中：F為M×n維基函數(shù)矩陣；R為訓(xùn)練點(diǎn)之間相關(guān)函數(shù)矩陣；r為u與各訓(xùn)練點(diǎn)之間的相關(guān)函數(shù)列向量；R與r可分別表示為：

此外，由式(10)可知u處的函數(shù)值z滿足：

Kriging 模型的詳細(xì)介紹可參考文獻(xiàn)[26 - 27]，且Kriging 模型的建模及運(yùn)算均可由MATLAB 軟件的DACE 工具箱實(shí)現(xiàn)。

根據(jù)顯式表達(dá)式(13)，可方便地推導(dǎo)出梯度的顯式表達(dá)式；由于梯度值的計(jì)算亦可由DACE工具箱直接計(jì)算，其表達(dá)式不再詳述。

實(shí)用中，利用DACE 工具箱計(jì)算Kriging 模型在u處的預(yù)測值、梯度值以及方差命令為“[y,dy,mse]=predictor(u,dmodel)”。

2.2.2 Kriging 模型的誤差估計(jì)

為考察Kriging 模型在u處的預(yù)測精確性，令：

由式(17)可得：

式中，Pr{·}表示概率。若取S(u)=50，則z的誤差小于6%的概率達(dá)到99.7%。

值得指出的是，在AK-MCS[29]方法中，盡管Kriging 模型是針對功能函數(shù)建立的，但是其誤差估計(jì)是針對功能函數(shù)的示性函數(shù)展開的[29-30]，因此僅適用于主動(dòng)學(xué)習(xí)Kriging 模型(AK)與Monte Carlo 模擬(MCS)相結(jié)合的方法中。相比較而言，文中的式(19)直接針對功能函數(shù)進(jìn)行誤差估計(jì)，為主動(dòng)學(xué)習(xí)Kriging 模型與其他類型可靠度分析方法的結(jié)合奠定了基礎(chǔ)。

2.2.3 Kriging 模型的更新

同理，若已得到第l次更新的Kriging 模型(U)及與之對應(yīng)的迭代點(diǎn)u(q+l)，可根據(jù)式(18)計(jì)算(U)在u(q+l)處的S(u(q+l))，并根據(jù)S(u(q+l))進(jìn)行Kriging 模型更新，即：

式中，P(l+1)、m(U)(l+1)為將u(q+l)加入訓(xùn)練點(diǎn)集重新更新后P、m(U)的對應(yīng)項(xiàng)。

不難發(fā)現(xiàn)，上述的Kriging 模型更新屬于基于主動(dòng)學(xué)習(xí)的Kriging 模型更新，S(·)則對應(yīng)著其中的學(xué)習(xí)函數(shù)。

2.3 改進(jìn)一次可靠度方法實(shí)現(xiàn)流程圖

綜合全局搜索階段與局部搜索階段，本文建議的改進(jìn)一次可靠度方法的流程圖如圖2 所示。

圖2 改進(jìn)一次可靠度方法流程圖Fig.2 Flow chart of the improved FORM

仔細(xì)考察不難發(fā)現(xiàn)，改進(jìn)一次可靠度方法盡管引入了Kriging 模型，但是Kriging 模型所需的初始訓(xùn)練點(diǎn)集以及訓(xùn)練點(diǎn)的增加過程均是確定的，因此該方法的計(jì)算結(jié)果亦是確定的，與文獻(xiàn)[23 - 24]的隨機(jī)性結(jié)果具有顯著區(qū)別。

3 算例分析

文中分別通過數(shù)值算例和工程算例驗(yàn)證建議方法的有效性，其中各算例均視為隱式功能函數(shù)情形。為考察改進(jìn)一次可靠度方法的精度與效率，將其與通用的一次可靠度方法[21](即第1 節(jié)介紹方法，記為FORM)以及兩種Kriging 模型與一次可靠度方法相結(jié)合的方法(即文獻(xiàn)[23]和文獻(xiàn)[24]的方法)進(jìn)行對比。文獻(xiàn)[23]采用LHS 方法抽取N1個(gè)樣本點(diǎn)，計(jì)算對應(yīng)的響應(yīng)值并構(gòu)建Kriging 模型，利用Kriging 模型替代功能函數(shù)并結(jié)合一次可靠度方法進(jìn)行可靠度計(jì)算；文獻(xiàn)[24]采用蒙特卡洛法隨機(jī)抽取N2個(gè)樣本點(diǎn)(N2≥n+1)，計(jì)算對應(yīng)的響應(yīng)值并構(gòu)建Kriging 模型，基于此模型計(jì)算梯度值并結(jié)合一次可靠度方法確定新的迭代點(diǎn)，并將迭代點(diǎn)的計(jì)算結(jié)果用于更新Kriging 模型。顯然，樣本點(diǎn)集的隨機(jī)性必然導(dǎo)致Kriging 模型的隨機(jī)性，因此，與FORM 不同，文獻(xiàn)[23]與文獻(xiàn)[24]的計(jì)算結(jié)果是隨機(jī)的。

由于文獻(xiàn)[23]和文獻(xiàn)[24]方法不適用于相關(guān)變量情形，因此將兩種方法中的一次可靠度方法代之以文獻(xiàn)[20]的通用一次可靠度方法，且分別簡稱為K-FORM-1 和K-FORM-2。由于文獻(xiàn)[23]中未給出隨機(jī)點(diǎn)數(shù)量的確定準(zhǔn)則，后文中的K-FORM-1均采用30 個(gè)樣本。為了算法比較的公平性，各方法迭代過程的收斂準(zhǔn)則統(tǒng)一采用||u(k+1)-u(k)||<10-6。

各方法精度比較時(shí)，以抽樣數(shù)為108的蒙特卡洛模擬方法(記為MCS)結(jié)果作為標(biāo)準(zhǔn)解，各方法的誤差為：

式中：Value 為對應(yīng)方法的可靠指標(biāo)；MCS 為蒙特卡洛模擬方法計(jì)算的可靠指標(biāo)。

各方法效率比較時(shí)，主要考察其調(diào)用功能函數(shù)的次數(shù)。此外，由于K-FORM-1 和K-FORM-2中初始選點(diǎn)的隨機(jī)性導(dǎo)致可靠度計(jì)算結(jié)果在精度與效率上產(chǎn)生隨機(jī)性。因此為了體現(xiàn)其隨機(jī)性，上述兩種方法的可靠度計(jì)算結(jié)果均為分別運(yùn)算5 次的均值。

算例1.數(shù)值算例

考察如下功能函數(shù)[2,31]：

式中：X1～N(0,0.5),X2～N(0,0.5)且變量相互獨(dú)立。根據(jù)k值的不同可分為2 種工況。

1)工況1：k=5

本文建議方法基于FORM 在全局搜索階段迭代5 次后，迭代點(diǎn)u(6)與u(5)之間的距離(||u(6)-u(5)||=0.0863<0.1)滿足設(shè)定精度，進(jìn)入局部搜索階段。在局部搜索階段基于改進(jìn)一次可靠度方法繼續(xù)迭代了7 次，其中調(diào)用功能函數(shù)計(jì)算迭代點(diǎn)函數(shù)值并以此更新Kriging 模型4 次，因此本文建議方法的函數(shù)調(diào)用次數(shù)為5×(1+2)+4=19。結(jié)果示于表1。此外，表1 亦給出了抽樣次數(shù)為108的MCS、FORM、K-FORM-1 和K-FORM-2 的計(jì)算結(jié)果；表2 給出了K-FORM-1 和K-FORM-2 各5 次運(yùn)算的結(jié)果，其中由于K-FORM-1 各次運(yùn)算效率相同，因此僅給出可靠指標(biāo)，K-FORM-2 各次的可靠指標(biāo)基本相同，因此僅給出功能函數(shù)調(diào)用次數(shù)。不難發(fā)現(xiàn)，建議方法與FORM 精度一致但效率最高；K-FORM-1 在與FORM 效率一致時(shí)，其精度是變化的，但總體而言，略低于FORM；K-FORM-2 與FORM 精度一致且效率高于FORM，但其效率是變化的。

表1 工況1 的可靠度計(jì)算結(jié)果Table 1 Reliability results for Case 1

表2 已有方法的隨機(jī)性分析Table 2 Stochastic analysis of existing methods

2)工況2：k=10

各類方法的計(jì)算結(jié)果如表3 所示，K-FORM-1和K-FORM-2 各5 次運(yùn)算的結(jié)果示于表2；其中，由于K-FORM-1 的5 次運(yùn)算中僅2 次收斂，因此表3 中K-FORM-1 的結(jié)果為這2 次運(yùn)算的平均值，而K-FORM-2 的5 次運(yùn)算均不收斂。由上述結(jié)果可知，對于此工況，建議方法與FORM 精度基本相當(dāng)，但效率顯著提升；而在工況1 中較FORM 更優(yōu)的K-FORM-2 不能給出收斂的結(jié)果；K-FORM-1 仍存在精度不高的問題，同時(shí)收斂性亦受到影響。

表3 工況2 的可靠度計(jì)算結(jié)果Table 3 Reliability results for Case 2

總體而言，算例1 的兩種工況對比驗(yàn)證了本文建議方法能夠在保證精度與FORM 一致的情況下顯著提高計(jì)算效率；且相比K-FORM-1 和K-FORM-2，既避免了結(jié)果的隨機(jī)性，尚能更好地兼顧精度、效率與收斂性。

算例2.抗力-效應(yīng)工程算例

已知某工程功能函數(shù)形式為[1-2]：

式中：永久荷載效應(yīng)G服從正態(tài)分布，μG= 5.3，σG=0.371；可變荷載效應(yīng)Q服從極值Ⅰ型分布，μQ= 7，σQ= 2.03；抗力R服從對數(shù)正態(tài)分布，μR=30.92，σR= 5.26。

各類方法的計(jì)算結(jié)果如表4 所示。由于K-FORM-1 和K-FORM-2 在計(jì)算結(jié)果的精度與效率上分別存在隨機(jī)性，因此將兩種方法各5 次運(yùn)算的結(jié)果示于表5。不難發(fā)現(xiàn)，建議方法與FORM精度一致但效率最高；K-FORM-1 相比FORM 而言計(jì)算效率有所提高，但其精度具有隨機(jī)性，且平均誤差大于FORM；K-FORM-2 與FORM 精度一致且效率高于FORM，但其效率具有隨機(jī)性。比較而言，建議方法既避免了結(jié)果的隨機(jī)性，尚能更好地兼顧精度、效率。

表4 算例2 的可靠度計(jì)算結(jié)果Table 4 Reliability results for Example 2

表5 已有方法的隨機(jī)性分析Table 5 Stochastic analysis of existing methods

算例3.柔性基礎(chǔ)沉降算例

考慮一個(gè)矩形基礎(chǔ)沉降問題，功能函數(shù)為[32]：

式中：基礎(chǔ)寬度B為30 m；基礎(chǔ)長度L為40 m；基礎(chǔ)嵌入深度D為3 m；地層厚度H為10 m；角數(shù)m為4；影響因素I1、I2和IF分別為0.073、0.089、0.95；均布荷載q0服從均值為280 kPa、標(biāo)準(zhǔn)差為40 kPa 的極值I 型分布，泊松比ν服從均值為0.25、標(biāo)準(zhǔn)差為0.08 的正態(tài)分布，彈性模量Es服從均值為80 MPa、標(biāo)準(zhǔn)差為2.5 MPa 的對數(shù)正態(tài)分布。

各類方法的計(jì)算結(jié)果如表6 所示，其中K-FORM-1 進(jìn)行5 次運(yùn)算的可靠指標(biāo)在[3.4018,3.4163]內(nèi)波動(dòng)，K-FORM-2 進(jìn)行5 次運(yùn)算的函數(shù)調(diào)用次數(shù)在[23, 27]內(nèi)波動(dòng)。比較而言，不難發(fā)現(xiàn)各方法精度基本相當(dāng)；建議方法仍然效率最高，K-FORM-2 優(yōu)于K-FORM-1，且兩者均優(yōu)于FORM。

表6 算例3 的可靠度計(jì)算結(jié)果Table 6 Reliability results for Example 3

算例4.矩形截面懸臂梁算例

考察如圖3 所示的矩形截面彈性懸臂梁的位移可靠度問題，功能函數(shù)為[21]：

圖3 矩形截面懸臂梁Fig.3 Cantilever beam with rectangular section

式中：L為梁的跨度；q為均布荷載；E為材料的彈性模量；b、h分別為截面的寬和高，假設(shè)上述變量中E、L為確定性變量，分別取值為2.6×105MPa、6000 mm。q服從均值為1000 kN/m、標(biāo)準(zhǔn)差為200 kN/m 的對數(shù)正態(tài)分布，h服從均值250 mm、標(biāo)準(zhǔn)差為37.5 mm 的對數(shù)正態(tài)分布，b服從均值為100 mm、標(biāo)準(zhǔn)差為20 mm 的正態(tài)分布。假設(shè)b和h為相關(guān)變量，相關(guān)系數(shù)為0.6。

各類方法的計(jì)算結(jié)果如表7 所示。其中，采用建議方法時(shí)，由式(4)確定的等效相關(guān)系數(shù)為0.61888，與文獻(xiàn)[33]直接求解二維積分方程得到的結(jié)果0.61887 非常吻合。

表7 算例4 的可靠度計(jì)算結(jié)果Table 7 Reliability results for Example 4

由表7 結(jié)果可知，建議方法與FORM 精度一致但效率最高；K-FORM-1 在效率上略高于FORM，但其精度低于FORM 且具有隨機(jī)性；K-FORM-2與FORM 精度一致但效率低于FORM，且其效率亦具有隨機(jī)性。

4 結(jié)論

為了進(jìn)一步提高一次可靠度方法在可靠度分析時(shí)的計(jì)算效率，本文將一次可靠度方法分為全局搜索階段和局部搜索階段，在全局搜索階段沿用通用的一次可靠度方法迭代過程，在局部搜索階段充分利用全局搜索階段的迭代結(jié)果建立初始Kriging 模型，并引入基于主動(dòng)學(xué)習(xí)的Kriging 模型更新，借助更新后的模型提供迭代點(diǎn)處的預(yù)測值和梯度值，發(fā)展了基于主動(dòng)學(xué)習(xí)Kriging 模型的改進(jìn)一次可靠度方法，文中分析結(jié)果表明：

(1) 建議方法能夠較好的處理具有隱式功能函數(shù)的結(jié)構(gòu)可靠度問題，且能夠在保證精度與通用一次可靠度方法一致的條件下顯著提高其計(jì)算效率；

(2) 相比已有的結(jié)合Kriging 模型和一次可靠度方法的組合方法而言，建議方法對于特定的結(jié)構(gòu)可靠度分析問題，其構(gòu)建Kriging 模型的樣本點(diǎn)集是唯一的，避免了可靠度分析結(jié)果的隨機(jī)性，且能更好地兼顧精度、效率與收斂性。