地震動演變功率譜模型參數的統計建模

范穎霏 姜云木 劉章軍

摘要：工程上常將規范反應譜擬合的確定性取值作為地震動演變功率譜參數，這雖然便于工程應用，但同時存在明顯弊端：一方面，參數并非從實測強震記錄中獲得，具有較強的經驗性;另一方面，確定性參數生成的地震動樣本過于規則，且工程特性單一，無法全面地反映地震動的隨機性。為了克服上述難題，從PEER的NGAGWest2地震動數據庫中篩選1 766條主軸方向的實測強震記錄，并根據《中國地震動參數區劃圖》建議的場地類別和聚類方法劃分為15組;隨后，識別實測強震記錄的演變功率譜參數，并結合K-S檢驗和BIC信息準則，確定每個參數最優概率模型;最后，根據演變功率譜參數統計建模結果，結合降維方法，生成了Ⅱ類場地與Ⅲ類場地的代表性時程。與《建筑抗震設計規范》中建議取值不同，研究參數來源于實測地震動，且具有隨機性，避免地震動樣本工程特性的單一化和規則化。

關鍵詞：參數識別; 全隨機性地震動; 實測記錄; 統計建模; 降維模擬

中圖分類號： P315.9;TU311.3????? 文獻標志碼：A?? 文章編號： 1000-0844（2024）01-0126-11

DOI：10.20000/j.1000-0844.20220130001

Statistical modeling of evolutionary power spectral density parameters of ground motions

Abstract：?In engineering practices， the deterministic values of the response spectrum in accordance with seismic design codes are generally applied as evolutionary power spectral density （EPSD） parameters of ground motions. Although this treatment is convenient for engineering applications， it comes with two evident disadvantages. On the one hand， parameters are not obtained from measured strong motion records， which makes them highly empirical. On the other hand， the ground motion samples generated by deterministic parameters are extremely regular and exhibit unitary engineering characteristics， which cannot fully reflect the randomness of ground motions. To overcome the abovementioned challenges， we selected 1 766 strong motion records in the main axis direction from PEER and divided them into 15 groups based on the site categories suggested by the Seismic Ground Motion Parameters Zonation Map of China. Subsequently， the EPSD parameters of strong motion records were identified. Combined with the K-S test and BIC information criterion， the optimal probability models for each parameter were determined. Finally， based on statistical modeling results of EPSD parameters， representative time histories of site categories II and III were generated through the dimension reduction method. In contrast to the recommended values from the Code for Seismic Design of Buildings， the parameters suggested in this study were derived from measured records and possessed significant randomness， thus avoiding the simplification and regularization of engineering characteristics of ground motion samples.

Keywords：parameter identification; fully stochastic ground motion; measured records; statistical modeling; dimension reduction simulation

0 引言

相對于震源機制等地震學因素和“震源-傳播途徑-局部場地”的物理模型，工程上更加關注地震動本身的工程特性。演變功率譜模型從地震動的幅值、持時、頻譜等工程特性出發，是地震工程中應用最廣泛的隨機地震動模型。一般來說，演變功率譜模型由強度調制函數和平穩功率譜兩部分組成，強度調制函數對平穩信號進行調幅;平穩功率譜描述局部場地作用在白噪聲過程下的能量分布規律。在確定了演變功率譜模型之后，便可用諧波合成法［1］模擬地震動加速度過程。

目前演變功率譜的常用參數是根據規范反應譜擬合得到的確定性取值，這雖然便于工程應用，但是忽略了地震震級、傳播途徑、場地條件等因素的不確定性，從而忽略了演變功率譜參數顯著的隨機性。這將導致生成的地震動樣本過于規則，且工程特性單一，不能反映地震動豐富的概率信息和細部概率結構。

對于參數取值的確定性研究，已有較多學者以抗震規范為基礎，取得了較為完善的成果［2-5］。對于參數的隨機性建模，丁艷瓊等［6］基于工程隨機地震動物理模型，對7 000條地震動進行了精細的分組與參數識別，并對其進行了統計建模;李杰等［7］同樣基于工程隨機地震動物理模型，對5 000條地震動進行了參數識別與統計建模。雖然隨機性參數與功率譜的概念似有矛盾，但卻更符合工程上對地震動認知與應用。

為了體現震級、距離以及場地因素對地震動造成的影響，本文對實測強震記錄進行了精細的分組。利用最小二乘法，識別演變功率譜的模型參數;并結合擬合優度檢驗和貝葉斯信息準則，確定每個參數的最優概率模型;最后引入隨機函數約束［8-10］，僅用6個隨機變量便可實現地震動過程在概率層面上的模擬。由于地震動代表性時程具有賦得概率，可進而與概率密度演化理論相結合［11-12］，對結構進行精細的動力響應和可靠度分析。

1 地震動工程隨機模型及人工合成方法

根據文獻［13］，非平穩地震動過程的演變功率譜密度函數可表示為：

SUg（ω，t;λUg）=f（t;λf）2S（ω;λS） （1）

式中：SUg（ω，t;λUg）為非平穩地震動加速度過程Ug（t）的單邊演變功率譜;f（t;λf）為強度調制函數;S（ω;λS）為單邊功率譜。

對于強度調制函數，采用工程上常用的Amin-Ang模型［14］。該模型優點是可以較為全面地刻畫地震動的波形;缺點是采用了分段形式，形式復雜且無法準確地反映地震動峰值到達時刻［15］。

式中：t1、t2和c分別代表非平穩地震動過程的平穩段起始時間、平穩段結束時間和衰減段的衰減速度。因此，Amin-Ang模型的參數向量可表示為λf=（t1，t2，c）。

對于單邊功率譜，采用工程上常用的Kanai-Tajimi模型［16］。該模型優點是從動力學角度推導，參數具有明確的物理意義;缺點是S（0;λS）≠0，不符合實際情況：

式中：ωg和ξg分別為場地土的卓越圓頻率和阻尼比;S0為白噪聲激勵的強度，可表示為［17］：

式中：Amax表示地震動峰值加速度;r為峰值因子;ωe為S0=1時S（ω;λS）與坐標軸圍成的面積。因此，Kanai-Tajimi模型的參數向量可以表示為λS=（ωg，ξg，r）。

這樣，演變功率譜密度函數的參數向量λUg便可表示為：

λUg=（λf λS）=（t1，t2，c，ωg，ξg，r） （5）

確定演變功率譜密度函數SUg（ω，t;λUg）之后，便可采用諧波合成法生成地震動過程。根據文獻［10］，非平穩地震動過程的源譜表達為：

式中：Δω為頻率步長;Δω=ωu/N;N為截斷項數;ωk=k×Δω（k=1，2，…，N）;ωu為截斷頻率。

式（6）中，Xk，Yk為一組零均值的標準正交隨機向量，滿足以下基本條件：

E［Xk］=E［Yk］=0，E［XjYk］=0

E［XjXk］=E［YjYk］=δjk （7）

式中：E［·］為數學期望;δjk為Kronecker-delta記號。

傳統的諧波合成法在本質上屬于Monte Carlo方法：一方面，生成的樣本數量巨大，極大地增加了結構動力反應分析的計算量;另一方面，Monte Carlo方法在本質上屬于隨機抽樣方法，所生成的樣本概率信息不完備，無法進行結構精細化動力反應分析和動力可靠性評價。為克服傳統Monte Carlo方法的上述挑戰，本文引入了隨機函數的降維思想［8-10］，將標準正交隨機變量集{Xk，Yk}（k=1，2，…，N）定義為基本隨機向量的正交函數形式：

2 實測強震記錄

本文從太平洋地震動工程研究中心（PEER）NGA-West2地震動數據庫中篩選了1 766條主軸方向的實測強震記錄，且嚴格遵循以下原則挑選：

（1） 斷層距離應大于10 km，以減少近斷層地震動的數量。

（2） 實測記錄的矩震級應大于5，以排除對結構影響較小的地震動。

在篩選的實測強震記錄中，Rjb（Joyner-Boore距離）范圍為0～540 km;vS30（地下30 m平均剪切波速）范圍為106.83～1 525.85 m/s;M（矩震級）范圍為：5.0～7.9。

《中國地震動參數區劃圖》［18］中建議了Ⅰ0～Ⅳ五種場地類別。本文結合文獻［9］，根據每條實測強震記錄的vS30對實測強震記錄進行分類。同時，為了進一步體現震級與距離因素對地震動工程特性造成的影響，本文以M和Rjb為指標，采用K均值聚類的方法，將每個場地的實測強震記錄另分為3組，可概括為：遠場大震、近場小震和近場大震［6］。由于小震無法傳播到較遠的距離，因此沒有遠場小震這一項。

Ⅱ類場地的聚類結果如圖1所示?？梢钥闯?，聚類方法的分組結果與《建筑抗震設計規范（GB 50011—2010）》［19］中的地震分組有著相似性。

表1給出了不同聚類分組和場地對應的實測強震記錄數量：

此外，實測強震記錄還需要四階Butterworth濾波處理以及1%～99%范圍的能量截取，以保證研究結果的可靠性［20］。

圖2和圖3分別給出了不同場地的反應譜均值對比圖和Ⅱ類場地不同分組的反應譜均值對比圖。需要說明的是，為了對比清晰，將實測強震記錄均調幅至0.2g·m/s2。由圖2可知，隨著場地土變軟，反應譜峰值逐漸向長周期移動，且長周期成分也逐漸增多。由圖3可知，不同分組的反應譜峰值和長周期成分有著顯著差異，充分證明了震級和距離對地震動工程特性的顯著影響。

3 演變功率譜參數識別

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

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

對于非平穩地震動過程，其隨時間變化的歸一化模型能量曲線P（t）為：

以Ii（t）為目標值，采用最佳平方逼近原則，便可識別ai（t）的參數向量λf，i：

根據帕薩瓦爾定理，信號在時域上的能量與頻域上的能量相等，可以得到：

令

于是：

同時結合式（3）和式（4）得：

由于ai（t）的峰值加速度Amax，i是已知的，將參數向量λf，i代入式（14）中便可得到Ei，再根據式（15b）便可得到與ai（t）對應的峰值因子ri。

對于平穩地震動功率譜的場地參數，采用擬合反應譜的方法識別。為簡便起見，在對場地參數進行識別時，將地震動視為等效平穩過程。Vanmarcke將隨機過程的反應譜定義為單質點體系反應峰值系數的平均值與反應方差的乘積［22］，因此，地震動反應譜與功率譜轉換公式可以定義為：

式中：（ω0;λf，i）為等效平穩過程的峰值因子;σ（ω0，ξ;λS，i）為等效平穩過程的反應方差。ω0與ξ分別為結構的固有圓頻率和阻尼比，在本文中ω0≥1.05 rad/s，ξ=0.05;Td，i為等效平穩過程的持續時間，即為強度超過峰值50%的持續時間，對于Amin-Ang模型，Td，i的表達式為：

以ai（t）前6 s反應譜αi（ω0）為目標值，根據式（16），采用最佳平方逼近原則，對參數向量λS，i進行識別：

這樣，便完成了與ai（t）對應的演變功率譜參數向量λUg，i的識別。類似地，可依次識別1 766組實測強震記錄的參數向量λUg，i。

以SFERN地震為典型實例，根據式（11）與式（16），反應譜和地震動歸一化能量的結果見圖4（a）

與圖4（b），可以看出，模型與實測記錄均擬合良好，驗證了本文采用的識別方法的有效性。參數識別結果如表2所列。

4 演變功率譜參數統計建模

為了體現演變功率譜參數的隨機性，本文選取了5種備選概率模型對每個參數的概率分布進行擬合。這5種備選概率模型分別為：對數正態分布（LOGN）、耿貝爾分布（GUM）、廣義極值分布（GEV）、韋布爾分布（WEI）、伽馬分布（GAM），它們的概率密度函數如表3所列。

以Ⅱ類場地的參數ωg為例，概率分布擬合結果如圖5?？梢钥闯?，難以直接判斷ωg的最優概率模型。

為此，本文引入了K-S檢驗。一般來說，當K-S檢驗為0時，表示接受該備選概率模型;為1時，表示拒絕該備選概率模型。進一步地，當有多個備選概率模型通過K-S檢驗時，則根據貝葉斯信息準則進一步篩選。貝葉斯信息準則可表示為［24］：

BIC=ln（n）u-2ln（L） （21）

式中：u為概率模型中參數的個數;L為模型最大似然函數;n為樣本容量。一般地，BIC值越小，代表模型擬合度越好。

結合K-S檢驗與BIC信息準則，每個演變功率譜參數的最優概率模型以及對應的模型參數如表4所列。

在表4中，當概率模型為LOGN時，Par1為均值，Par2為標準差;概率模型為GUM時，Par1為位置參數，Par2為形狀參數;概率模型為GEV時，Par1為形狀參數，Par2為尺度參數，Par3為位置參數;概率模型為WEI時，Par1為比例參數，Par2為形狀參數;概率模型為GAM時，Par1形狀參數，Par2為尺度參數。

為了更清晰地表現參數隨場地變化的趨勢，圖6展示了每個參數均值不同場地的變化趨勢圖。

圖6（a）可知，隨著土質變軟，場地土卓越圓頻率逐漸變小;而峰值因子和場地土阻尼比沒有明顯的變化。對于場地土阻尼比，文獻［6］指出，阻尼比不僅與場地土的軟硬有關，還與其厚薄、孔隙率、形狀、體積等因素有關。因此不能通過簡單的場地類別與分組反映其影響的強弱。對于峰值因子，文獻［3］指出，峰值因子對場地變化不敏感。

圖6（b）可知，隨著土質的變弱，地震動的平穩段起始時間和地震動衰減因子逐漸縮小，而地震動的平穩段持續時間逐漸延長。這說明隨著土質的變弱，地震動的總持時逐漸延長。

綜上，建議將場地土阻尼比和峰值因子近似做均值處理，而其他參數做隨機變量處理。

5 數值算例

5.1 降維模擬的實現

由前文可知，模擬地震動共需要6個隨機變量{t1，t2，c，ωg，Xk，Yk}。且根據式（8），Xk，Yk均是同一個基本隨機變量的函數形式。這樣，模擬地震動過程共需要5個基本隨機變量。

式中：l=1，…，nsel;F-1ωg、F-1t1、F-1t2及F-1c分別代表任意分組下ωg、t1、t2及c最優累計分布函數的反函數。這樣，便得到了具有隨機性的演變功率譜模型參數。

式中：l=1，…，nsel。

這樣，便得到了地震動過程的基本隨機變量集Θ={Θ1，Θ2，Θ3，Θ4，Θ5}，將基本隨機變量Θ1～Θ4，場地土阻尼比均值和峰值因子均值代入式（1）中，便可得到非平穩地震動的演變功率譜模型SUg（ω，t;λUg）;將基本隨機變量Θ5代入式（8）中，得到隨機變量{Xk，Yk}。再根據式（6），便可模擬地震動樣本。

5.2 數值模擬

在本文算例中，以Ⅱ類和Ⅲ類場地為例。根據表4可知，場地土阻尼比ξg和峰值因子r的取值如表5所列。

其他模擬參數為：樣本數量144;頻率截斷項數N=1 600;截斷頻率ωu=240 rad/s;頻率步長Δω=0.15 rad/s;模擬持時40 s ;時間步長Δt=0.01 s。

圖7分別給出了Ⅱ類和Ⅲ類場地的代表性樣本?？梢钥闯?，樣本在頻譜和持時上均具有明顯的隨機性，更符合實際地震動。

圖8分別給出了Ⅱ類和Ⅲ類場地的模擬值均值和標準差與各自目標值的對比結果?？梢钥闯?，模擬值與目標值對應效果良好，證明了降維方法的有效性。為了進一步證明降維方法的精確性，表6中分別給出了Ⅱ類場地144、377和987條樣本的模擬誤差，可以看出隨著樣本數量的增加，均值及標準差誤差逐漸變小，體現了降維方法的收斂性;987條樣本的模擬誤差小于5%，符合工程要求，體現了降維方法的工程適用性。

圖9分別給出了Ⅱ類和Ⅲ類場地的模擬反應譜均值與實測反應譜均值的對比?？梢钥闯?，實測反應譜均值在模擬反應譜均值的±1倍標準差范圍內，且與其擬合較為一致，證明了本文參數模型的正確性和有效性。

6 結論

本文基于實測強震記錄，對演變功率譜模型進行了統計建模，得到以下結論：

（1） 為了反映震級以及傳播距離對地震動參數的影響，對實測強震記錄進行精細化分組，通過對比各組反應譜均值可知，震級以及傳播距離等因素對地震動有直接且顯著的影響。

（2） 對演變功率譜參數識別與統計建?？芍?，震級、傳播距離等不確定因素均間接地體現在了演變功率譜參數上。因此，有必要考慮演變功率譜參數的隨機性。

（3） 本文的參數模型來源于實測強震記錄，且具有顯著的隨機性，克服了傳統經驗化和固定化的參數取值而導致地震動樣本工程特性單一的缺陷;同時，采用降維方法生成的樣本具有良好的效率和精度。因此，將本文的參數模型和降維方法相結合，可高效且精確地模擬出能夠全面反映地震動工程特性的樣本。

參考文獻（References）

［1］ SHINOZUKA M，JAN C M.Digital simulation of random processes and its applications［J］.Journal of Sound and Vibration，1972，25（1）：111-128.

［2］ 任文杰，趙利軍，穆蒙蒙.現行抗震規范下地震動隨機模型參數研究［J］.世界地震工程，2017，33（2）：33-38.

REN Wenjie，ZHAO Lijun，MU Mengmeng.Study on parameters of random earthquake ground motion model based on current seismic code［J］.World Earthquake Engineering，2017，33（2）：33-38.

［3］ 張治勇，孫柏濤，宋天舒.新抗震規范地震動功率譜模型參數的研究［J］.世界地震工程，2000，16（3）：33-38.

ZHANG Zhiyong，SUN Baitao，SONG Tianshu.Study on the parameters of seismic power spectrum model based on the new seismic code［J］.World Information on Earthquake Engineering，2000，16（3）：33-38.

［4］ 白國良，朱麗華.基于現行抗震規范的Kanai-Tajimi模型參數研究［J］.世界地震工程，2004，20（3）：114-118.

BAI Guoliang，ZHU Lihua.Study on the parameters of Kanai-Tajimi model based on the code （GB 50011-2001）［J］.World Information on Earthquake Engineering，2004，20（3）：114-118.

［5］ 歐進萍，朱荻濤.設計用隨機地震動的模型及其參數確定［J］.地震工程與工程振動，1991，11（3）：45-54.

OU Jinping，ZHU Ditao.Random earthquake ground motion model and its parameter determination used in aseismic design［J］.Earthquake Engineering and Engineering Vibration，1991，11（3）：45-54.

［6］ 丁艷瓊，李杰.工程隨機地震動物理模型的參數識別與統計建模［J］.中國科學：技術科學，2018，48（12）：1422-1432.

DING Yanqiong，LI Jie.Parameters identification and statistical modelling of physical stochastic model of seismic ground motion for engineering purposes［J］.Scientia Sinica （Technologica），2018，48（12）：1422-1432.

［7］ 李杰，王鼎.工程隨機地震動物理模型的參數統計與檢驗［J］.地震工程與工程振動，2013，33（4）：81-88.

LI Jie，WANG Ding.Parametric statistic and certification of physical stochastic model of seismic ground motion for engineering purposes［J］.Journal of Earthquake Engineering and Engineering Vibration，2013，33（4）：81-88.

［8］ LIU Z J.Dimension reduction of Karhunen-Loeve expansion for simulation of stochastic processes［J］.Journal of Sound and Vibration，2017，408：168-189.

［9］ 劉章軍，劉子心，阮鑫鑫，等.地震動隨機場的POD降維表達［J］.中國科學：技術科學，2019，49（5）：589-601.

LIU Zhangjun，LIU Zixin，RUAN Xinxin，et al.POD-based dimension reduction representation of stochastic ground motion fields［J］.Scientia Sinica （Technologica），2019，49（5）：589-601.

［10］ LIU Z J.Random function based spectral representation of stationary and non-stationary stochastic processes［J］.Probabilistic Engineering Mechanics，2016，45：115-126.

［11］ LI J，CHEN J B.Stochastic dynamics of structures［M］.Chichester，UK：John Wiley & Sons，Ltd，2009.

［12］ LI J.The principle of preservation of probability and the generalized density evolution equation［J］.Structural Safety，2008，30（1）：65-77.

［13］ 劉章軍，劉子心.基于規范反應譜的全非平穩地震動過程模擬［J］.振動工程學報，2017，30（3）：457-465.

LIU Zhangjun，LIU Zixin.Simulation of fully non-stationary ground motion based on seismic design response spectrum［J］.Journal of Vibration Engineering，2017，30（3）：457-465.

［14］ AMIN M，ANG A H S.Nonstationary stochastic models of earthquake motions［J］.Journal of the Engineering Mechanics Division，1968，94（2）：559-584.

［15］ BOGDANOFF J L，GOLDBERG J E，BERNARD M C.Response of a simple structure to a random earthquake-type disturbance［J］.Bulletin of the Seismological Society of America，1961，51（2）：293-310.

［16］ KANAI K，YOSHIZAWA S，SUZUKI T.An empirical formula for the spectrum of strong earthquake motions Ⅱ［J］.東京大學地震研究所彙報，1961，41：261-270.

［17］ CLOUGH R W，PENZIEN J.Dynamic of structures［M］.New York：McGraw-HILL Inc，1993.

［18］ 國家質量監督檢驗檢疫總局，中國國家標準化管理委員會.中國地震動參數區劃圖：GB 18306—2015［S］.北京：中國標準出版社，2016.

General Administration of Quality Supervision，Inspection and Quarantine of the People's Republic of China，Standardization Administration of the People's Republic of China.Seismic ground motion parameters zonation map of China：GB 18306—2015［S］.Beijing：Standards Press of China，2016.

［19］ 中華人民共和國住房和城鄉建設部，國家質量監督檢驗檢疫總局.建筑抗震設計規范：GB 50011—2010［S］.北京：中國建筑工業出版社，2010.

Ministry of Housing and Urban-Rural Development of the People's Republic of China，General Administration of Quality Supervision，Inspection and Quarantine of the People's Republic of China.Code for seismic design of buildings：GB 50011—2010［S］.Beijing：China Architecture & Building Press，2010.

［20］ ASCE.Minimum design loads for buildings and other structures （ASCE/SEI 7-10）［S］.Reston：American Society of Civil Engineers，2010.

［21］ VANMARCKE E H，LAI S S P.Strong-motion duration and RMS amplitude of earthquake records［J］.Bulletin of the Seismological Society of America，1980，70（4）：1293-1307.

［22］ VANMARCKE E H.Structural response to earthquakes[J].Developments in Geotechnical Engineering，1976，15：287-337.

［23］ KOTZ S，NADARAJAH S.Extreme value distributions：theory and applications［M］.Singapore：Upper Saddle River：Prentice Hall，2000.

[24] SCHWARZ G.Estimating the dimension of a model[J].The Annals of Statistics，1978，6（2）：461-464.

［25］ SHIELDS M D，ZHANG J X.The generalization of Latin hypercube sampling［J］.Reliability Engineering & System Safety，2016，148：96-108.