具有时变功率谱的非高斯随机过程的数值模拟

李锦华， 李建丰， 陈水生， 余维光， 李春祥

(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.