阮鑫鑫, 劉章軍, 姜云木

(1. 信陽師范大學(xué) 建筑與土木工程學(xué)院,河南 信陽 464000; 2. 武漢工程大學(xué) 土木工程與建筑學(xué)院, 武漢 430074)

強(qiáng)烈地震往往會(huì)導(dǎo)致工程結(jié)構(gòu)坍塌,造成重大人員傷亡和經(jīng)濟(jì)損失。為了保障工程結(jié)構(gòu)在地震作用下的安全性和可靠性,通常需要選用一定數(shù)量的地震動(dòng)輸入來進(jìn)行工程結(jié)構(gòu)的動(dòng)力時(shí)程反應(yīng)分析,因而選用合理的地震動(dòng)輸入時(shí)程至關(guān)重要。雖然目前已積累了豐富的強(qiáng)震記錄,但從中選取適合各種場地和結(jié)構(gòu)的地震動(dòng)輸入時(shí)程仍然較為困難[1]。因此,人工地震動(dòng)模擬成為地震工程領(lǐng)域的研究熱點(diǎn)之一。

地震動(dòng)是典型的隨機(jī)過程。自Housner首次提出使用隨機(jī)過程模擬地震動(dòng)以來,隨機(jī)地震動(dòng)模型的研究得到了長足的發(fā)展[2-5]。目前,隨機(jī)地震動(dòng)模型可分為地震學(xué)模型和工程學(xué)模型兩大類。由于地震學(xué)模型需要對(duì)地震源和地震波傳播路徑等物理機(jī)制了解,而工程學(xué)模型的應(yīng)用相對(duì)簡單。因此,本文主要關(guān)注工程學(xué)模型,此類模型一般有頻域方法和時(shí)域方法兩大主流。頻域方法可直接利用地震動(dòng)的頻域統(tǒng)計(jì)特性,具有理論完善、計(jì)算精度高等特點(diǎn),便于工程應(yīng)用。但通過平穩(wěn)功率譜與強(qiáng)度調(diào)制函數(shù)的地震動(dòng)難以體現(xiàn)地震動(dòng)的時(shí)頻全非平穩(wěn)性,而演變過程的時(shí)頻調(diào)制函數(shù)除了少數(shù)模型[6-7]外,都存在難以參數(shù)化估計(jì)的問題。隨著時(shí)頻非平穩(wěn)分析技術(shù)的發(fā)展,基于小波變換[8-9]、希爾伯特-黃變換[10-11]等時(shí)頻非平穩(wěn)地震動(dòng)模擬方法也逐漸興起并完善,具有良好的應(yīng)用前景。與頻域方法相比,時(shí)域方法無需復(fù)雜的頻譜分析,具有模擬所需存儲(chǔ)空間少、計(jì)算效率高的特點(diǎn),可直接在時(shí)域上生成大量的非平穩(wěn)地震動(dòng)。然而,無論是頻域方法還是時(shí)域方法,在進(jìn)行地震動(dòng)的隨機(jī)模擬時(shí),通常需要處理大量的隨機(jī)變量。例如,譜表示模型中隨機(jī)相位角的數(shù)量取決于諧波分量的數(shù)量,而時(shí)域過濾白噪聲模型[12-13]中包括的高斯隨機(jī)變量的數(shù)量取決于模擬時(shí)刻的數(shù)量,這導(dǎo)致模擬公式中需要成百上千的隨機(jī)變量。在這種情況下,通常采用Monte-Carlo方法進(jìn)行隨機(jī)抽樣,這導(dǎo)致了生成的地震動(dòng)樣本數(shù)量巨大以及概率信息不完備,難以應(yīng)用于工程結(jié)構(gòu)精細(xì)化的地震反應(yīng)和動(dòng)力可靠性分析[14]。為了減少非平穩(wěn)地震動(dòng)過程模擬時(shí)基本隨機(jī)變量的數(shù)量,Liu等[15-16]提出了基于譜表示以及基于Karhunen-Loeve分解的降維方法,實(shí)現(xiàn)了僅用1～2基本隨機(jī)變量即可描述原地震動(dòng)過程。另一方面,上述模型低估了地面運(yùn)動(dòng)的自然變異性,因?yàn)樵诖祟惸Ｐ椭?通常根據(jù)經(jīng)驗(yàn)將模型參數(shù)取為確定性值,而忽略了模型參數(shù)的不確定性[17],這將可能導(dǎo)致后續(xù)得到結(jié)構(gòu)動(dòng)力可靠度被高估。

基于上述研究進(jìn)展,本文擬在平穩(wěn)過程的時(shí)域表達(dá)理論基礎(chǔ)上,引入隨機(jī)函數(shù)的降維思想,建議一類全非平穩(wěn)地震動(dòng)過程的時(shí)域降維模型。為了進(jìn)一步捕獲地震動(dòng)時(shí)域降維模型中隨機(jī)參數(shù)的概率分布,基于實(shí)測(cè)強(qiáng)震記錄進(jìn)行隨機(jī)參數(shù)的樣本值識(shí)別,利用最大似然估計(jì)方法對(duì)幾種常用的概率分布進(jìn)行其參數(shù)擬合,并通過AIC(Akaike information criteria)準(zhǔn)則和K-S(Kolmogorov-Smirnov)檢驗(yàn)來確定最優(yōu)的概率分布。在地震動(dòng)過程的時(shí)域降維模擬時(shí),將時(shí)域降維表達(dá)中的基本隨機(jī)變量與模型參數(shù)中的隨機(jī)變量進(jìn)行統(tǒng)一處理,生成具有賦得概率的代表性點(diǎn)集,進(jìn)而得到具有相同賦得概率的地震動(dòng)代表性樣本集合,且代表性樣本能夠體現(xiàn)地震動(dòng)的自然變異性。同時(shí),地震動(dòng)過程的時(shí)域降維模型在本質(zhì)上能夠與概率密度演化理論相結(jié)合[18],進(jìn)而實(shí)現(xiàn)工程結(jié)構(gòu)的精細(xì)化動(dòng)力反應(yīng)與可靠度分析。

1 地震動(dòng)過程的時(shí)域降維模擬

1.1 平穩(wěn)過程時(shí)域表達(dá)的數(shù)學(xué)理論

任一零均值的實(shí)值平穩(wěn)過程X(t)可表示為如下的Fourier-Stieltjes積分形式[19]

(1)

式中,Z(ω)和U(ω)均為復(fù)的正交增量過程,它們的增量分別滿足如下條件

E[dZ(ω)dZ*(ω′)]=SX(ω)δωω′dω

(2)

E[dU(ω)dU*(ω′)]=δωω′dω

(3)

式中:E[·]為數(shù)學(xué)期望;“*”為取復(fù)共軛;δωω′為Kronecker符號(hào);SX(ω)為平穩(wěn)過程X(t)的雙邊功率譜密度函數(shù)。

在式(1)中,H(ω)是一個(gè)復(fù)值函數(shù),且滿足

|H(ω)|2=SX(ω)

(4)

式(1)即為平穩(wěn)過程X(t)的頻域表達(dá)形式。進(jìn)一步,若定義h(t)為H(ω)的Fourier逆變換,即

(5)

可以證明(見附錄A),平穩(wěn)過程X(t)的時(shí)域表達(dá)形式為

(6)

式中,W(τ)為一個(gè)實(shí)的正交增量過程,其增量滿足

E[dW(τ)dW(τ′)]=2πδττ′dτ

(7)

進(jìn)一步地,式(6)還可以寫成白噪聲表達(dá)的形式

(8)

(9)

式中,δ(t-t′)為Dirac函數(shù)。

式(8)表明,任意一個(gè)實(shí)值平穩(wěn)過程(二階矩過程)都可以表示為一個(gè)過濾白噪聲過程。而且,函數(shù)h(t)往往與線性時(shí)不變系統(tǒng)的脈沖響應(yīng)函數(shù)相關(guān)。因此,平穩(wěn)過程也可看作在白噪聲激勵(lì)下的線性時(shí)不變系統(tǒng)的反應(yīng)過程。

總之,式(1)和式(6)分別在頻域和時(shí)域上提供了平穩(wěn)過程的兩種等價(jià)表達(dá)形式。

1.2 地震動(dòng)過程的時(shí)域模型

式中:tk=kΔt,k=0,1,…,n;int(t/Δt)=n,Δt為時(shí)間離散步長。若T為地震動(dòng)持時(shí),且int(T/Δt)=N,則n≤N。

在微小時(shí)間段ti-1≤t≤ti內(nèi),假定h(t-τ)為一常數(shù),同時(shí)略去式(10)中的最后一項(xiàng)。于是,式(10)可以近似寫成

(11)

(12)

E[Vi]=0,E[ViVj]=δij

(13)

式中,δij為Kronecker符號(hào)。

對(duì)于全非平穩(wěn)地震動(dòng)加速度過程a(t),采用具有隨機(jī)參數(shù)的強(qiáng)度調(diào)制函數(shù)f(t),以及具有隨機(jī)和時(shí)變參數(shù)的函數(shù)h(t-τ)來描述時(shí)-頻全非平穩(wěn)特性。為此,根據(jù)式(12),全非平穩(wěn)地震動(dòng)加速度過程a(t)可表示為

(14)

其中,

(15)

在式(14)中,由于求和部分是一個(gè)具有時(shí)變頻率成分的標(biāo)準(zhǔn)化過程(單位方差過程)。因此,強(qiáng)度調(diào)制函數(shù)f(t)就完全控制了過程的強(qiáng)度非平穩(wěn)性,而脈沖響應(yīng)函數(shù)(濾波器)及其隨機(jī)時(shí)變參數(shù)則控制了過程的頻譜非平穩(wěn)性。由于強(qiáng)度非平穩(wěn)性和頻譜非平穩(wěn)性的完全分離描述,這為后續(xù)的參數(shù)識(shí)別提供了便利。

1.3 正交隨機(jī)變量集的降維表達(dá)

為了實(shí)現(xiàn)地震動(dòng)過程時(shí)域表達(dá)的降維模擬,根據(jù)隨機(jī)函數(shù)的降維思想,本文建議正交隨機(jī)變量集的隨機(jī)函數(shù)形式為三角正交基,即

(16)

式中:i,j=1,2,…,N/2;基本隨機(jī)變量Θ服從區(qū)間(0,2π]上的均勻分布;α為任意常數(shù),本文取α=π/3。

在式(16)中,i與j(i,j=1,2,…,N/2)之間存在某種確定性的一一映射關(guān)系。這種一一映射關(guān)系可采用MATLAB工具箱中的rand(‘state’,0)和randperm(N)函數(shù)來實(shí)現(xiàn),即它們的一一映射關(guān)系為i=temp(j)。需要指出的是,這一映射關(guān)系正是降維方法的一個(gè)充分條件。

2 地震動(dòng)時(shí)域模型的參數(shù)識(shí)別

2.1 模型的參數(shù)化

在式(14)中,強(qiáng)度調(diào)制函數(shù)f(t)采用三段式模型,即

(17)

式中:參數(shù)λf=(σmax,t1,t2,c,T);σmax為地震動(dòng)過程a(t)的標(biāo)準(zhǔn)差函數(shù)的最大值;t1和t2分別為強(qiáng)震平穩(wěn)段的起始和結(jié)束時(shí)刻;c為控制下降段衰減快慢的參數(shù);T為地震動(dòng)過程a(t)的持時(shí)。

對(duì)于平穩(wěn)地震動(dòng)的功率譜,采用胡聿賢-周錫元模型[20],即

(18)

式中:ωg和ξg分別為場地土的卓越圓頻率和阻尼比;ωh為控制地震動(dòng)低頻含量的參數(shù),可取ωh=2π/T;S0為地震動(dòng)的譜強(qiáng)度因子,由于式(15)的標(biāo)準(zhǔn)化,可取S0=1。

根據(jù)式(4)和式(18),可以給出復(fù)值函數(shù)H(ω)為

(19)

再根據(jù)式(5),并利用留數(shù)定理,可得脈沖響應(yīng)函數(shù)h(t)為

(21)

式中,α=ωh/ωg。

對(duì)于式(20),其脈沖響應(yīng)函數(shù)h(t-τ)的時(shí)變參數(shù)λh(τ)=[ωg(τ),ξg(τ)]。其中,ωg(τ)取為關(guān)于時(shí)間的指數(shù)函數(shù)[21],ξg(τ)取為常數(shù),即

(22)

式中,λω=(η0,γ)為場地土卓越圓頻率的待定參數(shù)。

綜上,地震動(dòng)時(shí)域模型的參數(shù)向量λa=(λf,λh)=(σmax,t1,t2,c,T,η0,γ,ξg)。

2.2 參數(shù)識(shí)別方法

對(duì)于強(qiáng)度調(diào)制函數(shù),其控制著地震動(dòng)過程的幅值與能量。因此,可用地震動(dòng)的累積能量來對(duì)參數(shù)λf進(jìn)行識(shí)別。對(duì)于第i條實(shí)測(cè)地震動(dòng)記錄ai(t),其累積能量可表示為[22]

(23)

式中,Ti為第i條實(shí)測(cè)強(qiáng)震記錄的持時(shí)。

對(duì)于非平穩(wěn)地震動(dòng)過程a(t),其累積能量的期望I(t;λf)為

(24)

根據(jù)式(23)和式(24),采用最佳平方逼近原則,即可對(duì)參數(shù)向量λf,i進(jìn)行識(shí)別

(25)

(26)

于是,根據(jù)式(27)便可對(duì)參數(shù)向量λω,i進(jìn)行識(shí)別

(27)

(28)

式中,Tu為反應(yīng)譜周期的上限,本文取Tu=6 s。

2.3 參數(shù)識(shí)別結(jié)果

在太平洋地震工程研究中心的NGA-West2地震動(dòng)數(shù)據(jù)庫中,按照如下原則篩選實(shí)測(cè)強(qiáng)震記錄:①斷層距離在10～100 km,以減少近場效應(yīng)的影響以及排除小強(qiáng)度的地震動(dòng);②實(shí)測(cè)強(qiáng)震記錄的矩震級(jí)應(yīng)大于6,以排除對(duì)結(jié)構(gòu)影響較小的地震動(dòng);③實(shí)測(cè)強(qiáng)震記錄的地表以下30 m內(nèi)平均剪切波速VS,30范圍在265～550 m/s[24],對(duì)應(yīng)于GB 18306—2015《中國地震動(dòng)參數(shù)區(qū)劃圖》[25]中的II類場地。

經(jīng)過篩選得到了31次地震的454條地震動(dòng)記錄。表1給出了本文篩選實(shí)測(cè)強(qiáng)震記錄的地震名稱、臺(tái)站數(shù)量、震級(jí)、斷層距范圍和VS,30范圍等信息。進(jìn)一步,對(duì)實(shí)測(cè)強(qiáng)震記錄進(jìn)行四階Butterworth濾波處理以及截取1%～99%的能量,以保證結(jié)果的可靠性。同時(shí),為了便于參數(shù)識(shí)別,將實(shí)測(cè)地震動(dòng)記錄的峰值加速度調(diào)整為200 cm/s2。

表1 選取的實(shí)測(cè)強(qiáng)震記錄基本信息Tab.1 Basic information of the measured strong motion records

以Taft Lincoln School臺(tái)站記錄的Kern County地震為例,采用上述的時(shí)域模型與參數(shù)識(shí)別方法對(duì)地震動(dòng)的調(diào)制函數(shù)參數(shù)和脈沖響應(yīng)函數(shù)參數(shù)進(jìn)行識(shí)別,識(shí)別結(jié)果如圖1。由圖1可知:對(duì)于累計(jì)能量曲線和累計(jì)向上穿零次數(shù),模型值與實(shí)測(cè)值擬合效果良好;對(duì)于加速度反應(yīng)譜,實(shí)測(cè)反應(yīng)譜也基本包括在樣本反應(yīng)譜之內(nèi),驗(yàn)證了本文參數(shù)識(shí)別方法的有效性。典型實(shí)例的參數(shù)擬合結(jié)果如表2所示。

圖1 典型強(qiáng)震記錄的參數(shù)識(shí)別Fig.1 Parameter identification of typical strong motion records

表2 典型強(qiáng)震記錄的參數(shù)識(shí)別結(jié)果

如2.1節(jié)所述,地震動(dòng)的時(shí)域模型中含有(σmax,t1,t2,c,T,η0,γ,ξg)共8個(gè)參數(shù)。對(duì)于所篩選的454條地震波,這8個(gè)參數(shù)的統(tǒng)計(jì)結(jié)果如表3所示。為便于分析,表中用ts=t2-t1代替了t2。由表3可知,參數(shù)σmax和ξg的變異性相對(duì)較小,其他參數(shù)的變異性均較大。同時(shí),考慮到持時(shí)T在模擬時(shí)可直接給定,因此下文僅對(duì)參數(shù)t1、ts、c、η0和γ等5個(gè)參數(shù)進(jìn)行概率密度擬合。此外,由表3可知卓越圓頻率函數(shù)的衰減指數(shù)γ的最小值為負(fù)值,這表明實(shí)測(cè)地震動(dòng)的主頻率并非總是衰減,也可能出現(xiàn)隨時(shí)間增大的情況。

表3 地震動(dòng)參數(shù)的統(tǒng)計(jì)結(jié)果Tab.3 Statistical results of ground motion parameters

下面分別采用7種概率分布模型對(duì)上述5個(gè)參數(shù)進(jìn)行擬合,即正態(tài)分布、對(duì)數(shù)正態(tài)分布、Gamma分布、Weibull分布、Gumbel分布、廣義極值分布和指數(shù)分布。通過最大似然估計(jì)方法估計(jì)分布參數(shù),并通過最小化AIC值[26]來確定最佳擬合模型。然后,再進(jìn)行K-S檢驗(yàn)[27],以判斷所獲得的模型是否能夠真實(shí)地表示地震動(dòng)參數(shù)的邊緣分布。K-S檢驗(yàn)的顯著性水平通常取5%。最后,得到參數(shù)t1、ts、c、η0和γ的最優(yōu)概率分布擬合結(jié)果如表4所示。

表4 地震動(dòng)參數(shù)的最優(yōu)概率分布擬合結(jié)果

圖2給出了參數(shù)c的最優(yōu)概率密度擬合結(jié)果。由圖2可知,擬合的概率模型能夠較好地描述參數(shù)c的分布特征。

圖2 參數(shù)c的最優(yōu)概率密度Fig.2 Optimal probability density of parameter c

3 數(shù)值算例

在本算例中,參數(shù)t1、ts、c、η0和γ視為隨機(jī)變量(概率分布見表4),則t2=t1+ts。為簡化,忽略參數(shù)之間的相關(guān)性。參數(shù)σmax和ξg取其均值,即σmax=68.98 cm/s2,ξg=0.38。模擬持時(shí)取其均值加一倍標(biāo)準(zhǔn)差,近似取為T=60 s。本算例的其他參數(shù)為:時(shí)間間隔Δt=0.01 s;代表性樣本數(shù)量為610。

圖3給出了本文方法生成的II類場地的代表性地震動(dòng)加速度樣本。從時(shí)域特性來看,3條代表性樣本的峰值加速度、強(qiáng)震平穩(wěn)段的起始和結(jié)束時(shí)刻以及下降段衰減快慢均有顯著的差異。從頻域特性來看,3條代表性樣本的頻率成分以及其變化特性也都不同。因此,模擬樣本在時(shí)域和頻域上均具有明顯的隨機(jī)性,反映了實(shí)測(cè)地震動(dòng)的自然變異性。

圖3 地震動(dòng)加速度代表性樣本Fig.3 Representative samples of ground motion acceleration

圖4分別給出了模擬的610條地震動(dòng)加速度代表性樣本集合的均值和標(biāo)準(zhǔn)差與其目標(biāo)值的對(duì)比結(jié)果?？梢?均值和標(biāo)準(zhǔn)差的模擬值與目標(biāo)值均對(duì)應(yīng)良好,驗(yàn)證了降維方法的有效性。值得說明的是,盡管式(18)的調(diào)制函數(shù)為三段式模型,但地震動(dòng)過程的標(biāo)準(zhǔn)差曲線仍是一個(gè)光滑的函數(shù),這是因?yàn)樵谡{(diào)制函數(shù)模型中參數(shù)t1、t2和c當(dāng)作了隨機(jī)變量。

圖4 地震動(dòng)加速度代表性樣本集合的均值和標(biāo)準(zhǔn)差Fig.4 Mean and standard deviation of representative samples set of ground motion acceleration

為了進(jìn)一步分析本文方法的模擬精度,表5分別給出了377、610和987條代表性樣本集合的模擬誤差,其中均值和標(biāo)準(zhǔn)差相對(duì)誤差的定義參見Liu等的研究?？梢?隨著代表性樣本數(shù)量的增加,均值和標(biāo)準(zhǔn)差相對(duì)誤差逐漸變小,體現(xiàn)了降維方法的收斂性;當(dāng)代表性樣本數(shù)量達(dá)到610條時(shí),模擬誤差均小于5%,滿足精度要求。

表5 地震動(dòng)過程的模擬誤差Tab.5 Simulation error of ground motion process 單位:%

圖5分別給出了本文方法模擬的地震動(dòng)加速度的反應(yīng)譜和傅里葉幅值譜與實(shí)測(cè)強(qiáng)震記錄的比較。可以看出,實(shí)測(cè)強(qiáng)震記錄的均值在模擬均值加減一倍標(biāo)準(zhǔn)差范圍內(nèi),且與模擬均值的擬合較為一致,驗(yàn)證了本文方法具有良好的工程特性。

圖5 模擬的加速度反應(yīng)譜及傅里葉幅值譜與實(shí)測(cè)記錄的比較Fig.5 Comparison of response spectrum and Fourier amplitude spectrum between simulated acceleration and measured records

4 結(jié) 論

從平穩(wěn)過程時(shí)域表達(dá)的基本理論出發(fā),導(dǎo)出了平穩(wěn)和非平穩(wěn)地震動(dòng)的時(shí)域降維表達(dá)?；趯?shí)測(cè)強(qiáng)震記錄,對(duì)地震動(dòng)時(shí)域模型的參數(shù)進(jìn)行了統(tǒng)計(jì)建模。通過算例驗(yàn)證了本文方法的有效性,主要結(jié)論如下:

(1)通過平穩(wěn)過程時(shí)域表達(dá)的推導(dǎo),闡明了平穩(wěn)過程在頻域和時(shí)域上的兩種等價(jià)表達(dá)形式。無論是頻域表達(dá)還是時(shí)域表達(dá),隨機(jī)過程均是通過一系列正交隨機(jī)變量所調(diào)制的確定性函數(shù)的線性疊加來表達(dá)的。

(2)通過正交隨機(jī)變量集的隨機(jī)函數(shù)構(gòu)造,實(shí)現(xiàn)了僅用一個(gè)基本隨機(jī)變量即可精細(xì)地表達(dá)時(shí)域模型的目的,極大地降低了地震動(dòng)過程的隨機(jī)度。

(3)建議了地震動(dòng)時(shí)域模型的參數(shù)識(shí)別方法,通過典型實(shí)測(cè)強(qiáng)震記錄驗(yàn)證了其有效性。根據(jù)所選的強(qiáng)震記錄,獲得了地震動(dòng)時(shí)域模型基本參數(shù)的概率密度函數(shù),完善了地震動(dòng)時(shí)域降維建模。

(4)本文方法生成的同一類場地的代表性地震動(dòng)加速度樣本在時(shí)域和頻域上均具有明顯的隨機(jī)性,能夠反映地震動(dòng)的自然變異性。

附錄A

式(6)和式(7)的證明以及原過程的相關(guān)函數(shù)推導(dǎo)如下。

根據(jù)式(5)可知

(A.1)

同時(shí),令

(A.2)

利用式(A.1)和式(A.2),式(1)可寫為

(A.3)

式(6)得證。

由式(A.2)可知

(A.4)

(A.5)

將式(A.5)乘以其共軛后求期望,可得

E[dW(τ)dW*(τ′)]=

(A.6)

式(7)得證。

于是,原過程的相關(guān)函數(shù)為

(A.7)