李锦华， 李春祥， 邓 莹， 蒋 磊

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

非高斯过程的数值模拟可以分为两类。第一类，根据任意指定的高阶特征统计参数(例如斜度和峰度)和目标功率谱模拟出非高斯过程。第二类，根据任意指定的边缘概率分布函数和目标功率谱模拟出非高斯过程。对于第一类， Seong等[1-4]做了许多的研究工作。他们都是先模拟生成潜在的高斯过程，然后通过建立的合适非线性转换关系间接地将潜在的高斯过程转换成具有目标功率谱、高阶特征统计参数的非高斯过程。线性滤波法是除谐波合成法外用于模拟随机过程的主要方法之一。其中，自回归(Autoregressive, AR)模型和自回归滑动平均(Autoregressive Moving Average, ARMA)模型因其计算量小、速度快，在工程中被广泛用于模拟高斯随机过程[5-7]。李春祥等[8]建立了基于AR和ARMA模型直接模拟单变量非高斯过程的数值算法。本文将这个发展的单变量非高斯过程AR和ARMA模型模拟算法扩展至多变量非高斯过程的数值模拟。

1 基于AR或ARMA模型的多变量非高斯模拟算法

基于AR(p) 模型的多变量非高斯随机过程的模拟公式可表示为

(1)

式中： [Ai]n×n为n行n列的自回归系数矩阵； [L]n×n为n行n列的系数矩阵； [Uk]n×1为n行1列的多变量非高斯过程； [Xk]n×1为n行具有均值为0、方差为1、斜度为SkX、 峰度为KX的非正态分布白噪声序列。模型系数矩阵仅与多变量非高斯随机程过的互相关函数有关。因此，通过模型系数矩阵的确定，可以实现具有相关性的多变量随机过程的模拟。为了进一步模拟具有非高斯特征的随机过程，需要建立模型输入与输出随机变量之间的高阶特征参数的转换关系。

于是，基于AR模型的两变量非高斯随机过程的显式可表示为

(2)

(3)

式中：Un,k为第n变量的非高斯随机过程的第k序列；Xn,k为第n行的非正态分布的白噪声的第k序列；an-n,i为自回归系数矩阵[Ai]n×n的第n行第n列位置上的元素；ln-n为系数矩阵[L]n×n的第n行第n列位置上的元素。

根据式(2)和式(3)的递推关系，可知两变量均是白噪声的线性化表示，则有

(4)

(5)

式中，ψ、γ参数a、l与模型系数、有关。根据式(4)，第一变量的非高斯随机过程U1,k的第二、三阶统计矩为

(6)

(7)

(8)

(9)

因此，第一变量非高斯过程的斜度可表示为

(10)

式中，η1,3为待定的转换系数，可通过数据拟合来确定。

对于第一变量的非高斯随机过程U1,k的第四阶统计矩为

(11)

同理， 根据Xn,k各个元素具有相互独立的性质，有

(12)

因此，第一变量非高斯过程的峰度可表示为

(13)

式中，η1,4，ξ1,4为待定的转换系数，可通过数据拟合来确定。

ELISA法需要的仪器和设备不多，试剂盒需要冷藏保存，接触到的标准品一般浓度较低甚至没有[5]。酶标仪可以携带，一次可以检测大批量的样品，适合于样品的初步筛选[6]。

同理，可知第二变量的非高斯过程的斜度、峰度也可表示为

SkU2=η2,3SkX

(14)

KU2=η2,4KX+ξ2,4

(15)

式中，η2,4，ξ2,4为待定的转换系数，可通过数据拟合来确定。

因此，通过AR(p)模型模拟的多变量非高斯过程的斜度、峰度与输入白噪声的斜度、峰度具有下列转换关系

SkUn=ηn,3SkX,n=1,2,3,…

(16)

KUn=ηn,4KX+ξn,4,n=1,2,3,…

(17)

式中，ηn,3，ηn,4，ξn,4为待定的转换系数，可以通过数据拟合来确定。

多变量非高斯随机过程的ARMA(p,q)模型可表示为

(18)

2 数值验证

为验证算法有效性，考虑变量间互相关进行两点非高斯脉动风压模拟。目标功率谱密度函数

(19)

由上述知，基于AR和ARMA模型的多变量非高斯模拟，其输入非正态分布白噪声斜度和峰度与输出非高斯斜度和峰度分别近似呈线性关系。因此，可通过任意选定几组非正态分布的白噪声，输入特定的AR和ARMA模型产生非高斯样本，估计样本的斜度和峰度，然后通过线性拟合获得待定的线性转换系数ηn,3，ηn,4，ξn,4。图1和图2分别为基于ARMA(28,2)模型，输入非正态分布白噪声的斜度与输出第一和二点非高斯脉动风压斜度的线性关系。通过对数据的线性拟合，知第一点与第二点的斜度线性转换系数分别为η1,3=0.235 2，η2,3=0.156 3。 图3和图4分别为第一和二点非高斯脉动风压峰度与输入非正态分布白噪声峰度的关系，通过数据的线性拟合知η1,4=0.074 48，ξ1,4=2.827；η2,4=0.032 99，ξ2,4=2.872。下面将基于ARMA(28,2)模型，根据拟合出的线性转换关系，进行了两点脉动风压的四种情况：高斯脉动风压Ⅰ、非高斯脉动风压Ⅱ、非高斯脉动风压Ⅲ和非高斯脉动风压Ⅳ的数值模拟分析(见表1)。

图1 基于ARMA(28, 2)模型的第一点脉动风压斜度转换关系的线性拟合Fig.1 Linear fitting of the first point fluctuating wind pressure skewness transformation relationship based on ARMA (28, 2)

图2 基于ARMA(28, 2)模型的第二点脉动风压斜度转换关系的线性拟合Fig.2 Linear fitting of the first second fluctuating wind pressure skewness transformation relationship based on ARMA (28, 2)

图3 基于ARMA(28, 2)模型的第一点脉动风压峰度转换关系的线性拟合Fig.3 Linear fitting of the first point fluctuating wind pressure kurtosis transformation relationship based on ARMA (28, 2)

图4 基于ARMA(28, 2)模型的第二点脉动风压峰度转换关系的线性拟合Fig.4 Linear fitting of the first second fluctuating wind pressure kurtosis transformation relationship based on ARMA (28, 2)

表1 基于ARMA(28, 2)的多变量脉动风压高阶特征值对比

图5为模拟产生的第一点与第二点高斯脉动风压Ⅰ的时程。图6为高斯脉动风压Ⅰ的功率谱与目标谱对比，基本相互吻合。第一点与第二点高斯脉动风压Ⅰ的自、互相关函数与目标自、互相关函数的对比如图7所示，也基本吻合。模拟产生的高斯脉动风压Ⅰ斜度、峰度如表1所示，与目标斜度、峰度也基本吻合。

图5 基于ARMA(28, 2)模型模拟的两点脉动风压ⅠFig.5 Simulated two point fluctuating wind pressure Ⅰ based on ARMA (28, 2)

图6 两点脉动风压Ⅰ功率谱与目标功率谱的对比Fig.6 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅰ

图7 两点脉动风压Ⅰ自、互相关函数与目标自、互相关函数的对比Fig.7 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅰ

图8为模拟产生的第一点与第二点非高斯脉动风压Ⅱ的时程，从图中可观察出其斜度较小。从图9中可观察出，图8的低斜度非高斯脉动风压Ⅱ的功率谱与目标谱较好地吻合，其自、互相关函数与目标自、互相关函数也很好地吻合(见图10)。从表1中还知，图8的低斜度非高斯脉动风压Ⅱ的斜度、峰度与目标斜度、峰度也基本吻合。

图8 基于ARMA(28, 2)模型模拟的两点脉动风压ⅡFig.8 Simulated two point fluctuating wind pressure Ⅱ based on ARMA (28, 2)

图9 两点脉动风压Ⅱ功率谱与目标功率谱的对比Fig.9 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅱ

图10 两点脉动风压Ⅱ自、互相关函数与目标自、互相关函数的对比Fig.10 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅱ

图11为模拟产生的第一和二点中斜度非高斯脉动风压Ⅲ的时程，斜度明显大于图8非高斯脉动风压斜度。中斜度风压Ⅲ的第一和二点时程的功率谱与目标值吻合良好(见图12)。图13(a)、图13(b)分别为第一和二点自相关函数与目标值的对比，图13(c)为两点互相关函数与目标值的对比，均吻合良好。第一和二点中斜度风压Ⅲ的斜度、峰度与目标值也吻合良好(见表1)。

图11 基于ARMA(28, 2)模型模拟的两点脉动风压ⅢFig.11 Simulated two point fluctuating wind pressure Ⅲ based on ARMA (28, 2)

图12 两点脉动风压Ⅲ功率谱与目标功率谱的对比Fig.12 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅲ

图13 两点脉动风压Ⅲ自、互相关函数与目标自、互相关函数的对比Fig.13 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅲ

为说明该方法能模拟高斜度多变量非高斯过程，进行了高斜度两点非高斯脉动风压模拟。图14为模拟产生的第一和二点高斜度非高斯脉动风压Ⅳ的时程。图15和图16分别给出高斜度非高斯脉动风压Ⅳ的第一和二点时程功率谱与目标功率谱，自、互相关函数与目标自、互相关函数，斜度、峰度与目标斜度、峰度的对比，均吻合良好。

图14 基于ARMA(28, 2)模型模拟的两点脉动风压ⅣFig.14 Simulated two point fluctuating wind pressure Ⅳ based on ARMA (28, 2)

图15 两点脉动风压Ⅳ功率谱与目标功率谱的对比Fig.15 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅳ

图16 两点脉动风压Ⅳ自、互相关函数与目标自、互相关函数的对比Fig.16 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅳ

3 结 论

在基于AR和ARMA的单变量非高斯随机过程模拟算法的基础上，考虑了多变量随机过程之间的相关性，将其扩展至多变量非高斯随机过程的数值模拟。数值分析表明：AR和ARMA模型的模拟算法能有效地模拟出低斜度、中斜度和高斜度多变量非高斯随机过程。

[ 1 ] SEONG S H, PETERKA J A. Computer simulation of non-Gaussian multiple wind press time series[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 72(1): 95-105.

[ 2 ] SEONG S H, PETERKA J A. Digital generation of surface-pressure fluctuations with spiky features[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1998, 73(2): 181-192.

[ 3 ] SURESH KUMAR K, STATHOPOULOS T. Computer simulation of fluctuating wind pressures on low building roofs[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 69/70/71(97): 485-495.

[ 4 ] GURLEY K, KAREEM A. Analysis, interpretation, modeling and simulation of unsteady wind and pressure data[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 69/70/71(97): 657-669.

[ 5 ] SAMARAS E, SHINOZUKA M, TSURUI A. ARMA representation of random process[J]. Journal of Engineering Mechanics, 1985, 111(3): 449-461.

[ 6 ] 李春祥, 谈雅雅, 李锦华. 基于ARMA模型模拟高架桥的脉动风速时程[J]. 振动与冲击, 2009, 28(6): 46-51.

LI Chunxiang, TAN Yaya, LI Jinhua. Simulation of fluctuating wind speed time series applied on overpass bridges with resorting to ARMA model[J]. Journal of Vibration and Shock, 2009, 28(6): 46-51.

[ 7 ] 李锦华, 李春祥. 土木工程随机风场数值模拟研究的进展[J]. 振动与冲击, 2008, 27(9): 116-125.

LI Jinhua, LI Chunxiang. Development of numerical simulations for stochastic wind fields in civil Engineering[J]. Journal of Vibration and Shock, 2008, 27(9): 116-125.

[ 8 ] LI Jinhua, LI Chunxiang. Simulation of non-Gaussian stochastic process with target power spectral density and lower-order moments[J]. Journal of Engineering Mechanics，2012, 138(5): 391-404.