基于廣義子集模擬和自適應Kriging模型的非線性隨機動力系統的時變可靠性分析

唐和生， 郭雪媛， 薛松濤

(1. 同濟大學 結構防災減災工程系， 上海 200092； 2. 同濟大學 土木工程防災國家重點實驗室， 上海 200092)

工程結構的可靠性是指結構系統在規定的時間間隔內完成其預定功能的概率，近年來受到國內外學者的廣泛關注。結構系統在服役期內會受到地震等具有強隨機性的動力作用，同時結構系統自身材料特性與幾何參數也存在固有不確定性。由于兩種不確定性的來源及表征形式的不同，并且在非線性系統中疊加原理不再成立，經典的隨機振動理論難以推廣到非線性動力系統中，因此非線性隨機動力系統在非平穩隨機激勵下的時變可靠性分析仍然面臨巨大挑戰。

根據首次超越破壞準則的基本理論，現有時變可靠性研究方法主要分為基于跨越率的方法和基于極值的方法。基于跨越率的方法的核心是計算特定時間間隔內，隨機過程超過給定閾值的“跨越事件”發生次數的期望與失效概率的關系。這一概念最早由Rice[1]提出，隨后利用關注響應量的各種統計數據，基于跨越事件發生服從Poisson過程假設的解析方法得到進一步發展。近年來，Andrieu-Renuad等[2]引入一階可靠度方法(FORM)，將跨越率的計算轉換為并聯的靜態問題提出一種PHI2方法，實現了跨越率的高效求解。Hu等[3]針對跨越事件相互獨立假設的不足，提出一種具有較高精度的聯合跨越率方法求解時變可靠性問題。張德權等[4]提出了一種考慮認知不確定的PHI2方法，得到結構在設計周期內的時變可靠性區間。雖然基于跨越率與FORM相結合的方法的準確性已得到提高，然而在實際工程中，動力系統的時變可靠性分析是一個高維的隱式問題，這一類方法并不適用。目前針對動力系統，主要采用發展成熟的功率譜密度函數方法[5-6]。然而這類方法往往需要通過等效線性化方法才能推廣到非線性系統中[7]，因此存在較大誤差。

基于極值的方法的關鍵在于識別極限狀態函數的極值，對極值的不確定性進行量化獲得極值分布，這類方法不受系統是否為線性的約束。在極值概率分布已知的情況下，時變可靠性問題可以轉化為一個時不變問題。在實際應用中，獲取極限狀態極值的概率特性通常是十分困難的，特別是對于復雜的工程應用，Monte Carlo模擬(MCS)在計算上非常昂貴。為了解決計算效率問題，Wang等[8]提出了一種嵌套極值響應曲面(NERS)方法，在可靠性分析中消除了隨機過程，構造了一種時間響應曲面，用于評價極限狀態函數的極值響應。其他基于代理模型的方法也被廣泛研究[9-10]，然而絕大多數代理模型都存在維數詛咒且不適用于極限狀態函數與隨機過程存在雙向嵌套關系的動力系統。概率密度演化方法[11-12]也被用來近似極值分布，該方法的計算精度較高，但對多變量情形，樣本點數量較多導致計算耗時。數值模擬方法的目的是通過一些改進的采樣技術，如重要性采樣[13]和子集模擬[14]來估計時變可靠性。Hu等[15]提出一種針對極值分布的抽樣方法，把時變可靠度問題轉換為時不變問題進行求解，大大提高了計算效率。然而遺憾的是，即使是高效的計算方法，為得到離散時間區間內的累積失效概率，極限狀態函數改變，需要對每一區間進行重復運行，造成效率的明顯下降。本文采用Li等[16]提出的廣義子集模擬方法(generalized subset simulation, GSS)，實現僅通過單次運行得到所有時間區間內的累積失效概率，解決非平穩隨機激勵造成的高維輸入不確定性。

傳統的非線性動力系統的時變可靠性研究往往將兩種不確定分離，即分別考慮隨機系統在確定地震激勵下的可靠性問題和確定性系統在隨機激勵下的可靠性問題，難以實現兩者的統一處理。基于全概率定理的方法可以同時考慮系統的隨機性與外部隨機激勵，從多角度處理不確定性來源對可靠性的影響，分析過程清晰且可以將廣泛研究的Kriging代理模型應用于構建系統隨機參數與累積失效概率之間的非線性關系。考慮到計算資源有限性，構建代理模型的一個關鍵問題就是如何通過更“好”地選取樣本點來提高代理模型的精度。為此，采用Echard等[17]提出的一種結合Kriging(adaptive Kriging, AK)和原始MCS(Monte Carlo simulation, MCS)的主動學習可靠性方法(AK-MCS)，主動地訓練代理模型精確逼近重要區域的累積失效概率。

為解決非線性隨機動力系統的時變可靠性分析中系統參數不確定性與非平穩隨機激勵造成的計算成本難以承受的問題，本文提出一種基于廣義子集模擬和自適應Kriging模型的非線性隨機動力系統時變可靠性分析的GSS-AK-MCS方法。該方法采用全概率定理對兩種不確定性來源融合，將基于極值的時變可靠性分析與自適應代理模型方法相結合，顯著改善不確定傳播分析的計算效率，并且通過兩個數值算例對所提出方法的有效性進行探討。

1 隨機動力系統時變可靠性中的全概率定理

非線性隨機動力系統受到非平穩隨機地震激勵時的運動方程

(1)

(2)

式中：Pf(0,T)為累積失效概率;FE表示失效事件;P(FE/x)為X一次樣本下的條件累積失效概率;fX(x)表示隨機變量X的聯合概率密度函數，且Ω為其支撐集。考慮系統失效由單一響應隨機過程Y控制，對于隨機變量的一次樣本實現，動力系統參數確定，由累積失效概率定義

P(FE/x)=Pr(?t∈[0,T],Y(x,t)>b)

(3)

一般來說，式(2)中的精確數值積分難以計算，因此建立條件累積失效概率P(FE/x)在隨機變量空間內的代理模型，并在代理模型基礎上通過Monte Carlo方法直接計算雙隨機時變可靠性。特別地，在處理非正態和相關的輸入隨機變量時，采用全概率定理不需要額外的假設從而保證了計算結果的準確性。

2 累積失效概率的廣義子集模擬算法

2.1 累積失效概率的極值理論表達

動力系統的失效定義為系統的某一與位移或速度響應相關的量(例如，絕對位移、相對位移、控制點的應力或應變)超過給定閾值b(確定性的常量)。這種基于時變失效概率的可靠性評估問題也被稱為首次超越問題，根據結構或系統的動力反應的臨界值，考慮雙側界限的動力系統時變可靠性問題(隨機過程的絕對值超過給定閾值)。

傳統的Poisson過程法基于任意兩次跨越事件與其發生的時間統計獨立，服從(無記憶)Poisson過程假設。然而對于低閾值水平和/或窄帶過程，Poisson過程法得到的是累積失效概率的保守值，將低估時變可靠性。為克服Poisson過程法的限制，Vanmarcke認為交叉事件是成群出現的且各群之間相互獨立，在雙態Markov過程假設的基礎上提出了一種累積失效概率的改進估計。這類基于跨越率的方法對于高度非線性的動力系統求解技術難度大，不僅精度低而且計算成本極高。

基于極值的方法關注于工程系統在特定時間內的極值性能，在時變可靠性分析中，定義響應在時間段[0,T]內的極值變量

(4)

式中，θ為隨機激勵基本變量的一次樣本。極值法將動力系統的時變可靠性問題轉換為對應等效極值分布的求解，則累積失效概率可表達為

Pf(0,T)=P(ye(T)>b)

(5)

針對工程中常見的非平穩隨機過程，Priestley提出一種演變譜過程模型，表示為

(6)

式中：A(t,ω)是t與ω的確定性調制函數；Z(ω)是一個正交增量過程。如果用A(t)代替A(t,ω)，則非均勻調制演變譜退化為工程中常用的均勻調制演變譜。由非平穩隨機過程的演變譜表示理論，導出了非平穩隨機過程模擬的一個譜表示方法，其樣本函數是由余弦級數公式計算產生的

(7)

式中：Δω=ωu/N，ωu表示計算截止頻率，超過ωu所對應功率譜能量可忽略不計;θk為在[0,2π]區間上均勻分布、相互獨立的隨機相位角。因此非平穩隨機過程激勵可轉換為一組隨機變量θ=[θ1,…,θk]表示，時變可靠性問題被轉化為與時間無關的一個隱式非線性高維可靠性問題。

2.2 廣義子集模擬算法

基于抽樣的方法由于與維數無關且未對極限狀態函數做任何假設，目前是解決高維可靠性問題的唯一可行方法。針對高維小失效概率情況，MCS需要大量的計算樣本才能獲得較準確的失效概率，而子集模擬方法對于此類問題非常有效。然而時變可靠性分析需要對不同時間段內可靠性評估或當極限狀態函數閾值發生改變時，子集模擬方法需要進行反復計算，大大增加了計算成本。因此采用Li等提出的GSS，進一步推導基于極值的累積失效概率計算公式，解決時變可靠性分析這一多失效模式問題。

首先對原始子集模擬進行簡要介紹。子集模擬方法的基本思想是通過自適應地引入m層嵌套的中間事件，滿足FE1?FE2?…?FEm=FE，將小失效概率轉換成一系列相對較大的失效概率的乘積。中間事件與目標失效事件FE的表達相似，定義為FEi={ye>bi,i=1, …m}(b=bm>…>b2>b1)。其中，m為中間事件的總數。因此原始子集模擬的目標失效概率可表示為

P(FE)=P(FEm|FEm-1)P(FEm-1)=…=

(8)

在隨機動力系統的時變可靠性分析中，多時間點的累積失效概率計算是一個多極限狀態問題，原始的子集模擬方法需要重復計算。為克服這一局限性，采用改進后的GSS算法定義多個失效事件的并聯事件驅動隨機樣本的生成，高效求解多時間點的累積失效概率。

基于極值的GSS算法中，將時間段[0,T]離散化，時變可靠性的關注時間點依次為[T1,…,Tnt]，由首次超越破壞的定義，其相對應的累積失效概率呈隨時間遞增，即P(FE1)

(9)

式中，下標i表示第i個關注時間點。所有統一中間事件同樣滿足嵌套關系

G1?G2?…?Gm=G

(10)

因此與原始子集模擬方法相似，GSS算法首先根據不確定參數概率分布通過MCS產生第一層隨機樣本。然后采用MCMC方法不斷產生各層統一中間事件相應的樣本。為自適應地確定中間事件，設定所有的中間條件概率為一常數p0，通常為兼顧準確性與計算效率取值p0∈[0.1,0.3]。最終僅通過一次GSS計算得到各關注時間點Tj的累積失效概率

(11)

3 基于AK-MCS算法的時變可靠性分析

Kriging模型作為一種半參數化的方差估計最小的隨機過程插值算法，具有局部和全局的統計特征，已廣泛應用于工程中的可靠性分析。由于累積失效概率計算成本高昂，而在工程應用中往往針對某一閾值規定臨界失效概率，因此臨界失效概率對應的隨機參數樣本附近區域的計算精度更為重要。采用AK-MCS方法，使代理模型不斷優化，提高計算精度的同時節約計算成本。具體計算步驟如下：

步驟2采用拉丁超立方抽樣產生N0個初始試驗設計點xS={x(1),…,x(N0)}，計算相應的失效概率歸一化指標gS={g(x(1)),…,g(x(N0))}，構成樣本集{xS,gS}。

步驟3基于樣本集應用MATLAB中的DACE工具箱構建一系列Kriging模型，時間點Ti處的代理模型為

(12)

其中h(x)為基函數向量，β為回歸系數向量，Z(x)是一個均值為0服從正態分布的平穩高斯過程，其協方差由相關函數計算得到。

步驟4生成大量的MC樣本S={s(1),…,s(NMC)}，根據構建的代理模型計算每一樣本的預測值及方差。

步驟5定義U函數作為學習函數

(13)

步驟6判斷是否達到收斂條件。當達到收斂準則miniU(s(i))>2，i=1,…,NMC時，代理模型更新終止。在優化后的一系列代理模型上，運用MCS計算相應的累積失效概率及時變可靠性。

根據全概率定理，非線性隨機動力系統的時變可靠性問題轉化為一個雙層嵌套問題。首先通過拉丁超立方抽樣生成系統隨機參數空間內樣本點；在內層中針對高維輸入不確定性，通過GSS算法計算非平穩隨機激勵下累積失效概率；而在外層中解決系統參數不確定性，構建一系列Kriging代理模型擬合系統隨機參數與累積失效概率之間的非線性關系，并通過主動學習技術對代理模型更新；最后在更新完成的Kriging模型基礎上，基于全概率定理求解隨機系統的時變可靠性。本文所提出的非線性隨機動力系統的時變可靠性分析方法流程如圖1所示。

圖1 隨機動力系統的時變可靠性分析流程圖Fig.1 Flowchart of time-variant reliability analysis forstochastic dynamic systems

4 算例分析

4.1 附加非線性黏滯阻尼器的單自由度消能減震結構

以設有黏滯阻尼器的單自由度消能減震結構系統(如圖2)為例，研究GSS-AK-MCS方法在非線性動力系統時變可靠性分析中的效率和精度。

圖2 單自由度消能減震結構系統計算模型Fig.2 Analytical model of the SDOF-VD structural system

該系統在地震激勵作用下的運動方程為

(14)

(15)

式中：cd為阻尼系數；α為速度指數；sgn(·)為符號函數。系統不確定參數cd與α均服從正態分布，結合小概率事件發生規律一般遵從“3σ”法則，各隨機變量分布及代理建模范圍如表1所示。

表1 非線性系統隨機變量

結構所受非平穩隨機激勵的調制函數為

(16)

式中，T=20 s，ta=2.5 s，tb=10 s，β=0.1，計算步長Δt=0.02 s，功率譜密度函數取Kanai-Tajimi模型

(17)

GSS-AK-MCS方法中在隨機變量[μ-3σ,μ+3σ]范圍內，采用拉丁超立方抽樣均勻抽取16個初始樣本點，構建初始Kriging模型。通過AK-MCS方法更新代理模型，新增樣本點為20個。作為比較GSS-Kriging-MCS方法選取隨機參數變量樣本點總數相同為36個。以T=20 s時為例，兩種方法樣本點分布及累積失效概率計算結果分別見圖3、4與表2。由上述結果可以得到如下結論：① 與傳統Kriging代理模型相比，GSS-AK-MCS方法采用的自適應代理模型更新技術可以使樣本更集中分布于關鍵區域；② 在隨機變量樣本數量相同的情況下，GSS-AK-MCS方法實現了在計算成本相同時提高計算精度。

圖3 b=0.055 m時GSS-Kriging-MCS樣本點設計Fig.3 Design of sample points of the GSS-Kriging-MCS forb=0.055 m

圖4 b=0.055 m時GSS-AK-MCS樣本點設計Fig.4 Design of sample points of the GSS-AK-MCS forb=0.055 m

表2 單自由度消能減震結構系統累積失效概率Tab.2 Cumulative probability of failure for the SDOF structure system

為進一步研究GSS-AK-MCS方法在不同頻譜特征非平穩隨機激勵下的性能，對表3所示各計算工況采用本文所提出方法與MCS方法結果進行對比。各工況對Kanai-Tajimi模型中的地基土卓越圓頻率ωg和基巖擾動白噪聲強度S0的設計，滿足非平穩隨機激勵所包含總能量相同。以工況3為例，在T=20 s時采用GSS-AK-MCS方法構建的累積失效概率代理模型如圖5。GSS-AK-MCS方法中隨機變量樣本點共36個，平均每個樣本點通過GSS方法進行28 325次時程分析，為MCS方法時程分析次數的(36×28 325)/(104×106)≈0.010%，具有相當高的計算效率。

圖5 b=0.055 m時累積失效概率代理模型Fig.5 Surrogate model of cumulative probability of failurefor b=0.055 m

基于已構建的一系列代理模型，得到各工況累積失效概率隨持時的變化曲線如圖6所示。對各計算工況，以MCS方法結果作為累積失效概率的準確值，MCS方法計算結果在圖6中由實線表示。對以上結果進行分析可以得出，本文所提出的GSS-AK-MCS方法基于時變可靠性的極值理論，無需引入對首次超越事件發生次數的假設，在不同頻譜特征的非平穩隨機激勵情況下，均保持良好的精度；自適應代理模型的技術的引入大大降低了非線性隨機動力系統時變可靠性分析的計算成本。

表3 計算工況

圖6 累積失效概率計算結果Fig.6 Results for cumulative probability of failure

4.2 非線性多自由度剪切框架

為驗證本文所提出方法的有效性，本算例采用基于Bouc-Wen模型的五層非線性剪切框架。該非線性剪切框架的結構參數為：每層質量mi=50 kg，線性剛度ki=1.5×105N/m，阻尼系數ci=1 000 N·s/m，i=1,2,…,5。假設非線性滯回分量符合Bouc-Wen模型，結構在地震作用下的運動方程為

(18)

式中：α為屈服后與屈服前的水平剛度之比;Δu=ui-ui-1為層間位移;β、γ和n為Bouc-Wen模型無量綱滯回參數。在本算例中，結構非線性參數α=0.5，n=1.25，不確定參數β與γ均服從正態分布，各隨機變量分布及代理建模范圍見表4。結構所受非平穩隨機激勵與上一算例相同。

表4 Bouc-Wen模型隨機變量

表5 T=20 s非線性框架累積失效概率

圖7和圖8表明本文提出的自適應代理模型方法給予預測條件累積失效概率在臨界失效概率的鄰近區域內的樣本點更高的權重，提高代理模型在關鍵區域的預測精度。通過構建累積失效概率的代理模型響應面，代替非線性結構在非平穩隨機過程激勵下響應的復雜時程分析求解。表5中的對比分析驗證了與傳統Kriging代理模型和MCS方法相比，GSS-AK-MCS方法具有良好的計算精度并且顯著提高了時變可靠性分析的計算效率。

圖7 b=0.04 m時GSS-AK-MCS樣本點設計Fig.7 Design of sample points of the GSS-AK-MCS forb=0.04 m

圖8 b=0.04 m時累積失效概率代理模型Fig.8 Surrogate model of cumulative probability of failurefor b=0.04 m

5 結 論

本文提出了非線性隨機動力系統時變可靠性的GSS-AK-MCS方法，推導了累積失效概率的計算公式并構建了累積失效概率的代理模型，以一個裝配非線性黏滯阻尼器的單自由度系統和一多自由度非線性剪切框架為例驗證了所提出方法的有效性，主要結論如下：

(1) 基于全概率定理可將非線性隨機系統的時變可靠性分析中兩種來源的不確定性分別處理。采用GSS方法解決隨機激勵的高維不確定性，通過Kriging代理模型擬合累積失效概率與系統隨機參數的非線性關系，避免了理論推導的繁瑣和困難。

(2) 本文所提出的GSS-AK-MCS方法作為一種本質為數值模擬的方法，對于小失效概率問題的可靠性求解，與傳統的MCS方法相比具有相當高的計算效率，并且可以通過一次計算解決時變可靠性中多時間點造成的多失效模式問題。

(3) GSS-AK-MCS方法在不同頻譜特征的非平穩隨機激勵情況下均具有良好的計算精度，且由于自適應代理模型更新技術的引入，在不增加計算成本的條件下提高系統隨機參數重要區域的可靠性精度。因此，本文所提出的GSS-AK-MCS方法為同時考慮系統參數不確定性和隨機地震動作用的非線性動力系統時變可靠性分析提供了準確高效的計算方法，在工程中復雜隨機結構的時變可靠性分析中具有廣闊應用前景。