蔣磊,李春祥,鄧瑩

(上海大學(xué)土木工程系,上海200444)

基于LPZ譜分析法的非高斯脈動(dòng)風(fēng)壓模擬

蔣磊,李春祥,鄧瑩

(上海大學(xué)土木工程系,上海200444)

提出一種基于線性預(yù)測和Z轉(zhuǎn)換(linear prediction and Z-transform,LPZ)譜分析法的非高斯脈動(dòng)風(fēng)壓模擬算法.運(yùn)用Johnson變換系統(tǒng)實(shí)現(xiàn)高斯隨機(jī)過程到非高斯白噪聲的轉(zhuǎn)換;再使用LPZ對(duì)其進(jìn)行數(shù)字濾波,進(jìn)而得到所需的非高斯脈動(dòng)風(fēng)壓.采用基于LPZ譜分析法對(duì)單變量非高斯隨機(jī)信號(hào)和非高斯脈動(dòng)風(fēng)壓進(jìn)行了數(shù)值模擬.通過對(duì)模擬的非高斯信號(hào)和脈動(dòng)風(fēng)壓的統(tǒng)計(jì)參數(shù)(峰度和偏度)與目標(biāo)統(tǒng)計(jì)參數(shù),以及模擬的功率譜與目標(biāo)功率譜進(jìn)行比較,驗(yàn)證了基于LPZ譜分析法的非高斯脈動(dòng)風(fēng)壓模擬算法的有效性.

非高斯脈動(dòng)風(fēng)壓模擬;線性預(yù)測和Z轉(zhuǎn)換譜分析法;Johnson變換系統(tǒng);數(shù)字濾波

非高斯過程的模擬方法可以分為兩類:第一類是根據(jù)任意指定的高階特征統(tǒng)計(jì)參數(shù)(例如斜度與峰度)及其目標(biāo)功率譜密度(power spectral density,PSD)函數(shù)來實(shí)現(xiàn);第二類是通過任意指定的邊緣概率分布函數(shù)(probability density function,PDF)及其目標(biāo)功率譜密度(PSD)函數(shù)來實(shí)現(xiàn).Gurley等[1]、Seong等[2]和Suresh等[3]等在20世紀(jì)90年代就開始了第一類非高斯過程的模擬研究.Seong等[4]采用指數(shù)峰值模型進(jìn)行了單變量和多變量非高斯風(fēng)壓時(shí)程的模擬.Suresh等[3,5]基于快速傅里葉變換(fast Fourier transform,FFT)技術(shù),采用參數(shù)較少的指數(shù)峰值模型模擬了一維單變量非高斯風(fēng)壓時(shí)程,并用于大跨低矮屋蓋的風(fēng)振分析.第一類模擬方法過于繁瑣,建立的轉(zhuǎn)換關(guān)系往往為非線性,且在建立過程中并沒有考慮PSD,因而轉(zhuǎn)換后的PSD可能出現(xiàn)負(fù)值.對(duì)于第二類非高斯模擬,主要有KLE(Karhunen-Loeve)法[6]和Hermite譜修正法[7-9]等.Yang等[9]通過Hermite PDF模型模擬出了多變量非高斯風(fēng)壓.Ma等[10]提出了基于增強(qiáng)現(xiàn)實(shí)(augmented reality,AR)模型與Johnson逆變換的風(fēng)壓模擬方法.第二類模擬通常需要采用迭代算法.然而,隨著迭代次數(shù)的增加,模擬數(shù)值不一定會(huì)逐漸逼近目標(biāo)值,且當(dāng)?shù)揭欢ù螖?shù)時(shí),反而可能出現(xiàn)發(fā)散現(xiàn)象.針對(duì)第一類模擬算法,本工作提出一種基于線性預(yù)測和Z轉(zhuǎn)換(linear prediction and Z-transform,LPZ)譜分析法的非高斯脈動(dòng)風(fēng)壓非迭代模擬算法,簡化了已有算法,提高了運(yùn)算效率與模擬精度.

1 LPZ譜分析法

根據(jù)隨機(jī)過程理論,表征隨機(jī)變量數(shù)學(xué)特征的統(tǒng)計(jì)參數(shù)如下:均值(一階)μ=w=表∫示隨機(jī)過程的擺動(dòng)中心,其中p(w)表示隨機(jī)過程w的概率密度函數(shù);方差(二階)表示隨機(jī)過程在某時(shí)刻對(duì)于其均值的偏離程度;偏度(三階)是定量描述獨(dú)立隨機(jī)過程與高斯分布偏離程度的基本量;峰表示分布曲線尖峭或平坦的程度.

高斯過程的偏度為0而峰度為3.偏度不為0或峰度不為3的隨機(jī)過程即為非高斯過程.偏度反映了以其平均值為中心的分布的不對(duì)稱程度,正偏度表示不對(duì)稱部分的分布更趨向正值,負(fù)偏度表示不對(duì)稱部分的分布更趨向負(fù)值.峰度能夠表明時(shí)域信號(hào)峰值的概率大小,描述了概率密度函數(shù)中間部分的尖銳度或平滑度以及尾部的窄度或?qū)挾?當(dāng)峰度K>3時(shí),表示分布呈高峰度,具有相對(duì)多的高峰值和低谷值,此時(shí)的非高斯過程稱為超高斯隨機(jī)過程(softening non-Gaussian processes).當(dāng)K<3時(shí),表示分布呈低峰度,具有相對(duì)較少的高峰值和低谷值,此時(shí)的非高斯過程稱為亞高斯隨機(jī)過程(hardening non-Gaussian processes).

把離散的隨機(jī)脈動(dòng)風(fēng)壓序列{wr|r=0,1,2,···,N-1}作為研究的隨機(jī)過程,則表征該隨機(jī)脈動(dòng)風(fēng)壓的數(shù)學(xué)統(tǒng)計(jì)參數(shù)分別為

LPZ譜分析法是FFT的一種替代和改進(jìn)[11].對(duì)于一個(gè)有限時(shí)間內(nèi)無限數(shù)據(jù)序列的離散信號(hào),運(yùn)用LPZ譜分析法可以避免因FFT的局限性和數(shù)據(jù)截?cái)嗨斐傻娜毕莺驼`差.盡管FFT被普遍認(rèn)為是信號(hào)從時(shí)域轉(zhuǎn)換到頻域的最終機(jī)制,但LPZ譜分析法具有更為優(yōu)良的頻率分辨率和靈敏度.此外,LPZ譜分析法可降噪且不損失分辨率,這是FFT算法無法達(dá)到的.正因?yàn)閮?yōu)良靈敏度和高分辨率,LPZ譜分析法替代FFT可獲得一些不容易獲得的頻譜信息,提高了模擬計(jì)算的準(zhǔn)確性.

線性預(yù)測(linear prediction,LP)理論表明,每個(gè)復(fù)值數(shù)據(jù)點(diǎn)yn可由如下數(shù)值點(diǎn)yn+m的線性組合來求得其近似結(jié)果:

式中,am為未知線性系數(shù),M為濾波長度,N為采樣點(diǎn)總數(shù).

由式(1)可以推導(dǎo)出無限的數(shù)據(jù)序列以及由Z轉(zhuǎn)換定義的相應(yīng)譜函數(shù)S(z):

式中,z=exp(-iωτ),1/τ為采樣率,ω/2π為頻率值.式(2)可以通過式(1)進(jìn)行簡化,得到LPZ譜分析法公式:

同時(shí),式(1)的PSD函數(shù)可表示為

2 Johnson變換系統(tǒng)

Johnson變換系統(tǒng)可將高斯過程轉(zhuǎn)變?yōu)閾碛胁煌群头宥鹊姆歉咚剐蛄?基于統(tǒng)計(jì)參數(shù)算法,Johnson變換系統(tǒng)產(chǎn)生的頻率曲線可用于計(jì)算已知四階統(tǒng)計(jì)參數(shù)的隨機(jī)分布.在這一系統(tǒng)中,主要有3種不同形式的頻率曲線:無界系統(tǒng)(unbound system),即有界系統(tǒng)(bounded system),即為高斯隨機(jī)序列,η′′為給定偏度與峰度的非高斯隨機(jī)序列,γ,δ,ξ和λ為常數(shù),可以通過擬合前四階參數(shù)得到.Hill等[12]給出了求解這些常數(shù)的公式.

3 基于LPZ譜分析法的非高斯脈動(dòng)風(fēng)壓模擬新算法

本工作采用LPZ數(shù)字濾波技術(shù)[13]和Johnson變換系統(tǒng),建立基于LPZ譜分析法的非高斯脈動(dòng)風(fēng)壓模擬新算法.本算法屬于第一類模擬,即根據(jù)任意指定的高階特征統(tǒng)計(jì)參數(shù)(如斜度與峰度)及其目標(biāo)功率譜密度(PSD)函數(shù)來實(shí)現(xiàn)非高斯隨機(jī)過程的模擬.

首先假設(shè){ηr}為一組高斯隨機(jī)序列,則這組高斯序列可通過有限脈沖響應(yīng)(finite impulse response,FIR)濾波器表示為

對(duì)式(5)的兩邊分別進(jìn)行LPZ轉(zhuǎn)換,得到

式中,Ak為高斯序列{ηr}的LPZ轉(zhuǎn)換形式,Wk為wr的LPZ轉(zhuǎn)換形式,Hk為系統(tǒng)的轉(zhuǎn)換方程(或頻率特性),具體等式關(guān)系為

其中z=exp(-iωτ),r,k=0,1,2,···,N-1.Hk的逆LPZ轉(zhuǎn)換形式,即為濾波因子h:

這二者之間的PSD函數(shù)存在如下轉(zhuǎn)化關(guān)系:

式中,Sk為所求的非高斯隨機(jī)序列PSD函數(shù),而Sη為非高斯白噪聲隨機(jī)序列的PSD函數(shù),為常數(shù).隨機(jī)生成的高斯序列{ηr},通過Johnson變換系統(tǒng)轉(zhuǎn)化成擁有一定偏度SSkη和峰度Kη的非高斯白噪聲輸入序列.最后輸出序列中的偏度SSk和峰度K與輸入序列中SSkη和Kη的關(guān)系為

綜上所述,為模擬指定偏度SSk、峰度K和PSD函數(shù)SSk的非高斯脈動(dòng)風(fēng)壓,整個(gè)模擬過程如下.

步驟1已知Sk,并令Sη=C(C為常數(shù)),使用式(9)求解濾波方程

步驟2使用式(8)求解濾波系數(shù)hr.

步驟3已知非高斯序列的峰度K和偏度SSk,使用式(10)和(11)求解原輸入序列,即求解非高斯白噪聲風(fēng)壓時(shí)間序列η的偏度SSkη和峰度Kη.

步驟4通過Johnson變換系統(tǒng),將隨機(jī)生成的高斯序列{ηr}轉(zhuǎn)換為具有指定偏度SSkη和峰度Kη的非高斯白噪聲隨機(jī)序列η.

步驟5已知η和hr,通過方程(5)對(duì)步驟4中求出的非高斯白噪聲隨機(jī)風(fēng)壓η進(jìn)行降噪處理,求出擁有指定偏度SSk、峰度K和PSD函數(shù)SSk的非高斯隨機(jī)風(fēng)壓.

4 數(shù)值驗(yàn)證

4.1 白噪聲非高斯信號(hào)模擬

首先,對(duì)一維單變量白噪聲非高斯信號(hào)進(jìn)行模擬,以驗(yàn)證上述方法的有效性.選圖1中的目標(biāo)PSD曲線為所需模擬的非高斯信號(hào)目標(biāo)功率譜密度函數(shù)Sk,該頻譜為200～500 Hz上的平直譜,均方根值σ約為5.96 g.由于fu=2 000 Hz,根據(jù)采樣定理確定fs= 5 120 Hz,N=4 096,Δf=1.25 Hz.模擬的目標(biāo)PSD函數(shù)為Sk,目標(biāo)峰度K=7.0,偏度SSk=-1.10.假設(shè)Sη=1,由式(10)得Hk=使用式(11)和(12),求得非高斯白噪聲隨機(jī)序列η的指定偏度SSkη=-2.33,峰度Kη=17.53.

圖1 濾波輸出后的非高斯白噪聲隨機(jī)序列PSD函數(shù)F ig.1 PSD function of the fi ltered non-Gaussian white noise randomsequences

通過Johnson變換系統(tǒng),使隨機(jī)生成的高斯序列{ηr}轉(zhuǎn)換為具有指定偏度SSkη和峰度Kη的非高斯白噪聲隨機(jī)序列η,如圖2所示.圖3為非高斯白噪聲隨機(jī)序列η的PSD函數(shù).

圖2 Johnson變換系統(tǒng)產(chǎn)生的非高斯白噪聲隨機(jī)序列F ig.2 Non-Gaussian white noise randomsequences generated by Johnson translator system

經(jīng)過FIR濾波器的降噪過濾,輸出的非高斯脈動(dòng)風(fēng)壓如圖4所示,其偏度和峰度分別為-1.06與7.08,與目標(biāo)偏度SSk=-1.10和峰度K=7.0相差很小.圖1為濾波輸出后的非高斯脈動(dòng)風(fēng)壓的PSD函數(shù)與目標(biāo)PSD函數(shù)的對(duì)比,發(fā)現(xiàn)二者吻合很好.值得指出的是,該信號(hào)模擬過程中的功率譜屬于分貝頻率坐標(biāo)設(shè)置下的形式,與線性頻率坐標(biāo)下的PSD不為負(fù)數(shù)并不沖突.因此,本工作中非高斯脈動(dòng)風(fēng)壓模擬將使用線性坐標(biāo).

圖3 Johnson變換系統(tǒng)產(chǎn)生的非高斯白噪聲隨機(jī)序列PSD函數(shù)F ig.3 PSD function of the non-Gaussian white noise randomsequences generated by Johnson translator system

圖4 濾波輸出后的非高斯白噪聲隨機(jī)序列Fig.4 Filtered non-Gaussian white noise randomsequences

4.2 非高斯脈動(dòng)風(fēng)壓模擬與算法對(duì)比

下面進(jìn)行一維單變量非高斯脈動(dòng)風(fēng)壓的模擬.使用李錦華等[14]推導(dǎo)出的脈動(dòng)風(fēng)壓功率譜方程模擬符合目標(biāo)脈動(dòng)風(fēng)壓功率譜、偏度和峰度的非高斯風(fēng)壓.Davenport水平脈動(dòng)風(fēng)速功率譜為

表1 相關(guān)模擬參數(shù)Tab le 1 Related simu lation parameters

圖5為通過Johnson變換系統(tǒng)轉(zhuǎn)換得到的非高斯白噪聲風(fēng)壓.圖6為非高斯白噪聲隨機(jī)風(fēng)壓的PSD函數(shù).經(jīng)過FIR濾波器的降噪過濾,輸出的非高斯脈動(dòng)風(fēng)壓如圖7所示,其偏度和峰度分別為-2.49與4.77,與目標(biāo)偏度SSk=-2.4和峰度K=4.8非常接近.圖8為濾波輸出后的非高斯脈動(dòng)風(fēng)壓的PSD函數(shù),與目標(biāo)PSD函數(shù)吻合良好.

圖5 Johnson變換系統(tǒng)產(chǎn)生的非高斯白噪聲脈動(dòng)風(fēng)壓Fig.5 Non-Gaussian white noise fluctuating w ind pressure generated by Johnson translator system

圖6 Johnson變換系統(tǒng)產(chǎn)生的非高斯白噪聲脈動(dòng)風(fēng)壓PSD函數(shù)Fig.6 PSD function of the non-Gaussian white noise fluctuating w ind pressure generated by Johnson translator system

5 結(jié)束語

本工作提出了一種基于LPZ譜分析法的非高斯脈動(dòng)風(fēng)壓模擬新算法,即運(yùn)用Johnson變換系統(tǒng)實(shí)行高斯隨機(jī)過程到非高斯白噪聲的轉(zhuǎn)換,再使用LPZ譜分析實(shí)現(xiàn)了非高斯白噪聲譜到目標(biāo)譜的轉(zhuǎn)化.模擬譜與目標(biāo)譜基本一致,偏度、峰度值與目標(biāo)值近似,故而所模擬的非高斯脈動(dòng)風(fēng)壓與目標(biāo)風(fēng)壓趨于一致.

圖7 濾波輸出后的非高斯白噪聲脈動(dòng)風(fēng)壓Fig.7 Filtered non-Gaussian white noise fluctuating w ind pressure

圖8 模擬非高斯白噪聲脈動(dòng)風(fēng)壓PSD函數(shù)與目標(biāo)PSD函數(shù)的比較Fig.8 Comparisons of the PSD function of the simulated and the target non-Gaussian white noise fluctuating w ind pressure

通過MATLAB編程,使用本算法分別對(duì)非高斯隨機(jī)信號(hào)和非高斯脈動(dòng)風(fēng)壓進(jìn)行了模擬.通過比較模擬非高斯脈動(dòng)風(fēng)壓的統(tǒng)計(jì)參數(shù)(峰度與偏度)與目標(biāo)統(tǒng)計(jì)參數(shù),以及模擬功率譜與目標(biāo)功率譜,驗(yàn)證了本算法的有效性.

[1]G URLEY K,K AREEMA.Analysis interpretation modeling and simu lation of unsteady w ind and pressure data[J].JournalofW ind Engineering and Industrial Aerodynamics,1997,69:657-669. [2]S EONG S H,P ETERKAJ A.Computer simulation of non-Gaussian mu ltiple w ind press time series[J].Journal ofW ind Engineering and Industrial Aerodynamics,1997,72(1):95-105.

[3]S URESH K K,S TATHOpOULOS T.Synthesis of non-Gaussian w ind pressure time series on low building roofs[J].Engineering Structures,1999,21(12):1086-1100.

[4]S EONG S H,P ETERKAJ A.D igital generation of surface-pressure fluctuations with spiky features[J].Journal of W ind Engineering and Industrial Aerodynamics,1998,73(2):181-192.

[5]S URESH K K,S TATHOpOULOS T.Computer simulation of fluctuating w ind pressures on low building roofs[J].Journal ofW ind Engineering and Industrial Aerodynamics,1997,69(1):485-495.

[6]V OˇRCHOVSK′Y M.Simulation of simply cross correlated randomfields by series expansion methods[J].Structural Safety,2008,30(4):337-373.

[7]Y ANG Q,T IAN Y.Amodel of probability density function of non-Gaussian w ind pressure with multiple samples[J].Journal of W ind Engineering and Industrial Aerodynamics,2015,140: 67-78.

[8]H UANG G,L UO Y,G URLEY K R,et al.Revisiting moment-based characterization for w ind pressures[J].Journal ofW ind Engineering and Industrial Aerodynamics,2016,151:158-168.

[9]Y ANG L,G URLEY K R.Effi cient stationary mu ltivariate non-Gaussian simulation based on a Hermite PDF model[J].Probabilistic Engineering Mechanics,2015,42(4):31-41.

[10]MAX L,X U F Y,K AREEMA,et al.Estimation of surface pressure extremes:hybrid data and simulation-based approach[J].Journal of Engineering Mechanics,2016,DOI: 10.1061/(ASCE)EM.1943-7889.0001127.

[11]T ANG J,N ORRIS J R.LPZ spectral analysis using linear prediction and the Z transform[J]. Journal of Chemical Physics,1986,84(9):5210-5211.

[12]H ILL ID,H ILL R,H OLDER R L.Fitting Johnson curves by moments[J].Applied Statistics, 1976,25:149-176.

[13]H U Y Z,T ONDER K.Simulation of 3-D randomrough surface by 2-D digital fi lter and Fourier analysis[J].International Journal of Machine Tools Manufacturing,1992,32(1/2):83-90.

[14]李錦華,李春祥.土木工程隨機(jī)風(fēng)場數(shù)值模擬研究的進(jìn)展[J].振動(dòng)與沖擊,2008,19(9):116-125, 187.

Simu lation of non-Gaussian fl uctuating w ind pressu re based on LPZ spectral analysis

JIANG Lei,LIChunxiang,DENG Ying
(Department of Civil Engineering,Shanghai University,Shanghai 200444,China)

This paper presents an algorithmbased on linear prediction and Z-transform(LPZ)spectralanalysis to simulate non-Gaussian fluctuating w ind pressure.The Gaussian process is converted into a non-Gaussian white noise process by Johnson translator system. The non-Gaussian white noise process is then fi ltered with LPZ spectral analysis to obtain fluctuating w ind pressure.The univariate non-Gaussian randomsignals and non-Gaussian fluctuating w ind pressure are simulated using the algorithm.By comparing the statistical parameters including skewness,kurtosis and power spectral density of the non-Gaussian white noise process and fl uctuating w ind pressure with their target values,it is confi rmed that the algorithmbased on LZP spectral analysis can eff ectively simulate non-Gaussian fluctuating w ind pressure.

non-Gaussian fl uctuating w ind pressure simu lation;linear prediction and Z-transform(LZP)spectral analysis;Johnson translator system;digital fi lter

TU 311

A

1007-2861(2017)04-0600-09

DO I:10.12066/j.issn.1007-2861.1852

2016-09-30

國家自然科學(xué)基金資助項(xiàng)目(51378304)

李春祥(1964—),男,教授,博士生導(dǎo)師,博士,研究方向?yàn)榻Y(jié)構(gòu)振動(dòng)控制、結(jié)構(gòu)風(fēng)工程、結(jié)構(gòu)抗震.

E-mail:Li-chunxiang@vip.sina.com