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

阮鑫鑫, 劉章軍, 姜云木

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

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

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

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

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

1.1 平穩過程時域表達的數學理論

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

(1)

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

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

(2)

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

(3)

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

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

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

(4)

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

(5)

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

(6)

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

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

(7)

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

(8)

(9)

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

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

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

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內,假定h(t-τ)為一常數,同時略去式(10)中的最后一項。于是,式(10)可以近似寫成

(11)

(12)

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

(13)

式中,δij為Kronecker符號。

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

(14)

其中,

(15)

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

1.3 正交隨機變量集的降維表達

為了實現地震動過程時域表達的降維模擬,根據隨機函數的降維思想,本文建議正交隨機變量集的隨機函數形式為三角正交基,即

(16)

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

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

2 地震動時域模型的參數識別

2.1 模型的參數化

在式(14)中,強度調制函數f(t)采用三段式模型,即

(17)

式中:參數λf=(σmax,t1,t2,c,T);σmax為地震動過程a(t)的標準差函數的最大值;t1和t2分別為強震平穩段的起始和結束時刻;c為控制下降段衰減快慢的參數;T為地震動過程a(t)的持時。

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

(18)

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

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

(19)

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

(21)

式中,α=ωh/ωg。

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

(22)

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

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

2.2 參數識別方法

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

(23)

式中,Ti為第i條實測強震記錄的持時。

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

(24)

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

(25)

(26)

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

(27)

(28)

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

2.3 參數識別結果

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

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

表1 選取的實測強震記錄基本信息Tab.1 Basic information of the measured strong motion records

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

圖1 典型強震記錄的參數識別Fig.1 Parameter identification of typical strong motion records

表2 典型強震記錄的參數識別結果

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

表3 地震動參數的統計結果Tab.3 Statistical results of ground motion parameters

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

表4 地震動參數的最優概率分布擬合結果

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

圖2 參數c的最優概率密度Fig.2 Optimal probability density of parameter c

3 數值算例

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

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

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

圖4分別給出了模擬的610條地震動加速度代表性樣本集合的均值和標準差與其目標值的對比結果。可見,均值和標準差的模擬值與目標值均對應良好,驗證了降維方法的有效性。值得說明的是,盡管式(18)的調制函數為三段式模型,但地震動過程的標準差曲線仍是一個光滑的函數,這是因為在調制函數模型中參數t1、t2和c當作了隨機變量。

圖4 地震動加速度代表性樣本集合的均值和標準差Fig.4 Mean and standard deviation of representative samples set of ground motion acceleration

為了進一步分析本文方法的模擬精度,表5分別給出了377、610和987條代表性樣本集合的模擬誤差,其中均值和標準差相對誤差的定義參見Liu等的研究。可見,隨著代表性樣本數量的增加,均值和標準差相對誤差逐漸變小,體現了降維方法的收斂性;當代表性樣本數量達到610條時,模擬誤差均小于5%,滿足精度要求。

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

圖5分別給出了本文方法模擬的地震動加速度的反應譜和傅里葉幅值譜與實測強震記錄的比較。可以看出,實測強震記錄的均值在模擬均值加減一倍標準差范圍內,且與模擬均值的擬合較為一致,驗證了本文方法具有良好的工程特性。

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

4 結 論

從平穩過程時域表達的基本理論出發,導出了平穩和非平穩地震動的時域降維表達。基于實測強震記錄,對地震動時域模型的參數進行了統計建模。通過算例驗證了本文方法的有效性,主要結論如下:

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

(2)通過正交隨機變量集的隨機函數構造,實現了僅用一個基本隨機變量即可精細地表達時域模型的目的,極大地降低了地震動過程的隨機度。

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

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

附錄A

式(6)和式(7)的證明以及原過程的相關函數推導如下。

根據式(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)得證。

于是,原過程的相關函數為

(A.7)