匹配多目標參數的地震動合成方法

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

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

在工程地震研究中如何得到場地相關的地震動一直受到研究者的廣泛關注。人工合成地震動可為地震動記錄匱乏地區的地震危險性分析和結構抗震分析、性能評估提供輸入。1947 年Housner[1]從工程研究的角度將地震動視作大小一定、隨機到達的脈沖的疊加來模擬地震動開始，國內外許多學者對地震動合成方法進行了研究，也提出了多種滿足部分地震動統計特征，如功率譜、傅氏譜、反應譜等的人造地震動方法。目前地震動的合成方法主要分為三種：一是以直接運用地震動三要素與震源、路徑和場地的關系為基礎的經驗方法，如比例法、經驗格林函數法[2-3]；二是采用隨機振動模型的地震動合成方法，如脈沖法、隨機點源法、隨機有限斷層法[4-6]；三是采用考慮強度和頻率非平穩的地震動合成方法，如各種時頻非平穩合成、相位譜、相位差譜、三角級數法等[7-12]。

在工程中應用最廣泛的則是擬合目標反應譜的地震動時程合成方法，以合成時程的反應譜與給定的目標反應譜相等或相近為控制條件。但是，僅以匹配目標譜為控制條件會一定程度上得到不唯一的合成地震動時程，而且這些合成地震動由于缺乏更多的約束條件，基于此進行時程分析可能會增加結構響應的離散性和不確定性。究其原因，反應譜雖然可以反映地震動幅值和頻譜的特性，但是對地震動三要素中的持時并沒有體現。作為地震動三要素之一的持時對于結構的響應有重要的影響[13-14]，Bommer 等[15]認為地震發生過程中的強震持續時間對于基礎材料和結構響應有著重要的影響，特別是在強度和剛度退化時。因此，地震危險性分析也應該包括對強震持時的估計。Mahin[16]也指出地震動持續時間對結構的彈性變形和耗能需求有顯著影響，尤其是對于相對薄弱的短周期結構。Cabanias 等[17]指出使用同時考慮振幅和持時的響應參數可以更好的估計結構的潛在損傷。

地震動的幅值、頻譜及持時等特征主要受到震源、路徑和場地條件影響和控制。研究表明不同地震構造區域的地震動的幅值、頻譜和持時存在區域性差異[18]，因此，在合成和模擬場地相關地震動時，應該通過約束幅值、頻譜和持時這些工程相關的主要參數，以期使得合成的地震動中包含地震動的區域特征。基于此，本研究的目的是建立同時匹配設定區域地震動目標反應譜和持時參數的地震動合成方法。一方面，這樣合成的地震動可以綜合的考慮到地震動三要素中的幅值、頻譜和持時特征；另一方面，增加地震動持時的合成控制條件可以在實現對目標反應譜匹配的基礎上，為研究輸入地震動持時對結構響應的影響提供基礎。

近年來機器學習方法被多個領域廣泛應用[19-21]。地震學是一個數據密集型領域，目前正在經歷著大數據界中體量、種類、速度的快速變化[22-24]，這些變化為基于數據驅動的機器學習在地震學中的應用提供了條件。Khosravikia 等[25]開發了一種人工神經網絡用于美國俄克拉荷馬州等地區的自然和誘發地震的地震動預測模型；Kotha 等[26]在對日本活動性淺源地震的地震動預測模型回歸分析的基礎上，利用無監督學習中的聚類技術對代表各個場地的經驗放大函數的估計結果進行了分類，進而提出了一種修正的數據驅動的場地分類方法；Alimoradi 等[27]開發了一種將一組稀疏的被稱作特征地震記錄的標準正交基向量進行線性組合合成地震動的方法，其組合系數是匹配由高斯過程回歸得到設定情景下的反應譜確定。這些研究表明，機器學習具有很好的數據分析處理能力，可以為解決復雜的地震動模擬和合成問題提供新的思路。

為了在地震動合成中考慮地震動幅值、頻譜和持時的特征[28]，將地震動數據與機器學習相結合，運用智能算法尋優求解實現。本文首先應用主成分分析[27,29-31]從區域地震動數據庫中提取具有本區域地震動特征的地震動母波，為了實現合成的地震動能夠同時考慮地震動幅值、頻譜和持時特征，選擇同時匹配目標反應譜和目標持時。采用多目標遺傳算法求解地震動母波的線性組合系數并且給出合成地震動的評價標準來控制合成效果，從最優解集中篩選出最合理的母波組合系數用于合成地震動。以2019 年中國四川長寧Ms6.0地震的主余震數據作為區域地震動數據庫，對本文方法進行了可行性和有效性的驗證，最后將本文方法同三角級數法進行了對比。

1 基于PCA 和MOP 的地震動合成方法

1.1 主成分分析方法

為了在合成地震動時考慮地震動區域特征，需要對大量原始地震動數據的特征進行提取，這實質上就是數據降維。數據降維主要有投影和流行學習兩種常用的方法。主成分分析(Principal Component Analysis, PCA)是基于投影的方式來實現線性數據降維和去除相關性的一種方法，其在數量地理學、數學建模和數理分析等學科領域中都有廣泛應用。本文將PCA 算法引入對地震動數據庫的降維處理中，提取表征數據庫主要特征成分的特征地震動母波，將這些特征地震動母波用于合成地震動時程。PCA 算法的核心是先識別出最接近數據的超平面，然后，將數據投影到該超平面之上，進而實現降低數據向量的維度并去除掉各個分量之間的相關性。目標超平面的識別可以轉化為如何得到投影矩陣。該問題是通過優化目標函數算法來實現的，這與其他機器學習算法的思路是一致的。

以一個包含N條地震動記錄的數據集為例來說明PCA 方法的應用原理。該地震動數據集可視為一個擁有N個維度的數據空間，則可以通過N個兩兩互相垂直坐標軸上的單位向量的線性組合來表示數據空間中的任意一條地震動記錄。將整個數據集投影到第i個的單位向量所在坐標軸上后就稱為該投影結果是第i個主成分，即第i條特征地震動母波，且這些用于投影的坐標軸排序是按照數據空間向本軸單獨投影后保留數據差異性的大小確定。具體對地震動數據集應用PCA 算法的降維流程如下：

1)計算N維實際地震動記錄數據集的均值向量，將所有向量減去各自的均值向量，即白化處理；

2)計算數據集的散度矩陣S(或協方差矩陣C)；

3)對散度矩陣進行特征值分解，利用奇異值分解(Singular Value Decomposition, SVD)，得到所有的特征值和特征向量；

4)將所得特征值按從大到小排序，保留最大一部分的K個特征值對應的特征向量，以它們為列，形成降維所需的投影矩陣；

5)將白化后的數據庫矩陣右乘投影矩陣，得到降維后全新的正交特征地震動數據，即為特征地震動母波。

對于選擇降維后的維度K(即特征地震動母波個數)，可以用SVD 分解時產生的散度矩陣S來表示。利用式(1)計算可得滿足誤差條件最小的K值，若其中t值取0.15，則代表了該PCA 算法保留了原數據85%的主要信息，從而可以確定該允許誤差下的提取的最少特征地震動母波數。具體t的取值可根據對誤差的要求確定，由此選擇不同的K值。

1.2 合成系數的控制條件

在得到特征地震動母波的基礎上，需要考慮合成地震動時各條母波的對應系數的取值問題。不同于傳統的合成地震動時程只考慮目標反應譜的匹配，本文從多參數的角度約束合成地震動時程。工程輸入中地震動的特征參數一般包括幅值、頻譜和持時。地震動的反應譜可以體現地震動的幅值和頻譜兩方面的特征，而地震動持時是結構地震反應分析中的重要參數。在結構地震反應的彈性階段，循環往復的地震作用可能引起構件的低周疲勞，而在結構進入非彈性階段后，持時的長短可能對結構反應、殘余變形和倒塌存在重要影響[32]。因此，需要引入地震動持時作為約束合成的地震動時程的控制條件。

1.2.1 合成地震動的評價函數

為了對持時進行約束，首先需要選擇適合目標地區的地震動持時預測方程，進而得到目標地震動持時，記為d5-95，新合成地震動的持時記為D5-95。持時項誤差絕對值為：

1.2.2 合成地震動的評價標準

為了對上述合成地震的兩項誤差給出一個合理的評價標準。對于地震動預測模型(GMPE)，均采用式(4)和式(5)中的評價標準來篩選出合成的地震動時程。

在評價標準中，以地震動反應譜和持時預測方程中給出的標準差σ 為評價條件，合成反應譜在各控制周期點與目標反應譜譜值之差的絕對值小于一倍的反應譜衰減關系中的標準差時，視為合成地震動反應譜與目標反應譜良好匹配；對于地震動持時也以一倍的標準差為控制值來判斷合成地震動持時的匹配效果。

1.3 確定合成系數

為了將PCA 算法提取出的特征地震動母波線性組合得到合成地震動，需要確定組合系數。由于目標評價函數增加，傳統上用于求解評價函數最優化解的單目標優化算法(SOP)已不再適用，而選擇多目標優化算法(MOP)中最常用的多目標遺傳算法來求解合成系數。其核心是協調各個目標函數之間的關系，從而得到使各子目標函數都盡可能達到最優的最優解集，這不同于單目標優化的解通常唯一，多目標優化問題的解是一組均衡解集。在應用算法得到一系列合成系數解集后，通過設定的評價函數指標篩選出最滿意的一組合成系數。

多目標遺傳算法中的帶精英策略的非支配排序遺傳算法(NSGA-II)是2000 年Deb 等[33]為改進NSGA 而提出的。由于改進后的二代算法引入了快速非支配排序法、精英策略、采用擁擠度和擁擠度比較算子使得算法的計算復雜度大大降低，計算結果在Pareto 域中均勻分布，維持了種群的多樣性。具體算法流程如下，并以流程框圖1說明該優化算法對于合成控制函數的優化過程：

圖1 多目標遺傳算法流程圖Fig. 1 Flow chart of multi-objective genetic algorithm

1)隨機產生設定規模N的初始種群，經非支配排序后通過遺傳算法的選擇、交叉、變異后得到第一代子代種群；

2)從第二代開始，將父代與子代種群合并，進行快速非支配排序，同時計算每個非支配層中的個體的擁擠度，根據非支配關系以及個體擁擠度選取合適的個體組成新的父代種群；

3)最后通過遺傳算法的基本操作產生新子代種群：循環上述3 個步驟，直到達到終止代數，輸出最優解集。

1.4 地震動合成整體流程

本文提出的合成方法是應用PCA 技術從區域地震動數據庫中提取地震動母波，在合成中考慮了對目標反應譜和顯著持時進行同時匹配，使用多目標遺傳算法確定出母波的線性組合系數得到最終的合成地震動。合成的整體流程見圖2。

圖2 地震動合成流程圖Fig. 2 Flow chart of ground motion synthesis method

2 方法驗證

2019 年6 月17 日在四川省東南部的宜賓市長寧縣發生了Ms6.0 級地震。在震后的幾日內，震源附近接連發生了多次的強余震，我國強震臺網觀測到了此次群震的大量地震動[34]。對于長寧地震這類數據特征明顯且記錄豐富的數據，應用PCA技術提取的特征母波可以在滿足設定的誤差條件下，代表整個主余震地震序列中的主要特征信息。通過對長寧地震動的合成和比較可以驗證本方法的可行性和有效性。

2.1 地震動數據庫

本文收集了此次長寧地震的主震及其4 級以上余震的記錄作為地震動數據庫，具體的地震信息見表1。數據庫包含9 次地震事件中，來自69 個臺站的286 條水平向地震記錄，并采用Boore 等[35]提出的方法對地震動數據進行了濾波和基線調整。

表1 長寧主余震信息Table 1 Changning mainshock-aftershock events

2.2 地震動母波的提取

首先按式(1)計算確定提取的地震動母波數目，為了使地震動母波能夠保留數據集中95%的主要信息，t取值為0.05，通過計算得到滿足誤差條件的最小K值為11。在數據庫中提取按照保留信息占比大小進行排名后的前11 條地震動母波。表2 中給出了提取的地震動母波對原始數據庫信息保留占比情況，圖3 給出了排名前八位的地震動母波時程和時頻特征，從圖3 中可以看出提取到的地震動母波與實際地震動記錄的時域和頻域非平穩特性相似，這一特征使地震動母波適合作為一組基向量用于地震動合成。

圖3 排名前8 的地震動母波Fig. 3 Top 8 mother wave ground motions

表2 提取的地震動母波保留信息占比Table 2 Proportion of information retained by extracted ground motion mother waves

2.3 地震動預測模型

盡管從實際地震動數據中提取而來的地震動母波包含了區域地震動的特征信息，但是直接采用提取的母波合成地震動并不能考慮震級、距離和場地等因素的約束，為了考慮這些區域地震相關因素對地震動幅值、頻譜和持時的影響，需要在合成地震動時程對地震動的工程相關參數加以約束，使得合成地震動能夠與區域實際地震動特征很好的匹配。

四川長寧地震發生于我國西部，因此需要選取適合于此地區的中國西部地區地震動模型[36]，以及基于我國強震動數據建立的地震動持時預測方程[37]。通過設定地震情景，利用反應譜和持時地震動模型分別計算得到目標反應譜和目標持時作為對長寧地震進行地震動合成的約束條件。選用的反應譜地震動模型和地震動持時模型的具體形式如下：

圖4(a)～圖4(d)中分別給出了長寧地震記錄在0.01 s、0.20 s、1.00 s 三個周期點處的反應譜值和對應地震動模型得到的衰減曲線的對比，以及實際地震動數據的5%～95%顯著持時同持時預測方程所得預測值的對比。從圖4 中可以看出，選取的兩個地震動模型可以較好預測本區域的實際地震動特征，也反映出了作為合成約束條件的目標反應譜和持時選取的合理性。

圖4 長寧地震地震動記錄反應譜和持時與地震動預測模型對比Fig. 4 Comparison between response spectra and durations of Changning earthquake and ground motion prediction models

2.4 合成目標地震動時程

為了驗證合成方法的可行性，采用本文方法合成一條地震動并將其與實際地震動數據庫中選取的一條地震動記錄進行匹配，檢驗匹配結果是否合理。任意選取長寧地震中的一個臺站51NXT，選取的地震動信息如下：震級5.4 級，震源深度10 km，震中距46.92 km。然后，基于本方法進行合成，表3 中給出了合成的地震動加速度幅值、反應譜、持時與實際地震動記錄的匹配效果。圖5中為合成的地震動與實際地震記錄的匹配效果和誤差分布圖。圖5(b)中橫縱軸分別表示合成地震動的反應譜和持時與地震動模型所得目標值之間的誤差值，具體計算采用式(2)和式(3)，圖中分布的每個散點均代表最優解集中的一組母波線性組合系數，正方形散點代表著在評判標準下最終確定的組合系數。圖5(a)、圖5(c)中的合成地震動即為利用該散點所對應組合系數將母波線性組合得到的結果。結合表3 和圖5 可以看出合成地震動的峰值加速度和反應譜都實現了同實際記錄的良好匹配；更重要的是，合成地震動的90%顯著持時與實際記錄完美匹配，將合成地震的持時控制到了區域持時預測模型的目標水平。實現了考慮區域地震動預測模型約束的地震動合成，可以為工程結構動力時程分析提供體現區域特征的地震動輸入。

圖5 合成地震動與實際地震動的比較Fig. 5 Comparison between synthetic ground motion and actual ground motions

表3 合成地震動幅值、頻譜和持時的匹配情況Table 3 Matching of amplitude, spectrum and duration of synthetic ground motion

3 與其他方法的比較

三角級數法是工程上廣泛應用的合成地震動的方法，為了將本文方法與其結果進行對比，采用三角級數法進行地震動合成。三角級數法合成中采用的地震動時程強度包絡函數如式(8)所示。地震動包絡函數參數取值根據霍俊榮[38]建立的回歸式(9)確定。

式中：c為衰減關系；t1和t2分別是控制平穩段首末時刻；t3是地震動持時。

式中：Ms為地震震級；D為震中距；ts為地震平穩段持續時間，t2=t1+ts。

合成地震動所需的設定地震信息以及持時和反應譜匹配情況見表4，圖6 和圖7 分別給出了兩種方法合成的地震動時程、持時曲線、多目標優化最優解集的誤差分布、反應譜的譜形匹配誤差。在同時滿足對目標反應譜和持時的允許誤差的解集中，篩選出與反應譜匹配最佳的一組合成系數，并以此系數合成本文方法的地震動。

圖6 Ms=5, R=10 km 時兩種方法合成結果對比Fig. 6 Comparison of synthetic results of two methods for Ms=5 and R=10 km

圖7 Ms=6, R=10 km 時兩種方法合成結果對比Fig. 7 Comparison of synthetic results of two methods for Ms=6 and R=10 km

表4 設定地震下兩種方法持時和反應譜匹配情況Table 4 Matching of duration and response spectrum of two synthesis methods under twoearthquake scenarios

對比兩種方法合成的地震動，從圖6(c)、圖7(c)兩次設定地震下本文方法的最優解集誤差分布上來看，均存在較多組的系數組合可以使得合成地震動滿足預先設定的評價標準。在圖6(c)、圖7(c)中也可以看出單一追求同反應譜的匹配會導致合成地震動的持時與地震動持時預測方程得到的目標持時偏差過大的情況出現，不能在合成地震動中很好地約束工程抗震中所關注的地震動持時，這對于將合成的地震動作為輸入研究結構響應會造成影響。這一結果在圖6(b)、圖6(d)、圖7(b)、圖7(d)本文方法和三角級數法的合成結果對比中也可以看出。三角級數法考慮通過與反應譜的匹配可以實現對地震動幅值和頻譜的約束，但不能對地震動持時參數進行約束。而應用本文的地震動合成方法在約束條件和評價標準的控制下，可以很好的將區域地震動預測方程所體現的區域地震動特征在合成結果中體現。

4 結論與討論

本文提出一種基于機器學習合成考慮區域地震動幅值、頻譜和持時特征的地震動的方法，基于長寧地震的主余震數據對方法進行了驗證，并與三角級數方法進行了對比。本方法采用了從區域實際地震動中提取的特征地震動母波，保留了原始地震動的時頻特征，合成的地震動考慮了區域地震動的持時特征。總結本文可得如下結論：

(1)將機器學習中的主成分分析應用于從大量實際地震動數據中提取地震動特征成分是可行的、有效的。以降維提取得到的少量特征地震動母波來代表整個數據庫的整體特征，反映地震動的震源、路徑和場地的綜合效應。從母波時程和時頻分析中可以看出，特征地震動母波與實際地震動的時頻非平穩特征一致，可以將其用于合成包含本區域數據特征的目標地震動。

(2)不同于一般方法僅將目標反應譜作為匹配條件，本文方法可以實現對地震動幅值、頻譜和持時進行同時約束，考慮了地震動的時頻特征。對地震動三要素的綜合考慮也能夠更好得到符合實際地震動特征的合成結果。

(3)與三角級數法對比表明，僅匹配反應譜可能會得到與實際地震動持時相差很大的合成地震動，這會對結構動力響應結果引入較大的不確定性。