FK 法合成地震动的频带范围研究＊

曹泽林 陶夏新 陶正如 王可意

1） 中国河北邯郸 056038 河北工程大学土木工程学院

2） 中国哈尔滨 150090 哈尔滨工业大学土木工程学院

3） 中国哈尔滨 150080 中国地震局工程力学研究所

引言

地震动合成广泛用于估计未来大地震产生的近场地震动.几十年来，多种理论和方法相继发展以用于合成地震动（Douglas，Aochi，2008）.常用的地震动合成方法可以分为数值方法、随机合成法、经验格林函数法、混合方法等.其中，数值方法受计算资源的限制，网格离散尺寸无法足够精细，难以保证合成地震动高频段的精度；随机合成的地震动在低频段又不够可靠.对于强震动观测记录的可靠频段和工程结构抗震分析所需地震动的低频、高频范围，宽频带（0.1—20 Hz）地震动成为当前的研究热点.目前，宽频带地震动合成常使用混合方法（Frankel，2009；Graves，Pitarka，2010），分别在低频段、高频段通过数值模拟、随机合成获得整个断层引起的地震动，再将两部分滤波、叠加得到宽频带地震动.其中，数值方法可以模拟三分量地震动，随机合成法只能提供一个无明确方位的水平分量，两种方法的叠加有不协调之处.此外，叠加所用交叉频率的选取具有主观性，两种方法对含义相同参数的表达也有不一致现象（孙晓丹，陶夏新，2012）.因此，需要研究计算精度高、逻辑严密的多分量宽频带地震动合成方法.

震源模型和格林函数是各种地震动合成方法的两个关键内容.Hartzell （1978）提出的有限断层震源模型广泛用于地震动合成，但震源描述形式各不相同.例如：经验格林函数法通过大震和小震的震源机制“调制”来表达震源（Irikura，Miyake，2011）；随机有限断层法用震源谱模型表达震源（Motazedian，Atkinson，2005）；包括本文方法在内的多种方法采用运动学震源模型，即，用有限断层面上的错动量分布、破裂速度、各子源时间函数和上升时间及滑动角等震源参数通过较为简单的数学表达描述颇为复杂的震源破裂过程.格林函数的计算是近场地震动合成的关键一环，也是不同地震动合成方法的核心区别.例如：经验格林函数法采用包含了真实场地效应和路径效应的小震（或余震）记录作为大震的格林函数（Irikura，Miyake，2011）；随机有限断层法用有限带宽白噪声描述随机格林函数（Motazedian，Atkinson，2005）.这些格林函数主要用于高频地震动合成，在低频段具有明显的偏差，且数值格林函数用于高频地震动合成时需要极大的计算资源（Graves，Pitarka，2010）.Wang （1999）、Zhu 和Rivera （2002）等提出的格林函数算法解决了水平成层介质内源引起的三维位移响应问题，在地震动合成研究中受到关注（Hartzellet al，2005；Kielinget al，2014；Sunet al，2015）.这类地震动合成方法具有严密的理论基础，可以避免混合方法存在的逻辑缺陷，优势在于能够表达地壳介质内地震波的复杂传播过程，难点在于震源建模，有研究称之为基于物理的方法（physics-based）（Sunet al，2015）或基于震源的方法（source-based）（Hartzellet al，2011）.从工程应用的角度看，这类方法所用的水平成层地壳速度模型是与当前对地壳结构的认知水平相匹配的，既可以表达地壳的主要分层结构，也具有较高的计算效率.

1 合成地震动的FK 法

1.1 FK 法的计算原理

根据位移表示定理，点源引起的地表一点的地震动位移可以表示为（Aki，Richards，2002）

1.2 本文采用的算例

为便于分析比较，本文采用图1 所示的一个算例.首先，设定MW6.5 直立走滑地震，断层长度和宽度分别为32 km 和16 km，上缘埋深为5 km.整个断层面划分为16×8 个2 km×2 km的子源.依据震源参数定标律（姜伟等，2017），平均错动量取为59 cm.有限断层震源模型的总地震矩为6.31×1018N·m，由于采用均匀错动分布，每个子源的地震矩相等，均为4.93×1016N·m.有限断层震源模型的一个子源的地震矩相当于一次MW5.1 地震的地震矩.本文将图1 所示子源视作点源，以便于后文分析.然后，设置距离断层中心地表投影30 km 的基岩场地计算点A，用于地震动合成和比较分析.地壳速度结构模型采用四川芦山地区的模型（Haoet al，2013）（表1），原因在于该地区研究基础良好，地壳速度模型的浅处有厚度分别为3 km 和5 km 的两个薄层，便于分析地壳速度结构对破裂时间和传播时间的影响.

图1 本文采用的算例Fig. 1 Calculation example adopted in this study

表1 本文采用的地壳速度结构模型Table 1 Crustal velocity structure model adopted in this study

1.3 FK 法对地壳速度结构的表达

FK 法采用的频率波数域格林函数表达了水平成层地壳模型的内源引起的响应，具体可参见相关文献（Zhu，Rivera，2002；Kennett，2009），本文不再赘述.除了格林函数，地壳速度结构还影响着震源参数.本小节讨论FK 法中地壳速度结构对子源地震矩分配和两个时滞的影响.

在常用的地震动合成方法中，以子源的错动量为权将整个断层的地震矩分配给各子源，表达式为

式（2）中的破裂时间表达子源破裂的触发时刻，一般由破裂起始点到子源的距离除以破裂速度得到，即

式中，Lij为破裂起始点到第ij个子源的距离，vr为破裂速度.破裂速度常按下式确定

式中，β为破裂速度与剪切波速vS的比值，常取为0.8.

[1]向邱.品管圈活动在呼吸重症监护病房的应用与效果[J].护理管理杂志,2013,13(2):104-105.

地壳介质剪切波速通常随深度变化，对破裂速度的分布和各子源的破裂时间有一定的影响.对于表1 的地壳速度模型和图1 的震源模型（此时断层埋深另取为0 km，以突出近地表波速复杂变化的影响），采用同一破裂速度2.7 km/s 和深度相关破裂速度，计算得到的破裂时间的分布如图2 所示.由图可见，考虑地壳介质剪切波速变化时，浅部破裂时间变大，深部破裂时间变小，差值最大达到0.5 s.这说明有必要考虑地壳结构对破裂速度的影响.

图2 不同破裂速度计算的破裂时间.星号表示破裂起始点（a） 同一破裂速度；（b） 深度相关破裂速度；（c） 深度相关破裂速度的分布；（d） 图（b）相对图（a）的差值Fig. 2 Rupture time calculated by different rupture velocity vr（a） The same vr ；（b） Depth-related vr；（c） Distribution of depth-related vr； （d） Difference between Fig.（b） and Fig.（a）

式（2）中的传播时间一般由子源到地表点的距离除以剪切波速得到，即

式中，Rij为第ij个子源到地表点的距离.

地壳介质中的剪切波速随深度发生变化，这会影响地震波传播速度，影响特点需要仔细分析.将震源设置于不同深度，对于不同震中距的地表点，采用同一剪切波速3.4 km/s 和表1 中的深度相关剪切波速计算S 波的传播时间，结果如图3 所示.由图可见，对于深度相关剪切波速，传播时间的等值线有多处转折且浅部和深部的增加趋势不同.图中的粗折线是两类初至波的分界线，左侧是震源向上传播的直达波先到达地表，右侧是震源向下传播的体波反射波先到达地表，粗折线水平段的位置正好对应地壳结构中速度不连续界面的深度，说明了地壳结构对传播时间的重要影响.考虑地壳介质剪切波速变化时，浅部传播时间变大，深部传播时间变小，差距最大达到3 s.考虑地壳结构对传播时间的影响也是有必要的，这对于FK 法合成地震动的可靠性和有效频带具有重要意义.

图3 不同剪切波速计算的S 波的传播时间（a） 同一剪切波速；（b） 深度相关剪切波速；（c） 图（b） 相对图（a） 的差值Fig. 3 Propagation time of S wave calculated by differentvS（a） The same vS；（b） Depth-related vS ；（c） Difference between Fig.（b） and Fig.（a）

2 格林函数计算的影响

考察式（4）可知，格林函数的计算式是严谨的解析表达式，对频率、波数的积分域均是从0 到无穷.在实际数值计算中，频率和波数均截止到某最大值，因此FK 法被称为一种半解析方法.本节分析频率和波数积分范围的截断处理对地震波计算的影响，并考察频率波数域格林函数是否能传播宽频带地震波.

2.1 关于频率的积分

一个离散的时程能够表达的最高频率称为奈奎斯特（Nyquist）频率.根据采样定理，奈奎斯特频率fN由下式确定，即

在格林函数计算中，Nt和dt的取值要考虑震源至地表点的距离范围和地表点的地震动持时，也要保证fN和df都满足宽频带计算的要求.在FK 法中，子源格林函数的零时刻是子源破裂的起始时刻.一般情况下，格林函数的总时长T都大于20 s.根据式（11），df小于0.05 Hz，完全能够表达格林函数的谱幅值随频率的复杂变化.

根据式（10），时间步距dt确定了格林函数能表达的最高频率，对FK 法合成地震动的频带具有显著影响.下面分析fN对合成地震动频带的影响.取Nt为4 096，dt为0.01，0.02，0.05，0.1，0.5 s，对应的fN为50，25，10，5，1 Hz.采用图1 中的点源，分别计算相应的格林函数，然后合成点A的断层垂直（fault normal，缩写为FN）分量的地震动.不同fN合成地震动的加速度幅值谱（Fourier amplitude spectrum，缩写为FAS）的比较如图4 所示.图中，左侧是上升时间为0 s 的结果，表达单位冲激函数源的影响；右侧是上升时间为1.5 s 的结果，表达考虑震源破裂过程的影响.为了更直观地表达合成地震动的频谱特征，本文采用了Konno 和Ohmachi （1998）的方法对幅值谱进行了光滑处理.由图4 可见，频率波数域格林函数可以表达fN以下的频率成分.当dt足够小，fN足够大时，格林函数能够传播宽频带地震波.在接近于fN时，格林函数一定程度偏小，这是由于数值计算的精度损失.对于0.1 Hz 以下的频段，fN较小的结果有一定的偏差，这是由于dt较大时，计算精度较差.对比图4a 和4b，震源破裂过程的时长（即上升时间）影响合成地震动的谱幅值水平，但在可靠的频率范围内没有中间频段缺失现象.dt除了确定格林函数的最高频率，还影响FK 法的计算量.表2 说明不同dt时FK 法合成地震动的计算量.可见，当dt从0.02 s 变为0.01 s 时，地震动合成的计算时间和数据存储空间都变为原来的2 倍（本文计算所用的台式计算机配置为CPU i5-10400、运行内存16 G）.对于本文算例，计算量差别不大.但对于汶川MW7.9 地震之类的大地震，除了距离范围大，断层面也需要划分为更多的子源，此时计算量的差异就十分显著了.因此，考虑到数值计算效率，dt取0.02 s，可以满足宽频带地震动合成的需要.

表2 不同dt 时FK 法合成地震动的计算量Table 2 Calculation cost of ground motion simulation using FK approach for different dt

图4 上升时间为0 s （a）和1.5 s （b） 时奈奎斯特频率对合成地震动加速度幅值谱FAS 的影响Fig. 4 Influence of Nyquist frequency on FAS of synthetic ground motion with rise time of 0 s （a） and 1.5 s （b）

2.2 关于波数的积分

FK 法的格林函数在波数域求解震源引起的响应，很好地解决了数值稳定问题（Zhu，Rivera，2002）.式（4）中关于波数的积分，数值计算截止到最大波数.零频最大波数kmax和波数间隔dk是提前设定的.理论上讲，kmax越大、dk越小，越有利于格林函数的计算精度，但dk变小时计算量会显著增加.以表1 中的地壳速度结构模型为例，表3 说明了不同dt、dk和kmax时FK 法格林函数的计算量.当dt从0.02 s 变为0.01 s 时，格林函数计算时间和数据存储空间分别为原来的4 倍和2 倍.下面分析kmax和dk对合成地震动频带和计算量的影响，除了确定参数的建议值，还可以了解参数的合理取值范围.

表3 不同dt，dk 和kmax 时FK 法格林函数的计算量Table 3 Calculation cost of Green’s function in FK approach with different dt, dk and kmax

kmax和dk的单位是1/Xmax，Xmax是点源埋深与点源至地表点的震中距两者中的较大值.对于图1 中的点源，Xmax为30 km.取dk为0.1，kmax取为2，5，10，15，20，分别计算相应的格林函数，然后合成点A的FN 分量地震动.不同kmax合成地震动的FAS 和位移时程的比较如图5a 所示.从图中可见，kmax过小时，0.1 Hz 以下和接近fN的频段的FAS 有一定偏差，其它频率的差距不明显.kmax越大，位移计算结果越精确.kmax对格林函数计算量的影响小，主要影响计算精度与稳定性.kmax取为10 以上时可满足计算精度要求，为避免特殊问题，本文建议kmax取为15.

取kmax为15，dk分别取为0.05，0.1，0.2，0.3，0.5.采用图1 中的点源，分别计算相应的格林函数，然后合成点A的FN 分量地震动.不同dk合成地震动的FAS 和位移时程的比较如图5b 所示.从图中可见，对于1 Hz 以下频段，dk较大时，合成地震动的偏差随频率减低而增加；在1—10 Hz 频段，dk的影响很小.从fN以上频段的FAS 可以发现，FAS 受dk的显著影响并不是单调变化的，这是因为在式（4）中，关于波数的积分因子随波数有明显的振荡现象，dk从小变大时，积分点并非单调变化.dk应足够小，以保证能够描述积分因子随波数的复杂变化，dk大于0.2 时，位移时程出现计算错误及失稳现象.表3 说明，当dk从0.1 变为0.05，0.3 时，格林函数的计算时间分别为原来的2 倍和0.33 倍.对于汶川MW7.9 地震之类的大地震，格林函数的计算量差距会十分显著.考虑到计算精度和计算效率，本文建议dk取为0.1，更高要求时可取为0.05.

图5 最大波数（a）和波数间隔（b）对合成地震动加速度幅值谱FAS 和位移时程的影响Fig. 5 Influence of the maximum wavenumber （a） and the wavenumber interval （b）on FAS and displacement time history of synthetic ground motion

2.3 地壳速度结构影响一例

地壳速度结构对格林函数的幅值、波形具有控制作用.在地壳速度结构中，品质因子是描述地壳介质耗能衰减的参数.下面以品质因子为例，分析其对格林函数频带的影响特点.对于表1的地壳速度结构模型，分别将品质因子取为原值的0.5，1.0，2.0 倍.采用图1 中的点源，分别计算相应的格林函数，然后合成点A的FN 分量地震动.不同品质因子合成地震动的FAS 的比较如图6 所示.由图可见，品质因子主要影响1 Hz 以上的高频地震动，影响程度随频率增加.这一分析表明，地壳结构影响不同频段谱幅值的总体水平，但不影响FK 法合成地震动可靠的频带宽度.

图6 品质因子对合成地震动加速度幅值谱FAS 的影响Fig. 6 Influence of quality factor on FAS of synthetic ground motion

3 震源破裂过程的影响

根据式（3），为了表达震源破裂过程释放地震矩需要一定时间、地震矩释放速率随时间的变化，地震动合成需要卷积格林函数与震源时间函数.震源时间函数（STF）描述地震矩率（或错动速率）的时间历程，目前，还很难对每一个子源都估计一个不同的震源时间函数.一般做法是采用相同的函数形式，对每一个子源估计各不相同的上升时间和地震矩.此外，破裂速度控制各子源破裂触发时间，也是一个重要参数.本节讨论震源时间函数、上升时间、破裂速度对合成地震动频带的影响，考察标准主要包括两方面：合成地震动的频带宽度是否覆盖0.1—20 Hz；在整个宽频带范围内，这些参数（或模型）是否会造成某些频段地震动幅值过小或过大.

3.1 震源时间函数的影响

迄今为止，地震反演和正演研究中已经有多种STF 解析函数表达式，包括简单的数学函数和具有震源动力学含义的函数.根据来源和函数形式，震源时间函数可以分为简单类型、三角函数类型、Brune 类型和Yoffe 类型.本文参考前期研究（Caoet al，2019），选取三角形STF、Liu STF （Liuet al，2006）、Hartzell STF （Hartzellet al，2007）、Brune STF （Brune，1970）和三角形规则化的Yoffe STF （Tintiet al，2005）等五种典型模型，分析震源时间函数对合成地震动频带的影响.图7 比较了五种模型表达的错动速率及其幅值谱.采用图1 中的点源，分别合成点A的FN 分量地震动.图8 比较了五种震源时间函数模型合成地震动的幅值谱.其中，幅值比较是以具有f−2高频谱衰减速率的Brune STF 的结果为参考值.

图7 五种震源时间函数的错动速率（a）及其加速度幅值谱FAS （b）Fig. 7 Slip rate （a） and its FAS （b） for five source time functions

图8 震源时间函数对合成地震动加速度幅值谱FAS 的影响Fig. 8 Influence of source time function on FAS of synthetic ground motion

一般来讲，STF 的函数形式对散射地震波的表达存在几方面问题（Dregeret al，2007），包括错动速率的对称性、错动速率FAS 中间频段幅值很小的谱洞、谱高频衰减速率，一般从这几方面分析震源时间函数的模型及其影响.在震源破裂过程中，错动速率通常先快速增加到最大值然后逐渐变小（Day，1982）.这说明STF 的错动速率随时间的变化不是对称的，STF 的错动速率应在错动早期具有一个有限的峰值.这描述了震源破裂过程释放能量随时间的变化，主要影响合成地震动的波形.STF 的幅值谱中较大的振荡称为谱洞，会导致频率相关地震动参数出现较大偏差.例如，三角形STF 合成地震动的FAS 在1 Hz 和2 Hz 处有明显低估，Yoffe STF 使相应FAS 在3 Hz 处亦有明显低估，正好对应STF 的FAS 的谱洞位置.震源破裂过程释放地震波的谱幅值随频率的变化应当是连续的，不应有剧烈振荡.因此，无谱洞的STF 更适用于描述震源破裂过程和地震动合成.此外，STF 应具有f-2高频谱衰减速率，衰减过快或过慢会低估或高估地震动，都不利于合成宽频带地震动，例如，Liu STF 在中频段明显较大的幅值会高估这个频段的地震动.由式（3）可知，对于FK 法的震源模型，上升时间含义明确的Hartzell STF 比Brune STF 更适用于描述子源错动的时间变化，也是本文推荐的最优模型.

3.2 上升时间的影响

上升时间表征破裂面上一点从破裂开始直至最终达到静态位错所历经的时间.根据相关研究（Brune，1970；Beresnev，2002），上升时间的频域等效参数是震源谱模型中的拐角频率.对于Hartzell STF，上升时间分别取为0，0.5，1.0，2.0，5.0 s，表示震源错动时间变化对合成地震动频带的影响.采用图1 中的点源，分别合成点A的FN 分量地震动.不同上升时间合成地震动的FAS 的比较如图9a 所示.从图中可见，上升时间是控制地震动幅值的最重要因素，比格林函数的控制作用要大很多.随着上升时间的增大，地震动谱幅值显著降低.这说明描述错动速率的上升时间控制了震源破裂释放高频地震波的能量，错动越快释放的高频地震动越强.结合图4 和图9 可知，上升时间显著影响合成地震动的谱幅值水平，但在可靠的频率范围内没有中间频段缺失现象.

图9 上升时间（a）和破裂速度（b）对合成地震动加速度幅值谱FAS 的影响Fig. 9 Influence of rise time （a） and rupture velocity （b） on FAS of synthetic ground motion

3.3 破裂速度的影响

根据式（2），震源时间函数和上升时间描述一个子源的错动过程，各子源的错动过程按一定时滞叠加组成整个震源破裂面的破裂过程.下面研究子源叠加是否影响合成地震动的频带.为此，分别将破裂速度取为剪切波速的0.6，0.7，0.8，0.9，1.0 倍，计算破裂速度和破裂时间.采用Hartzell STF 和1.5 s 上升时间，采用图1 中的整个有限断层震源模型，合成点A的FN 分量地震动.不同破裂速度合成地震动的FAS 的比较，如图9b 所示.从图中可见，在0.1—0.5 Hz 范围内，FAS 随破裂速度增加.这与预期相符，原因在于破裂速度越大，各子源地震动的时滞越小，合成地震动越强.在1 Hz 以上高频段，破裂速度对FAS 的影响不显著.在0.5—1.0 Hz 范围内，破裂速度的影响规律与低频段相反，分析发现，这是由于图1 中均匀错动分布引起的有限断层模型人为周期性.对于不均匀错动分布，人为周期性会大大减小.在实际地震动合成中，破裂速度对合成地震动可靠的频带宽度的影响较小.

4 讨论与结论

FK 法是基于有限断层震源模型的地震动合成方法，采用频率波数域格林函数表达地震波在地壳介质内的复杂传播过程，采用震源时间函数和上升时间描述子源的破裂过程，具有严密的理论基础.FK 法充分表达了地壳速度结构对格林函数、子源地震矩、破裂时间和传播时间的影响，更符合地震波传播规律.格林函数的计算分析表明，FK 法的格林函数具有表达宽频带地震波传播的能力.在计算格林函数时，时间步距、波数间隔以及零频最大波数应取值合理，以保证合成地震动的频带宽度、计算精度和计算效率，本文建议dt取0.02 s、dk取0.1、kmax取15、Nt按需取值.本文以品质因子为例，分析了地壳速度结构对合成地震动频带的影响.研究表明，地壳结构影响地震动谱幅值的总体水平，但可靠的频带宽度不变.震源破裂过程的影响研究表明，FK 法的震源破裂过程能够辐射宽频带地震波，这需要选取合适的震源时间函数、对上升时间和破裂速度进行合理约束.震源时间函数对震源辐射地震波的频谱成分具有控制性作用，应满足错动速率对称性、谱洞、谱衰减速率等方面的要求，本文建议选用Hartzell STF 描述子源错动过程.上升时间和破裂速度影响合成地震动的谱幅值的总体水平，在可靠的频率范围内没有中间频段缺失现象.子源地震动的叠加应注意避免人为周期性的影响.总之，在合理约束震源上升时间和破裂速度、选用震源时间函数和地壳速度结构模型时，FK 法合成地震动能够表达宽频带范围，在多分量宽频带地震动合成中有良好的应用前景.