基于平穩非高斯結構響應前四階矩的首次穿越概率計算

張龍文， 盧朝輝,2， 趙衍剛,2

(1.中南大學 土木工程學院，長沙 410075； 2.中南大學 高速鐵路建造技術國家工程實驗室，長沙 410075)

在工程應用中，首次穿越概率用以描述結構在動力荷載作用下的動力可靠性問題。對此類問題，結構失效以其動力反應(如控制點的應力，或控制點、控制層的位移、速度、加速度等)首次穿越臨界值水平為標志[1]。臨界界限可以是屈服或極限應力、應變或位移。國內外眾多學者對此問題進行了大量的研究[2-7]，但是到目前為止，首次穿越問題甚至是最簡單的單自由度線性振子在高斯白噪聲激勵下的首次穿越概率都沒有精確的解析解[8-11]。在大多數情況下，假設結構響應是高斯的，以及結構響應的穿越服從Poisson分布。

當線性結構的激勵被定義為高斯過程模型時，結構的響應也是高斯的。在這種情況下，從概率與統計的角度來講，結構的響應可被均值和標準差完全定義。然而，當結構行為是非線性，或者激勵是非高斯過程，亦或兩者兼而有之，結構的響應將不是高斯過程，因此利用高斯模型計算首次穿越概率是不準確的。此時，對于非高斯過程的描述，需要用到更高階的統計矩。

非高斯過程的描述，可以通過一系列的分布方法[12-13]找到簡單的轉換函數進行非高斯過程與標準高斯過程之間的轉換。基于轉換的思路，許多學者提出了多種非線性模型用以描述非高斯過程。例如，基于結構反應矩(偏態系數和峰度系數等)，利用Gram-Charlier和Edeworth級數可以得到反應的全概率分布，進而變換為標準高斯過程[14]。但Gram-Charlier和Edeworth級數可能產生多個模型甚至負的概率密度函數和超越率。Winterstein提出了基于Hermite 矩模型的Winterstein (1988)多項式，利用結構響應的前四階矩將非高斯結構響應變換為標準高斯過程[15]。基于Wintertein(1988)的多項式，He等[16]考慮了穿越的群超效應和初始條件計算了平穩非高斯結構響應的首次穿越概率。然而，基于Winterstein (1988)多項式計算結果與Monte-Carlo模擬計算結果差別顯著，尤其是硬化反應情況中。該計算結果產生差異的原因在于轉換精度不足。另外，在一些相關文獻中[17-19]也發現，Winterstein (1988)多項式在軟化和硬化的非高斯過程中并不精確。因此，對于非高斯結構響應的首次穿越概率的計算，提高轉化模型的精度是必要的。

本文針對Winterstein (1988)多項式轉換的不足，采用精度更能滿足要求的Winterstein修正模型Winterstein(1994)[20]，以及Ding和Chen模型[21]，再利用考慮群超效應和初始條件的Poisson過程模型[22]，建立軟化和硬化非高斯結構響應的首次穿越概率解析表達式。以線性與非線性單自由度振動系統為例，對本文修正的方法、已有方法以及Monte Carlo模擬計算結果進行比較分析，驗證該方法在首次穿越概率計算中的有效性與準確性。

1 平穩非高斯隨機過程的首次穿越概率

在結構安全分析中，首次穿越概率Pf(T), 按照超越數服從泊松分布的假定表示為

Pf(T)=1-exp[-E[N+(T)]]

(1)

式中:E[N+(T)]是在時間段[0,T]內X(t)穿過界限的超越數。對于工程中感興趣的平穩情況，有E[N+(T)]=v+T,其中v+代表由Rice公式計算的平均超越率。它可以表示為

(2a)

(2b)

(2c)

式中:ω0是振動系統的無阻尼固有頻率。

如果結構響應是窄帶反應過程或者超越界限較低，此時穿越傾向于成群發生，需要對平均超越率進行修正[22-23]。當同時考慮平穩開始的初始條件及群超尺度〈cs〉時，式(1)可以被進一步改進為[24]

(3a)

其中

(3b)

(3c)

(3d)

Pf(0)=1-Φ[u(x)]

(3e)

式中:Pf(0)是T=0時的瞬時失效概率；RU(τ)是標準高斯過程U(t)的自相關函數，可以表示為

(3f)

2 基于修正模型的首次穿越概率計算

由于Winterstein(1988)多項式轉化精度的不足，本節首先介紹了Winterstein(1994)軟化模型以及Ding和Chen硬化模型，并給出了相應的等效高斯分位數。接著，對兩種模型的精度進行分析，討論了模型的精度。最后，說明了兩種模型在軟化與硬化非高斯結構響應的首次穿越概率計算過程。

2.1 Winterstein(1994)多項式及其等效高斯分位數

Winterstein等[25]提出了修正的Winterstein(1994)多項式

(4a)

式中:μX和σX分別為隨機過程X(t)的均值和標準差；U(t)是標準高斯過程；c3、c4以及κ為多項式系數，可通過偏態系數和峰度系數計算得到

(4b)

(4c)

(4d)

(4e)

對于修正的Winterstein(1994)多項式的等效高斯分位數可以表示為

(5a)

(5b)

2.2 改進的硬化模型及其等效高斯分位數

Ding等基于正交展開的方法提出了一個更為匹配的硬化模型(α4X<3)。該硬化模型可以表示為

X(t)=μX+σXg(U)

(6a)

式中

(6b)

j=b3/(3b4),k=(b2-b3α3X-b4α4X)/(3b4)

(6c)

式中:b2，b3，b4為模型系數，表示為

φ=[1-0.06(3-α4X)]1/3

(6d)

該模型的等效高斯分位數為

(7a)

式中，

放療的毒性反應較強，治療期間應采取以下護理措施：①照射野皮膚護理。照射野皮膚會產生發癢、紅斑、脫皮等癥狀，或發生放射性皮炎，不宜抓撓，避免冷熱刺激，不宜反復清洗，放射性皮炎患者可在患處涂抹冰片滑石粉；②骨髓抑制的處理。長時間大面積放療，會引起骨髓造血功能損傷，白細胞下降，病情嚴重時，應停止放療，應用升白細胞藥物，同時做好保護性隔離措施；③食道反應的護理。患者存在吞咽困難的問題，治療期間宜進流質/半流質食物，忌食粘性或帶骨頭、魚刺的食物，吞咽疼痛患者，可在餐前10 min口服適量黏膜表面麻醉劑，以緩解疼痛。

m1=-b4α3X,m2=b2-b3α3X-b4α4X+3b4,

(7b)

它的適用范圍為(1.35α3X)2+1.25≤α4X<3。因此，該模型適用于硬化非高斯過程。

2.3 模型精度分析

2.3.1 Winterstein(1994)多項式精度調查

由于Winterstein(1994)以及Winterstein(1988)的軟化模型為一元三次多項式，在此根據Fleishman[27]的矩匹配方法，計算多項式系數的準確值，用以說明Winterstein(1994)以及Winterstein(1988)的軟化模型的精度問題。對于一般的一元三次多項式可以表達為

(8)

在X(t)的前四階矩已知的情況下，基于矩匹配方法，多項式系數a1,a2,a3,和a4與前四階矩的關系，可以表達為[28]

a1+a3=0

(9a)

(9b)

(9c)

(9d)

簡化等式(9a)～(9d)，可以根據等式(10a)～(10f)計算a2與a4，表達如下

(10a)

3A1A3+3A4=α4X

(10b)

式中，

(10c)

(10d)

(10e)

(10f)

當a2與a4計算后，a1與a3可以通過下式計算，

(11)

圖1表示在不同偏態系數α3X下，隨著峰度系數α4X變化的Winterstein(1988)與Winterstein(1994)多項式系數以及準確值。從圖中的計算結果表明：

(1) Winterstein(1988)模型計算的各系數與準確值的差異大，且隨著偏態系數α3X的增大而變大。

(2) Winterstein(1994)模型計算的各系數改進了Winterstein(1988)模型的不足，能夠與準確值吻合，特別是a2與a4系數與準確值基本重合。另外，從圖1(f)中可以看出，當偏態系數α3X=2.0時，隨著峰度系數α4X的增大，Winterstein(1994)模型各系數的精度開始下降。因此，雖然Winterstein(1994)模型優于Winterstein(1988)模型，但更適用于α3X<2.0的情況。從等式(4b)～(4e)可以看出α3X的正負值只改變c3的正負號,因此，當α3X取負值時，同樣可以得到圖1的類似結果。

(a) α3X=0

(b) α3X=0.4

(c) α3X=0.8

(d) α3X=1.2

(e) α3X=1.6

2.3.2 Ding和Chen模型精度調查

對于等式(6a)所示Ding和Chen模型的逆函數形式可以表達為

(12)

對等式(12)取期望，可得，

(13)

該修正的硬化模型保證了標準高斯過程U(t)均值恒為0，改進了Winterstein(1988)模型中U(t)均值不為0的缺陷。另外，該模型的模型系數在文獻[21]中通過最小二乘法擬合得到，驗證了模型的精度滿足要求。為了進一步分析Ding和Chen模型的精度，通過計算該模型的偏態系數與峰度系數，并與它們的目標值(α3X,α4X)進行比較分析。該模型的前四階中心矩可以通過下式計算

(14a)

(14b)

式中φ(·)為標準正態分布的概率密度函數，表達為

(15)

則該模型計算的偏態系數γ3與峰度系數γ4可以表達為

(16)

為了便于分析Ding和Chen模型計算的偏態系數γ3與峰度系數γ4與目標值的誤差，將偏態系數與峰度系數的誤差值分別記為δα3X與δα4X。它們可以表達為

δα3X=γ3-α3X,δα4X=γ4-α4X

(17)

通過以上方法，計算得到偏態系數與峰度系數的誤差如圖2與圖3所示。圖2，圖3分別為偏態系數與峰度系數誤差結果。從圖中計算結果表明：

(1) 從圖2可以看出，Ding和Chen硬化模型偏度系數計算值與目標值誤差基本控制在-0.1與0.1之間。當峰度系數α4X=1.5～2時，離散性較大。當峰度系數α4X=2～3時離散性變小，且誤差趨近于0。

(2) 從圖3可以看出，Ding和Chen硬化模型峰度系數計算值與目標值誤差控制在-0.05與0.05之間。當偏態系數α3X=0時，模型峰度系數計算值與目標值的誤差最小。

2.4 基于修正方法的首次穿越概率

根據Winterstein(1994)以及Ding和Chen模型，基于改進方法計算平穩非高斯結構響應的首次穿越概率的步驟如下：

(1) 計算平穩非高斯過程的前四階矩。

(2) 確定轉換模型及其系數；如果峰度系數α4X>3，選用Winterstein(1994)模型，再根據式(4b)和(4e)計算轉換模型系數；如果峰度系數α4X<3，選用Ding和Chen模型，根據式(6b)～(6d)計算轉換模型系數。

(3) 如果峰度系數α4X>3，根據式(5a)～(5b)計算等效高斯分位數；如果峰度系數α4X<3，根據式(7a)～(7b)計算等效高斯分位數；確定結構響應的超越界限或超越水平x，將平穩非高斯結構響應映射為標準高斯過程。

(4) 最后，根據基于平穩高斯結構響應的Poisson模型即式(3a)～(3f)計算平穩非高斯結構響應的首次穿越概率。

如圖2和3所示。

圖2 Ding和Chen硬化模型偏態系數計算值與目標值誤差

Fig.2 bias between the estimations based on Ding and Chen model and the target values for skewness

圖3 Ding和Chen硬化模型峰度系數計算值與目標值誤差

Fig.3 bias between the estimations based on Ding and Chen model and the target values for kurtosis

3 算例分析

3.1 轉換模型的x-u變換

為了進一步調查Winterstein(1994)以及Ding和Chen模型的轉換精度，選取兩種分布：① Gamma分布(峰度系數=9>3)；② 截尾正態分布(峰度系數=2.794<3)進行分析。兩種分布的參數及其四階矩列于表1中。圖4和圖5分別對應了Gamma和截尾正態分布兩種情況的x-u變換。在圖4給出了利用Rosenblatt變換，Winterstein(1994)，Winterstein(1988) 多項式變換計算結果。圖5給出了利用Rosenblatt變

換，Ding和Chen模型，Winterstein(1988)多項式變換計算結果。由于Rosenblatt變換[29]完全保留了邊緣分布的信息，在此可以作為精確結果，用以作為已有轉換模型比較的標準。從圖4以及圖5結果說明如下：

(1) 在圖4中，隨著xs的逐漸增大，使用Winterstein(1988)多項式計算的結果與Rosenblatt變換計算結果差異也逐漸變大。

(2) 在圖5中，當xs的絕對值越來越大時，使用Winterstien(1988)多項式計算的結果與Rosenblatt變換計算結果差異也變大。

(3) 從圖4及圖5分別可以看出，Winterstien(1994)以及Ding和Chen計算結果基本與Rosenblatt變換結果重合，說明了Winterstien(1994)模型以及Ding和Chen模型轉換的精度及有效性。

表1 概率分布及其前四階矩

圖4 Gamma分布的x-u變換

圖5 截尾正態分布的x-u變換

3.2 平穩非高斯過程的首次穿越概率計算

選取線性與非線性單自由度系統為例，利用已有的He和Zhao模型，本文修正方法以及Monte Carlo模擬計算結構響應為軟化反應(α4X>3)與硬化反應(α4X<3)的首次穿越概率，并進行對比分析。

例1結構響應為軟化反應的首次穿越概率

二次力函數F(t)激勵線性單自由度系統(SDOF)如圖6所示。力函數F(t)=α1U(t)+α2U2(t)，其中U(t)代表一個平穩零均值的標準高斯過程，α1與α2是常數。

圖6 二次力函數F(t)激勵線性單自由度系統

Fig.6 A single-degree-of-freedom (SDOF) structure excited by a quadratic forcing functionF(t)

考慮阻尼比ζ=0.1和0.2兩種情況，根據文獻[30] 的方法，對于α2U2(t)激勵部分的結構響應前四階矩計算于表2中。從表2中可以看出在ζ=0.1和0.2兩種情況下計算的結構響應為軟化反應(α4X>3)。對應于表2的前四階矩，在界限x=3σX情況下，采用本文修正方法(α4X>3的情況)計算的首次穿越概率如圖7所示。同時，圖中給出了由Monte Carlo模擬[31](10 000個樣本)，He和Zhao方法和傳統高斯模型計算結果。

表2 平穩結構響應X(t)的前四階矩

圖7(a)和圖7(b)計算結果表明：

(1) 傳統高斯模型計算結果與Monte Carlo模擬結果有很大的差異，具有很大的保守性。

(2) 雖然He和Zhao模型相對于傳統高斯模型有很大的提高，但是與Monte Carlo計算結果還有一些差距。這些誤差產生的原因在于Winterstein多項式轉換的精度不夠。

(3) 本文修正的方法比較He和Zhao模型，以及傳統高斯模型有很大的提高，在整個計算的時間段內均能與Monte Carlo模擬結果更為吻合。因此，對于軟化反應，本文修正方法更為適用于首次穿越概率的計算。

(a) ζ=0.1

(b) ζ=0.2

例2結構響應為硬化反應的首次穿越概率

考慮平穩高斯白噪聲激勵Duffing振子的平穩響應。對于單邊譜密度為1/π時，該振子的運動方程為

(18)

式中:c>0為阻尼系數;ε為控制系統的非線性參數。

Duffing振子的平穩概率密度函數f(X)有精確解[32]，可以表達為

(19a)

式中，

(19b)

式中:K1/4(·)是修正的Bessel函數1/4階。

該概率密度函數f(X)說明X(t)是非高斯的。根據式(19a)和(19b)，計算得到位移反應的前四階矩，結果及其相關參數列于表3中。從表3中可以看出該結構響應的峰度系數小于3，所以該結構響應為硬化反應。

表3結構參數及其位移響應的前四階矩

Tab.3Thestructuralparametersandthefirstfourmomentsofdisplacement

式(18)的相關參數位移響應的前四階矩μXσXα3Xα4Xω0=1rad/s,ε=0.1,c=0.201.3466702.47436

考慮界限水平x=2σX，利用本文修正方法、He和Zhao方法、傳統高斯模型以及Monte Carlo模擬(10 000個樣本)計算在[0,T]時段內的首次穿越概率。計算結果如圖8所示。從圖8中說明了結構在硬化非高斯結構響應下，本文修正的方法與Monte Carlo模擬結果擬合最好，而其余兩種方法均有不同的差異。因此，本文修正的方法也適用于硬化反應的首次穿越概率計算。

圖8 界限水平為x=2σX的首次穿越概率

4 結 論

本文基于平穩高斯結構響應的Poisson模型修正了一個平穩非高斯結構響應(軟化與硬化)的首次穿越概率的解析方法。在軟化與硬化的非高斯結構響應中，基于結構響應的前四階矩，分別采用Winterstein(1994)模型，Ding和Chen模型的等效高斯分位數及非高斯結構響應的界限水平將平穩非高斯結構響應映射為標準高斯過程。接著，利用考慮初始條件與群超效應的平穩高斯結構響應的Poisson模型計算首次穿越概率。

通過分析Winterstein(1994)軟化模型、Ding和Chen硬化模型的轉換精度及適用情況，驗證了模型運用于平穩非高斯結構響應的首次穿越概率計算的準確性與適用性。通過線性與非線性單自由度系統算例分析，說明了本文修正的方法計算結果比其他已有解析方法能更好地與Monte Carlo模擬結果吻合,同時進一步驗證了本文修正的方法更適用于軟化與硬化平穩非高斯結構響應的首次穿越概率計算。

[1] 李桂青，李秋勝. 工程結構時變可靠度理論及其應用[M]. 北京：科學出版社，2001.

[2] RICE S O. Mathematical analysis of random noise[J]. Bell System Technical Journal, 1944, 23(3): 282-332.

[3] RICE S O. Mathematical analysis of random noise. Part III: Statistical properties of random noise currents[J]. Bell System Technical Journal, 1945, 24(1): 46-156.

[4] COLEMAN J J. Reliability of aircraft structures in resisting chance failure [J]. Operations Research, 1959, 7(5): 639-945.

[5] 劉佩. 基于貝葉斯理論的結構動力可靠度更新方法與分析[J]. 振動與沖擊，2015，34(12)：29-34.

LIU Pei. Structural dynamic reliability updating method based on Bayesian theorem [J]. Journal of Vibration and Shock, 2015, 34(12): 29-34.

[6] VANMARCKE E H. On the distribution of the first passage time for normal stationary process [J]. Journal of Applied Mechanics, 1975, 42(1): 215-220.

[7] 陳建兵，李杰. 非線性隨機地震響應的概率密度演化分析[J]. 武漢理工大學學報，2010，32(9)：6-10.

CHEN Jianbing, LI Jie. Probability density evolution analysis for stochastic seismic response of nonlinear structures [J]. Journal of Wuhan University of Technology, 2010, 32(9):6-10.

[8] CRANDALL S H. First-crossing probabilities of the linear oscillator [J]. Journal of Sound and Vibration, 1970, 12(3): 285-299.

[9] NAESS A. Approximate first-passage and extremes of narrow-band Gaussian and non-Gaussian random vibrations [J]. Journal of Sound and Vibration, 1990, 138(3): 365-380.

[10] BARBATO M, CONTE J P. Structural reliability applications of nonstationary spectral characteristics [J]. Journal of Engineering Mechanics, 2011, 137(5): 371-382.

[11] GHAZIZADEH S, BARBATO M, TUBALDI E. New analytical solution of the first-passage reliability problem for linear oscillators [J]. Journal of Engineering Mechanics, 2012, 138(6):695-706.

[12] GRIGORIU M. Crossings of non-Gaussian translation processes [J]. Journal of Engineering Mechanics, 1984, 110(4): 610-620.

[13] OCHI M. Non-Gaussian random processes in ocean engineering [J]. Probabilistic Engineering Mechanics, 1986, 1(1): 28-39.

[14] JOHNSOM N L, KOTZ S. Continuous univariate distribution-1 [M]. Boston: Houghton Mifflin Company, 1970.

[15] WINTERSTEIN S R. Nonlinear vibration models for extremes and fatigue [J]. Journal of Engineering Mechanics, 1988, 114(10): 1772-1790.

[16] HE J, ZHAO Y G. First passage times of stationary non-Gaussian structural responses [J]. Computers and Structures, 2007, 85(7/8): 431-436.

[17] ZHAO Y G, LU Z H. Fourth-moment standardization for Structural Reliability Assessment [J]. Journal of Structural Engineering, 2007, 133(7): 916-924.

[18] CHOI M, SWEETMAN B. The Hermite moment model for highly skewed response with application to tension leg [J]. Journal of Offshore Mechanics and Arctic Engineering, 2010, 132(2): 1-8.

[19] HUANG M F, LOU W J, CHAN C M, et al. Peak distribution and peak factors of wind-induced pressure processes on tall buildings [J]. Journal of Engineering Mechanics, 2013, 139 (12): 1744-1756.

[20] WINTERSTEIN S R, UDE T C, KLEIVEN G. Springing and Slow-Drift Responses: Predicted Extremes and Fatigue vs. Simulation [C]//BOSS-94. Cambrige: Massachusetts Institute of Technology, 1994: 1-15.

[21] DING J, CHEN X. Moment-based translation model for hardening non-Gaussian response processes [J]. Journal of Engineering Mechanics, 2016, 142(2): 1-7.

[22] DITLEVSEN O. Duration of visit to critical set by Gaussian process [J]. Probabilistic Engineering Mechanics, 1986, 1(2): 82-93.

[23] HE J. Approximate method for estimating extreme value responses of nonlinear stochastic dynamic systems [J]. Journal of Engineering Mechanics, 2015, 141(7): 1-9.

[24] LANGLEY R S. A first passage approximation for normal stationary random processes [J]. Journal of Sound Vibration, 1988, 122(2): 261-275.

[25] WINTERSTEIN S R, KASHEF T. Moment-based load and response models with wind engineering applications [J]. Journal of Solar Energy Engineering, 2000, 122(3): 122-128.

[26] WINTERSTEIN S R, MACKENZIE C A. Extremes of nonlinear vibration: comparing models based on moments, L-moments, and maximum entropy [J]. Journal of Offshore Mechanics Arctic Engineering, 2011, 135(2): 185-195.

[27] FLEISHMAN A L. A method for simulation non-normal distributions [J]. Psychometrika, 1978, 43(4):521-532.

[28] ZHAO Y G, LU Z H. Cubic normal distribution and its significance in structural reliability [J]. Structural Engineering and Mechanics, 2008, 28(3):263-280.

[29] HOHENBICHLER M, RACKWITZ R. Zon-normal dependent vectors in structural safety[J]. Journal of Engineering Mechanics Division, 1981, 107(6): 1227-1238.

[30] NAESS A. The response statistics of nonlinear, second order transformations to Gaussian loads [J]. Journal of Sound and Vibration, 1987, 115(1): 103-129.

[31] BAYER V, BUCHER C. Importance sampling for first passage problems of nonlinear structures [J]. Probabilistic Engineering Mechanics, 1999, 14(1/2): 27-32.

[32] SONG T T, GRIGORIU M. Random vibration of mechanical and structural systems [M]. Englewood Cliffs, NJ: Prentice Hall, 1993.