胡進(jìn)軍，王中偉，張 輝，靳超越，胡 磊

(中國地震局工程力學(xué)研究所地震局地震工程與工程振動重點實驗室，黑龍江，哈爾濱 150080)

在結(jié)構(gòu)抗震設(shè)計中，重要建筑和不規(guī)則建筑通常需要通過動力時程分析得到結(jié)構(gòu)的非線性反應(yīng)特性，選擇不同的地震動作為結(jié)構(gòu)非線性動力分析的輸入荷載會得到不同的結(jié)果，而地震動輸入結(jié)構(gòu)前通常需要調(diào)整其反應(yīng)譜與設(shè)計反應(yīng)譜一致[1]，因此，如何獲得合理、真實且匹配目標(biāo)反應(yīng)譜的地震動時程，是減少結(jié)構(gòu)響應(yīng)不確定性的重要問題。目前相關(guān)地震動模擬方法主要有三角級數(shù)法[2?4]、經(jīng)驗格林函數(shù)法[5?6]、隨機(jī)有限斷層法[7? 8]和寬頻帶模擬方法[9]等，這些方法通常根據(jù)與震源相關(guān)的頻域參數(shù)合成人工地震動，難以全面考慮區(qū)域地震構(gòu)造背景、地殼介質(zhì)和場地特征等引起的差異[10?11]。以三角級數(shù)法為例，由于受強(qiáng)度包線的控制，該方法生成的地震動時程具有很強(qiáng)的相似性，不能針對不同區(qū)域合成考慮區(qū)域性差異的地震動。雖然目前考慮場地項合成地震動有傳遞函數(shù)法和四分之一波長法[12]，但這要求詳細(xì)、準(zhǔn)確的區(qū)域地質(zhì)資料。由于目前地震動模擬很難考慮地震發(fā)生全過程的所有因素，且合成的地震動缺乏區(qū)域性差異，這些方法構(gòu)造的人工地震動與真實地震動之間有較大的差別，在用于抗震設(shè)計或評估時會增大輸入地震動的不確定性。

基于實際強(qiáng)震記錄合成反應(yīng)譜匹配的地震動主要有頻域和時域方法。頻域方法[13? 15]首先利用一組特定形式的函數(shù)分解地震動時程，再通過調(diào)整這組函數(shù)的幅值使其反應(yīng)譜與目標(biāo)反應(yīng)譜一致，最后根據(jù)新的幅值重組得到新地震動時程，重復(fù)以上步驟最終實現(xiàn)反應(yīng)譜匹配的目的。與頻域方法相比，時域方法的匹配精度更高，且合成地震動相較原始地震動具有改變程度小的優(yōu)點[16]，因此時域上反應(yīng)譜匹配方法受到了廣泛的研究。KAUL[17]首先提出了時域上反應(yīng)譜匹配的方法，之后LILHANAND 和TSENG[18]將這種反應(yīng)譜匹配方法應(yīng)用于匹配多阻尼反應(yīng)譜，但是其采用的小波不能保持地震動的非平穩(wěn)性。ABRAHAMSON[19]借鑒前者的算法提出了新的小波函數(shù)，雖然這種調(diào)整保留了原始地震動的非平穩(wěn)性，但未考慮其對速度和位移時程的影響，需要通過基線校正來消除速度和位移時程的漂移。HANCOCK 等[20]通過調(diào)整小波函數(shù)消除了速度和位移時程的漂移，由于其小波形式中包含了基線校正，導(dǎo)致其無法得到解析解，并影響計算速度和效率。AL ATIK和ABRAHAMSON[16]提出了一種不產(chǎn)生速度和位移時程漂移的調(diào)整小波函數(shù)，從而使反應(yīng)譜匹配過程中能夠得到小波系數(shù)的解析解。張郁山和趙鳳新[21]構(gòu)造了一種增量位移小波函數(shù)，通過微分、時移以及線性調(diào)幅得到加速度小波函數(shù)。ADEKRISTI 等[22]利用Broyden 修正方法來降低反應(yīng)譜匹配中由于不同小波之間相互影響所產(chǎn)生的非線性問題。HONG 和HUANG[23]每次只匹配單個周期的反應(yīng)譜值，避免了由多個小波函數(shù)引起的非線性問題。ZENGIN 和ABRAHAMSON[24]將反應(yīng)譜匹配方法應(yīng)用于匹配瞬時功率譜，并得到能夠同時匹配反應(yīng)譜和瞬時功率譜的地震動，從而為合成近場脈沖型地震動提供了新方法。雖然這些時域上地震動合成方法的原理基本相同，但是用于調(diào)整反應(yīng)譜的小波函數(shù)在改進(jìn)，解決反應(yīng)譜匹配中非線性問題的方法也在改進(jìn)。

目前頻域和時域合成地震動的方法均有缺陷。頻域方法在調(diào)整函數(shù)幅值時改變了地震動的非平穩(wěn)特征[16]，且反應(yīng)譜匹配的收斂性較差。時域方法雖然能保留地震動的頻譜特征，但由于引入反應(yīng)譜峰值時間不變的假設(shè)來簡化小波幅值的求解過程，將非線性問題用線性方程組近似求解，從而忽略了小波的相互作用以及地震動與其反應(yīng)譜之間復(fù)雜的關(guān)系，導(dǎo)致這類方法存在迭代效率低、收斂性較差和匹配精度低等問題。此外，由于地震的產(chǎn)生受震源效應(yīng)、傳播路徑和場地效應(yīng)的影響[25?26]，且地震動存在隨機(jī)不確定性[26]，頻域和時域方法都無法合成出合理考慮區(qū)域特征的地震動。

相較于頻域合成的地震動，時域方法合成的地震動的非線性動力反應(yīng)與實際地震動的結(jié)果更為相似[27]，因此，本文基于時域方法合成地震動。為了實現(xiàn)考慮區(qū)域特征、降低不確定性且匹配反應(yīng)譜的目標(biāo)，本文提出了一套完整且高效的地震動合成新方法。該方法首先基于主成分分析算法得到包含目標(biāo)區(qū)域地震動特征的種子地震動，再基于粒子群算法優(yōu)化時域上反應(yīng)譜匹配方法，用此方法調(diào)整種子地震動的反應(yīng)譜，最終得到可模擬目標(biāo)區(qū)域地震的合成地震動。為了驗證該方法的可行性、合理性和有效性，本文基于中國四川地區(qū)的地震動數(shù)據(jù)，給出了該方法的實際應(yīng)用案例與分析結(jié)果。

1 方法

1.1 主成分分析提取種子地震動

主 成 分 分 析[28](Principal components analysis,PCA)是一種降低數(shù)據(jù)維度的算法，該算法通常用于去除高維數(shù)據(jù)的噪聲和提取重要特征[29]。高維數(shù)據(jù)在某一維度上特征的重要程度用數(shù)據(jù)投影到該維度后的方差衡量，方差越大則數(shù)據(jù)的特征越重要。PCA 算法將m維數(shù)據(jù)映射到n維上(m>n)，不僅要讓數(shù)據(jù)在這n維上對應(yīng)的方差最大，且這些維度之間相互正交，為此該算法首先找出使數(shù)據(jù)投影的方差最大的坐標(biāo)軸，然后從與之前所有坐標(biāo)軸正交的超平面中找出使數(shù)據(jù)投影的方差最大的坐標(biāo)軸，以此類推，最終找出n個坐標(biāo)軸并完成數(shù)據(jù)降維。將這一方法應(yīng)用于地震工程，能從一組地震動數(shù)據(jù)中得到包含原始數(shù)據(jù)特征的種子地震動。

主成分分析方法得到種子地震動的具體步驟如下：

1)將m條地震動數(shù)據(jù)按行組成矩陣X。

2)對矩陣X按行進(jìn)行零均值化處理。

3)計算矩陣X的協(xié)方差矩陣C，并求出矩陣C的特征值與特征向量：

4)將矩陣C的特征值與特征向量按特征值從大到小排列，取前n個特征向量按順序組成矩陣P，則矩陣P右乘矩陣X即可得到新地震動數(shù)據(jù)—矩陣Y。

5)矩陣Y中有n條種子地震動，種子地震動的主成分貢獻(xiàn)率(Contributing Rate of Principal Component, CRP)越大，其包含的原始數(shù)據(jù)的區(qū)域地震動特征[30](eigenquakes)越明顯，主成分貢獻(xiàn)率的計算公式如下：

其中，λi為矩陣C的第i個特征值(降序)。

時域上反應(yīng)譜匹配方法的基本思路是在地震動的加速度時程上添加一組幅值很小的小波，使其反應(yīng)譜值與目標(biāo)反應(yīng)譜值之間的差值減小，通過迭代計算實現(xiàn)反應(yīng)譜匹配；求解小波幅值是一個非線性問題，為了用線性方程解決該問題，假設(shè)在地震動加速度時程上添加小波后，自振周期為Ti的單自由度結(jié)構(gòu)的加速度達(dá)到峰值的時間ti不變。

反應(yīng)譜匹配的方法如下：

1)計算周期為Ti時的反應(yīng)譜差值，公式如下：

式中：n為周期的個數(shù)；Qi為目標(biāo)反應(yīng)譜值；Ri為地震動a(t)的反應(yīng)譜值；Pi由結(jié)構(gòu)加速度峰值確定，當(dāng)加速度峰值大于零時，Pi=1，否則Pi=?1。

2)在a(t)上疊加小波函數(shù)[16]fj(t) 以消除反應(yīng)譜差值Ej，此時Ri的改變量為cij：

式中：hi(t)為單自由度結(jié)構(gòu)的加速度脈沖響應(yīng)函數(shù)[16]；ti為加速度達(dá)到峰值的時間。

3)當(dāng)周期為Ti時，令添加的所有小波對反應(yīng)譜值的總改變量等于反應(yīng)譜差值Ei，則：

1.2 反應(yīng)譜匹配方法

其中，bj為小波函數(shù)fj(t)的幅值。

4)根據(jù)式(5)可求解出所有小波的幅值，則新地震動a1(t) 的計算公式如下：

其中，γ 為在0 到1 之間的縮放系數(shù)[16]。

用a1(t)代替a(t)，重復(fù)以上步驟，直到反應(yīng)譜之間的差值小于設(shè)定差值為止。

1.3 基于粒子群算法的優(yōu)化原理

目前時域上反應(yīng)譜匹配方法并沒有對式(6)中γ 的取值做出明確的說明，因此現(xiàn)有方法普遍存在迭代效率低、反應(yīng)譜匹配精度低的問題。本文采用機(jī)器學(xué)習(xí)中的粒子群算法[31](Particle swarm optimization, PSO)來尋找每次迭代過程中γ 的最優(yōu)解，由于PSO 算法能夠共享所有粒子的信息來調(diào)整粒子在尋優(yōu)過程中的速度和位置，采用該算法能快速找到問題的最優(yōu)解，從而提高反應(yīng)譜匹配方法的迭代效率與匹配精度。

粒子群算法的原理如下：

1)參數(shù)設(shè)置與種群初始化。

設(shè)置相關(guān)參數(shù)，例如：粒子的數(shù)目、搜索范圍和速度范圍以及迭代次數(shù)等，根據(jù)設(shè)置的參數(shù)隨機(jī)確定粒子的位置和速度。

2)計算極值，確定最優(yōu)解。

根據(jù)定義的適應(yīng)度函數(shù)計算粒子的適應(yīng)度值，找出每個粒子的個體極值，再從所有個體極值中找出群體極值作為最優(yōu)解。

3)更新粒子的速度和位置，公式如下：

式中：i=1, 2, ···,N，N為選取的粒子總數(shù)；t為當(dāng)前迭代步數(shù)；vi為粒子的速度；xi為粒子的當(dāng)前位置；w為慣性因子(非負(fù))，用于調(diào)節(jié)算法的尋優(yōu)能力；r1和r2是0 到1 之間的隨機(jī)數(shù)且相互獨立；c1、c2是學(xué)習(xí)因子；pbesti為個體最優(yōu)解；gbest為全局最優(yōu)解。

本文的慣性因子采用線性遞減權(quán)值策略，公式如下：

式中：g為迭代次數(shù)；w(g)為第g次迭代的慣性因子；Gk為最大迭代次數(shù)；wini為初始慣性權(quán)值；wend為最終慣性權(quán)值。

4)判斷迭代次數(shù)與最優(yōu)解是否符合要求。

當(dāng)達(dá)到最大迭代次數(shù)或者最優(yōu)解滿足設(shè)定值時，迭代結(jié)束，否則返回步驟2。

在上述原理中，單個粒子的優(yōu)劣要通過適應(yīng)度函數(shù)來評價。本文構(gòu)造了四種誤差函數(shù),其中任意一種均可作為適應(yīng)度函數(shù)，從而滿足不同情況下對反應(yīng)譜匹配精度的要求。四種誤差函數(shù)的計算公式如表1 所示。

表1 誤差函數(shù)的定義Table 1 Definition of error function

式中：N為選取周期點的個數(shù)；Ti為第i個自振周期；Sac(Ti)為周期為Ti時新地震動的反應(yīng)譜值，Sat(Ti)為周期為Ti時目標(biāo)反應(yīng)譜值。

本文應(yīng)用粒子群算法求反應(yīng)譜匹配過程中γ 最優(yōu)解的流程圖如圖1。

圖1 粒子群算法求最優(yōu)解γ 的流程圖Fig. 1 Flow chart of particle swarm optimization for solving optimal solution γ

2 地震動合成整體思路

由上述方法可知本文合成地震動的整體思路分為兩步：首先從目標(biāo)區(qū)域選取一組區(qū)域地震動，用PCA 算法從中得到包含區(qū)域地震動特征的種子地震動；再用基于PSO 算法優(yōu)化的時域內(nèi)反應(yīng)譜匹配方法調(diào)整種子地震動，通過疊加小幅值的小波調(diào)整其反應(yīng)譜與目標(biāo)反應(yīng)譜一致，最終得到合成地震動。整體流程圖如圖2 所示。

圖2 本文的地震動合成流程圖Fig. 2 The flow chart of ground motion synthesis

PCA 算法得到種子地震動的優(yōu)點：1)降低了地震動的不確定性，地震動在時域和頻域上的不確定性會導(dǎo)致結(jié)構(gòu)地震響應(yīng)也具有不確定性[32]，而基于最大方差理論的PCA 算法能從一組地震動中得到不確定性更低的種子地震動；2)具有區(qū)域特征，由于PCA 算法能找到一組數(shù)據(jù)的特征并得到包含此特征的新數(shù)據(jù)，因此，如果選擇的一組地震動具有區(qū)域特征，則種子地震動也具有區(qū)域特征。此外，從目標(biāo)區(qū)域的地震動數(shù)據(jù)庫中挑選地震動應(yīng)遵循震級、震中距和場地條件等與設(shè)定條件基本相符的原則，則用PCA 算法得到的種子地震動也滿足這三個條件的約束。正是由于種子地震動具有上述兩個優(yōu)點且滿足一定的約束條件，用它來模擬設(shè)定地震在目標(biāo)區(qū)域可能產(chǎn)生的地震動是合理的。

由于種子地震動的反應(yīng)譜一般與設(shè)定區(qū)域的目標(biāo)反應(yīng)譜不匹配，因此，還需要用優(yōu)化的反應(yīng)譜匹配方法調(diào)整種子地震動的反應(yīng)譜與目標(biāo)反應(yīng)譜一致。由于該方法不會對種子地震動的時域和頻域特征產(chǎn)生較大的改變，因此，最后合成的地震動可用于模擬設(shè)定地震在目標(biāo)區(qū)域產(chǎn)生的地震動。

AL ATIK 和ABRAHAMSON[16]提出的反應(yīng)譜匹配的改進(jìn)方法為了保證迭代的收斂性，在單次匹配過程中使匹配周期從短周期逐漸擴(kuò)大到全周期，且迭代過程中未優(yōu)化縮放系數(shù)γ 的取值，因此，該方法的迭代效率和匹配精度相對較低。用PSO 算法優(yōu)化該方法后，求解γ 的最優(yōu)值能提高迭代的收斂性，同時，全周期匹配也避免了多次疊加小波降低匹配精度的問題，因此，本文的優(yōu)化方法提高了迭代的收斂性、效率和匹配精度。

3 方法的驗證

為了驗證本文提出的基于實際記錄的地震動時域合成方法的可行性，以四川地區(qū)的地震動為基礎(chǔ)進(jìn)行目標(biāo)區(qū)域的地震動構(gòu)建。本文的反應(yīng)譜匹配方法僅考慮單一阻尼，下文算例中所有反應(yīng)譜對應(yīng)的阻尼比均為5%。

3.1 基于PCA 方法獲取種子地震動

假定需要合成目標(biāo)區(qū)域為四川地區(qū)、震級5 級、震中距30 km、場地條件為土層的地震動，為了體現(xiàn)區(qū)域地震動特征，選取34 條四川雅安余震的地震動記錄，以此合成雅安地區(qū)的種子地震動，這些地震動記錄的震級為4.5～5.4、震中距為25 km～35 km，詳細(xì)的地震信息如表2 所示。

表2 區(qū)域地震動信息Table 2 Regional ground motion information

通過PCA 算法的分析計算，可依據(jù)主成分貢獻(xiàn)率從大到小選取對應(yīng)的種子地震動，本文選取對應(yīng)CRP最大的4 條地震動作為種子地震動，依次將其記為地震動1～4，其加速度時程圖如圖3 所示。

圖3 種子地震動的加速度時程Fig. 3 The acceleration time history of the seed ground motion

3.2 優(yōu)化反應(yīng)譜匹配方法合成地震動

為了使目標(biāo)反應(yīng)譜具有四川地區(qū)的地震特性，本文選擇了基于實際數(shù)據(jù)得到的川滇地區(qū)地震動預(yù)測模型[33]，根據(jù)設(shè)定的震級和震中距即可得到目標(biāo)反應(yīng)譜(控制周期點見參考文獻(xiàn)[33]，不含PGA)。如圖4 所示，灰色曲線是目標(biāo)反應(yīng)譜，其他曲線是4 條種子地震動的反應(yīng)譜。

圖4 種子地震動反應(yīng)譜與目標(biāo)反應(yīng)譜Fig. 4 Seed ground motion response spectrum and target response spectrum

為了保證收斂性，種子地震動需要適當(dāng)?shù)恼{(diào)幅和補(bǔ)零[16]。由于迭代效率的提高，本文設(shè)置迭代次數(shù)為20 次即可使反應(yīng)譜差值降到很低的水平，實際應(yīng)用中也可以預(yù)設(shè)合適的匹配精度。為了找出γ 的最優(yōu)解，PSO 算法的粒子數(shù)和迭代次數(shù)不宜過低，本文設(shè)定粒子數(shù)為20 個、迭代次數(shù)為5 次。為了使合成地震動的反應(yīng)譜與目標(biāo)譜之間的差值更小，本文選擇反應(yīng)譜匹配差值的總體平均值(ERR1mean)作為粒子群算法的適應(yīng)度函數(shù)。由種子地震動調(diào)整得到的4 條合成地震動的加速度時程如圖5 所示。

圖5 合成地震動的加速度時程Fig. 5 The acceleration history of the synthetic ground motion

圖7 是4 條合成地震動在控制周期點上ERR1和ERR2的曲線圖，可以看出在所有周期上ERR1不超過0.045 cm/s2，ERR2不超過0.3%，因此優(yōu)化的反應(yīng)譜匹配方法能得到匹配精度很高的合成地震動，至此合成地震動方法的可行性得到了驗證。

圖6 匹配誤差隨迭代次數(shù)增加的變化曲線Fig. 6 The variation curve of matching error with the increase of iteration number

3.3 合成地震動的結(jié)果分析

從圖7 可以看出，合成地震動的反應(yīng)譜與目標(biāo)譜基本一致，但在匹配反應(yīng)譜時種子地震動的時程需要調(diào)幅和疊加小波函數(shù)，因此要從時域和頻域上分析此過程是否改變了地震動的區(qū)域特征。由于種子地震動調(diào)幅為初始地震動的過程不會改變區(qū)域特征，且小波函數(shù)均疊加在初始地震動上，因此分析比較初始地震動與合成地震動即可確定后者是否仍具有區(qū)域特征，為此選擇種子地震動1 對應(yīng)的初始地震動和合成地震動為本節(jié)的分析對象。

圖7 合成地震動的匹配誤差Fig. 7 Matching error of the synthetic ground motion

如圖8 所示，初始地震動和合成地震動的加速度和速度時程非常相似，兩者的位移時程基本相同，且最大位移相差不大；從圖9 可以看出初始地震動和合成地震動的歸一化Arias 強(qiáng)度曲線基本一致，且兩者的顯著持時D70和D90相差不超過0.024 s，故兩者的能量分布是一致的，因此這兩條地震動在時域上的特征基本沒有差別。

圖8 初始地震動和合成地震動的時程Fig. 8 Time history of initial and synthetic ground motion

圖9 初始地震動和合成地震動的歸一化Arias 強(qiáng)度Fig. 9 Normalized Arias intensities for initial and synthetic ground motion

如圖10 所示，兩條地震動的傅里葉幅值譜基本一致。由于地震動是非平穩(wěn)信號，因此還需要用短時傅里葉變換分析地震動的時頻特征。圖11分別是兩條地震動的時頻圖，與圖11(a)相比，圖11(b)在2 s～5 s 處存在2 Hz～8 Hz 的低頻成分減少、在3 s 附近存在大于10 Hz 的高頻成分增加的現(xiàn)象，這是由于種子地震動在不同周期下達(dá)到加速度峰值的時間ti基本在2 s～5 s，使得疊加小波后這一段時程的改變較大，但從整體上可以看出兩條地震動的時頻特征基本相同，因此兩者在頻域上的特征基本沒有差別。

圖10 初始地震動和合成地震動的傅里葉幅值譜Fig. 10 Fourier amplitude spectra of initial and synthetic ground motion

圖11 初始地震動和合成地震動的時頻圖Fig. 11 Time-frequency diagrams of initial and synthetic ground motion

本節(jié)從時域和頻域上對合成地震動和初始地震動進(jìn)行比較，得出兩者的相關(guān)特征基本相同的結(jié)論，而初始地震動保留了種子地震動中的區(qū)域特征，因此合成地震動也具有設(shè)定區(qū)域的地震動特征，同時其反應(yīng)譜與目標(biāo)譜一致，所以用本文合成方法得到的地震動模擬目標(biāo)區(qū)域可能產(chǎn)生的地震是合理的。

4 結(jié)論與討論

本文基于實際地震動記錄并結(jié)合機(jī)器學(xué)習(xí)算法得到具有區(qū)域特征且與目標(biāo)譜一致的合成地震動，本文的主要結(jié)論如下：

(1)提出了一種時域上合成地震動的方法，該方法先由PCA 算法獲得包含目標(biāo)區(qū)域地震動特征的種子地震動，再根據(jù)PSO 算法優(yōu)化的反應(yīng)譜匹配方法調(diào)整種子地震動，使其反應(yīng)譜與目標(biāo)譜一致，最終得到可用于模擬目標(biāo)區(qū)域地震的地震動。

(2)與現(xiàn)有的方法相比，用PCA 算法降維處理地震動數(shù)據(jù)時不僅降低了地震動的不確定性，同時保留了地震動的區(qū)域特征，因此相對于選擇實際地震動記錄或是人工地震動，用種子地震動模擬目標(biāo)區(qū)域地震更為合理。

(3)用PSO 算法對目前已有的時域上反應(yīng)譜匹配方法優(yōu)化，從而在每次反應(yīng)譜匹配過程中獲取用于約束小波函數(shù)的最優(yōu)縮放系數(shù)γ，此方法大幅提高了反應(yīng)譜匹配過程中的迭代效率和匹配精度，并使地震動高精度匹配目標(biāo)譜成為可能。

由于獲取實際地震動特征需要目標(biāo)區(qū)域?qū)嶋H地震動記錄，對于強(qiáng)震數(shù)據(jù)不豐富的地區(qū)可從地震地質(zhì)背景相似地區(qū)的地震動記錄中提取種子地震動。本文考慮區(qū)域?qū)嶋H地震動特征的構(gòu)建方法為解決面向區(qū)域的地震動場構(gòu)建提供了參考。