马完君　卢建旗　唐丽华　谢志南

摘要：基于地震动能量传播速度为群速度这一概念，回顾了相位微分和窄带滤波两种频率相依波包时延计算方法及其在非平稳地震动相位拟合中的应用，归纳给出了一种非平稳地震动生成方法，通过模拟2条典型的海底与陆地地震动验证了方法的有效性，并将该方法应用于区域非平稳地震动场模拟中。结果表明：利用该方法所得模拟地震动记录平均加速度反应谱与原记录加速度反应谱有较高的一致性，而且能很好地再现地震的非平稳特性。

关键词：非平稳地震动;波包时迟;等效群速度;地震动模拟

中图分类号：P315.914 文献标志码：A 文章编号：1000-0666（2019）04-0523-08

0 引言

加速度地震动是非平稳时间过程，不同频率成分地震动组合得到窄带波包所携带主要能量的到达时间是不一致的。一般而言，高频成分的主要能量先到，低频成分的主要能量后到。结构弹塑性地震反应分析结果表明，地震动非平稳特性对结构反应具有重要影响（Wang et al，2012;张郁山，赵凤新，2014）。《超限高层建筑工程抗震设防专项审查技术要点》（中华人民共和国住房和城乡建设部，2015）要求高度超过200 m的各类建筑结构需要进行大震弹塑性时程分析（王亚勇，2017）。因此，有必要研究可计入其非平稳特征的地震动模拟方法，为重要工程结构抗震设计以及大中型城市地震灾害情景构建提供合理的地震动输入。

地震动模拟方法有2类：地震学方法和工程方法。地震学方法一般应用于模拟区域地震动场，可大致区分为：确定性方法、随机方法和混合方法。确定性方法受限于地下介质资料不足，目前大多用于1 Hz或以下频带范围地震动的模拟（Roten et al，2011）。随机方法产生的相位受限于高斯白噪声是平稳随机过程，模拟结果难以反映地震动的非平稳特性（Boore，2003，2009;Dabaghi et al，2018）。混合方法在如何实现高低频地震动的合理衔接以保证地震动非平稳特征的一致性上仍需进一步研究（孙晓丹等，2011;方小丹等，2014）。工程方法一般应用于模拟工程结构场点位置处的地震动，但它一般不考虑震源以及地震动传播的物理过程。为此，廖振鹏和魏颖（1988）提出了一种耦合工程方法与地震学方法的非平稳地震动模拟方法，基于能量传播速度为群速度这一概念，建立了可用于统计区域范围强地震动相位规律的等效频散模型，进而提供了一个比相位差谱更容易理解的物理框架，以给出统计平均意义上符合区域地下介质波传播特征的强地震动相位生成技术。

本文是廖振鹏和魏颖（1988）研究工作的延拓。首先回顾了窄带波包能量行进速度——群速度与波包时延这2个概念，然后对比分析了已有的波包时延计算方法，采用Boore（2003）给出的相位微分法计算波包时延并建立等效群速度模型。基于该模型阐述了地震动非平稳相位的计算及地震动模拟方法。以典型非平稳地震动为例，通过分析模拟地震动和实际地震动对应的反应谱和RotD50残差分布，验证了该方法的合理性。在此基础上，探讨了其在区域非平稳地震动场模拟中的应用。

1 窄带波包时延与群速度

1.1 波包时延的定义

基于傅里叶变换，可将地震动记录分解为若干简谐波的叠加。考虑震源附近地震动记录的傅氏变换如下：

式中：Re[z]表示复数z的实部;exp表示自然对数为底的指数函数;ts（ω）表示震源附近处频率为ω的地震分量起始时间。相应地，可将从震源附近辐射至距离R处地震动记录表示为：

式中：HR（ω）涵盖了几何扩散与传播介质对幅值的耗散等效应;R/cphase（ω）表示频率为ω的地震动分量以相速度cphase（ω）传播至R处所需要的时间;cphase（ω）=ω/k（ω），k（ω）表示空间波数。由于实际介质存在阻尼，介质中不同频率成分波的相位传播速度不一致。进一步考虑以ω0为中心频率的窄频带（ω0-Δω，ω0+Δω）地震动分量的叠加，即：

假定在对Δω取值较小情形下，有：

式中：HR（ω0）表示HR（ω0）的绝对值;ζ=（t-R/cg）Δω;cg=dω/dkω=ω0表示群速度。取HR（ω0）=π，Δω=π/10，ω0=4π，k0=π，cg=3 km/s，R分别取为45 km和60 km，图1给出上述窄频带地震动分量组合的传播示意图，从图中可知，窄频带地震动分量组合表现为波包形式，其幅值包络由（sinζ/ζ）函数所控制，由于这一地震动的能量主要包含在波包内部，可以看出波包传播的速度为群速度即能量传播速度。进一步分析式（5）可知，包络函数峰值所在位置对应于t-R/cg=0处，给定地震动起始时间，将包络函数峰值所在时间与地震动起始时间差值定义为包络时延时间。

1.2 波包时延计算方法

从上述定义可知，最简单的波包时延计算方法可通过对地震动记录进行窄带滤波得到（廖振鹏，魏颖，1988），滤波后记录的峰值所在位置即为波包时延（方法1）。以1990年加利福尼亚州高地5.6级地震中的海底地震台S3E水平加速度记录（Boore，2003）为例，选取中心频率f=ω/2 π=0.5，1.0，2.5，5.0，10.0 Hz对该记录在频域内用宽度为0.2 Hz的矩形窗滤波，然后返回至时域。图2分别给出了S3E原记录与窄带滤波后所得记录，结果表明：窄带滤波后的记录包含若干幅值差异较小甚至相当的波包，直接从滤波记录中提取包络时延存在不唯一的问题。另外这一方法需要使用多次傅里叶逆变换，极大地增加了計算量。针对这一问题，廖振鹏和金星（1995）首先建立了波包时延与记录傅里叶变换后所得相位导数之间的关系（方法2），即：

式（8）相当于在局部范围内对相位进行了解卷绕，然而相位解卷绕缺乏严格的数学基础（Shatilo，1992）。Boore（2003）建议了一种基于相位导数的波包时延计算方法（方法3）：

式中：Im[z]表示复数z的虚部;H（f）为原记录h（t）的傅里叶变换;G（f）为原记录与时间t的乘积所得记录t×h（t）的傅里叶变换。

针对S3E记录，图3给出了上述3种方法的计算结果，从图中可以看出，3种方法计算所得结果均反映了S3E记录中呈现的高频地震动分量先到、低频地震动分量后到的特征。第一种方法由于选取中心频率点较少，计算结果较为稀疏，后2种方法所得计算结果大致相同。但由于Boore（2003）的方法避免了相位解卷绕以及多次傅里叶逆变换，因此本文利用这一方法提取记录波包时延。

2 等效群速度模型与地震动拟合

基于波包时延的概念，在震源辐射出来的波包中，将以最快速度到达记录台站的波包所经历的时间t0作为地震动起始时间，得到不同中心频率成分地震动波包行进的等效群速度为：

式中：R是震中距;t0=R/Um，Um是等效群速度的最大值;t（f）为相对于t0的波包时延。考虑到地震动受震源物理过程、波在深部地下介质和浅层场地中传播过程的综合影响，且一般而言传播距离不等于震中距，称U（f）为等效群速度。因此基于地震动记录提取等效群速度，做如下处理：首先对地震动记录进行基线校正和高通滤波处理（Boore，2003），因为极低频地震动记录有可能是失真的。然后根据99%能量持时提取记录的强震段td=t2-t1，即根据I（t1）=0.001，I（t2）=0.991来确定，其中能量积分为：

式中：t1=t0。再由式（9）计算波包时延，这样保证波包时延计算用到的记录起始时间为最快波包到达台站的时间，然后利用式（10）计算等效群速度。考虑到地震动传播具有一定的随机性，可将等效群速度表示为：

式中：U（f）表示等效群速度随频率变化的平均趋势;η是一个具有零均值的随机过程。考虑如下U（f）拟合公式（廖振鹏，金星，1995）：

式中：a0，a1，a2系数可通过最小二乘或其他优化方法求解。一旦求得U（f），η的样本η便可以确定，但是得到的样本均值可能不是零，所以要对η进行微调：

同时还需要调整式（13）的求得的系数，即ai=ai（1+ε/fm），最后基于修正后η（f）样本确定η这一均匀分布的取值范围，本文取为95%置信区间-1.96S[η（f）]，1.96S[η（f）]，S[η（f）]为η（f）样本的标准差。图4中给出了基于S3E记录求取对应等效群速度模型的流程图，最终得到的参数分别为a0=1.110 9，a1=3.332 1，a2=

-1.002 9，S[η（f）]=0.285 7。基于群速度模型，建立地震动记录的相位信息，进而给出模拟地震动。具体步骤包括：①基于式（12）计算出等效群速度样本;②利用cphase（ω）=f∫f0[1/U（f）]df，计算得到相速度，在此基础上计算相位φ（f）=-2πfR/cphase（f）+2πfR/Um;③结合基于模拟或实际记录得到的地震动幅值谱，基于傅里叶逆变换给出模拟地震动。

3 基于等效群速度模型的地震动模拟

以上述S3E海底地震动以及1986年11月14日7级地震为例，利用震中距为79 km的Smart-1中心台站三分量观测地震动，基于上述方法给出符合记录非平稳特征的模拟地震动。结果表明：模拟地震动在统计平均意义上与实测记录的非平稳特征基本一致。最后讨论了该方法在区域地震动场模拟中的应用。

3.1 典型海域與陆域场点非平稳地震动模拟

在上文给出的S3E记录的一个相位拟合样本的基础上，在-1.96S[η（f）]，1.96S[η（f）]区间对η（f）进行随机采样，得到50组不同的相位样本，然后将相位与原记录傅里叶幅值谱结合，基于傅里叶逆变换得到模拟地震动。图5给出了S3E原记录与50组模拟地震动中的一个样本，以及50组模拟地震动记录计算所得加速度反应谱的平均值与原记录加速度反应谱的对比图，结果表明：模拟地震动记录不仅在时间过程上较好地反映了原记录的非平稳特性，且模拟地震动记录所得加速度反应谱的平均值与原记录加速度反应谱基本吻合。与之类似，本文将上述方法用于拟合Smart-1中心台地震动三分量观测记录，对每一个地震分量给出50次模拟结果。表1给出对应于三分量记录的群速度模型，地震动模拟结果如图6所示。此外本文还将水平两方向的50次模拟地震动进行了合成，求出其RotD50与原始的记录的

RotD50的残差（Boore，2010），结果如图7所示，本文建立的地震动模拟方法在保证反应谱拟合精度的前提下，很好地再现了地震的非平稳特性。

4 结论与讨论

本文回顾了相位微分和窄带滤波2种频率相依波包时延计算方法，阐明了窄带波包能量行进速度——群速度与波包时延二者之间的关系，基于波包时延建立等效群速度模型及其在非平稳地震动相位拟合中的应用，并建议根据Boore（2003）给出的相位微分法计算波包时延。在此基础上，建立了一种非平稳地震动生成方法，讨论了该方法在场点地震动模拟中的应用，利用2条典型的海底与陆地地震动记录模拟验证了这一方法的有效性，给出了该方法在区域非平稳地震动场模拟中的推广措施。该方法可进一步应用于历史地震与设定地震作用下非平稳地震动场的模拟，为重要结构弹塑性地震响应分析、大型城市地震灾害评估与预测提供合理的地震动输入。后续研究工作还需深入分析断层破裂效应以及不同类别场地或特殊场地如盆地效应对地震动场的影响。

本文的研究工作与Kumari等（2018）、廖振鹏和魏颖（1988）的结论基本一致，不同的是本文仅考虑点源，并且等效群速度模型计入了距离的影响，今后的研究有必要将二者结合，考虑震源在断层面上的破裂过程，还可进一步探讨剔除震源与场地影响的等效群速度模型构建方法。

參考文献：

方小丹，魏琏，周靖.2014.长周期结构地震反应的特点与反应谱[J].建筑结构学报，35（3）：6-23.

廖振鹏，金星.1995.强地震动相位的一个随机模型[J].地震学报，17（3）：353-361.

廖振鹏，魏颖.1988.设计地震加速度图的合成[J].地震工程与工程振动，8（1）：12-28.

孙晓丹，陶夏新，刘陶钧.2011.近场地震动估计中高，低频地震动的叠加[J].自然灾害学报，20（4）：84-89.

王亚勇.2017.结构时程分析输入地震动准则和输出结果解读[J].建筑结构，47（11）：6-11.

张郁山，赵凤新.2014.基于希尔伯特变换的非平稳地震动模拟方法的验证[J].地震学报，36（4）：686-697.

Boore D M.2003.Phase derivatives and simulation of strong ground motions[J].Bulletin of the Seismological Society of America，93（3）：1132-1143.

Boore D M.2009.Comparing stochastic point-source and finite-source ground-motion simulations：SMSIM and EXSIM[J].Bulletin of the Seismological Society of America，99（6）：3202-3216.

Boore D M.2010.Orientation-independent，nongeometric-mean measures of seismic intensity from two horizontal components of motion[J].Bulletin of the Seismological Society of America，100（4）：1830-1835.

Dabaghi M，Der Kiureghian A.2018.Simulation of orthogonal horizontal components of near‐fault ground motion for specified earthquake source and site characteristics.Earthquake Engineering & Structural Dynamics，47（6）：1369-1393.

Kumari N，Gupta I D，Sharma M L.2018.Synthesizing nonstationary earthquake ground motion via empirically simulated equivalent group velocity dispersion curves for Western Himalayan Region[J].Bulletin of the Seismological Society of America，108（6）：3469-3487.

Razafindrakoto H N，Bradley B A，Graves R W.2018.Broadband Ground‐Motion Simulation of the 2011 MW6.2 Christchurch，New Zealand，Earthquake[J].Bulletin of the Seismological Society of America，108（4）：2130-2147.

Roten D，Olsen K B，Pechmann J C，et al.2011.3D simulations of M7 earthquakes on the Wasatch fault，Utah，Part I：Long-period（0-1 Hz）ground motion[J].Bulletin of the Seismological Society of America，101（5）：2045-2063.

Shatilo A P.1992.Seismic phase unwrapping：methods，results，problems[J].Geophysical prospecting，40（2）：211-225.

Wang J，Fan L，Qian S，et al.2002.Simulations of nonstationary frequency content and its importance to seismic assessment of structures[J].Earthquake Engineering & Structural Dynamics，31（4）：993-1005.

中華人民共和国住房和城乡建设部.2015.超限高层建筑工程抗震设防专项审查技术要点[S].

A Method for Non-stationary Ground Motion Simulation Basedon the Equivalent Group Velocity Model

MA Wanjun1，LU Jianqi1，TANG Lihua2，XIE Zhinan1

（1.Key Laboratory of Earthquake Engineering and Engineering Vibration，Institute of Engineering Mechanics，China Earthquake Administration，Harbin 150080，Heilongjiang， China）（2.Earthquake Agency of Xinjiang Uygur Autonomous Region，Urumqi，830011，Xinjiang， China）

Abstract

Based on the concept that the energy propagation velocity of ground motion is group velocity，this paper briefly reviews the phase differential and narrow-band filtering methods for calculating frequency-dependent wave packet time delay and their applications in non-stationary ground motion phase simulation.Then a non-stationary ground motion simulation method is established.The effectiveness of this method is verified by simulating two typical ocean-bottom and land ground motions.Finally，the application of this method in regional non-stationary ground motion field simulation is discussed.The results show that the average acceleration response spectrum of the simulated ground motion recorded by this method is consistent with the original recorded acceleration response spectrum.In addition，the non-stationary characteristics of the earthquake can be well reproduced.

Keywords：nonstationary ground motion;envelope delay;equivalent group velocity;ground motion simulation