基于ar和arma時(shí)間序列的地震反應(yīng)分析

基于ar和arma時(shí)間序列的地震反應(yīng)分析

1采用速度隨機(jī)過(guò)程隨機(jī)振動(dòng)模擬方法是中國(guó)最早和最常用的模擬方法是波形合成法。該方法的想法是將隨機(jī)過(guò)程表示為許多隨機(jī)相位的正和負(fù)和合行的總和。該方法是由Shinozuka在20世紀(jì)70年代初首先提出的,現(xiàn)在已比較成熟,它適用于模擬具有任意形狀譜密度的平穩(wěn)隨機(jī)過(guò)程,但該模型具有大量的三角級(jí)數(shù)的運(yùn)算,計(jì)算機(jī)的用時(shí)是相當(dāng)大的。Shinozuka和同濟(jì)大學(xué)將快速傅里葉變換(FFT)引入合成模型,從而大大提高了計(jì)算速度,但該方法對(duì)于模擬二維及二維以上隨機(jī)過(guò)程已不可能。自回歸滑動(dòng)平均法(Auto-RegressionandMovingAverage縮寫ARMA)是20世紀(jì)80年代得到較系統(tǒng)地發(fā)展的一種用于隨機(jī)過(guò)程的數(shù)字模擬方法,該方法首先用于數(shù)字信號(hào)處理方面,Spanos和Mignolet首先將其用于風(fēng)、海浪和地震等隨機(jī)過(guò)程的模擬,本文采用該方法對(duì)隨機(jī)地震的數(shù)值模擬作了進(jìn)一步較為系統(tǒng)的研究。2arma基本理論由于用AR、MA和ARMA模型來(lái)模擬隨機(jī)地震動(dòng)的方法較新,盡管在國(guó)際上已有多篇文獻(xiàn)報(bào)道過(guò),但在國(guó)內(nèi)仍然是空白,至今還未見到有關(guān)這方面的任何文獻(xiàn)報(bào)道,因此在這里對(duì)它的基本理論作一簡(jiǎn)要介紹,并對(duì)本文所采用的方法作一說(shuō)明。本文采用的隨機(jī)地震動(dòng)的模擬方法有兩種,一種是AR模型,一種是ARMA模型。求解ARMA模型有多種方法,本文中的ARMA模型是先計(jì)算高階AR模型的參數(shù),再用此求解低價(jià)ARMA模型的參數(shù)。當(dāng)然,從MA模型也可以發(fā)展到ARMA模型,本文對(duì)此沒(méi)有涉及,本文還采用了文獻(xiàn)給出的方法求解ARMA模型。2.1隨機(jī)地震動(dòng)功率譜密度函數(shù)將隨機(jī)地震動(dòng)看作一個(gè)變量的隨機(jī)過(guò)程。對(duì)于多維和多變量的隨機(jī)過(guò)程,其原理是相似的,只不過(guò)方程中的各個(gè)系數(shù)應(yīng)是n×n的矩陣。對(duì)于AR模型,m階的自回歸隨機(jī)過(guò)程?y,可表示成AR(m),其第r個(gè)樣本?yr可由先前m個(gè)時(shí)刻的響應(yīng)值與同一時(shí)刻的激勵(lì)值按下式遞推得出:?yr=-m∑k=1?ak?yr-k+?b0wr(1)式中:?ak(k=1,?,m),?b0為AR模型的系數(shù),通常?b0=1;?yr-k(k=1,?,m)為先前的m個(gè)時(shí)刻的響應(yīng)值;?yr為當(dāng)前時(shí)刻的響應(yīng)值;wr為當(dāng)前時(shí)刻的白噪聲序列值。對(duì)式(1)兩邊進(jìn)行z變換,得傳遞函數(shù)為?Η(z)=?b01+m∑k=1?akz-k(2)因此,估計(jì)的譜密度為S?y?y(w)=?Η(z)*?Η*Τ(z)=|Η(z)|2Ζ=ejw(3)式中,“*”表示共軛,T表示轉(zhuǎn)置。AR方法的實(shí)質(zhì)就是選擇適當(dāng)?shù)南禂?shù)?ak,?b0,使得由上式確定的譜密度與目標(biāo)譜密度Syy(w)接近,其近似程度可由下式度量:εAR=12ωb[JF(Ζ]ωb-ωbtr[Η-*(ejωt)Syy(ω)Η-Τ(ejωt)dω][JF)](4)式中:ωb為截止圓頻率,滿足Nyquist條件ωb=πΤ?Τ為采樣周期,本文地震的截止圓頻率取50π,則T為0.02s;tr是運(yùn)算符,tr(A)表示對(duì)A角元素的和或A的特征值的和。由εAR關(guān)于ak為最小的條件,可以導(dǎo)出以下m個(gè)線性方程組:R?y?y(l)+m∑k=1?akR?y?y(k-1)=0,l=1?2,?,m(5)由上式可解出AR模型系數(shù)?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)上式中,自相關(guān)函數(shù)可由功率譜密度函數(shù)的傅里葉逆變換求得。因S?y?y(ω)是偶函數(shù),時(shí)間lT的自相關(guān)函數(shù)R?y?y(l)的值由下式求出:R?y?y(l)=[JF(Ζ]ωb-ωbSyy(ω)coslωΤdω[JF)](7)在實(shí)際計(jì)算中,R?y?y(l)的值可由離散傅立葉逆變換求得;對(duì)于某些功率譜密度函數(shù),也可由連續(xù)傅立葉逆變換直接求出,本文兩種方法均采用。AR模型的系數(shù)?b0可由AR過(guò)程與目標(biāo)過(guò)程的總能量相等條件得出: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模型的全部系數(shù),然后用式(1)求出時(shí)間序列?y,它是滿足目標(biāo)功率譜密度的隨機(jī)過(guò)程。本文的隨機(jī)地震動(dòng)功率譜密度函數(shù)采用了Kanai-Tajimi譜和Penzien譜,其表達(dá)式分別見式(37)和式(38),各采用AR模型分別生成了3條人工地震波。3條人工地震波的差別僅參數(shù)ωg不同,ωg依次為15rad/s,10rad/s,5rad/s,其余參數(shù)均相同,并進(jìn)行了譜密度和自相關(guān)性檢驗(yàn),表明它們均是滿足已知功率譜的人工地震波。本次模擬僅局限于平穩(wěn)隨機(jī)過(guò)程的模擬,實(shí)際上,對(duì)于非平穩(wěn)隨機(jī)地震動(dòng)的模擬,還應(yīng)將加速度時(shí)程乘以與時(shí)間有關(guān)的地震動(dòng)的振幅包絡(luò)曲線。文獻(xiàn)給出了存在有限階AR模型的條件:[JF(Ζ]ωb-ωblog[Syy(ω)]dω>-∞[JF)](9)實(shí)際上前面提到的Kanai-Tajimi譜和Penzien譜均滿足上述條件,因此存在有限階AR模型。文獻(xiàn)還給出了一個(gè)不能滿足上述條件的P-M波浪譜,提出將該譜轉(zhuǎn)化為滿足條件式(9)的有限譜,就可以滿意地用AR模型來(lái)擬合。本文通過(guò)計(jì)算發(fā)現(xiàn),自相關(guān)系數(shù)的計(jì)算對(duì)于計(jì)算AR模型的系數(shù)是相當(dāng)重要的。自相關(guān)系數(shù)的計(jì)算不精確常常會(huì)導(dǎo)致解的不可靠性,因此,作者建議對(duì)于那些能夠直接用連續(xù)傅里葉變換積分出的功率譜,應(yīng)優(yōu)先采用直接積分,不能積出的可采用離散傅立葉變換等數(shù)值積分方法。2.2傳遞函數(shù)的解析MA模型本文沒(méi)有專門研究,但該模型與ARMA模型聯(lián)系密切,故在這里略作介紹。對(duì)于一個(gè)m階MA模型響應(yīng)的第r個(gè)樣本yr,可由m個(gè)先前時(shí)刻與m個(gè)以后時(shí)刻的激勵(lì)值按下式算出:yr=m∑l=-mblwr-l(10)同理,對(duì)兩邊進(jìn)行z變換,得傳遞函數(shù)ΗΜA(z)=m∑l=-mblz-l(11)MR模型的譜矩陣可以由傳遞函數(shù)得出:SΜAyy(ω)=Η*ΜA(ejωΤ)ΗΤΜA(ejωΤ)(12)式中:“*”表示復(fù)共軛;T表示轉(zhuǎn)置。用滑動(dòng)平均模擬隨機(jī)過(guò)程,通常要對(duì)目標(biāo)譜密度矩陣Syy(ω)進(jìn)行Cholesky分解,其過(guò)程較為復(fù)雜,但對(duì)一個(gè)變量的情況,僅需進(jìn)行開方即可:Q(ω)=√Syy(ω)(13)bl=12ωb[JF(Ζ]ωb-ωbQ(ω)ejωtdω[JF)](14)2.3arma參數(shù)ak+bl的pom方法自回歸滑動(dòng)平均法是AR法和MA法的結(jié)合與推廣。對(duì)于一個(gè)(p,q)階的自回歸滑動(dòng)平均隨機(jī)過(guò)程ˉy,表示為ARMA(p,q),其第r個(gè)樣本ˉyr可由先前p個(gè)時(shí)刻的響應(yīng)值和先前q個(gè)時(shí)刻的激勵(lì)值來(lái)表示,即ˉyr=-p∑k=1akˉyr-k+q∑l=0blwr-1(15)式中:ak(k=1,…,p),bl(l=0,1,…,q)——ARMA模型的系數(shù),通常a0=1;ˉyr-k(k=1,?,p)——先前的p個(gè)時(shí)刻的響應(yīng)值;wr——當(dāng)前時(shí)刻的白噪聲序列值;ˉyr——當(dāng)前時(shí)刻的響應(yīng)列值。對(duì)該式兩邊z變換可得傳遞函數(shù)的z變換形式Η(z)=p∑l=0blwr-l1+p∑k=1akz-k(16)其功率譜密度Sˉyˉy(ω)為Sˉyˉy(ω)=Η(jωΤ)Η*Τ(jωΤ)=|Η(jωΤ)|2(17)式中,上式“*”表示復(fù)共軛,上標(biāo)T表示轉(zhuǎn)置確定ARMA模型的系數(shù)ak(k=1,2,…,p),bl(l=0,1,…,q)有多種方法,文獻(xiàn)給出了ACM(Auto-Cross-CorrelationMathing)和POM(PowerOrderMacthing)兩種方法,本文僅討論高階AR模型與低價(jià)ARMA模型匹配來(lái)確定ARMA參數(shù)ak,bl的POM方法。ACM方法是根據(jù)對(duì)于同一個(gè)時(shí)間序列的隨機(jī)過(guò)程,其ARMA模型的響應(yīng)自相關(guān)與響應(yīng)-激勵(lì)互相關(guān)應(yīng)匹配,可得如下方程: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)對(duì)式(15)兩邊同乘以ˉyr-i,并取數(shù)學(xué)期望,得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的匹配關(guān)系,得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,并取數(shù)學(xué)期望,得E[ˉyrwr-i]=-E[p∑k=1akˉyr-kwr-i]+q∑l=0[blwr-lwr-i](24)當(dāng)i=0時(shí)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)當(dāng)i≠0時(shí)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系統(tǒng)的參數(shù),式中互相關(guān)系數(shù)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系統(tǒng)和ARMA系統(tǒng)的傳遞函數(shù)等效,即各邊z的各次冪的系數(shù)的值相等來(lái)確定,即?b01+m∑s=1?asz-S=q∑l=0blz-l1+p∑k=1akz-k(33)展開上式,由z相同冪次對(duì)應(yīng)的系數(shù)應(yīng)相等,得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]可以很容易解出,當(dāng)p=q時(shí),無(wú)式(35),當(dāng)p>q時(shí),先用式(35)求出ak(k=q+1,q+2,…,p),并用式(34)求出ak(k=1,2,…,q),至此ARMA模型的所有參數(shù)就已確定。由式(15)求出滿足給定功率譜密度的隨機(jī)過(guò)程的時(shí)間序列yˉ。由于低階ARMA模型是從高階的AR模型而來(lái),其精度依賴于AR的精度,一般較AR模型的精度低。3功率譜密度模型地震動(dòng)的觀測(cè)資料表明,地震動(dòng)采用隨機(jī)模型可能更適宜些,在隨機(jī)地震分析中較多采用金井(Kanai)(1957)和田治見(Tajimi)(1960)提出的過(guò)濾白噪聲的加速度功率譜密度,其表達(dá)式如下:S(ω)=1+4ξg2(ω/ωg)2[1-(ω/ωg)2]2+4ξg(ω/ωg)2G0(39)式中:ωg——場(chǎng)地土的卓越頻率;ξg——場(chǎng)地土的阻尼比;G0——基巖白噪聲的譜強(qiáng)度。金井建議取ωg=15.6rad/s和ξg=0.6表示堅(jiān)硬土層,ωg和ξg的取值根據(jù)場(chǎng)地土的軟弱情況來(lái)確定。該模型是將場(chǎng)地土看成單自由度線性過(guò)濾器,考慮了場(chǎng)地土對(duì)地震動(dòng)頻率特性的影響,但得到的加速度功率譜函數(shù)仍然具有如下缺點(diǎn):(1)在ω=0處出現(xiàn)明顯的奇異點(diǎn),得到的地面加速度、速度和位移為無(wú)限值,這不符合地震常識(shí);(2)過(guò)份夸大了長(zhǎng)周期即在鄰近ω=ωg處的地震能量,過(guò)份削弱了高頻處的能量。Penzien(1975)為克服Kanai-Tajimi譜的缺點(diǎn),將場(chǎng)地土看作雙自由度的濾波器,將上述模型再經(jīng)過(guò)一次濾波,減消弱其極低頻率分量的能量,其表達(dá)式如下:S0(ω)=1+4ξg2(ω/ωg)2[1-(ω/ωg)2]2+4ξg(ω/ωg)2(ωf2-ω2)2+4ω2ωf2ξf2G0(40)當(dāng)式中ωg=15rad/s,ξg=0.55,ωf=3rad/s,ξf=0.6和G0=0.1378時(shí),其人工加速度的時(shí)程同ElCentro波接近式(40)中ωg,ξg含義同式(39)。頻率ωg和阻尼比ξg是為了給出所需要的過(guò)濾特性而選擇的,一般取ωf=0.1～0.2ωg,ξf=ξg。Penzien功率譜密度模型修正了Kanai-Tajimi功率譜密度模型的低頻部分,其余頻段幾乎不變,使之與實(shí)際地震動(dòng)更加吻合。Kanai-Tajimi和Penzien加速度功率譜的自相關(guān)函數(shù)可以用留數(shù)定理顯式積出,當(dāng)已知參數(shù)為上述給定的參數(shù)時(shí),它們的自相關(guān)函分別是: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)為了反映地面土層的影響,可適當(dāng)調(diào)整參數(shù)ωg,其它參數(shù)保持不變,當(dāng)ωg分別取10rad/s和5rad/s時(shí),其自相關(guān)函數(shù)分別為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)由于其總的階數(shù)較AR模型大為降低,因此運(yùn)算次數(shù)更少,時(shí)間更短,這對(duì)于生成大規(guī)模的時(shí)間序列如長(zhǎng)時(shí)間或多變量隨機(jī)過(guò)程是相當(dāng)快速而有效的。在前面AR模型的基礎(chǔ)上,用ARMA模型模擬的隨機(jī)地震動(dòng),自相關(guān)檢驗(yàn)及譜檢驗(yàn)表明,ARMA模型仍然是相當(dāng)有效的,如AR模型的階數(shù)取為59階,與之匹配的ARMA模型的階數(shù)取p=q=4即ARMA(4,4),其精度比AR(59)略為下降,但花的時(shí)間僅為AR(59)的1/7。文獻(xiàn)還提出另外一種直接的用ARMA模型進(jìn)行譜估計(jì)的修正尤力-沃克方法,由于該方法較為復(fù)雜,這里就給出用該方法進(jìn)行時(shí)間序列生成時(shí)所用的ARMA模型的系數(shù),見表2,該方法ARMA模型系數(shù)不須從AR模型轉(zhuǎn)化來(lái),計(jì)算效率高,但該方法要想得到較好的功率譜估計(jì)效果,必須采用更高階數(shù),本文采用20階,其譜檢驗(yàn)如圖2。用表1或表2中提供的數(shù)據(jù),就可以生成滿足相應(yīng)功率譜的隨機(jī)加速度時(shí)程曲線。從圖1可和2可以看出,本文的隨機(jī)地震動(dòng)的模擬是成功的,生成的加速度時(shí)間序列是符合給定加速度功率譜的,按照抗震設(shè)防要求,調(diào)整幅度用于工程應(yīng)用,但是在實(shí)際應(yīng)用時(shí)應(yīng)去除最初的幾秒時(shí)程,以保證較高的精度要求。續(xù)表1AR模型系數(shù)作者認(rèn)為,今后用ARMA生成人工地震波的二點(diǎn)改進(jìn)方向,一是由于不少國(guó)家包括中國(guó)、歐洲的抗震規(guī)范規(guī)定的結(jié)構(gòu)抗震設(shè)計(jì)均是按