李錦華， 李春祥， 鄧 瑩， 蔣 磊

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

非高斯過程的數(shù)值模擬可以分為兩類。第一類，根據(jù)任意指定的高階特征統(tǒng)計參數(shù)(例如斜度和峰度)和目標(biāo)功率譜模擬出非高斯過程。第二類，根據(jù)任意指定的邊緣概率分布函數(shù)和目標(biāo)功率譜模擬出非高斯過程。對于第一類， Seong等[1-4]做了許多的研究工作。他們都是先模擬生成潛在的高斯過程，然后通過建立的合適非線性轉(zhuǎn)換關(guān)系間接地將潛在的高斯過程轉(zhuǎn)換成具有目標(biāo)功率譜、高階特征統(tǒng)計參數(shù)的非高斯過程。線性濾波法是除諧波合成法外用于模擬隨機過程的主要方法之一。其中，自回歸(Autoregressive, AR)模型和自回歸滑動平均(Autoregressive Moving Average, ARMA)模型因其計算量小、速度快，在工程中被廣泛用于模擬高斯隨機過程[5-7]。李春祥等[8]建立了基于AR和ARMA模型直接模擬單變量非高斯過程的數(shù)值算法。本文將這個發(fā)展的單變量非高斯過程AR和ARMA模型模擬算法擴展至多變量非高斯過程的數(shù)值模擬。

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

基于AR(p) 模型的多變量非高斯隨機過程的模擬公式可表示為

(1)

式中： [Ai]n×n為n行n列的自回歸系數(shù)矩陣； [L]n×n為n行n列的系數(shù)矩陣； [Uk]n×1為n行1列的多變量非高斯過程； [Xk]n×1為n行具有均值為0、方差為1、斜度為SkX、 峰度為KX的非正態(tài)分布白噪聲序列。模型系數(shù)矩陣僅與多變量非高斯隨機程過的互相關(guān)函數(shù)有關(guān)。因此，通過模型系數(shù)矩陣的確定，可以實現(xiàn)具有相關(guān)性的多變量隨機過程的模擬。為了進一步模擬具有非高斯特征的隨機過程，需要建立模型輸入與輸出隨機變量之間的高階特征參數(shù)的轉(zhuǎn)換關(guān)系。

于是，基于AR模型的兩變量非高斯隨機過程的顯式可表示為

(2)

(3)

式中：Un,k為第n變量的非高斯隨機過程的第k序列；Xn,k為第n行的非正態(tài)分布的白噪聲的第k序列；an-n,i為自回歸系數(shù)矩陣[Ai]n×n的第n行第n列位置上的元素；ln-n為系數(shù)矩陣[L]n×n的第n行第n列位置上的元素。

根據(jù)式(2)和式(3)的遞推關(guān)系，可知兩變量均是白噪聲的線性化表示，則有

(4)

(5)

式中，ψ、γ參數(shù)a、l與模型系數(shù)、有關(guān)。根據(jù)式(4)，第一變量的非高斯隨機過程U1,k的第二、三階統(tǒng)計矩為

(6)

(7)

(8)

(9)

因此，第一變量非高斯過程的斜度可表示為

(10)

式中，η1,3為待定的轉(zhuǎn)換系數(shù)，可通過數(shù)據(jù)擬合來確定。

對于第一變量的非高斯隨機過程U1,k的第四階統(tǒng)計矩為

(11)

同理， 根據(jù)Xn,k各個元素具有相互獨立的性質(zhì)，有

(12)

因此，第一變量非高斯過程的峰度可表示為

(13)

式中，η1,4，ξ1,4為待定的轉(zhuǎn)換系數(shù)，可通過數(shù)據(jù)擬合來確定。

ELISA法需要的儀器和設(shè)備不多，試劑盒需要冷藏保存，接觸到的標(biāo)準(zhǔn)品一般濃度較低甚至沒有[5]。酶標(biāo)儀可以攜帶，一次可以檢測大批量的樣品，適合于樣品的初步篩選[6]。

同理，可知第二變量的非高斯過程的斜度、峰度也可表示為

SkU2=η2,3SkX

(14)

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

(15)

式中，η2,4，ξ2,4為待定的轉(zhuǎn)換系數(shù)，可通過數(shù)據(jù)擬合來確定。

因此，通過AR(p)模型模擬的多變量非高斯過程的斜度、峰度與輸入白噪聲的斜度、峰度具有下列轉(zhuǎn)換關(guān)系

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

(16)

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

(17)

式中，ηn,3，ηn,4，ξn,4為待定的轉(zhuǎn)換系數(shù)，可以通過數(shù)據(jù)擬合來確定。

多變量非高斯隨機過程的ARMA(p,q)模型可表示為

(18)

2 數(shù)值驗證

為驗證算法有效性，考慮變量間互相關(guān)進行兩點非高斯脈動風(fēng)壓模擬。目標(biāo)功率譜密度函數(shù)

(19)

由上述知，基于AR和ARMA模型的多變量非高斯模擬，其輸入非正態(tài)分布白噪聲斜度和峰度與輸出非高斯斜度和峰度分別近似呈線性關(guān)系。因此，可通過任意選定幾組非正態(tài)分布的白噪聲，輸入特定的AR和ARMA模型產(chǎn)生非高斯樣本，估計樣本的斜度和峰度，然后通過線性擬合獲得待定的線性轉(zhuǎn)換系數(shù)ηn,3，ηn,4，ξn,4。圖1和圖2分別為基于ARMA(28,2)模型，輸入非正態(tài)分布白噪聲的斜度與輸出第一和二點非高斯脈動風(fēng)壓斜度的線性關(guān)系。通過對數(shù)據(jù)的線性擬合，知第一點與第二點的斜度線性轉(zhuǎn)換系數(shù)分別為η1,3=0.235 2，η2,3=0.156 3。 圖3和圖4分別為第一和二點非高斯脈動風(fēng)壓峰度與輸入非正態(tài)分布白噪聲峰度的關(guān)系，通過數(shù)據(jù)的線性擬合知η1,4=0.074 48，ξ1,4=2.827；η2,4=0.032 99，ξ2,4=2.872。下面將基于ARMA(28,2)模型，根據(jù)擬合出的線性轉(zhuǎn)換關(guān)系，進行了兩點脈動風(fēng)壓的四種情況：高斯脈動風(fēng)壓Ⅰ、非高斯脈動風(fēng)壓Ⅱ、非高斯脈動風(fēng)壓Ⅲ和非高斯脈動風(fēng)壓Ⅳ的數(shù)值模擬分析(見表1)。

圖1 基于ARMA(28, 2)模型的第一點脈動風(fēng)壓斜度轉(zhuǎn)換關(guān)系的線性擬合Fig.1 Linear fitting of the first point fluctuating wind pressure skewness transformation relationship based on ARMA (28, 2)

圖2 基于ARMA(28, 2)模型的第二點脈動風(fēng)壓斜度轉(zhuǎn)換關(guān)系的線性擬合Fig.2 Linear fitting of the first second fluctuating wind pressure skewness transformation relationship based on ARMA (28, 2)

圖3 基于ARMA(28, 2)模型的第一點脈動風(fēng)壓峰度轉(zhuǎn)換關(guān)系的線性擬合Fig.3 Linear fitting of the first point fluctuating wind pressure kurtosis transformation relationship based on ARMA (28, 2)

圖4 基于ARMA(28, 2)模型的第二點脈動風(fēng)壓峰度轉(zhuǎn)換關(guān)系的線性擬合Fig.4 Linear fitting of the first second fluctuating wind pressure kurtosis transformation relationship based on ARMA (28, 2)

表1 基于ARMA(28, 2)的多變量脈動風(fēng)壓高階特征值對比

圖5為模擬產(chǎn)生的第一點與第二點高斯脈動風(fēng)壓Ⅰ的時程。圖6為高斯脈動風(fēng)壓Ⅰ的功率譜與目標(biāo)譜對比，基本相互吻合。第一點與第二點高斯脈動風(fēng)壓Ⅰ的自、互相關(guān)函數(shù)與目標(biāo)自、互相關(guān)函數(shù)的對比如圖7所示，也基本吻合。模擬產(chǎn)生的高斯脈動風(fēng)壓Ⅰ斜度、峰度如表1所示，與目標(biāo)斜度、峰度也基本吻合。

圖5 基于ARMA(28, 2)模型模擬的兩點脈動風(fēng)壓ⅠFig.5 Simulated two point fluctuating wind pressure Ⅰ based on ARMA (28, 2)

圖6 兩點脈動風(fēng)壓Ⅰ功率譜與目標(biāo)功率譜的對比Fig.6 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅰ

圖7 兩點脈動風(fēng)壓Ⅰ自、互相關(guān)函數(shù)與目標(biāo)自、互相關(guān)函數(shù)的對比Fig.7 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅰ

圖8為模擬產(chǎn)生的第一點與第二點非高斯脈動風(fēng)壓Ⅱ的時程，從圖中可觀察出其斜度較小。從圖9中可觀察出，圖8的低斜度非高斯脈動風(fēng)壓Ⅱ的功率譜與目標(biāo)譜較好地吻合，其自、互相關(guān)函數(shù)與目標(biāo)自、互相關(guān)函數(shù)也很好地吻合(見圖10)。從表1中還知，圖8的低斜度非高斯脈動風(fēng)壓Ⅱ的斜度、峰度與目標(biāo)斜度、峰度也基本吻合。

圖8 基于ARMA(28, 2)模型模擬的兩點脈動風(fēng)壓ⅡFig.8 Simulated two point fluctuating wind pressure Ⅱ based on ARMA (28, 2)

圖9 兩點脈動風(fēng)壓Ⅱ功率譜與目標(biāo)功率譜的對比Fig.9 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅱ

圖10 兩點脈動風(fēng)壓Ⅱ自、互相關(guān)函數(shù)與目標(biāo)自、互相關(guān)函數(shù)的對比Fig.10 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅱ

圖11為模擬產(chǎn)生的第一和二點中斜度非高斯脈動風(fēng)壓Ⅲ的時程，斜度明顯大于圖8非高斯脈動風(fēng)壓斜度。中斜度風(fēng)壓Ⅲ的第一和二點時程的功率譜與目標(biāo)值吻合良好(見圖12)。圖13(a)、圖13(b)分別為第一和二點自相關(guān)函數(shù)與目標(biāo)值的對比，圖13(c)為兩點互相關(guān)函數(shù)與目標(biāo)值的對比，均吻合良好。第一和二點中斜度風(fēng)壓Ⅲ的斜度、峰度與目標(biāo)值也吻合良好(見表1)。

圖11 基于ARMA(28, 2)模型模擬的兩點脈動風(fēng)壓ⅢFig.11 Simulated two point fluctuating wind pressure Ⅲ based on ARMA (28, 2)

圖12 兩點脈動風(fēng)壓Ⅲ功率譜與目標(biāo)功率譜的對比Fig.12 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅲ

圖13 兩點脈動風(fēng)壓Ⅲ自、互相關(guān)函數(shù)與目標(biāo)自、互相關(guān)函數(shù)的對比Fig.13 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅲ

為說明該方法能模擬高斜度多變量非高斯過程，進行了高斜度兩點非高斯脈動風(fēng)壓模擬。圖14為模擬產(chǎn)生的第一和二點高斜度非高斯脈動風(fēng)壓Ⅳ的時程。圖15和圖16分別給出高斜度非高斯脈動風(fēng)壓Ⅳ的第一和二點時程功率譜與目標(biāo)功率譜，自、互相關(guān)函數(shù)與目標(biāo)自、互相關(guān)函數(shù)，斜度、峰度與目標(biāo)斜度、峰度的對比，均吻合良好。

圖14 基于ARMA(28, 2)模型模擬的兩點脈動風(fēng)壓ⅣFig.14 Simulated two point fluctuating wind pressure Ⅳ based on ARMA (28, 2)

圖15 兩點脈動風(fēng)壓Ⅳ功率譜與目標(biāo)功率譜的對比Fig.15 Power spectral density and prescribed value of two point fluctuating wind pressure Ⅳ

圖16 兩點脈動風(fēng)壓Ⅳ自、互相關(guān)函數(shù)與目標(biāo)自、互相關(guān)函數(shù)的對比Fig.16 Comparison of correlation functions and their prescribed values for two point fluctuating wind pressure Ⅳ

3 結(jié) 論

在基于AR和ARMA的單變量非高斯隨機過程模擬算法的基礎(chǔ)上，考慮了多變量隨機過程之間的相關(guān)性，將其擴展至多變量非高斯隨機過程的數(shù)值模擬。數(shù)值分析表明：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模型模擬高架橋的脈動風(fēng)速時程[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 ] 李錦華, 李春祥. 土木工程隨機風(fēng)場數(shù)值模擬研究的進展[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.