基于實測強震記錄的主余震向量過程降維模擬

劉子心，姜云木，劉章軍，阮鑫鑫

（1.防災科技學院中國地震局建筑物破壞機理與防御重點實驗室，河北三河 065201；2.武漢工程大學土木工程與建筑學院，湖北武漢 430074）

引言

研究表明，地震動不僅具有較強的隨機性，同時具有明顯的序列性，即在一次主震后往往會伴隨著多次余震的發生，例如2008年汶川地震，在主震后共發生了13 718次余震，且余震最高震級達6.5級。由此可見，余震不僅數量較多，且部分震級較大，成為結構二次破壞的主要因素。然而，由于主余震型地震動實測記錄數量有限，且需要對實測強震記錄進行一定處理才能滿足結構抗震分析的需求，因此，采用人工合成方法生成主余震型地震動時程近年來受到了廣泛關注。

為了構造主余震型地震動，近年來學者們開展了對主余震參數關系的研究。歐進萍等［1］根據49 組主余震震級資料，線性回歸了主余震間震級的經驗公式；周昱辰等［2］利用蘆山的實測主余震強震記錄，研究了4 種地震預警震級估計參數隨震級與時間的變化趨勢；Zhang 等［3］以線性回歸的方式，較為全面地研究了主余震間在持時、幅值和頻譜方面的關系。上述研究為人工合成主余震型地震動奠定了良好基礎。一般來說，主余震型地震動的人工合成方法大致可分為確定性方法與隨機方法兩大類。對于確定性方法，Hatzigeorgiou等［4］提出了采用重復法構造主余震序列，該方法假設主余震的特性一致，但這明顯不符合實際狀況；Li 等［5］提出了隨機組合法構造主余震序列，即從主震和余震記錄庫中分別隨機挑選出一條地震動記錄組合成主余震序列，然而，該方法得到的主余震序列在本質上不具有隨機性。對于隨機方法，Hu等［6］利用調幅過濾白噪聲方法生成主震過程，再根據分支余震序列法生成余震。然而，該方法屬于Monte Carlo 方法，無法在概率層面上精確描述主余震過程；姜云木等［7］基于Copula 函數，建立了主余震參數的相關結構，并引入正交隨機變量的隨機函數方法，從而實現了給定主震參數條件下余震的高效模擬。然而，該研究采用了時-頻全非平穩演變功率譜模型，由于該模型在時頻域相互耦合，因此在識別主余震參數時不可避免地引入了誤差。同時，該研究將主余震型地震動視為2 個相互獨立的隨機過程，并沒有考慮主余震間的相干性，需要進行2 次模擬才能得到主余震型代表性時程。綜上可見，亟需發展一種能夠在概率層面上一體化高效模擬主余震型地震動的方法。

為此，本研究將主余震型地震動視作一個隨機向量過程。首先，根據實測強震記錄，建立了主余震相干函數模型；其次，根據地震動的能量曲線與反應譜，對主余震的演變功率譜模型參數進行識別，并給出了對應不同場地的所有參數取值；最后，結合基于本征正交分解（POD）的隨機函數方法［8-9］，實現了主余震向量過程的降維模擬。數值算例驗證了主余震相干函數模型的正確性以及參數識別方法的有效性。值得說明的是，采用隨機函數-降維模擬方法生成的每一條主余震型代表性時程均具有賦得概率，且代表性時程集合構成一個完備的概率集，因此可與概率密度演化方法［10-11］相結合，實現復雜工程結構在主余震作用下的精細化動力響應與可靠性分析。

1 主余震實測強震記錄的選取

文中從美國太平洋地震工程研究中心（PEER）的NGA-West2數據庫和中國國家地震科學數據共享中心（CSMNC）中篩選了20 次地震的488 組主余震實測記錄，基本信息如表1 所示。主余震記錄的篩選原則如下［12］：

表1 選取的主余震實測記錄基本信息［7］Table 1 The basic information of the measured main-aftershock records selected in this paper

（1）主震記錄與其對應的余震記錄必須來自同一臺站；

（2）僅選取與主震對應且震級最大的余震作為研究對象；

（3）斷層距離應大于10 km，以減少近場效應的影響；

（4）主震與余震的震級均應大于4，以排除對結構影響較小的地震動。

PEER 數據庫提供了VS,30作為場地劃分的依據，同時，文獻［13］給出了《中國地震動參數區劃圖》（GB 18306-2015）［14］中5 種場地類別與VS,30的對應關系。然而，與PEER 數據庫不同，CSMNC 數據庫僅以Rock或Soil區分場地。根據《建筑抗震設計規范》（GB 50011-2010）［15］中對場地的描述，文獻［16］建議Rock對應Ⅱ類場地，Soil對應Ⅲ類場地。綜上，文中選用的主余震實測記錄的場地類別與VS,30及地震數量的對應關系如表2所示。

表2 對應于不同場地類別的VS,30 范圍及實測記錄數量［13］Table 2 The range of VS,30 and number of measured motion records corresponding to different site classifications

為保證實測強震記錄能夠滿足分析需求，還需對其進行基線校正和四階Butterworth 濾波處理，以及在1%～99%的能量范圍內進行截取。

2 主余震的相干性分析

眾所周知，局部場地對地震動的頻譜特性有較大影響，而同一臺站記錄的主震和余震發生于同一局部場地，由此推知，主余震在頻譜特性上可能存在著某種相干性。為此，根據篩選的488 組實測強震記錄，對主震與余震之間的相干性進行統計分析。為簡便起，將主震和余震均視為平穩過程，則（遲滯）相干函數定義為：

式中：Sma(ω)為主震和余震過程的互功率譜密度函數；Sm(ω)與Sa(ω)分別為主震和余震過程的自功率譜密度函數。

對于任意一組主余震實測強震記錄，采用MATLAB 工具箱自帶的“cpsd”函數計算主余震的互功率譜，以及“pwelch”函數分別計算主震和余震的自功率譜，利用式（1）即可得到每組實測主余震的相干函數。圖1給出了II類和III類場地類別中典型主余震實測記錄的相干函數曲線。

圖1 II類和III類場地類別中典型主余震實測記錄相干函數Fig.1 Typical coherence function of the measured main-aftershock records for soil site classifications II and III

從圖1可知，在工程常用的頻率范圍(0～300) rad/s內，主余震的相干性總體上隨著頻率的增加呈現遞減趨勢。為此，可采用前三階傅里葉級數模型對主余震的（遲滯）相干函數進行擬合，即：

式中，λγ=(A,B1,C1,B2,C2,B3,C3,D)為傅里葉級數模型的參數向量。

在計算主余震的相干函數時，將每一組實測主余震記錄視作一個隨機過程，分別計算488組實測主余震記錄的相干函數，進而得到它們的均值相干函數曲線，并根據最佳平方逼近原則，以γma(ω)為目標值，對傅里葉級數模型的參數向量λγ進行識別，即：

式中，截斷頻率ωu=300 rad/s。傅里葉級數模型參數向量λγ的識別結果如表3所示。

表3 傅里葉級數模型的參數取值Table 3 Parameter values of the Fourier series model

進一步，傅里葉級數模型與主余震實測記錄（遲滯）相干函數的擬合結果如圖2所示。

圖2 相干函數模型與實測記錄的擬合結果Fig.2 Fitting result between the coherence function model and the measured records

需要說明的是，文中并未對相干函數進行場地類別分組處理，這是由于每一組實測主余震記錄均來自同一臺站，因此場地差異性對主余震相干性的影響可以忽略。

3 主余震過程的模型參數識別

3.1 非平穩地震動過程的演變功率譜模型

為簡便起見，文中采用強度調制的演變功率譜密度函數形式［17］：

式中：SJ(ω,t;λJ)為非平穩地震動過程的單邊演變功率譜密度函數；q(t;λq)為強度調制函數；S(ω;λS)為平穩地震動過程的單邊功率譜密度函數。

對于強度調制函數，采用Amin-Ang提出的三段式模型，該模型反映了地震動過程的上升段、平穩段和衰減段，其表達式為［18］：

式中，λq=(t1,t2,α)為強度調制函數q(t;λq)的參數向量。對于平穩地震動過程的單邊功率譜密度函數，采用經典的Clough-Penzien模型［19］：

式中

式中：ωg和ξg分別為場地土的卓越圓頻率和阻尼比；ωf和ξf分別為基巖的卓越圓頻率和阻尼比。一般地，可取ωf=0.1ωg，ξf=ξg。Amax為地震動峰值加速度的均值，為峰值因子。λS=為平穩地震動過程功率譜S(ω;λS)的參數向量。

由此，非平穩地震動過程的演變功率譜SJ(ω,t;λJ)的參數向量為

3.2 演變功率譜模型的參數識別

對于第i條實測強震記錄ai(t)，其隨時間變化的歸一化能量曲線Ii(t)可表示為［20］：

式中，Ti為第i條實測強震記錄的持時。

對于非平穩地震動過程，其隨時間變化的歸一化模型能量曲線PJ(t;λJ)為：

根據式（9）與式（10），以Ii(t)為目標值，采用最佳平方逼近原則，便可對參數向量λq,i進行識別：

由此，即可得到與第i條實測強震記錄相對應的強度調制函數參數向量λq,i=(t1,i,t2,i,αi)。

為了對平穩地震動過程功率譜S(ω;λS)的參數向量進行識別，首先，根據帕薩瓦爾定理，即信號在時域與頻域上的總能量相等，可得：

同時，令

于是，由式（12）和式（13）得到：

對于式（13），結合式（7）推導可知：

由于第i條實測強震記錄ai(t)的峰值加速度Amax,i是已知的，這樣，將所得的參數向量λq,i代入式（14）和式（15）中，即可得到與實測強震記錄ai(t)所對應的峰值因子。

對于平穩地震動功率譜的場地參數，可采用擬合反應譜的方法來識別。為簡便起見，在對場地參數進行識別時，將地震動視為等效平穩過程。

Vanmarcke將隨機過程的反應譜定義為單質點體系反應峰值系數的平均值與反應均方差的乘積［21］。由此，對于第i條實測強震記錄ai(t)，其反應譜與功率譜的轉換公式為：

式中：r(ω0;λq,i)為等效平穩過程峰值系數的平均值；σ(ω0,ξ;λS,i)為等效平穩過程的反應均方差。ω0與ξ分別為結構的固有圓頻率和阻尼比，在文中ω0≥1.05 rad/s，ξ=0.05；Td，i為等效平穩過程的持續時間，即強度超過峰值50%的持續時間，對于三段式強度調制模型，Td，i的表達式為：

最后，以實測強震記錄ai(t)的前6 s反應譜Sa,i(ω0)為目標值，根據式（16），采用最佳平方逼近原則，對參數向量λS,i進行識別：

由此可以得到功率譜模型的場地參數ωg,i和ξg,i。這樣，非平穩地震動過程的演變功率譜SJ(ω,t;λJ,i)的參數向量λJ,i識別完成。

進一步，文中引入了R2（決定系數）作為參數識別效果的衡量標準。任意一組待擬合數據y的R2可定義為［22］：

式中：M為待擬合數據y的總數據點數；分別為待擬合數據的均值和擬合數據。一般地，若R2越趨近于1，則擬合效果越好。

最后，對于非平穩主余震型向量過程J(t)=(J1,J2)=(m,a)，元素J1和J2分別代表主震過程和余震過程，將主震參數和余震參數分別代入式（4）中，即可得到主震過程和余震過程的演變功率譜密度函數。

現以TCU104 臺站記錄的CHICHI 主余震為例，根據式（11）與式（20），反應譜和地震動歸一化能量的擬合結果分別如圖3 與圖4 所示。從圖3 和圖4 可以看出，無論是主震還是余震過程，模型與實測記錄均擬合良好。

圖3 反應譜擬合結果Fig.3 Fitting result of response spectrum

圖4 地震動歸-化能量擬合結果Fig.4 Fitting result of ground motion normalized energy

上述主余震典型實例的參數識別結果如表4所示。

表4 TCU104臺站記錄CHICHI主余震參數識別結果以及誤差Table 4 Parameter identification results and error of CHICHI main aftershocks recorded by TCU104 station

可以看出，主余震型地震動的反應譜和歸一化能量曲線的R2均趨近于1，說明了文中參數識別方法的有效性。

類似地，利用主余震實測強震記錄，即可對主余震的演變功率譜模型參數向量λJ進行識別。根據不同場地類別的主余震實測強震記錄，應用最佳平方逼近對主余震向量過程的演變功率譜模型進行參數識別，并求得每個參數的均值，將其作為主余震向量過程的演變功率譜模型參數的建議取值，具體如表5所示。

表5 不同場地類別主余震演變功率譜參數建議取值Table 5 Recommended values of main-aftershock evolutionary power spectrum parameters for different soil site classifications

同樣，對于峰值加速度和峰值因子，計算488 組主余震實測強震記錄的峰值加速度比值和峰值因子比值，并取均值。具體結果如表6所示。

需要說明的是，由于I0和IV類場地實測強震記錄數量過少，因此沒有對這兩類場地進行分析。通過表5 和表6可以看出，相較于余震，主震持時更長，能量更高，場地土卓越頻率更小；不同場地相比，隨著土質變軟，持時更長，而場地土卓越頻率更小。

根據表5和表6，便可實現對主余震演變功率譜的建模。

表6 主余震幅值參數比值Table 6 Main-aftershock amplitude parameter ratio

4 基于POD的主余震向量過程降維模擬

4.1 基于正交隨機變量的POD表達

對于主余震型地震動過程，可將其視作一個零均值的1D-2V（一維雙變量）隨機向量過程，因此可采用本征正交分解（POD）方法進行模擬。非平穩主余震向量過程J(t)=(J1,J2)=(m,a)的演變功率譜密度矩陣SJ(ω,t)為：

根據文獻［9］，主余震的演變功率譜矩陣SJ(ω,t)可分解為如下形式：

式中：上角標*和T分別表示復共軛和矩陣轉置；D(ω,t)為對角矩陣，即D(ω,t)=；γ(ω)為主余震的相干函數矩陣：

主余震的相干函數矩陣γ(ω)是一個非負定的Hermitian 矩陣，其元素可由式（2）得到。對相干函數矩陣γ(ω)進行特征分解，即可得到：

式中：特征向量矩陣Ψ(ω)=[ψ1(ω),ψ2(ω)]；特征值對角矩陣Λ(ω)=diag[Λ1(ω),Λ2(ω)]；I為2×2 階的單位矩陣。由于存在Λi(-ω)=Λi(ω)以及ψi(-ω)=，因此可令

式中，i為虛數單位。根據文獻［13］，主余震向量過程J(t)的第r個分量過程Jr(t)的POD公式可表達為：

式中：N為頻率截斷項數；Δω為頻率步長；ωk=kΔω(k=1,2,...,N)；Rik和Iik為零均值的正交隨機變量，滿足如下基本條件：

4.2 正交隨機變量的降維表達

需要說明的是，隨機向量過程模擬的POD 方法在本質上屬于傳統Monte Carlo 方法，然而，傳統Monte Carlo方法通常需要進行大量的隨機抽樣才能達到令人滿意的模擬精度，一方面，由于生成樣本數量巨大，極大地增加了結構動力反應分析的計算量；另一方面，Monte Carlo 方法在本質上屬于隨機抽樣方法，導致所生成的樣本概率信息不完備，無法進行結構精細化動力反應分析和動力可靠性評價。為克服傳統Monte Carlo方法的上述挑戰，文中引入基于隨機函數的降維思想，將正交隨機變量{Rik,Iik}定義為如下的基本隨機變量的函數形式［13］：

式中，i，s=1，2；k，l=1，2，…，N；Θi(i=1,2)為在區間(0,2π]上均勻分布且相互獨立的基本隨機變量。

顯然，式（29）所定義的隨機函數形式完全滿足式（28）的基本條件。這樣，通過降維模擬方法，實現了僅需2個基本隨機變量即可精細化模擬主余震向量過程，且生成的代表性時程具有賦得概率，所有代表性時程集合構成一個完備的概率集，為與概率密度演化方法結合實現復雜結構精細化動力反應分析與動力可靠性評價奠定了基礎。

5 模擬步驟及數值算例

5.1 模擬步驟

主余震向量過程的具體模擬步驟如下：

（1）確定相干函數矩陣γ(ω)。該矩陣的非對角元素根據第2節建議的主余震相干函數模型來計算。

（2）構造主余震型的演變功率譜模型。根據第3 節建議的參數取值，代入式（4）、式（5）和式（6）中，即可確定主余震型的演變功率譜模型。

（4）確定正交隨機變量{Rik,Iik}。將上一步驟中生成的Θi(i=1,2)的代表性點集代入到式（29）中，此處，需要注意的是，為實現主余震向量過程的模擬，應進行確定性的映射變換(s,l)→(i,k)。這一確定性的一一映射過程可以通過MATLAB 工具箱中的rand（‘state’，0）和temp=randperm(2×N)函數實現。這樣，即可唯一確定正交隨機變量{Rik,Iik}。

（5）生成主余震型地震動的代表性時程。將主余震的相干函數矩陣γ(ω)、演變功率譜模型以及正交隨機變量{Rik,Iik}代入到式（27）中，即可生成主余震型地震動的代表性時程，且每條代表性時程的賦得概率與初始代表性點的賦得概率相同，均為Pl。

5.2 數值算例

在本算例中，以II 類場地為例，參數取值為：主震模擬持時25 s，余震模擬持時20 s，時間步長0.01 s，截斷頻率，頻率步長，代表性時程數量144，主震峰值加速度，主震峰值因子3，根據表6，余震峰值加速度，余震峰值因子2.75。其余參數按照表5取值。

圖5為采用降維模擬方法生成的主余震向量過程模擬結果，其中，圖5（a）為主余震向量過程的代表性時程，圖5（b）和圖5（c）分別為主余震向量過程的均值、標準差模擬值與目標值的比較結果。為進一步說明降維模擬方法的精確性，分別計算了144 和987 條代表性時程的均值與標準差誤差，結果列于表7 中。從圖5（a）可以看出，生成的代表性時程具有明顯的主余震序列型地震動的特征，即主震和余震在持時、頻譜以及幅值等方面都顯示出明顯的差異性。從圖5（b）、圖5（c）和表7 可以看出，主余震向量過程的均值和標準差均與對應的目標值擬合良好，其中均值誤差為0，標準差誤差隨著代表性時程數量增多而減小，驗證了降維模擬方法具有較好的收斂性。同時，在僅有144條代表性時程時標準差誤差即小于5%，能夠滿足工程需求，驗證了降維方法具有較好的精確性。

圖5 II類場地主余震向量過程模擬結果Fig.5 Simulation results of main-aftershock vector process for site classification II

表7 主余震向量過程的模擬誤差Table 7 Simulation error of main-aftershock vector process

圖6（a）為Ⅰ1、Ⅱ、Ⅲ三類場地（144×3 組）代表性時程的主余震間相干性與488 組實測強震記錄的對比結果；圖6（b）及圖6（c）為Ⅱ類場地144 組代表性時程的主余震反應譜模擬值與實測強震記錄的對比結果。其中，對于反應譜，需要將主震實測記錄調幅至200 cm/s2，余震實測記錄調幅至60.61 cm/s2。

從圖6 可以看出，主余震代表性時程間具有良好的相干性，且與實測記錄擬合一致，證明了文中所建議的相干函數模型的正確性。同時，主震和余震的反應譜模擬值均與對應的實測記錄擬合良好，驗證了文中所提出的參數識別方法的有效性。

圖6 主余震相干性和反應譜的模擬值與目標值比較Fig.6 Comparison of the simulated values and the corresponding target values of the coherence and response spectrum of the main-aftershocks

6 結論

文中從實測強震記錄出發，對主余震間的相干性和演變功率譜模型的參數關系進行了系統分析，并通過引入隨機函數思想，實現了基于POD方法的主余震向量過程降維模擬。文中得出的主要結論如下：

（1）實測強震記錄顯示出主余震間具有明顯的相干性，且隨著頻率的增大相干性總體上呈現遞減趨勢。據此，文中建議了一種用于模擬主余震相干性的傅里葉級數模型，并通過與實測記錄對比驗證了該相干函數模型的有效性。

（2）文中對主余震的演變功率譜模型參數進行了識別，并給出了對應于I1、II、III類場地的主余震演變功率譜模型參數的建議取值。結果表明，相較于余震，主震的持時更長，能量更高，卓越頻率更低。

（3）將主余震視作一個向量過程，實現了主余震序列型地震動的一體化高效模擬。通過該方法，可一次性生成主震與余震的代表性時程集合，而不必將主余震作為2個隨機過程分別進行模擬，這樣不僅體現了主余震間的相干性，同時可大幅提高計算效率。

（4）基于隨機函數的降維方法僅需數百條代表性時程即可在全概率信息上反映主余震向量過程的概率特性，這為結合概率密度演化理論進行主余震作用下工程結構的精細化動力反應分析與可靠性評價奠定了基礎。

致謝：感謝中國地震局工程力學研究所為本研究提供數據支持！