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

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

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

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

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

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

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

1 區域強地震動數據

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

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

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

2 區域地震動信息的提取

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

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

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

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

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

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

若有一組如下形式的數據，應用主成分分析的具體步驟為：

首先，對數據進行標準化處理：

其次，計算相關系數矩陣：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 地震動母波權重系數的確定

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

3.1 粒子群算法(PSO)

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

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

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

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

式中： ω為慣性因子，當取值較大時尋優能力強；C1和 C2為加速度常數； Pid為個體極值；Pgd為群體極值； Xid為粒子當前的位置； Vid粒子的速度；Maxgen 是迭代的次數。

4)設置迭代次數或者最小誤差。

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

3.2 權重系數求解

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

4 目標地震動的合成

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

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

通過粒子群算法尋優計算出的權重系數值如表4 所示。

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

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

表2 PSO 算法參數Table 2 Parameter of PSO

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

表4 地震動母波的權重系數Table 4 Weight coefficient of the mother wave of the ground motion

圖8 迭代次數與誤差之間的關系Fig.8 The relation between the number of iterations and the error

圖9 合成的地震動的反應譜與目標譜的比較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 結論

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

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

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

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

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