引言

夏垫断裂位于北京东南偏南的廊坊市附近，是该地区的主要活动断层，具有发生强震的危险，而这些地震会对北京市及其周边的人口稠密地区构成威胁（黄秀铭等，1991；焦青，邱泽华，2006；Zhang，Wang，2017）。已有研究表明夏垫断裂为大厂隐伏盆地的边界控制性断裂，是由从古近纪伸展变形到新近纪以来压剪构造环境下的应变动态变化所致，现在仍处于活动状态（冉勇康等，1997；徐锡伟等，2005；吴晶等，2007）。该断裂走向北东—北东东，呈上陡下缓的趋势，构造性质属于正断层，从深部到浅部倾角逐渐减小（江娃利等，2000；张先康等，2002；赵金仁等，2004；何付兵等，2013；雍凡等，2014；李赫等，2020））。夏垫断裂南北向长约120 km，向西南方向延伸至文安，可能与霸县—束鹿—邯郸断裂带发生相互作用，这些断裂与夏垫断裂一起构成夏垫断裂带（杨晓平等，2012；冉志杰等，2013；李赫等，2020）。

在夏垫断裂上开展的钻井和探槽开挖研究结果（张晚霞等，1995）显示夏垫断裂呈东南盘下降、西北盘抬升的正断倾滑性质，兼有右旋走滑分量。高景华等（2007）通过浅层地震勘探成功地探测到夏垫断裂的南段及其上断点，为隐伏夏垫断裂的地质构造、区域活动性特点等的进一步研究提供了可靠的依据。长期研究显示夏垫断裂两侧的灾害状况差异极大，产生的差异源于地层结构、沉积物高度和反射波能量等方面（高振寰，1979；Sun et al，1998；赵金仁等，2004；刘保金等，2009；毛昌伟等，2010；杨晓平等，2012），因此，从纵向和横向上来看震源区均处于地壳结构差异明显的接触带上，该接触带进而也可能成为孕育大地震的有利区域（Li et al，2012；Yu et al，2019）。同样，力学研究结果也表明，夏垫断裂是未来可能发生强震的关键之地（Stein et al，1997；Harris，2000；沈正康等，2004；Robinson，Zhou，2005；Liu et al，2007；Karakostas et al，2013）。

关于夏垫断裂未来活动性方面的研究主要有：高景华等（2008）经由探测试验得到夏垫断裂在1679年三河—平谷地震震源西南方向约3 km处的活动性减弱；肖修来等（2017）基于G-R关系利用最小二乘法对华北地区和三河—平谷地区的地震复发周期进行时间扫描计算，计算结果表明两地区的地震活动性自1976年唐山地震以来趋于稳定；余中元等（2020）采用时间相依的布朗模型、随机特征滑动模型和通用模型对断层破裂源上特征型地震活动的时间相依特征开展研究，其结果也表明夏垫断裂在未来三十年内发生强震活动的可能性较低。虽然强震发生的可能性较低意味着夏垫断裂的活动性较弱，但强震复发间隔和同震位移均有很大的不确定性，所以夏垫断裂的地震危险性仍然有必要进一步分析。

为了评估夏垫断裂潜在的地震危险性，诸多科研人员针对夏垫断裂开展了多方面的研究，包括其复杂的断裂行为、对大地震分布的控制以及与潜在地震风险评估的内在关系，还有该断裂的构造演化历史、地壳上地幔速度结构、地球化学特征以及断裂形成和变形机制（高清武，李霓，1998；刘博研等，2007；刘保金等，2009；雍凡等，2014；张培震等，2014；Wang et al，2014）。与此同时，夏垫断裂活动构造的几何结构特征和运动学、其最新活动期、全新世的古地震活动、晚第四纪的构造变形和分段以及强震模拟和地震危险性预测方面的研究也在展开（冉勇康等，1997；张先康等，2002；邓起东等，2003；徐锡伟等，2005；张培震等，2014；Wang et al，2014）。此外，前人还通过二十世纪九十年代分布在夏垫断裂带上的东六合屯村、潘各庄村和齐辛庄村的沟渠进行了地震危险性分析参数的研究（冉勇康等，1997；江娃利等，2000）。这些研究对于夏垫断裂所在地区的地震危险性分析和未来地震活动预测都具有重要意义。

本文拟在上述研究的基础上，利用多震源、多破裂方式对夏垫断裂开展地震动模拟，分析通州城区和北京中心城区可能发生的地震动灾害，为首都圈地震危险性分析提供理论支撑。

1.   方法

随机有限断层法广泛应用于地震动预测模型的构建中，其基本思想是将一个断层分成N个可视为点源的子断层，最终按照一定的时间延迟将各子断层的地震动模拟结果叠加起来即可得到整个大断层的模拟地震动效应（Beresnev，Atkinson，1997）。针对传统的随机有限断层法的模拟结果仅在高频部分可靠，Zhou和Chang （2019）提出了NNSIM （non-uniform stress and non-uniform window function simulation method）随机有限断层法，该方法与传统的随机有限断层法的不同主要在于：一是用不同源时间函数替换每个点源特定的时间函数，从而提高了低频模拟部分结果的可靠性；二是加入变化的应力降，以减少单个应力降在模拟过程中所产生的影响，使最终的模拟结果更加接近于实际情况。

NNSIM随机有限断层法沿用传统随机有限断层法的叠加思想，先把整个断层分为N个子断层，并将每个子断层看作一个点源，按照特定破裂模式和速度下的时间顺序，分别叠加子断层的模拟结果，从而得到有限断层破裂的的地震动效应。

式中：n_l和n_w分别表示断层在长和宽两个方向上划分的子断层个数，∆t_ij为第ij个子断层点源到观测点的相对延时，a_ij （t）是在第$ ij $个子断层上用随机点源法计算得到的地震加速度，a（t）为大断层在观测点的模拟加速度。计算第ij个子源产生的地震动a_ij （t）时，需要使用非均匀拐角频率（Zhou，Chang，2019），第ij个子源的角频率为

式中，∆σ_ij和M_0ij分别表示第ij个子断层的应力降和地震矩，其中下标ij表示各子源的地震矩和应力降不相同。Δσ_ij由Andrews （1980）的位错与应力降之间的数学关系式计算得到。NNSIM方法的详情请参阅之前相关研究（周红，2018；Zhou，Chang，2019；周红等，2021）。

2.   多震源设定

我们在前人研究（高孟潭等，2002；潘波等，2006，2009；刘博研等，2007；付长华等，2012；巴振宁等，2022）的基础上，将此次研究区域选取为夏垫断裂中间部分，该部分总长度约为80 km，宽度约为25 km。选取M_W7.9，M_W7.8和M_W7.5等三个震级，每个震级均分别按照多震源、多破裂的方式进行模拟。Somerville等（1999）指出震级大于M7.5时进行地震动模拟应当设置两个凹凸体，本文依此均设置两个凹凸体进行模拟。各震级的破裂模式按照Somerville等（1999）的震源设定方式来确定。具体设定方法如下：① 凹凸体为两个；② 两个凹凸体面积占总破裂面积的22%；③ 凹凸体破裂地震矩占总破裂地震矩的44%；④ 宏观震中处于主凹凸体之外；⑤ 主凹凸体面积占总破裂面积的17.5%。按照上面的震源设定方式，将夏垫断层均分成80个子断层，每个断层为5 km×5 km。按照凹凸体面积比例，得出主凹凸体的大小为350 km²，对应14个子断层，次凹凸体大小为90 km²，对应4个子断层，凹凸体位错为平均位错的两倍，其余62个子断层为背景位错。设置主凹凸体与次凹凸体在80个子断层上充分排列组合，且主、次凹凸体所对应的子断层不能重合，也不位于最边缘，按照这种分布方式，用排列组合的方法计算得到满足条件的破裂方式共有180种。本研究中断层模型兼顾走滑与倾滑的性质（向宏发等，1988；张先康等，2002），断层中每个凹凸体的埋深受倾角的影响不尽相同，设置顶点的埋深为1 km。

3.   夏垫断裂M_W7.8和M_W7.9地震的强地震动分析

夏垫断裂上发生于1679年的三河—平谷地震应该是弱M_W7.9或强M_W7.8的一次巨震。因此我们有必要评估如果夏垫断裂上再次发生这样的巨震会引起怎样的强地面运动，又会导致何等破坏。

如何评估夏垫断裂上未来可能发生的M_W7.9和M_W7.8地震的地震动？由于任何一次特定破裂模式的震源所激发的强地面运动不足以代表某一震级地震产生的地震动，而仅能代表特定破裂模式本身产生的地震动，因此对于某个震级引起的地震动变化，需要引用多破裂模式总体讨论该震级地震所产生的地震动变化。利用Zhou和Chang （2019）提出的NNSIM随机有限断层法，按照Somerville等（1999）的震源设定方式，最大凹凸体中心位于宏观震中之外，我们限定最大凹凸体在宏观震中的西南，这样对于夏垫断裂上可能的M_W7.9和M_W7.8复发地震的地震动，满足该限制条件的设定破裂模式共180种。下面我们将讨论这180种破裂模式的总体效果。

对于地表每个场点，分别模拟M_W7.9和M_W7.8地震的每个破裂模式产生的地震动。根据地震矩、应力降等参数的计算可以得到每个场点对应的振幅谱，基于该振幅谱的反傅里叶变换得到场点的地震动时程，然而由于所需相位谱是随机产生的，因此每个场点的地震动时程都需要经过多次计算才能得到比较可信的峰值。考虑到计算量，我们选取5次计算。这样，同一个震级下每个场点得到180×5个模拟结果，即900个加速度和速度记录。取每条加速度峰值（PGA）和速度峰值（PGV）或者反应谱为研究对象，讨论其空间分布特征。这里以PGA为例，每个场点900个PGA，分别取这900个PGA值的50%，60%，70%和80%分布概率所对应的PGA值即为某一震级50%，60%，70%和80%分位数的PGA值，其中50%分位数的PGA代表该场点的地震动相对较小，80%分位数的PGA代表该场点的地震动比较高，该震级大部分地震激发的地震动小于该值。

图1为M_W7.8震级的四个分位数所对应的PGA分布，402 cm/s²等值线对应于Ⅸ烈度，分数位从50%增至60%，70%和80%时，Ⅸ烈度范围逐渐加大，相对应的面积分别为6.92万km²，7.23万km²，7.67万km²和11.6万km²。从图1c和1d可见：作为北京副中心的通州城区，其PGA接近831 cm/s²，对应于Ⅹ烈度圈；北京处于195 cm/s²线上，靠近Ⅺ烈度圈。图2为对应M_W7.9地震50%，60%，70%和80%分位数所对应的PGA空间分布，可见：对应的Ⅸ烈度圈内的面积从10.3万km²依次增加至10.7万km²，11.2万km²，最后增至11.9万km²；Ⅺ烈度圈（图2中内红线）也有明显的增加；通州城区的PGA达到了1 000 cm/s²，在Ⅹ烈度圈里。

图  1  M_W7.8地震多破裂模式下50%，60%，70%和80%分位数的地震动峰值（PGA）分布

Figure  1.  Distribution of 50%，60%，70% and 80% quantile PGAs at M_W7.8 by multiple rupture patterns

下载: 全尺寸图片 幻灯片

图  2  M_W7.9地震多破裂模式下50%，60%，70%和80% 分位数的PGA分布

Figure  2.  Distribution of 50%，60%，70% and 80% quantile PGAs at M_W7.9 by multiple rupture patterns

下载: 全尺寸图片 幻灯片

图  2  M_W7.9地震多破裂模式下50%，60%，70%和80%分位数的PGA分布

Figure  2.  Distribution of 50%，60%，70% and 80% quantile PGAs at M_W7.9 by multiple rupture patterns

下载: 全尺寸图片 幻灯片

4.   夏垫断裂M_W7.5地震多破裂模式下的强地震动场分析

4.1   夏垫断裂多断层多破裂模式下的地震动分布

夏垫断裂真实的走向至今还不确定，多数研究成果显示夏垫断裂的总体走向处于30°—60°之间。为了使本文的讨论结果更具有代表性和普遍性，我们在多破裂模式下总体讨论M_W7.5地震产生的地震动变化，再在模拟中加入走向因素，针对走向为30°，45°，55°的三条断层，研究其强地面运动的分布特征。对于每组走向断层，断层长度仍然设定为破裂长度为80 km，宽度为25 km，子断层大小为5 km×5 km。对于多震源设定部分，我们取消凹凸体位置的限制，再按照排列组合的计算方法，得出每条断层发生破裂的可能方式共有294种。

三个走向下的断层在294种破裂方式下模拟的地震动也按50%，60%，70%和80%分位数筛选，对应的地震动分布如图3—5所示。从图中显然可见地震动分布呈现如下特征：① 50%，60%，70%，80%分位数地震动均环绕断层呈近似椭圆分布；② 80%分位数地震动对应的Ⅸ度烈度圈较50%分位数地震动对应的Ⅺ度烈度圈明显增大，即分位数越大，地震动所对应的Ⅸ度烈度圈的范围越大；③ 对于通州城区而言，在断层走向为30°，45°和55°三种情况下，30°走向断层距离Ⅸ度烈度圈最近；④ 同样分位数PGA在三种断层走向下，30°走向在北京中心城区处产生的地震动峰值最大；⑤ 不管哪个走向，均在70%和80%分位数时才会产生较为明显的大于750 cm/s²的加速度，这说明分位数增高时在同一走向条件下，PGA值会随之增大。

图  3  断层走向为30°时50%，60%，70%和80%分位数下的地震动PGA分布

Figure  3.  PGA distribution of ground motion at 50%，60%，70% and 80% quartile due to the fault with strike 30°

下载: 全尺寸图片 幻灯片

图  4  断层走向为45°时50%，60%，70%和80%分位数下的地震动PGA分布

Figure  4.  PGA distribution of ground motion at 50%，60% ，70% and 80% quartile due to the fault with strike 45°

下载: 全尺寸图片 幻灯片

图  5  断层走向为55°时50%，60%，70%和80%分位数下的地震动PGA分布

Figure  5.  PGA distribution of ground motion at 50%，60%，70%和80% quartile due to the fault with strike 55°

下载: 全尺寸图片 幻灯片

将三个走向断层的加速度合在一起进行综合评价，仍采用分位数方法，得到图6所示的50%，60%，70%和80%分位数的PGA分布图。相较于图3—5，变化主要表现在：① 等值线变宽，如402 cm/s²对应的Ⅸ烈度圈等值线有明显的加宽；② 由于断层南端不动，而北端发生移动，因此等值线形状由围绕断层较均匀的椭圆变成一头细一头粗的梨形，即等值线形状由于不同走向断层的加入而发生了变化，且变得更为复杂；③ 合成后的等值线值与走向为30°时的结果相差明显，说明30°走向时断层对地震动的影响占比较大；④ 此时通州城区位于402 cm/s² PGA等值线附近，即靠近Ⅸ烈度圈。

图  6  三个走向的加速度合成后50%，60%，70%和80%分位数的下地震动PGA分布

Figure  6.  Distribution of PGAs of the synthesis records at 50%，60%，70% and 80% quartiles for the faults with three strikes

下载: 全尺寸图片 幻灯片

4.2   夏垫断裂M_W7.5仿真地震破裂模式下地震动分布

Wang等^①针对夏垫断裂利用滑动摩擦准则模拟了M_W7.5地震的仿真破裂过程，本文将其计算结果作为发震震源，依然采用NNSIM随机有限断层法，按照第二节中的震源设定方式开展了M_W7.5的地表地震动模拟。图7为模拟的加速度峰值及其反应谱分布，图8为速度峰值及其对应的反应谱分布，图中断层（黄色线框）为非规则曲面，其几何形态更接近实际断层面。对于加速度（图7），我们选择长周期（T＝3 s）、中周期（T＝1 s）和短周期（T＝0.2 s）的反应谱作为代表以呈现加速度反应谱（pseudo-acceleration response spectrum，即 PSA）的空间分布，由图7显然可见：① 随着周期变短，PSA最大值减小，短周期（T＝0.2 s）的PSA最大值超过了1 000 cm/s²，长周期（T＝3 s）的PSA最大值小于100 cm/s²；② 在距离断层垂直距离20 km范围内，短周期（T＝0.2 s）的PSA均超过了1 000 cm/s²，距离断层垂直距离50 km范围内短周期（T＝0.2 s）的PSA均超过了400 cm/s²。图7d的PGA分布显示，仿真破裂产生了两个峰值区域（最内红圈），PGA超过750 cm/s²，Ⅸ烈度圈（第二红圈，对应402 cm/s²）的形态随内红圈变换，呈现为两头大中间小的形态；通州在Ⅸ烈度圈内，PGA达到500 cm/s²；北京中心的PGA达到160 cm/s²。

图  7  模拟M_W7.5加速度峰值及其对应的不同周期T的反应谱分布

Figure  7.  The distribution of simulated acceleration peaks and response spectra with different period T for M_W7.5

下载: 全尺寸图片 幻灯片

图  8  模拟M_W7.5速度峰值及其对应的不同周期T的反应谱分布

Figure  8.  Distribution of the simulated velocity peaks and response spectra with different period T for M_W7.5

下载: 全尺寸图片 幻灯片

图8为模拟的M_W7.5速度峰值及其对应的速度反应谱（pseudo-velocity response spectrum，即PSV）分布，可见：长周期（T＝3 s）、中周期（T＝1 s）和短周期（T＝0.2 s）的速度反应谱PSV大致为围绕断层变化的椭圆，断层南端PSV变化大于断层北端，表现为等值线变化剧烈；中周期（T＝1 s）的PSV最大值大于长周期（T＝3 s）和短周期（T＝0.2 s）的PSV最大值，其最大值大于100 cm/s；PGV的空间变化最为复杂，有明显的上盘效应，表现为PGV等值线向渤海方向延伸；通州的PGV约为25 cm/s，北京中心处的PGV小于10 cm/s。图7和图8显示对于模拟的M_W7.5地震，其地面地震动分布复杂，Ⅸ烈度圈（图7第二红线圈）大致距离断层40—50 km。

5.   讨论与结论

夏垫断裂是影响北京地区及首都圈地震安全的最重要的断裂带，1679年曾经发生过M_W8.0巨震，对北京地区产生了极为严重的破坏。本文通过全破裂模式设定震源，模拟了夏垫断裂发生强震在首都圈产生的地面运动及分布特征，可以为首都圈及其附近区域的地震预测及防震减灾工作提供参考，主要结论如下：

1） 多震级、多破裂模式的地震动模拟提供了震区各场点的地震动变化范围，较一次地震的模拟结果更可信，其模拟结果显示：夏垫断裂发生M_W≥7.5地震时，地面运动Ⅸ烈度圈的范围较大，首都圈可能遭受严重的破坏，而通州城区位于Ⅸ烈度圈内，受破坏程度可能更为严重；夏垫断裂发生M_W＜7.9地震时，北京中心城区位于Ⅷ烈度圈外，发生M_W7.9时，其位于Ⅷ烈度圈内，遭受严重破坏的可能性明显加大；

2） 基于摩擦准则提供的仿真断层破裂过程震源，可以提供更复杂的震源变化。基于仿真震源与设定震源的模拟结果量级一致，两结果可以相互印证、互为补充。