引言

1679年发生的三河平谷M8大地震是目前首都圈有记录的震级最大的地表破裂型地振动，地震震中烈度高达Ⅺ级。其震中位置位于河北省夏垫镇，震源方位为38.23°，倾角为82.54°，滑动角为−156.08°，地震产生的断裂断层从东流河屯附近开始，向东蔓延穿过夏垫镇北部，一直到东兴庄附近，地震极震区面积达1 800 km²，影响到整个华北地区，造成了巨大的人民生命与财产损失。其中破坏最严重的是三河、平谷地区，余震也波及到房山、延庆、昌黎、定县、新城、天津等20余县市（彭一民等，1981；向宏发等，1988；徐锡伟等，2000；刘保金等，2009；Wang et al，2014）。

地震区的地质构造情况也是构建地震模型的重要部分（Borcherdt，1991），近年来，三河—平谷的地震构造也成为许多研究的主题。地质结构研究表明，三河—平谷地震发生在北京东部的华北俯冲带与燕山隆起带交会的平原，存在断裂带的复合，地质构造和沉积形态复杂，呈现出强烈的新构造运动，并且在此地区附近小震活动频繁，地震震中呈带状分布，震源深度通常比较浅（田优平等，2014）。三河—平谷地震的发震断层上因有新生代沉积覆盖，不能完全在地表出露，这可能是造成此次地震震级较大的原因之一（彭菲等，2020）。震区上地壳存在部分不连续，倾斜不一致的反射叠层，展示出了此处地质构造的复杂性。从剖面表现出的特征来看，此处断裂构造体系存在有规律的脆-韧转换，容易积累应变，从而导致地震发生（刘保金等，2009）。

前人对京津唐地区地壳结构与强震发生关系的研究表明，三河—平谷地震与唐山地震有着非常相似的构造背景，均发生在速度变化较快的位置（孙若昧等，1993；张先康等，2002；于湘伟等，2010）。Huang等（2021）提出的拟波形地震层析方法基于北京及周边地区的多个高精度宽频数字地震台站所积累的丰富的高精度地震资料，利用数值求解波动方程反演了首都圈地区地壳的三维纵波、横波速度和泊松比结构，得出这两次大地震的震源区域有很大的相似性，为上述结论提供了新的证据。张先康等（2002）的研究表明，三河—平谷震区位于一个横波速度高、低速的过渡带，其速度结构与附近区域有着明显差异，在水平和垂直方向都表现出较强的非均匀性，于是推断其存在由上地幔物质上涌随后冷却而形成的岩脉或岩体，岩浆活动引起的局部应力分布差异可能导致了地壳深部断裂的形成，造成了三河—平谷大地震的发生。Huang等（2021）也持同样观点，其发现高速区域可能与三河—平谷地震震源下方发现的局部强反射能量群相对应。由此也可看出，地壳非均匀性是影响地震发生的重要因素。大地震的发生不仅是一个简单的力学问题，与地壳的物理化学性质也密切相关。

对历史地震震级的确定方法因地震发生年代及地震记录条件的不同而不同。对有台站记录的近代地震，依然存在地震参数不准确的问题。为了修正这些地震参数，确定地震震级，有些学者提出了相关衰减方程（张方浩等，2016；Lin，Wang，2020；Fu et al，2022）、回归公式（陈培善，刘家森，1975）、地震密度计算（王健，2001）、经验地震动预测方程（Peng et al，2014）和高频GPS重建地震波（Elósegui et al，2006；Wang et al，2007；Shi et al，2010；Yin et al，2010；Bock et al，2011；Zhang et al，2012；方荣新等，2013；Hung，Rau，2013）等方法。而年代久远的历史地震只能靠当时的文字记录来推断地震强度，进而估算历史地震震级。按前人的研究方法主要从四个方面估算历史地震震级，即震中区烈度、等烈度线面积、地表破裂长度和地震影响面积。1679年三河—平谷地震震级就是依据 《乾隆三河县志》 、 《民国23年平谷县志》 等相关史料记载，利用地震烈度对照方法进行评估的（李善邦，1957）。

之后也不断有学者进行过三河—平谷地震震级相关的研究（秦四清等，2016；殷娜等，2021），发现依据现有的史料记载来评估其所处的地震区烈度，进而估算的震级，尚存有一定的争议。其次，由于存在第四系强覆盖等的影响，对于此次地震的同震破裂长度以及端部变形特征尚未有确切定义，大都估计可确认地表断裂长度在10余千米（向宏发等，1988；徐锡伟等，2000；殷娜等，2021），通过研究极震区地表断裂长度等相关参数可以看出，这显然与全球范围内震级和破裂长度之间的传统统计关系不符（Wells，Coppersmith，1994），即并不能与M_W8.0震级有较好的对应，并小于震级与断层破裂长度的统计关系给出的长度，这点说明三河—平谷地震在震级与断层破裂长度统计关系的散点端，这是三河—平谷地震的一个特殊性，仅该点就足以说明该地震非常值得研究。杨百存等（2017）通过对获取的地震数据进行部分校正，认为三河—平谷地震为M_W7.8。由此可见，其发震构造和机制仍需深入研究，因此本文拟在前人研究的基础上，用多震级、多破裂方式模拟三河—平谷地震的强地面运动，开展对其震级的研究，以期达到确定三河—平谷地震真实震级的目的。

1.   方法

随机有限断层法是一种很常用的地震动模拟方法，在模拟过程中可以将研究的有限断层划分为多个子断层，每个子断层可看作一个点源，考虑时间延迟后叠加每个小点源的模拟结果便可以得到整个断层的地震动模拟结果（Beresnev，Atkinson，1997）。但是传统的随机有限断层法模拟结果大多在高频段地震动可靠，基于此，Zhou和Chang （2019）提出的随机有限断层法（non-uniform stress and non-uniform window function simulation method，缩写为NNSIM）主要作了两个修改：一是将每个点源的特定时间函数替换为不同源时间函数，以提高低频结果的可靠性；二是引入非均匀的应力降，来降低模拟结果对单个应力降的依赖性，使得最终结果与实际记录更加贴合。

NNSIM随机有限断层法也沿用了传统随机有限断层法的累积思想，将大断层分为N个可视为点源的子断层，在破裂延迟和传播延迟的作用下，将每个子断层的模拟结果叠加得到整个断层的地震动效应，即

式中，nl和nw分别表示整个断层长度和宽度范围内分出的子断层的个数，Δt_ij为第ij个子断层点源到观测点的相对延时，a_ij（t）为在第$ ij $个子断层上用随机点源法计算得到的地震加速度，a（t）为大断层在观测点的模拟加速度。Zhou和Chang （2019）详细说明了时窗函数在低频中的作用，在此我们使用这种方法，使模拟的地震动效应中的低频部分与实际记录更加接近，从而提高全频段模拟结果的可靠性。此处的拐角频率我们使用非均匀拐角频率（Zhou，Chang，2019），则第ij个子块的角频率为

式中Δσ_ij和M_0ij分别表示第ij个子断层的应力降和地震矩。Δσ_ij由Andrews （1980）推导出的位错与应力降之间的数学关系式计算得到。

传统的随机有限断层法不断得到优化，其中一个创新就是引入了动拐角频率（Motazedian，Atkinson，2005），它是对静拐角频率乘以不同时刻破裂子断层的个数开三次方而得到，由此使得构成各子单元的震源谱的拐角频率不再相同。式（2）的物理意义明显与动拐角频率不同，但同样使得每个子断层有不同的拐角频率，且对于断层划分尺度无依赖性。关于NNSIM随机有限断层法的具体公式、变量定义等详细内容，可见之前相关研究（周红，2018；Zhou，Chang，2019；周红等，2021），作者在这三篇文章中分别用该方法对2017年九寨沟M_S7.0地震、2013年芦山M_S7.0地震和2021年云南漾濞M_S6.4地震进行了地震动模拟，且模拟结果与实际结果重合度高，说明了此方法的有效性。现将此方法应用于1679年三河—平谷古地震的研究中。

2.   多震源的设定

通常大家都认为三河—平谷地震是一次M_W8.0巨震，但是否真的如此，本文将采用多震级多破裂震源模型地震动模拟的结果对比历史烈度，达到确定真实地震震源的目标。

三河—平谷地震发震断层为夏垫断裂。从震级与断层长度的经验公式我们可以看出，若三河—平谷地震的震级为M8，则其断层长度应具有一个比较大的规模（董瑞树等，1993；孙银涛等，2016），但是在诸多研究中计算得到的或借鉴采用的三河—平谷地震的断层长度多数小于震级与断层长度的统计关系，也体现出三河—平谷地震的实际震级应当小于M8。Wang等（2014）结合模拟退火算法和高斯—牛顿算法，反演得出三河—平谷地震断层长度约为80 km，在其它研究中三河—平谷的研究断层长度取63—200 km不等，且在60—90 km居多（向宏发等，1988；高孟潭等，2002；潘波等，2006，2009；刘博研等，2007；付长华等，2012，2015；何付兵等，2013；朱耿尚，2014；巴振宁等，2022）。同时有研究指出，震级并不一定与断层长度完全呈正比关系，更多取决于锁固段体积、剪切弹性模量与应力降或应变增量（秦四清等，2010；杨百存等，2020；靳天伟等，2021）。于是我们在前人研究的基础上，将此次研究区域取为夏垫断裂中间部分，总长度约为80 km，宽度为25 km。震级的选择我们按0.1级增加量，将M_W7.5至M_W8.0分成6个震级，每个震级都按照多震源、多破裂方式进行模拟。Somerville等（1999）研究中指出M7.5以上地震进行模拟时应采用两个凹凸体，我们考虑到三河—平谷地震虽然是历史地震，但依然是华北地区著名的大地震，所以在此次研究中我们设置了两个凹凸体进行模拟。每个震级的破裂模式的确定，按照Somerville等（1999）获得的震源设定方式给定。设定方式原则如下：凹凸体个数2.6个；凹凸体面积占总破裂面积的22%；凹凸体破裂地震矩占总破裂地震矩的44%；宏观震中在主凹凸体之外；主凹凸体面积占总破裂面积的17.5%。依据如上原则，我们把夏垫断层分成80个5 km×5 km的子断层，其中沿走向为16个，沿倾向为5个。主凹凸体大小为350 km²，占14个子断层，次凹凸体大小为90 km²，占4个子断层，凹凸体位错为平均位错的两倍，剩下的62个子断层为背景位错。主凹凸体和次凹凸体在80个子断层上充分排列组合，其中主凹凸体占14个，次凹凸体占4个，且主、次凹凸体不能重合，也不能在计算区域边缘，按照这种分布形式，借助排列组合的计算方法，得到满足条件的破裂模式共180个。基于每种破裂方式的位错分布，计算出相应的应力降分布、破裂速度分布、上升时间分布及不同的震源函数，即每种位错分布可以导出复杂性地震震源。先前探槽研究显示三河—平谷地震为倾滑型地震，且走滑分量不可忽略（向宏发等，1988；张先康等，2002），所以本研究中凹凸体模型兼顾走滑与倾滑的性质，取其两个分量的平均，即45°滑动方向。图1给出了随机选择的四个震源的破裂参数分布，其中地震矩分布（图1a）可由式（3）计算得到，即

图  1  震源1，73，102和180的破裂模式所对应的强地面运动模拟结果

图（c），（d）分别为为断层破裂面每个破裂子元的破裂延时和破裂初时

Figure  1.  Simulation results of strong ground motion corresponding to rupture modes of seismic sources 1，73，102 and 180

Fig.（c） and Fig.（d） are the rupture delay and rupture initiation time for each rupture subelement of the fault rupture surface

下载: 全尺寸图片 幻灯片

式中，$\mu $为介质剪切模量，A为破裂面面积，$ {D}_{ij}$表示每个单元的位错。应力降分布（图1b）按Ripperger和Mai （2004）给定的式（4）计算，即

式中，${\boldsymbol{k}} $为空间波数矢量，D为波数域位错，Δσ_∥和Δσ_⊥分别为沿走向和倾向的应力降，K_∥和K_⊥分别为走向和倾向的刚度函数。图1c为断层破裂面每个破裂子元的破裂延时，它由方程T_ij＝1/f_0ij得到。图1d为断层破裂面每个破裂子元的初始破裂时间，由t_ij＝r_ij/v计算得到，这里r_ij为子单元到破裂起始点的距离，v为破裂速度。从四幅图中可以看出，此断层主要有两个破裂点，分别在长度约20 km处和约60 km处，其中20 km处的震源破裂更加剧烈。断层的初始破裂时间集中在0.02—0.05 s之间，最大持续时长可达22 s，产生的地震矩有7×10²⁵ N·m，最大应力降在8 MPa附近。这四个分布是本文所使用的NNSIM模拟方法需要的震源参数。基于以上震源参数的设定开展多破裂模式下强地面运动模拟。

3.   历史烈度的对比

三河—平谷地震有详细的史料记载，可据此进行烈度评定。本文进行的历史地震烈度对比采用的是中国地震烈度标度，烈度级别共有12个，分别用罗马数字Ⅰ到Ⅻ表示。Bakun和Wentworth （1999）提出了一种直接利用烈度数据点来确定历史地震震级和震中的方法，并且得到了一些应用（Parsons et al，2000；Hinzen，Oemisch，2001；Bakun et al，2002，2003；Quadros et al，2019）。图2蓝线为近震源历史地震烈度分布曲线，棕红线为夏垫断裂带现今地震活动的地震密度，红点为三河—平谷地震的宏观震中，其位于夏垫潘各庄一带（孟宪梁等，1983）。地震活动密集分布空间展布显示夏垫断裂走向约为N50°E，从密集程度可以看出，该断裂现今地震活动仍较频繁。“地震密集值”是王健（2001）提出的参量，其全称为“震中分布分震级网格点密集值”。对于给定的震级和时间段，其能够合理地反映出地震震中的相对密集和稀疏所形成的图案，从而体现出地震的数量和震中位置的空间集中程度，在划定的研究区域内，计算出所有网格点的密集值并绘制等值线，就可以得到给定时间段内特定震级的地震活动性图像。具体计算方法在其研究（王健，2001）中有详细介绍。

图  2  1679年三河—平谷地震等震线与地震密集分布图

Figure  2.  Isoseismic lines and earthquake density distribution map of the 1679 Sanhe-Pinggu earthquake

下载: 全尺寸图片 幻灯片

为了确定三河—平谷地震的震级，首先我们把历史地震烈度曲线Ⅸ，Ⅹ，Ⅺ依据中国地震烈度表GB/T 17742−2020 （中国地震局，2020）转化成PGA和PGV峰值，即Ⅸ对应402 cm/s²，Ⅹ对应831 cm/s²，Ⅺ对应1730 cm/s²，每条曲线取约50个点。Somerville等（1999）认为最大凹凸体中心在宏观震中之外，为此我们限定最大凹凸体在宏观震中的西南。由震源设定部分的计算得出满足该限制条件的设定破裂模式共180个，利用NNSIM方法（Zhou，Chang，2019）计算比较参考场点的加速度、速度以及PGA和PGV。模拟及转化PGA和PGV的差为：

式中，si为模拟数据，hi为从历史事件中获得的数据。由式（5）可以看出，R和S同样反映参考场点的加速度、速度的PGA和PGV与历史烈度的差异。图3表示M_W7.5到M_W8.0五个震级下，每个震级的180种破裂模式对应的R值，其中R值越小，代表其模拟的PGA，PGV与历史烈度差异越小。图3显示M_W7.8模拟PGA，PGV与历史烈度相差最小，其次是M_W7.9，而M_W8.0相差较大。其中R最小值在M_W7.8的模式26。

图  3  R随震级和破裂模式的变化

Figure  3.  Variation of R with seismic magnitude and rupture mode

下载: 全尺寸图片 幻灯片

S与R值仅表明参考场点上模拟结果与转化PGA和PGV总的差异，并不能全面反映设定地震的空间分布特征，哪个R值对应的震源激发的地震动空间分布与历史烈度一致或相近，仍需要对M_W7.8，M_W7.9和M_W8.0模拟地震动进行空间展布并作详细对比分析。

针对M_W7.8，M_W7.9和M_W8.0各自对应的80种破裂方式作强地面运动模拟，分别绘制Ⅸ， Ⅹ和Ⅺ烈度曲线，分析发现其中四种破裂模式（图4）产生的烈度空间分布形态与史料烈度形态接近（图5）：图5a绿线为M_W7.8地震四种破裂模式下的模拟烈度，其中Ⅸ， Ⅹ烈度曲线与历史烈度曲线最接近，但是M_W7.8地震四种破裂模式几乎不产生Ⅺ烈度，这说明M_W7.8地震在集震区不足以激发PGA大于1730 cm/s²的地震动；而M_W7.9地震Ⅸ， Ⅹ烈度范围超过历史烈度圈相应范围，Ⅺ烈度圈的空间范围接近历史烈度圈的大小，但中心偏南（图5b）；而M_W8.0地震激发的地震烈度的范围都远大于历史烈度圈，Ⅺ， Ⅹ烈度与历史Ⅹ，Ⅸ烈度大小相似，Ⅸ烈度远大于历史Ⅸ烈度的范围（图5c），这点说明M_W8.0地震产生的破裂远强于历史资料的记载。图5的对比结果显示，历史上的三河—平谷地震应该小于M_W8.0，大于M_W7.8，接近M_W7.9，应该是一个M_W7.9−或M_W7.8＋ ＋的一次强巨震。

图  4  模式17，21，25和29四种破裂模式下的地震矩分布（左）和应力降分布（右）

Figure  4.  Distribution of seismic moment （left） and stress drop （right） in four rupture modes for modes 17，21，25 and 29

下载: 全尺寸图片 幻灯片

图  5  M_W7.8 （a），M_W7.9 （b）和M_W8.0 （c）地震在四种破裂方式下地震烈度模拟结果

绿色等值线上的数字表示PGA，单位cm/s²

Figure  5.  Simulated earthquake intensity under four rupture modes for M_W7.8 （a），M_W7.9 （b） and M_W8.0 （c） earthquake

The numbers on green contour lines stand for PGA with unit cm/s²

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文根据多震源设定方式建立起基于地震距、应力降、断层破裂面每个破裂子元的初始破裂时间和破裂初始持时的震源破裂模型，对1679年三河—平谷大地震进行强地面运动模拟，分析其地震动破裂特征。其次，为了确定其具体震级，我们把历史地震烈度曲线依据中国地震烈度表进行转化，得到该地区烈度曲线对应的PGA和PGV峰值。利用NNSIM方法对参考场点的相关参数进行计算比较，确定M_W7.8，M_W7.9和M_W8.0地震对应的参考场点的模拟结果与转化的PGA，PGV差异较小。基于以上结果，再用多震级多破裂方式对三个震级进行地震动模拟，并作详细对比分析，认为三河—平谷地震震级达不到M_W8.0，但大于M_W7.8，是一次M_W7.9−或M_W7.8＋ ＋的强巨震。

本次研究有以下几个特点：

1） 多震级、多破裂方式的地震动模拟，可提供震区各场点的地震动变化范围，较一次地震的模拟结果更可信；

2） 通过对比三个震级模拟产生的地震动烈度分布，并与历史烈度对比，确定M_W7.9与其一致性最高，也证明模拟中所使用的模拟方法是合适的。

研究中建立的三河—平谷地震动模拟，主要特征有以下几点：

1） 三河—平谷地震时主要有两个破裂点，并且呈现一大一小的趋势；

2） 该地震破裂面破裂持时在15—25 s内，说明破裂迅速，符合大地震的特征；

3） 主要震源破裂点在长度35 km，宽度15 km的研究区域附近，位于（117°E，40°N），此处的受灾也最为严重。

以上分析可以看出三河—平谷大地震的震级达不到M_W8.0，应为M_W7.9−或M_W7.8＋ ＋，这可为后续三河—平谷相关的研究提供地震动强度方面的参考。

审稿专家在本文撰写过程中提出了宝贵的意见和建议，作者在此表示感谢。