基于ar和arma时间序列的地震反应分析

基于ar和arma时间序列的地震反应分析

1采用速度随机过程随机振动模拟方法是中国最早和最常用的模拟方法是波形合成法。该方法的想法是将随机过程表示为许多随机相位的正和负和合行的总和。该方法是由Shinozuka在20世纪70年代初首先提出的,现在已比较成熟,它适用于模拟具有任意形状谱密度的平稳随机过程,但该模型具有大量的三角级数的运算,计算机的用时是相当大的。Shinozuka和同济大学将快速傅里叶变换(FFT)引入合成模型,从而大大提高了计算速度,但该方法对于模拟二维及二维以上随机过程已不可能。自回归滑动平均法(Auto-RegressionandMovingAverage缩写ARMA)是20世纪80年代得到较系统地发展的一种用于随机过程的数字模拟方法,该方法首先用于数字信号处理方面,Spanos和Mignolet首先将其用于风、海浪和地震等随机过程的模拟,本文采用该方法对随机地震的数值模拟作了进一步较为系统的研究。2arma基本理论由于用AR、MA和ARMA模型来模拟随机地震动的方法较新,尽管在国际上已有多篇文献报道过,但在国内仍然是空白,至今还未见到有关这方面的任何文献报道,因此在这里对它的基本理论作一简要介绍,并对本文所采用的方法作一说明。本文采用的随机地震动的模拟方法有两种,一种是AR模型,一种是ARMA模型。求解ARMA模型有多种方法,本文中的ARMA模型是先计算高阶AR模型的参数,再用此求解低价ARMA模型的参数。当然,从MA模型也可以发展到ARMA模型,本文对此没有涉及,本文还采用了文献给出的方法求解ARMA模型。2.1随机地震动功率谱密度函数将随机地震动看作一个变量的随机过程。对于多维和多变量的随机过程,其原理是相似的,只不过方程中的各个系数应是n×n的矩阵。对于AR模型,m阶的自回归随机过程ˆy,可表示成AR(m),其第r个样本ˆyr可由先前m个时刻的响应值与同一时刻的激励值按下式递推得出:ˆyr=-m∑k=1ˆakˆyr-k+ˆb0wr(1)式中:ˆak(k=1,⋯,m),ˆb0为AR模型的系数,通常ˆb0=1;ˆyr-k(k=1,⋯,m)为先前的m个时刻的响应值;ˆyr为当前时刻的响应值;wr为当前时刻的白噪声序列值。对式(1)两边进行z变换,得传递函数为ˆΗ(z)=ˆb01+m∑k=1ˆakz-k(2)因此,估计的谱密度为Sˆyˆy(w)=ˆΗ(z)*ˆΗ*Τ(z)=|Η(z)|2Ζ=ejw(3)式中,“*”表示共轭,T表示转置。AR方法的实质就是选择适当的系数ˆak,ˆb0,使得由上式确定的谱密度与目标谱密度Syy(w)接近,其近似程度可由下式度量:εAR=12ωb[JF(Ζ]ωb-ωbtr[Η-*(ejωt)Syy(ω)Η-Τ(ejωt)dω][JF)](4)式中:ωb为截止圆频率,满足Nyquist条件ωb=πΤ‚Τ为采样周期,本文地震的截止圆频率取50π,则T为0.02s;tr是运算符,tr(A)表示对A角元素的和或A的特征值的和。由εAR关于ak为最小的条件,可以导出以下m个线性方程组:Rˆyˆy(l)+m∑k=1ˆakRˆyˆy(k-1)=0,l=1‚2,⋯,m(5)由上式可解出AR模型系数ˆak。为方便编程,将上式写成矩阵形式:l=1,⋯,mk=1⋮m[l,1,⋯,m一般元素Rˆyˆy(k-l)]{ˆa1ˆa2⋯ˆam}={-Rˆyˆy(1)-Rˆyˆy(2)⋯-Rˆyˆy(m)}(6)上式中,自相关函数可由功率谱密度函数的傅里叶逆变换求得。因Sˆyˆy(ω)是偶函数,时间lT的自相关函数Rˆyˆy(l)的值由下式求出:Rˆyˆy(l)=[JF(Ζ]ωb-ωbSyy(ω)coslωΤdω[JF)](7)在实际计算中,Rˆyˆy(l)的值可由离散傅立叶逆变换求得;对于某些功率谱密度函数,也可由连续傅立叶逆变换直接求出,本文两种方法均采用。AR模型的系数ˆb0可由AR过程与目标过程的总能量相等条件得出:b20=12ωbm∑k=0ˆakRˆyˆy(k)=12ωb(Rˆyˆy(0)+m∑k=1ˆakRˆyˆy)(8)至此,由式(6)和(8)就可以确定AR模型的全部系数,然后用式(1)求出时间序列ˆy,它是满足目标功率谱密度的随机过程。本文的随机地震动功率谱密度函数采用了Kanai-Tajimi谱和Penzien谱,其表达式分别见式(37)和式(38),各采用AR模型分别生成了3条人工地震波。3条人工地震波的差别仅参数ωg不同,ωg依次为15rad/s,10rad/s,5rad/s,其余参数均相同,并进行了谱密度和自相关性检验,表明它们均是满足已知功率谱的人工地震波。本次模拟仅局限于平稳随机过程的模拟,实际上,对于非平稳随机地震动的模拟,还应将加速度时程乘以与时间有关的地震动的振幅包络曲线。文献给出了存在有限阶AR模型的条件:[JF(Ζ]ωb-ωblog[Syy(ω)]dω>-∞[JF)](9)实际上前面提到的Kanai-Tajimi谱和Penzien谱均满足上述条件,因此存在有限阶AR模型。文献还给出了一个不能满足上述条件的P-M波浪谱,提出将该谱转化为满足条件式(9)的有限谱,就可以满意地用AR模型来拟合。本文通过计算发现,自相关系数的计算对于计算AR模型的系数是相当重要的。自相关系数的计算不精确常常会导致解的不可靠性,因此,作者建议对于那些能够直接用连续傅里叶变换积分出的功率谱,应优先采用直接积分,不能积出的可采用离散傅立叶变换等数值积分方法。2.2传递函数的解析MA模型本文没有专门研究,但该模型与ARMA模型联系密切,故在这里略作介绍。对于一个m阶MA模型响应的第r个样本yr,可由m个先前时刻与m个以后时刻的激励值按下式算出:yr=m∑l=-mblwr-l(10)同理,对两边进行z变换,得传递函数ΗΜA(z)=m∑l=-mblz-l(11)MR模型的谱矩阵可以由传递函数得出:SΜAyy(ω)=Η*ΜA(ejωΤ)ΗΤΜA(ejωΤ)(12)式中:“*”表示复共轭;T表示转置。用滑动平均模拟随机过程,通常要对目标谱密度矩阵Syy(ω)进行Cholesky分解,其过程较为复杂,但对一个变量的情况,仅需进行开方即可:Q(ω)=√Syy(ω)(13)bl=12ωb[JF(Ζ]ωb-ωbQ(ω)ejωtdω[JF)](14)2.3arma参数ak+bl的pom方法自回归滑动平均法是AR法和MA法的结合与推广。对于一个(p,q)阶的自回归滑动平均随机过程ˉy,表示为ARMA(p,q),其第r个样本ˉyr可由先前p个时刻的响应值和先前q个时刻的激励值来表示,即ˉyr=-p∑k=1akˉyr-k+q∑l=0blwr-1(15)式中:ak(k=1,…,p),bl(l=0,1,…,q)——ARMA模型的系数,通常a0=1;ˉyr-k(k=1,⋯,p)——先前的p个时刻的响应值;wr——当前时刻的白噪声序列值;ˉyr——当前时刻的响应列值。对该式两边z变换可得传递函数的z变换形式Η(z)=p∑l=0blwr-l1+p∑k=1akz-k(16)其功率谱密度Sˉyˉy(ω)为Sˉyˉy(ω)=Η(jωΤ)Η*Τ(jωΤ)=|Η(jωΤ)|2(17)式中,上式“*”表示复共轭,上标T表示转置确定ARMA模型的系数ak(k=1,2,…,p),bl(l=0,1,…,q)有多种方法,文献给出了ACM(Auto-Cross-CorrelationMathing)和POM(PowerOrderMacthing)两种方法,本文仅讨论高阶AR模型与低价ARMA模型匹配来确定ARMA参数ak,bl的POM方法。ACM方法是根据对于同一个时间序列的随机过程,其ARMA模型的响应自相关与响应-激励互相关应匹配,可得如下方程:Rˉyˉy(k)=Rˆyˆy(k),l=0,1,⋯,p(18)Rˉyω(-l)=Rˆyw(-l),l=0,1,⋯,q(19)ˉyr=ˆyr(20)对式(15)两边同乘以ˉyr-i,并取数学期望,得E[ˉyrˉyr-i]=-E[ˉyr-ip∑k=1akˉyr-k]+E[ˉyr-iq∑l=0blwr-l](21)Rˉyˉy(i)=-p∑k=1akRˆyˉy(k-i)+q∑l=0blRˉyw(i-l),i=1,2,⋯,p(22)由与AR的匹配关系,得p∑k=1Rˆyˆy(k-i)ak-q∑l=0Rˆyw(i-l)bl=-Rˆyˆy(i),i=1,2,⋯,p(23)同理将式(20)两边乘以wr-i,并取数学期望,得E[ˉyrwr-i]=-E[p∑k=1akˉyr-kwr-i]+q∑l=0[blwr-lwr-i](24)当i=0时2ωbˆb0=-p∑k=1akRˉyw(k)+2ωbb0(25)2ωbˆb0=-p∑k=1akRˆyw(k)+2ωbb0(26)b0=ˆb0(27)当i≠0时Rˆyw(-i)=-p∑k=1akRˉyw(k-i)+2ωbbi(28)即p∑k=1Rˆyw(k-i)ak-2ωbbi=-Rˆyw(-i),i=1,2,⋯,q(29)由式(25),(26),(27)和(28),可以确定ARMA系统的参数,式中互相关系数Rw(i)可由AR模型按下式确定:Rˆyw(l)=0,l>0(30)Rˆyw(0)=2ωbˆb0,l=0(31)Rˆyw(-l)=min(m,l)∑k=1ˆakRˆyw(k-l),l>0(32)本文采用另一种方法,即POM法。该方法是由AR系统和ARMA系统的传递函数等效,即各边z的各次幂的系数的值相等来确定,即ˆb01+m∑s=1ˆasz-S=q∑l=0blz-l1+p∑k=1akz-k(33)展开上式,由z相同幂次对应的系数应相等,得ak=1ˆb0[k-1∑l=0blˆak-l+ˆbk],k=1,2,⋯,q(34)ak=1ˆb0q∑l=0blˆak-l,k=q+1,⋯,p(35)q∑l=0blˆak-l=-bˆak,k=p+1,⋯,p+q(36)b0=ˆb(37)该一系列方程的解法如下,首先,将式(36)写成矩阵形式,[ˆap+1-1ˆaΡ+1-2⋯ˆap+1-qˆap+2-1ˆap+2-2⋯ˆap+2-q⋯⋯⋯⋯ˆap+q-1ˆap+q-2⋯ˆap+q-q]{b1b2⋯bq}=-b0{ˆap+1ˆap+2⋯ˆap+q}(38)则[b1,b2,…,bq]可以很容易解出,当p=q时,无式(35),当p>q时,先用式(35)求出ak(k=q+1,q+2,…,p),并用式(34)求出ak(k=1,2,…,q),至此ARMA模型的所有参数就已确定。由式(15)求出满足给定功率谱密度的随机过程的时间序列y¯。由于低阶ARMA模型是从高阶的AR模型而来,其精度依赖于AR的精度,一般较AR模型的精度低。3功率谱密度模型地震动的观测资料表明,地震动采用随机模型可能更适宜些,在随机地震分析中较多采用金井(Kanai)(1957)和田治见(Tajimi)(1960)提出的过滤白噪声的加速度功率谱密度,其表达式如下:S(ω)=1+4ξg2(ω/ωg)2[1-(ω/ωg)2]2+4ξg(ω/ωg)2G0(39)式中:ωg——场地土的卓越频率;ξg——场地土的阻尼比;G0——基岩白噪声的谱强度。金井建议取ωg=15.6rad/s和ξg=0.6表示坚硬土层,ωg和ξg的取值根据场地土的软弱情况来确定。该模型是将场地土看成单自由度线性过滤器,考虑了场地土对地震动频率特性的影响,但得到的加速度功率谱函数仍然具有如下缺点:(1)在ω=0处出现明显的奇异点,得到的地面加速度、速度和位移为无限值,这不符合地震常识;(2)过份夸大了长周期即在邻近ω=ωg处的地震能量,过份削弱了高频处的能量。Penzien(1975)为克服Kanai-Tajimi谱的缺点,将场地土看作双自由度的滤波器,将上述模型再经过一次滤波,减消弱其极低频率分量的能量,其表达式如下:S0(ω)=1+4ξg2(ω/ωg)2[1-(ω/ωg)2]2+4ξg(ω/ωg)2(ωf2-ω2)2+4ω2ωf2ξf2G0(40)当式中ωg=15rad/s,ξg=0.55,ωf=3rad/s,ξf=0.6和G0=0.1378时,其人工加速度的时程同ElCentro波接近式(40)中ωg,ξg含义同式(39)。频率ωg和阻尼比ξg是为了给出所需要的过滤特性而选择的,一般取ωf=0.1～0.2ωg,ξf=ξg。Penzien功率谱密度模型修正了Kanai-Tajimi功率谱密度模型的低频部分,其余频段几乎不变,使之与实际地震动更加吻合。Kanai-Tajimi和Penzien加速度功率谱的自相关函数可以用留数定理显式积出,当已知参数为上述给定的参数时,它们的自相关函分别是:Kanai-Tajimi谱R1Κ(t)=12π[e-11.12t(70.21cos(10.06t)-7.37sin(10.06t)](41)Penzien谱R1Ρ(t)=12π[9.69e-11.12tcos(10.06t)-0.797e-11.12tsin(10.06t)-0.531e-1.8tcos(2.4t)-1.264e-1.8tsin(2.4t)](42)为了反映地面土层的影响,可适当调整参数ωg,其它参数保持不变,当ωg分别取10rad/s和5rad/s时,其自相关函数分别为R2Κ(t)=12π[e-7.416t(46.809cos(6.708t)-4.917sin(6.708t)(43)R2Ρ(t)=-12π[-6.529e-7.416tcos(6.708t)+0.345e-7.416tsin(6.708t)+0.600e-1.8tcos(2.4t)+1.278e-1.8tsin(2.4t)](44)R3Κ(t)=-12π[e-3.708t(-23.405cos(3.354t)+2.459sin(3.354t)](45)R3Ρ(t)=12π[3.506e-3.708tcos(3.354t)+0.507e-3.708tsin(3.354t)-0.9478e-1.8tcos(2.4t)-1.494e-1.8tsin(2.4t)](46)由于其总的阶数较AR模型大为降低,因此运算次数更少,时间更短,这对于生成大规模的时间序列如长时间或多变量随机过程是相当快速而有效的。在前面AR模型的基础上,用ARMA模型模拟的随机地震动,自相关检验及谱检验表明,ARMA模型仍然是相当有效的,如AR模型的阶数取为59阶,与之匹配的ARMA模型的阶数取p=q=4即ARMA(4,4),其精度比AR(59)略为下降,但花的时间仅为AR(59)的1/7。文献还提出另外一种直接的用ARMA模型进行谱估计的修正尤力-沃克方法,由于该方法较为复杂,这里就给出用该方法进行时间序列生成时所用的ARMA模型的系数,见表2,该方法ARMA模型系数不须从AR模型转化来,计算效率高,但该方法要想得到较好的功率谱估计效果,必须采用更高阶数,本文采用20阶,其谱检验如图2。用表1或表2中提供的数据,就可以生成满足相应功率谱的随机加速度时程曲线。从图1可和2可以看出,本文的随机地震动的模拟是成功的,生成的加速度时间序列是符合给定加速度功率谱的,按照抗震设防要求,调整幅度用于工程应用,但是在实际应用时应去除最初的几秒时程,以保证较高的精度要求。续表1AR模型系数作者认为,今后用ARMA生成人工地震波的二点改进方向,一是由于不少国家包括中国、欧洲的抗震规范规定的结构抗震设计均是按