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

(中國(guó)地震局工程力學(xué)研究所地震局地震工程與工程振動(dòng)重點(diǎn)實(shí)驗(yàn)室，黑龍江，哈爾濱 150080)

隨著城市化進(jìn)程的加速，城市的規(guī)模越來(lái)越大，并且逐漸向地震區(qū)延伸，龐大復(fù)雜的城市基礎(chǔ)設(shè)施增加了潛在的地震危險(xiǎn)性。地震危險(xiǎn)性分析和評(píng)估一直是地震工程研究的關(guān)鍵問題[1 ? 3]，合理、準(zhǔn)確的地震動(dòng)輸入是進(jìn)行震前的地震風(fēng)險(xiǎn)評(píng)估和震后的地震災(zāi)害損失評(píng)估的前提。然而由于地震受到諸多因素的影響，震源、路徑以及場(chǎng)地的未知性和不確定性一直是模擬和合成地震動(dòng)中難以解決的問題。

確定工程場(chǎng)地地震動(dòng)的方法有多種，基于概率的地震危險(xiǎn)性分析(PSHA)方法是中常用的一種，其從概率的角度定量描述了地震作用[4 ? 5]?；诘卣饎?dòng)預(yù)測(cè)模型(GMM)[6 ? 7]，利用PSHA 方法可以建立一套具有不同超越概率的危險(xiǎn)譜，作為確定設(shè)防水準(zhǔn)的基礎(chǔ)。但是僅反應(yīng)譜還不能滿足結(jié)構(gòu)的抗震設(shè)計(jì)和性能評(píng)估所需，峰值加速度(PGA)以及譜加速度(Sa)也常被用來(lái)表征地震動(dòng)的強(qiáng)度指標(biāo)[8 ? 9]，一些其他的地震動(dòng)參數(shù)如峰值速度(PGV)、峰值位移(PGD)、有效峰值加速度(EPV)、累積絕對(duì)速度(CAV)也在逐漸引起人們的重視。但是地震動(dòng)信號(hào)是一種非平穩(wěn)的隨機(jī)過(guò)程，難以僅通過(guò)幅值或者能量參數(shù)就能夠完全包含其復(fù)雜過(guò)程對(duì)結(jié)構(gòu)的影響。因此地震動(dòng)時(shí)程是結(jié)構(gòu)非線性動(dòng)力時(shí)程分析的必須，這個(gè)所需的時(shí)程可以通過(guò)人造或者對(duì)天然的地震動(dòng)調(diào)幅得到[10]。除了基于PSHA 確定輸入地震動(dòng)的方法以外，基于物理過(guò)程的地震動(dòng)模擬也是常用的方法，其可以考慮地震斷層的破裂和傳播過(guò)程，具有物理機(jī)制[11]，但是這種確定性的模擬方法依賴于震源模型、速度結(jié)構(gòu)模型以及采用的模擬方法，需要準(zhǔn)確的震源、路徑和場(chǎng)地的模型，而目前工程上很難快速準(zhǔn)確確定這些參數(shù)，因此，模擬的結(jié)果存在很大的不確定性。

近年來(lái)全球強(qiáng)震臺(tái)站逐漸增多，區(qū)域強(qiáng)震數(shù)據(jù)日益豐富，我國(guó)的強(qiáng)震數(shù)據(jù)記錄的數(shù)目已經(jīng)超過(guò)了3 萬(wàn)條，而全球的強(qiáng)震數(shù)據(jù)已超過(guò)數(shù)十萬(wàn)條[12 ? 13]。這為研究基于區(qū)域強(qiáng)震數(shù)據(jù)和數(shù)值模擬方法結(jié)合構(gòu)建符合本區(qū)域地震構(gòu)造、路徑和場(chǎng)地特征的地震動(dòng)提供了可能。但是面對(duì)數(shù)萬(wàn)計(jì)的地震動(dòng)數(shù)據(jù)，如果采用人工的方法處理和挖掘?qū)⒎浅?fù)雜和繁瑣，在以往基于地震動(dòng)數(shù)據(jù)的統(tǒng)計(jì)研究中，研究者一般先要篩選一部分?jǐn)?shù)據(jù)，再針對(duì)選取的地震動(dòng)數(shù)據(jù)進(jìn)行分析，以減少計(jì)算量。隨著計(jì)算機(jī)技術(shù)的發(fā)展以及機(jī)器學(xué)習(xí)理論的逐漸成熟，應(yīng)用機(jī)器學(xué)習(xí)方法對(duì)海量地震動(dòng)數(shù)據(jù)進(jìn)行處理和挖掘成為可能。

為了探討應(yīng)用機(jī)器學(xué)習(xí)中的智能算法基于實(shí)際地震動(dòng)合成目標(biāo)地震動(dòng)，本文采用機(jī)器學(xué)習(xí)方法中的主成分分析(PCA)算法，從目標(biāo)區(qū)域地震動(dòng)數(shù)據(jù)庫(kù)中提取包含區(qū)域特征信息的地震動(dòng)母波，同時(shí)基于目標(biāo)地區(qū)的GMM 得到給定場(chǎng)地的加速度反應(yīng)譜，基于特征母波和設(shè)計(jì)譜構(gòu)建包含本地地震動(dòng)特征的地震動(dòng)時(shí)程。計(jì)算過(guò)程中，為了改進(jìn)地震動(dòng)時(shí)程的合成算法并提高計(jì)算效率，本文利用粒子群算法(PSO)快速找到母波地震動(dòng)的權(quán)重系數(shù)，使得合成的地震動(dòng)加速度反應(yīng)譜與GMM 得到的目標(biāo)譜誤差最小，最終通過(guò)母波線性組合得到目標(biāo)地震動(dòng)時(shí)程。為了闡述本文方法的可行性和合理性，本文結(jié)合我國(guó)西部地區(qū)的中強(qiáng)強(qiáng)震數(shù)據(jù)開展研究。

1 區(qū)域強(qiáng)地震動(dòng)數(shù)據(jù)

我國(guó)西部地區(qū)7 級(jí)以下的中強(qiáng)震數(shù)據(jù)相對(duì)豐富，但是7 級(jí)～8 級(jí)的大震事件仍然比較缺乏，在預(yù)測(cè)模型回歸時(shí)缺少大震強(qiáng)震數(shù)據(jù)的約束，為了減少模型的不確定性。本文采用了2007 年?2019 年間四川以及周邊省份的震級(jí)5.0 級(jí)～7.0 級(jí)的中強(qiáng)震地震動(dòng)數(shù)據(jù)，以及基于此數(shù)據(jù)建立的中強(qiáng)震地震動(dòng)預(yù)測(cè)模型。詳細(xì)的地震信息如表1 所示，震中和臺(tái)站位置見圖1。數(shù)據(jù)庫(kù)中包含了21 次地震中174 個(gè)臺(tái)站的4551 條水平向地震動(dòng)記錄。對(duì)原始地震動(dòng)記錄進(jìn)行了濾波和基線調(diào)整[14 ? 15]。

表1 選取的西部地區(qū)的中強(qiáng)震Table 1 Selected earthquake events in west region of China

圖1 臺(tái)站和地震震中分布圖Fig.1 Map of stations and earthquake epicenters

2 區(qū)域地震動(dòng)信息的提取

研究表明：地震動(dòng)由于受到地震構(gòu)造、地殼結(jié)構(gòu)和場(chǎng)地條件的影響，不同區(qū)域地震動(dòng)可能具有不同的特征[15 ? 17]。在模擬設(shè)定區(qū)域的地震動(dòng)時(shí)需要考慮本區(qū)域的實(shí)際地震動(dòng)的特征信息，因此，需要采用合理的方法從實(shí)際地震動(dòng)中提取區(qū)域地震動(dòng)的特征信息。

在數(shù)據(jù)挖掘和機(jī)器學(xué)習(xí)中，數(shù)據(jù)一般被表示為向量，與之類似，也可以把一條地震動(dòng)記錄視為1 列向量，那么n 條地震動(dòng)記錄就可視為n 列向量進(jìn)而組合得到如下矩陣：

A=[ ?→α1?→α2?→α3··· ?→αn]m×n(1)

這樣得到的合成地震動(dòng)具備了原地震動(dòng)所有的形狀特征以及隨機(jī)性特性。由于地震動(dòng)數(shù)據(jù)庫(kù)的記錄數(shù)目較大，為了提交計(jì)算效率可以對(duì)原始的地震動(dòng)矩陣A 進(jìn)行降維簡(jiǎn)化。在線性代數(shù)中，一個(gè)內(nèi)積空間的正交基是元素兩兩正交的基。在二維平面中，任意的二維向量都可以通過(guò)一組二維的正交基表示出來(lái)：

在三維平面中，任意的三維向量都可以通過(guò)一組三維的正交基表示：

當(dāng)把一條條地震動(dòng)向量視為一個(gè)個(gè)列向量時(shí)，那么肯定也存在一組正交基能夠表示任意一條地震動(dòng)所構(gòu)成列向量。主成分分析[18](Principal component analysis，PCA)方法正好可以滿足這方面的要求。它可以把數(shù)據(jù)降維，找出一組符合條件的正交基用于計(jì)算任意一條地震動(dòng)記錄。PCA 算法是一種對(duì)高維數(shù)據(jù)降維的方法，并將高維數(shù)據(jù)中重要的特征保留，去除噪聲和不重要的特征。

若有一組如下形式的數(shù)據(jù)，應(yīng)用主成分分析的具體步驟為：

首先，對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理：

其次，計(jì)算相關(guān)系數(shù)矩陣：

接著，用雅克比方法求解出相關(guān)系數(shù)矩陣R的特征值λ 和特征向量，這里的特征向量就是正交基。

最后，選擇重要的主成分，根據(jù)方差解釋率即：

這里的方差解釋率也稱主成分貢獻(xiàn)率，用于判斷單個(gè)主成分所包含的原始數(shù)據(jù)信息的多少，方差解釋率越大，所包含的原始信息越多因此本文基于此方法，從目標(biāo)區(qū)域原始地震動(dòng)數(shù)據(jù)庫(kù)中提取含有本地地震動(dòng)特征信息的母波。

2.1 基于PCA 算法的地震動(dòng)母波提取

以數(shù)據(jù)為驅(qū)動(dòng)提取母波的方法的具體流程如圖2 所示。首先從區(qū)域原始數(shù)據(jù)庫(kù)中應(yīng)用PCA 提取一組標(biāo)準(zhǔn)的正交基向量，并要求這些提取出的正交基向量能夠表征地震動(dòng)時(shí)程序列的主要成分。

圖2 主成分分析算法提取地震動(dòng)母波的流程圖Fig.2 The flow chart of ground motion mother wave extraction by PCA

應(yīng)用PCA 算法提取的地震動(dòng)主要成分的正交基在本文稱為地震動(dòng)母波，提取的地震動(dòng)母波和原始地震動(dòng)具有相同的時(shí)間采樣頻率。因此，合成的地震動(dòng)可以由提取的地震動(dòng)母波線性組合而成。

式中： ki為系數(shù)； ui為提取的地震動(dòng)母波； n為提取的地震動(dòng)母波的個(gè)數(shù)。地震動(dòng)母波是數(shù)據(jù)矩陣組成的特征向量，然后根據(jù)特征值大小進(jìn)行排序。

采用上述方法，可以從原始地震動(dòng)數(shù)據(jù)庫(kù)中提取n 條地震動(dòng)母波，圖3 給出了提取的4 條母波，從圖3 中給出的地震動(dòng)時(shí)程和傅里葉頻譜特性可以看出，提取的母波與實(shí)際地震動(dòng)記錄特征非常接近。

2.2 地震動(dòng)母波的合理性驗(yàn)證

為了驗(yàn)證基于PCA 算法提取的地震動(dòng)母波合成地震動(dòng)的合理性，本文以原始地震動(dòng)數(shù)據(jù)庫(kù)中的50 條近場(chǎng)(震中距R<30 km)數(shù)據(jù)為例提取母波，并進(jìn)行合成和驗(yàn)證。選取近場(chǎng)數(shù)據(jù)進(jìn)行母波提取和驗(yàn)證的原因是由于近場(chǎng)地震動(dòng)的特征更顯著、更加復(fù)雜，更具有代表性意義。

圖3 從原始地震動(dòng)數(shù)據(jù)庫(kù)中提取的4 條地震動(dòng)母波及其傅里葉頻譜Fig.3 Four ground motion mother waves and their Fourier spectra extracted from the original ground motion database

為了使得提取的母波能夠表征原始地震動(dòng)數(shù)據(jù)庫(kù)的特征，首先需要引入主成分累積貢獻(xiàn)率的概念，主成分累積貢獻(xiàn)率是選擇有效主成分的重要依據(jù)，它是主成分的方差在所考察的隨機(jī)變量的總方差中所占的比例；再引入累積方差解釋率概念，即多個(gè)主成分方差所占的比例之和，它是通過(guò)主成分貢獻(xiàn)率之和求得。當(dāng)累積方差解釋率比較高時(shí)，能夠較好的代表數(shù)據(jù)庫(kù)的特征。

為了確保能夠充分的提取地震動(dòng)數(shù)據(jù)庫(kù)的特征，本文累積方差解釋率取值為95%，即當(dāng)累積方差解釋率達(dá)到95%時(shí)，提取的含有本地地震動(dòng)特征信息的母波能夠很好表征原始地震動(dòng)數(shù)據(jù)庫(kù)的特征。圖4 給出了累積方差解釋率和母波地震動(dòng)數(shù)量的關(guān)系，圖中的拐點(diǎn)就是累積方差解釋率取值為95%的點(diǎn)。根據(jù)圖4 分析結(jié)果，當(dāng)滿足累積方差解釋率為95%時(shí)，提取的國(guó)內(nèi)近場(chǎng)地震動(dòng)數(shù)據(jù)的母波數(shù)目為19 條。

圖4 累積方差解釋率和母波地震動(dòng)數(shù)量的關(guān)系Fig.4 The relationship between the interpretation rate of cumulative variance and the number of ground motions of the mother wave

為了驗(yàn)證提取的母波的合理性，以提取的19 條地震動(dòng)母波來(lái)合成近場(chǎng)數(shù)據(jù)庫(kù)中的地震動(dòng)。首先任意選取近場(chǎng)地震動(dòng)數(shù)據(jù)庫(kù)中的一條記錄，計(jì)算該條地震動(dòng)記錄的反應(yīng)譜，對(duì)提取的19 條地震動(dòng)母波進(jìn)行線性組合，可使得組合的新的地震動(dòng)反應(yīng)譜與之前選取的地震動(dòng)反應(yīng)譜誤差最小，即可得到一條新的合成的地震動(dòng)。圖5 給出了實(shí)際地震動(dòng)與合成地震動(dòng)的反應(yīng)譜和地震動(dòng)時(shí)程的比較。從圖5 的反應(yīng)譜和地震動(dòng)時(shí)程的比較中可以發(fā)現(xiàn)，PCA 算法提取的地震動(dòng)母波能夠很好的合成原始地震動(dòng)數(shù)據(jù)庫(kù)中的任意一條地震動(dòng)記錄，因此，PCA 提取的地震母波能夠很好的表征原始地震動(dòng)數(shù)據(jù)庫(kù)的特征。

圖5 實(shí)際近場(chǎng)地震動(dòng)和應(yīng)用PCA 算法提取的母波合成的地震動(dòng)比較Fig.5 Comparison between the actual near-field ground motion and the synthetic ground motion by the mother waves extracted by PCA

3 地震動(dòng)母波權(quán)重系數(shù)的確定

需要求解出方程的最優(yōu)解，因此，引入粒子群算法。

3.1 粒子群算法(PSO)

粒子群算法(Particle Swarm Optimization, PSO)是由Kennedy 等[19]和Stefan 等[20]可以用于求解最優(yōu)化問題，能夠有效地實(shí)現(xiàn)計(jì)算機(jī)智能搜索和優(yōu)化。該方法所求出的解是全局最優(yōu)解而不是局部最優(yōu)，它能夠找出滿足條件的一組 ki使得 S最小，具體的要點(diǎn)如下，流程見圖6 所示。

1)參數(shù)的初始化。設(shè)置初始化參數(shù)，如：自變量 ki初始值，最大迭代次數(shù)，粒子的最大速度，粒子群的規(guī)模以及整個(gè)搜索空間。

2)個(gè)體極值以及全局最優(yōu)解。個(gè)體極值為每個(gè)粒子找到的最優(yōu)解，從這些最優(yōu)解找到一個(gè)全局值，叫做本次全局最優(yōu)解。與歷史全局最優(yōu)比較，進(jìn)行更新。

3)更新速度和位置公式，即式(15)。

式中： ω為慣性因子，當(dāng)取值較大時(shí)尋優(yōu)能力強(qiáng)；C1和 C2為加速度常數(shù)； Pid為個(gè)體極值；Pgd為群體極值； Xid為粒子當(dāng)前的位置； Vid粒子的速度；Maxgen 是迭代的次數(shù)。

4)設(shè)置迭代次數(shù)或者最小誤差。

圖6 粒子群算法求解權(quán)值ki 流程圖Fig.6 Flow chart of PSO algorithm to solve weight ki

3.2 權(quán)重系數(shù)求解

為了使得地震動(dòng)母波線性組合得到的新的地震動(dòng)的反應(yīng)譜與目標(biāo)反應(yīng)譜誤差最小，圖7 給出了基于PSO 算法[21 ? 22]求解權(quán)重系數(shù)，以及基于地震動(dòng)預(yù)測(cè)模型合成目標(biāo)地震動(dòng)的流程圖。首先，選取本地震動(dòng)數(shù)據(jù)庫(kù)區(qū)域合適的地震動(dòng)預(yù)測(cè)模型，應(yīng)用PCA 算法提取地震動(dòng)母波，通過(guò)地震動(dòng)預(yù)測(cè)模型[23]得到的反應(yīng)譜與組合地震動(dòng)母波得到的新的地震動(dòng)的反應(yīng)譜匹配，再用PSO 算法快速求解權(quán)重系數(shù)。PSO 算法的具體參數(shù)參考了文獻(xiàn)[24]，如表2 所示。

4 目標(biāo)地震動(dòng)的合成

圖7 應(yīng)用PCA 和PSO 算法合成地震動(dòng)時(shí)程的流程Fig.7 Flow chart of simulation ground motion time history by PCA and PSO

為了驗(yàn)證本文提出的方法的可行性，分別對(duì)中國(guó)西部地區(qū)的四個(gè)設(shè)定地震場(chǎng)景下的不同場(chǎng)點(diǎn)進(jìn)行地震動(dòng)合成。設(shè)定震級(jí)、斷層距以及場(chǎng)地條件如表3 所示。將設(shè)定震級(jí)、距離以及場(chǎng)地參數(shù)輸入到本區(qū)域的地震動(dòng)預(yù)測(cè)模型中，本文采用了文獻(xiàn)[25]基于四川地區(qū)的中強(qiáng)震數(shù)據(jù)建立的地震動(dòng)預(yù)測(cè)模型，與本文的研究區(qū)域一致。然后基于此模型對(duì)設(shè)定場(chǎng)點(diǎn)的地震動(dòng)反應(yīng)譜進(jìn)行估計(jì)，通過(guò)組合母波得到的新的地震動(dòng)時(shí)程并計(jì)算其反應(yīng)譜，當(dāng)計(jì)算的反應(yīng)譜與地震動(dòng)預(yù)測(cè)模型反應(yīng)譜一致時(shí)，則得到最終的地震動(dòng)時(shí)程，這是一個(gè)迭代過(guò)程。

通過(guò)粒子群算法尋優(yōu)計(jì)算出的權(quán)重系數(shù)值如表4 所示。

圖8 中給出了迭代次數(shù)和誤差S 之間的關(guān)系，從圖8 中可以看出在迭代到50 次時(shí)誤差都已收斂，因此，針對(duì)本次地震動(dòng)的計(jì)算模擬，可以取迭代次數(shù)為50。

圖9 給出了合成的地震動(dòng)時(shí)程的反應(yīng)譜與預(yù)測(cè)模型得到的目標(biāo)反應(yīng)譜的比較，圖中給出的分別是不同場(chǎng)點(diǎn)(R=10，30 km)和不同震級(jí)(M=5.5，6)的比較。從圖9 中能夠看出，通過(guò)PSO 智能算法求解出的地震動(dòng)的反應(yīng)譜能夠較好的匹配地震動(dòng)預(yù)測(cè)方程得到的目標(biāo)反應(yīng)譜。圖10 給出了最終合成的地震動(dòng)時(shí)程，從圖中可以看出合成的地震動(dòng)與實(shí)際地震動(dòng)非常接近，具有隨機(jī)性和非平穩(wěn)性，包含了區(qū)域地震動(dòng)的特征。因此，地震動(dòng)母波的線性組合能夠得到地震動(dòng)數(shù)據(jù)庫(kù)中的任意地震動(dòng)數(shù)據(jù)，合成的目標(biāo)地震動(dòng)既匹配了目標(biāo)譜，有能夠很好地代表本區(qū)域?qū)嶋H地震動(dòng)的特征。

表2 PSO 算法參數(shù)Table 2 Parameter of PSO

表3 設(shè)定地震信息和計(jì)算信息Table 3 Scenario earthquake and calculation information

表4 地震動(dòng)母波的權(quán)重系數(shù)Table 4 Weight coefficient of the mother wave of the ground motion

圖8 迭代次數(shù)與誤差之間的關(guān)系Fig.8 The relation between the number of iterations and the error

圖9 合成的地震動(dòng)的反應(yīng)譜與目標(biāo)譜的比較Fig.9 Comparison between the response spectra of the synthesized ground motion and the object spectra

圖10 機(jī)器學(xué)習(xí)方法合成的地震動(dòng)時(shí)程Fig.10 Time history of ground motion synthesized by machine learning method

5 結(jié)論

為了研究考慮區(qū)域地震動(dòng)特征信息的地震動(dòng)合成方法，本文引入了機(jī)器學(xué)習(xí)中PCA 算法，從地震動(dòng)數(shù)據(jù)庫(kù)中提取有效的地震動(dòng)母波信息，結(jié)合目標(biāo)區(qū)域的地震動(dòng)預(yù)測(cè)模型給出的特定場(chǎng)點(diǎn)的地震動(dòng)反應(yīng)譜，通過(guò)PSO 算法求解組合地震動(dòng)母波的權(quán)重系數(shù)，使得合成反應(yīng)譜與目標(biāo)譜誤差的最小，最終由地震動(dòng)母波線性疊加得到目標(biāo)地震動(dòng)時(shí)程。通過(guò)上述研究可以得到如下結(jié)論：

(1)應(yīng)用PCA 算法能夠從地震動(dòng)數(shù)據(jù)庫(kù)中能夠提取出代表性地震動(dòng)母波，地震動(dòng)母波能夠合理表征地震動(dòng)數(shù)據(jù)庫(kù)的特性。

(2)應(yīng)用PSO 算法能夠快速高效求解地震動(dòng)母波的組合權(quán)重，PSO 智能算法避免了應(yīng)用窮舉法求解權(quán)值，提升了計(jì)算速度。

(3)通過(guò)PCA 和PSO 智能算法，結(jié)合本區(qū)域?qū)嶋H地震動(dòng)和預(yù)測(cè)模型來(lái)合成新的地震動(dòng)時(shí)程，能夠合理的包含區(qū)域?qū)嶋H地震動(dòng)的特性，能夠匹配目標(biāo)地震動(dòng)的頻譜特征。

本文提出的方法考慮了區(qū)域?qū)嶋H地震動(dòng)的特征，使得合成的地震動(dòng)時(shí)程既包含了時(shí)程上的區(qū)域特征，又匹配了目標(biāo)譜，滿足了譜型上的一致性。采用PCA 和PSO 智能算法，提高了計(jì)算效率，滿足了地震動(dòng)合成時(shí)效性的需求，因此，可為未來(lái)面向工程的抗震性能評(píng)估提供合理的地震動(dòng)時(shí)空分布場(chǎng)。