具有時變功率譜的非高斯隨機過程的數值模擬

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

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

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

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

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

1 具有時變功率譜的隨機過程譜表示

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

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

(1)

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

(2)

非平穩隨機過程的均值為：

(3)

其相關函數為：

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

(4)

當τ=0時，

(5)

(6)

由平穩隨機過程譜表示的數值模擬公式，經過一系列的證明可以得到非平穩隨機過程譜表示的數值模擬公式為[13]：

(7)

或者寫成：

(8)

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

(9)

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

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

(10)

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

(11)

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

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

(12)

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

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

(13)

2 非高斯隨機過程的靜態轉換

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

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

(14)

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

(15)

(16)

(17)

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

對于非高斯隨機過程的第二類數值模擬通常采用無記憶非線性平移方式的靜態轉換

(18)

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

當具有目標功率譜的標準化高斯隨機過程，經過式(18)非線性平移前后，生成的標準化非高斯隨機過程的功率譜必然發生非線性變化。因此，標準化非高斯隨機過程的目標功率譜不能直接作為高斯隨機過程目標功率譜進行模擬。為此，需要建立標準化非高斯隨機過程的目標功率譜與高斯隨機過程目標功率譜的轉化關系。根據相關函數與功率譜之間的相互轉換關系，非高斯隨機過程與高斯隨機過程的目標功率譜之間的非線性關系，可通過目標相關函數之間的非線性關系進行表達。根據式(18)的非線性平移，標準化的非高斯隨機過程的相關函數可表示為

(19)

式中：ψ為具有相關系數ρx的任意兩個時刻非高斯隨機變量的聯合概率密度函數；φ為具有相關系數ρf的任意兩個時刻高斯隨機變量的聯合概率密度函數。

(20)

而相關函數與相關系數的轉化關系為

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

(21)

根據式(19)～式(21)，可建立非高斯隨機過程與高斯隨機過程的相關系數或相關函數之間的轉換關系。

3 具有時變功率譜特征的非高斯風荷載模擬

為了驗證上述方法的有效性，本節將以具有非平穩非高斯特征的下擊暴流脈動風速的模擬作為數值算例。在數值模擬過程中，脈動風速的目標非平穩特征主要表現在功率譜隨時間變化的時變功率譜，該時變功率譜可采用kaimal非均勻調制函數

(22)

調制kaimal譜，

(23)

圖1 下擊暴流時變平均風速Fig.1 Time-varying average velocity of a downburst

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

脈動風速標準化的目標邊緣概率分布函數考慮為對數分布函數，即

(24)

該對數分布函數的參數考慮為。

由非線性平移的靜態轉換可知，在生成具有目標非高斯特征的非平穩隨機過程之前，首先需要模擬潛在的非平穩高斯隨機過程。根據式(18)～式(21)、式(24)，可建立非平穩非高斯隨機過程與潛在的非平穩高斯隨機過程之間的相關系數轉換關系，如3圖所示。

圖3 高斯與非高斯隨機過程相關系數之間的轉換關系Fig.3 The transformation relationship of the correlation coefficients between Gaussian and non-Gaussian random processes

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

圖4 模擬的非平穩高斯隨機過程Fig.5 The simulated non-stationary Gaussian stochastic process

圖5 非線性平移生成的非平穩非高斯脈動風速Fig.5 The generated non-stationary non-Gaussian fluctuating wind velocity through nonlinear translation

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

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

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

圖8 脈動風速的瞬時相關函數與目標相關函數的對比Fig.8 Instantaneous correlation functions of the fluctuating wind velocity with regard to the corresponding targets

圖9 脈動風速的瞬時概率密度函數與目標函數的對比Fig.9 Probability density functions of the fluctuating wind velocity with regard to the corresponding targets

4 結 論

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

[ 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 ] 李璟, 韓大建. 大跨度屋蓋結構非高斯風壓場的一種模擬方法[J]. 工程力學, 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模型的非平穩非高斯隨機過程的數值模擬[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.