地震動過程的時域降維建模

阮鑫鑫, 劉章軍, 姜云木

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

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

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

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

1 地震動過程的時域降維模擬

1.1 平穩(wěn)過程時域表達(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符號;SX(ω)為平穩(wěn)過程X(t)的雙邊功率譜密度函數(shù)。

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

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

(4)

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

(5)

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

(6)

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

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

(7)

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

(8)

(9)

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

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

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

1.2 地震動過程的時域模型

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

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

(11)

(12)

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

(13)

式中,δij為Kronecker符號。

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

(14)

其中,

(15)

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

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

為了實(shí)現(xiàn)地震動過程時域表達(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)系正是降維方法的一個充分條件。

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

2.1 模型的參數(shù)化

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

(17)

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

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

(18)

式中:ωg和ξg分別為場地土的卓越圓頻率和阻尼比;ωh為控制地震動低頻含量的參數(shù),可取ωh=2π/T;S0為地震動的譜強(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。

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

(22)

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

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

2.2 參數(shù)識別方法

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

(23)

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

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

(24)

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

(25)

(26)

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

(27)

(28)

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

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

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

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

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

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

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

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

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

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

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

表4 地震動參數(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。模擬持時取其均值加一倍標(biāo)準(zhǔn)差,近似取為T=60 s。本算例的其他參數(shù)為:時間間隔Δt=0.01 s;代表性樣本數(shù)量為610。

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

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

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

圖4 地震動加速度代表性樣本集合的均值和標(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)差相對誤差的定義參見Liu等的研究?？梢?隨著代表性樣本數(shù)量的增加,均值和標(biāo)準(zhǔn)差相對誤差逐漸變小,體現(xiàn)了降維方法的收斂性;當(dāng)代表性樣本數(shù)量達(dá)到610條時,模擬誤差均小于5%,滿足精度要求。

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

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

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

4 結(jié) 論

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

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

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

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

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

附錄A

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

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

(A.1)

同時,令

(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)