基于PCA 及PSO 智能算法的地震動合成方法—以中國西部中強地震為例

胡進軍，張 輝，靳超越，王中偉，胡 磊

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

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

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

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

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

1 區(qū)域強地震動數(shù)據(jù)

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

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

圖1 臺站和地震震中分布圖Fig.1 Map of stations and earthquake epicenters

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

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

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

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

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

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

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

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

首先，對數(shù)據(jù)進行標準化處理：

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

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

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

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

2.1 基于PCA 算法的地震動母波提取

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

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

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

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

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

2.2 地震動母波的合理性驗證

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

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

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

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

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

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

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

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

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

3.1 粒子群算法(PSO)

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

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

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

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

式中： ω為慣性因子，當取值較大時尋優(yōu)能力強；C1和 C2為加速度常數(shù)； Pid為個體極值；Pgd為群體極值； Xid為粒子當前的位置； 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ù)求解

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

4 目標地震動的合成

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

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

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

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

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

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

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

表4 地震動母波的權(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 合成的地震動的反應(yīng)譜與目標譜的比較Fig.9 Comparison between the response spectra of the synthesized ground motion and the object spectra

圖10 機器學習方法合成的地震動時程Fig.10 Time history of ground motion synthesized by machine learning method

5 結(jié)論

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

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

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

(3)通過PCA 和PSO 智能算法，結(jié)合本區(qū)域?qū)嶋H地震動和預(yù)測模型來合成新的地震動時程，能夠合理的包含區(qū)域?qū)嶋H地震動的特性，能夠匹配目標地震動的頻譜特征。

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