具有時變功率譜的非高斯隨機(jī)過程的數(shù)值模擬

李錦華， 李建豐， 陳水生， 余維光， 李春祥

(1.華東交通大學(xué) 土木建筑學(xué)院，南昌 330013； 2.上海大學(xué) 土木工程系，上海 200444)

對于隨機(jī)過程的模擬，蒙特卡洛隨機(jī)模擬技術(shù)使用最廣泛，能夠模擬產(chǎn)生具有目標(biāo)特征的隨機(jī)過程，包括一維或多維、單變量或多變量、平穩(wěn)或非平穩(wěn)、高斯或非高斯的隨機(jī)過程。為了使相關(guān)領(lǐng)域的研究能夠更符合實際情況，非高斯隨機(jī)過程的數(shù)值模擬越來越受到關(guān)注，特別是非平穩(wěn)非高斯隨機(jī)過程。目前，非高斯隨機(jī)過程的數(shù)值模擬可以分為兩類[1]：①根據(jù)指定的特征統(tǒng)計參數(shù)(例如均值、方差、偏度與峰度)和目標(biāo)功率譜密度(PSD)函數(shù)模擬產(chǎn)生非高斯隨機(jī)過程；②根據(jù)指定的邊緣概率密度函數(shù)(PDF)和目標(biāo)PSD 函數(shù)模擬產(chǎn)生非高斯隨機(jī)過程。

對于“①”非高斯隨機(jī)過程的模擬，Seong等[2]采用指數(shù)峰值模型進(jìn)行了單變量和多變量平穩(wěn)非高斯風(fēng)壓時程的模擬。Kumar等[3]基于FFT技術(shù)，采用了參數(shù)較少的指數(shù)峰值模型模擬了一維單變量平穩(wěn)非高斯風(fēng)壓時程，并用于大跨低矮屋蓋的風(fēng)振分析。但這類方法需要對峰值模型參數(shù)進(jìn)行不斷地優(yōu)化。Gurley等[4]提出了新的靜態(tài)轉(zhuǎn)換法，但該方法生成的單樣本平穩(wěn)非高斯過程的偏度和峰度與指定的偏度和峰度并不是十分吻合，需要對多樣本的偏度和峰度分別求均值才能較好地與目標(biāo)值吻合。在國內(nèi)，李璟等[5]采用三次多項式表達(dá)了非高斯隨機(jī)過程和潛在的高斯隨機(jī)過程之間的轉(zhuǎn)換關(guān)系，進(jìn)行了平穩(wěn)非高斯風(fēng)壓的模擬。Li等[6]通過建立AR和ARMA模型輸入與輸出統(tǒng)計特征之間的關(guān)系，建立了基于AR和ARMA模型的平穩(wěn)非高斯隨機(jī)過程模擬，并通過平穩(wěn)非高斯脈動風(fēng)荷載的模擬驗證了模擬的有效性。

根據(jù)指定的邊緣概率密度函數(shù)(PDF)和目標(biāo)PSD函數(shù)模擬產(chǎn)生非高斯隨機(jī)過程，即“②”非高斯隨機(jī)過程的模擬，主要方法有：采用反復(fù)迭代的譜方法[7]，但該方法為了使得模擬計算獲得的邊緣PDF 、PSD 或相關(guān)函數(shù)較好地吻合目標(biāo)函數(shù)，需要反復(fù)迭代，主要是進(jìn)行經(jīng)驗性的不斷調(diào)整，缺乏收斂的必然性依據(jù)[8]；比較有效的非迭代法是Grigoriu[9]提出基于非線性平移的非高斯隨機(jī)過程非迭代模擬算法。目前，國內(nèi)外主要研究是針對平穩(wěn)非高斯隨機(jī)過程進(jìn)行模擬。對于非平穩(wěn)非高斯隨機(jī)過程的模擬，不僅需要模擬具有相關(guān)函數(shù)或功率譜密度隨時間變化的特征，而且需要呈現(xiàn)出非高斯分布的特征。為此，李錦華等[10]進(jìn)一步探討了基于時變AR模型的第一類非平穩(wěn)非高斯隨機(jī)過程的有效模擬。本文將基于Grigoriu的非迭代方法進(jìn)一步進(jìn)行具有時變功率譜的“②”非高斯隨機(jī)過程的有效模擬。

1 具有時變功率譜的隨機(jī)過程譜表示

在平穩(wěn)高斯隨機(jī)過程的數(shù)值模擬中，較為成熟的方法主要有譜表示法和線性濾波法。譜表示法的模擬精度較高，但模擬效率較低，而線性濾波法由于計算量小、速度快，被廣泛用于工程中的隨機(jī)過程模擬。李錦華等基于線形濾波法中AR、ARMA模型實現(xiàn)了非高斯平穩(wěn)隨機(jī)過程的有效模擬，以及通過AR模型進(jìn)一步考慮時變特征，建立非平穩(wěn)高斯隨機(jī)過程的時變AR模型，并將此模型擴(kuò)展到第一類非平穩(wěn)非高斯隨機(jī)過程的模擬。此外，Li等[11]還通過模擬精度較高的譜表示法進(jìn)行了非平穩(wěn)脈動風(fēng)速的模擬。由于線性濾波法精度較低，因此本文針對具有目標(biāo)時變功率譜和目標(biāo)概率密度分布特征的非平穩(wěn)非高斯隨機(jī)過程模擬，將首先采用隨機(jī)過程的譜表示法來實現(xiàn)具有目標(biāo)時變功率譜特征的非平穩(wěn)高斯隨機(jī)過程的模擬。

根據(jù)Priestley非平穩(wěn)隨機(jī)過程進(jìn)化譜表示理論[12]，f(t)是一單變量，一維均值為零的非平穩(wěn)隨機(jī)過程，可用如下積分形式表示：

(1)

式中：A(t,ω)是非均勻調(diào)制函數(shù)；{Z(ω)}是增量正交的譜過程，且滿足

(2)

非平穩(wěn)隨機(jī)過程的均值為：

(3)

其相關(guān)函數(shù)為：

Rff(t,τ)=E[f*(t)f(t+τ)]=

(4)

當(dāng)τ=0時，

(5)

(6)

由平穩(wěn)隨機(jī)過程譜表示的數(shù)值模擬公式，經(jīng)過一系列的證明可以得到非平穩(wěn)隨機(jī)過程譜表示的數(shù)值模擬公式為[13]：

(7)

或者寫成：

(8)

因此，對于一維n變量的零均值非平穩(wěn)隨機(jī)過程{f(t)}=[f1(t)f2(t)…fn(t)]T，其譜密度矩陣是時間t和圓周頻率ω的函數(shù)，即：

(9)

其中，Γ(ω,t)為不同變量之間的相干函數(shù)。在每一時刻t時，譜密度矩陣S(ω,t)進(jìn)行Cholesky分解

S(ω,t)=H(ω,t)HT*(ω,t)

(10)

式中:H(ω,t)為下三角矩陣，HT*(ω,t)是其復(fù)共軛轉(zhuǎn)置矩陣。

(11)

考慮一般情況下G0(ω,t)為復(fù)數(shù)矩陣，因此H(ω,t)通常也是復(fù)數(shù)矩陣，其對角線元素為非負(fù)實數(shù)，非對角線元素為復(fù)數(shù)。矩陣H(ω,t)中的元素可以表示為：

Hjk(ω,t)=|Hjk(ω,t)|ejθkj(ω,t)

(12)

式中：θjk(ω,t)=tan-1{Im[Hjk(ω,t)]/Re[Hjk(ω,t)]}為Hjk(ω,t)的幅角。

因此，一維n變量的零均值非平穩(wěn)隨機(jī)過程可以通過下列譜表示來模擬

(13)

2 非高斯隨機(jī)過程的靜態(tài)轉(zhuǎn)換

目前，非高斯隨機(jī)過程的模擬主要是通過潛在高斯隨機(jī)過程的靜態(tài)轉(zhuǎn)換。對于非高斯隨機(jī)過程y(t)的第一類數(shù)值模擬可通過潛在高斯隨機(jī)過程f(t)進(jìn)行下列方式的靜態(tài)轉(zhuǎn)換

α[f(t)+h3(f2(t)-1)+h4(f3(t)-3f(t))]

(14)

式中：f(t)為零均值，方差為1的標(biāo)準(zhǔn)化高斯隨機(jī)過程；x(t)為零均值，方差為1的標(biāo)準(zhǔn)化非高斯隨機(jī)過程；my，σy為分別為原非高斯隨機(jī)過程y(t)的時變均值和時變均方差；α,h3,h4系數(shù)分別為

(15)

(16)

(17)

式中：γ3，γ4分別為非高斯隨機(jī)過程y(t)的目標(biāo)斜度和峰度。

對于非高斯隨機(jī)過程的第二類數(shù)值模擬通常采用無記憶非線性平移方式的靜態(tài)轉(zhuǎn)換

(18)

式中：φ為潛在高斯隨機(jī)過程f(t)的邊緣概率分布函數(shù)；F為非高斯隨機(jī)過程x(t)的邊緣概率分布函數(shù)。

當(dāng)具有目標(biāo)功率譜的標(biāo)準(zhǔn)化高斯隨機(jī)過程，經(jīng)過式(18)非線性平移前后，生成的標(biāo)準(zhǔn)化非高斯隨機(jī)過程的功率譜必然發(fā)生非線性變化。因此，標(biāo)準(zhǔn)化非高斯隨機(jī)過程的目標(biāo)功率譜不能直接作為高斯隨機(jī)過程目標(biāo)功率譜進(jìn)行模擬。為此，需要建立標(biāo)準(zhǔn)化非高斯隨機(jī)過程的目標(biāo)功率譜與高斯隨機(jī)過程目標(biāo)功率譜的轉(zhuǎn)化關(guān)系。根據(jù)相關(guān)函數(shù)與功率譜之間的相互轉(zhuǎn)換關(guān)系，非高斯隨機(jī)過程與高斯隨機(jī)過程的目標(biāo)功率譜之間的非線性關(guān)系，可通過目標(biāo)相關(guān)函數(shù)之間的非線性關(guān)系進(jìn)行表達(dá)。根據(jù)式(18)的非線性平移，標(biāo)準(zhǔn)化的非高斯隨機(jī)過程的相關(guān)函數(shù)可表示為

(19)

式中：ψ為具有相關(guān)系數(shù)ρx的任意兩個時刻非高斯隨機(jī)變量的聯(lián)合概率密度函數(shù)；φ為具有相關(guān)系數(shù)ρf的任意兩個時刻高斯隨機(jī)變量的聯(lián)合概率密度函數(shù)。

(20)

而相關(guān)函數(shù)與相關(guān)系數(shù)的轉(zhuǎn)化關(guān)系為

R(t,τ)=ρ(t,τ)·σ(t)2

(21)

根據(jù)式(19)～式(21)，可建立非高斯隨機(jī)過程與高斯隨機(jī)過程的相關(guān)系數(shù)或相關(guān)函數(shù)之間的轉(zhuǎn)換關(guān)系。

3 具有時變功率譜特征的非高斯風(fēng)荷載模擬

為了驗證上述方法的有效性，本節(jié)將以具有非平穩(wěn)非高斯特征的下?lián)舯┝髅}動風(fēng)速的模擬作為數(shù)值算例。在數(shù)值模擬過程中，脈動風(fēng)速的目標(biāo)非平穩(wěn)特征主要表現(xiàn)在功率譜隨時間變化的時變功率譜，該時變功率譜可采用kaimal非均勻調(diào)制函數(shù)

(22)

調(diào)制kaimal譜，

(23)

圖1 下?lián)舯┝鲿r變平均風(fēng)速Fig.1 Time-varying average velocity of a downburst

圖2 調(diào)制后的kaimal時變譜Fig.2 The modulated Kaimal time-varying power spectrum

脈動風(fēng)速標(biāo)準(zhǔn)化的目標(biāo)邊緣概率分布函數(shù)考慮為對數(shù)分布函數(shù)，即

(24)

該對數(shù)分布函數(shù)的參數(shù)考慮為。

由非線性平移的靜態(tài)轉(zhuǎn)換可知，在生成具有目標(biāo)非高斯特征的非平穩(wěn)隨機(jī)過程之前，首先需要模擬潛在的非平穩(wěn)高斯隨機(jī)過程。根據(jù)式(18)～式(21)、式(24)，可建立非平穩(wěn)非高斯隨機(jī)過程與潛在的非平穩(wěn)高斯隨機(jī)過程之間的相關(guān)系數(shù)轉(zhuǎn)換關(guān)系，如3圖所示。

圖3 高斯與非高斯隨機(jī)過程相關(guān)系數(shù)之間的轉(zhuǎn)換關(guān)系Fig.3 The transformation relationship of the correlation coefficients between Gaussian and non-Gaussian random processes

基于譜表示法生成潛在的非平穩(wěn)高斯隨機(jī)過程，如圖4所示，并將其經(jīng)過非線性平移變換生成的非平穩(wěn)非高斯脈動風(fēng)速，如圖5所示。對多組脈動風(fēng)速樣本進(jìn)行功率譜估計，其估計譜(如圖6所示)明顯具有如圖2所示的目標(biāo)時變功率譜的時變特征。任意時刻估計譜與目標(biāo)譜的對比，如圖7所示，在該圖中四個任意時刻的估計譜均與目標(biāo)譜較好的吻合，且相應(yīng)的相關(guān)函數(shù)也與目標(biāo)相吻合，如圖8所示，說明模擬的脈動風(fēng)速具有目標(biāo)特征的非平穩(wěn)特性。

圖4 模擬的非平穩(wěn)高斯隨機(jī)過程Fig.5 The simulated non-stationary Gaussian stochastic process

圖5 非線性平移生成的非平穩(wěn)非高斯脈動風(fēng)速Fig.5 The generated non-stationary non-Gaussian fluctuating wind velocity through nonlinear translation

圖6 脈動風(fēng)速樣本時變譜Fig.6 The time-varying power spectrum of the fluctuating wind velocity

為了進(jìn)一步說明模擬的有效性，通過模擬生成的多組脈動風(fēng)速樣本進(jìn)行了脈動風(fēng)速樣本的瞬時概率密度函數(shù)與目標(biāo)函數(shù)的對比如圖9所示。從圖中可以看出，任意時刻的概率密度函數(shù)均不相同，說明了概率密度函數(shù)具有時變性，這是由于目標(biāo)功率譜具有時變性造成方差隨時間變化的原因。從對比中可以觀察到，任意時刻脈動風(fēng)速樣本的概率密度函數(shù)均與目標(biāo)對數(shù)分布函數(shù)相互吻合。因此，模擬的脈動風(fēng)速樣本不僅具有目標(biāo)非平穩(wěn)特征而且還具有目標(biāo)非高斯特征，說明了非平穩(wěn)非高斯隨機(jī)過程模擬方法的有效性。值得指出的是：本文直接采用了精度較好的譜表示模擬生成單變量非平穩(wěn)高斯樣本，并沒有考慮計算效率；而對于譜表示模擬多變量的非平穩(wěn)高斯樣本則需進(jìn)一步考慮計算效率。

圖7 脈動風(fēng)速瞬時估計譜與目標(biāo)譜的對比Fig.7 Instantaneous power spectrums of the fluctuating wind velocity with regard to the corresponding targets

圖8 脈動風(fēng)速的瞬時相關(guān)函數(shù)與目標(biāo)相關(guān)函數(shù)的對比Fig.8 Instantaneous correlation functions of the fluctuating wind velocity with regard to the corresponding targets

圖9 脈動風(fēng)速的瞬時概率密度函數(shù)與目標(biāo)函數(shù)的對比Fig.9 Probability density functions of the fluctuating wind velocity with regard to the corresponding targets

4 結(jié) 論

為了有效地模擬具有目標(biāo)時變功率譜和目標(biāo)概率密度函數(shù)的非平穩(wěn)非高斯隨機(jī)過程，本文提出了基于目標(biāo)時變功率譜和目標(biāo)非高斯概率密度函數(shù)，通過建立非高斯與高斯隨機(jī)過程之間相互轉(zhuǎn)換的非線性平移關(guān)系，以及非線性平移前后高斯與非高斯隨機(jī)過程的功率譜或相關(guān)函數(shù)的轉(zhuǎn)換關(guān)系，將非平穩(wěn)非高斯隨機(jī)過程轉(zhuǎn)化為非平穩(wěn)高斯隨機(jī)過程的模擬。而非平穩(wěn)高斯隨機(jī)過程可通過譜表示進(jìn)行有效的模擬。為了驗證該方法的有效性，文中進(jìn)行了具有目標(biāo)非平穩(wěn)非高斯特征的脈動風(fēng)速模擬。模擬結(jié)果表明：模擬生成的脈動風(fēng)速樣本的功率譜具有時變特征，且瞬時功率譜和相關(guān)函數(shù)均與目標(biāo)相吻合，說明模擬的脈動風(fēng)速具有目標(biāo)特征的非平穩(wěn)特性；任意時刻脈動風(fēng)速樣本的概率密度函數(shù)與目標(biāo)函數(shù)的相互吻合，說明模擬的脈動風(fēng)速具有目標(biāo)特征的非高斯特性。因此，模擬的脈動風(fēng)速樣本不僅具有目標(biāo)非平穩(wěn)特征而且還具有目標(biāo)非高斯特征，說明了該非平穩(wěn)非高斯隨機(jī)過程模擬方法的有效性。

[ 1 ] DEODATIS G, MICALETTI R C. Simulation of highly skewed non-gaussian stochastic processes [J]. J. Engrg. Mech., ASCE, 2001, 127(12): 1284-1295.

[ 2 ] SEONG S H, PETERKA J A. Digital generation of surface-pressure fluctuations with spiky features [J]. J. Wind Eng. Ind. Aerodyn., 1998, 73(2): 181-192.

[ 3 ] KUMAR K S, STATHOPOULOS T. Synthesis of non-gaussian wind pressure time series on low building roofs [J]. Eng. Struct., 1999, 21: 1086-1100.

[ 4 ] GURLEY K, KAREEM A. Analysis interpretation modeling and simulation of unsteady wind and pressure data [J]. Wind Eng. Ind. Aerodyn., 1997, 69/70/71: 657-669.

[ 5 ] 李璟, 韓大建. 大跨度屋蓋結(jié)構(gòu)非高斯風(fēng)壓場的一種模擬方法[J]. 工程力學(xué), 2009, 26(5): 80-87.

LI Jing, HAN Dajian. A simulation method for non-gaussian wind pressure field of large-span roof structures [J]. Engineering Mechanics, 2009, 26(5): 80-87.

[ 6 ] LI J H, LI C X. Simulation of non-gaussian stochastic process with target power spectral density and lower-order moments [J]. Journal of Engineering Mechanics ASCE, 2012, 138(5): 391-404.

[ 7 ] DEODATIS G, POPESCU R, PREVOST J H. Simulation of homogenous non-gaussian stochastic vector fields [J]. Prob. Eng. Mech, 1998, 13(1): 1-13.

[ 8 ] PUIG B, AKIAN J L. Non-gaussian simulation using hermite polynomials expansion and maximum entropy principle [J]. Prob. Eng. Mech., 2004, 19(4): 293-305.

[ 9 ] GRIGORIU M. Simulation of stationary non-gaussian translation processes [J]. J. Eng. Mech., ASCE, 1998, 124 (2): 121-126.

[10] 李錦華, 陳水生, 吳春鵬, 等. 基于時變AR模型的非平穩(wěn)非高斯隨機(jī)過程的數(shù)值模擬[J]. 振動與沖擊, 2015,34(17): 142-146.

LI Jinhua, CHEN Shuisheng WU Chunpeng, et al. Simulation of non-stationary non-Gaussian stochastic process based on time-varying AR model[J]. Journal of Vibration and Shock, 2015, 34(17): 142-146.

[11] LI J H, LI C X, HE L, et al. Extended modulating functions for simulation of wind velocities with weak and strong nonstationarity[J]. Renewable Energy,2015,83: 384-397.

[12] PRIESTLEY M B. Power spectral analysis of non-stationary random processes [J]. Journal of Sound and Vibration, 1967, 6: 86-97.

[13] LIANG J, CHAUDHURI S R, SHINOZUKA M. Simulation of non-stationary stochastic processes by spectral representation [J]. Journal of Engineering Mechanics ASCE, 2007, 133 (6): 616-627.