0 引言

地震灾害的根源主要是断层突然运动产生地表位移，激发地震波造成地表强烈振动，即强地面运动。强地面运动可直接导致建、构筑物的倒塌、破坏，也可引起大量的崩塌、滑坡、滚石、地裂缝、砂土液化等，从而间接导致基础设施等遭受严重破坏。强地面运动研究的基础是强震记录数据，但是由于地震的随机性和不确定性，强震记录往往集中在某一些地区，地理环境和地质条件的不同使得许多数据的应用价值受到限制。为了在强震数据匮乏地区获得满足特定地理地质条件的地震数据，需要开展强地面运动模拟研究工作。

广西是我国少震、弱震地区，其所在的华南地块是中国大陆一个相对稳定的构造单元，地震活动性较弱(黄强强等，2022)。虽然广西地震活动性水平较低，但是在历史上却发生过华南内陆地区自有地震记载以来最大的地震——1936年灵山M6¾地震。1936年4月1日9时31分，广西灵山县平山圩东罗阳山发生M6¾强烈地震，震中位于灵山县东北约20km的罗阳山西北麓一带，震中烈度达Ⅸ度或Ⅸ度强(李冰溯等，2018)。该地震造成101人死亡，263人受伤，8000余间房屋倒塌破坏(陈恩民等，1984；李伟琦，1989；任镇寰等，1996)。对于缺乏仪器数据的历史地震，进行强地面运动模拟是评估该地区潜在地震危害性的重要手段。作为广西乃至华南内陆有记录以来最大的历史地震，开展灵山M6¾地震的强地面运动模拟，进而补充本地区强震记录数据库，具有重要的科学价值和工程意义。但是由于历史久远和仪器记录匮乏等原因，自发震至21世纪初70余年的时间里，灵山地震的研究程度一直较低(李细光等，2017b)。近年来，随着技术手段进步，相关学者对灵山地震开展了一系列研究，获得了一批突破性成果，已基本厘清此次地震的震级、震中位置、极震区等级范围和发震构造等关键信息(周本刚等，2008；何军等，2012；李细光等，2017a、2017b；李冰溯等，2018；李细光等，2018)。因此，针对灵山地震开展强地面运动模拟的条件已经成熟。

本研究使用随机振动有限断层法对1936年广西灵山M6¾地震进行了强地面运动模拟，在收集地震地质资料的基础上，确定发震断层的走向、倾角、长度、宽度、震级、应力降等参数，选取适当的传播路径参数和场地效应参数，将上述参数输入随机有限断层模型进行强地面运动模拟，并将模拟结果与地震烈度野外调查数据和地震动衰减关系进行对比，验证模拟结果的可靠性。本研究最终获得1936年广西灵山M6¾地震区域及极震区峰值加速度(PGA)、峰值速度(PGV)及烈度的空间分布，相关结果可为工程抗震设计、政府决策和地震灾害风险评估提供依据。

1 方法与原理 1.1 随机振动有限断层方法

随机振动法是由Hanks等(1981)和Boore(1983)于20世纪80年代提出的一种用随机振动来模拟高频地震波的方法。这种方法的本质是把远场地震动看作有限带宽、有限持时的高斯白噪声，并把影响地震动参数的一系列因素(包括震源、路径、场地)在频率域里表达为简单的乘积形式，然后通过窗函数调整噪声的傅里叶谱，最后通过傅里叶反变换得到地震动参数的时程曲线(Boore，1983、2003)。对于一次大地震，合成地震动的有效方法是基于小震模拟。在随机有限断层模型中，断层破裂面被分成N个大小相等的矩形子源，每个子源即为一个点源，每个子源满足ω平方衰减规律，通过设定震源模式和破裂传播速度，可以得到子源破裂的时间顺序，然后根据子源与场地的几何关系，计算每个子源对场地的影响。所有子源在观测点引起的地震动在时域中以适当的延迟时间叠加，便可获得场地的地震动时程a(t)，即

$ a(t)=\sum\limits_{i=1}^{N_L} \sum\limits_{j=1}^{N_W} a_{i j}\left(t+\Delta t_{i j}\right) $

(1)

其中，N_L、N_W分别为断层沿长度和宽度划分的子断层数，i和j分别代表沿断层走向和倾向方向的断层指标，Δt_ij为子源到场地的滞后时间，a_ij(t)为子源利用随机点源方法得到的在观测点的地震动时程。针对有限断层模型采用部分子断层多次破裂的方式来保证地震矩守恒的不合理之处，Motazedian等(2005)提出了“动力学拐角频率”的概念，认为子断层拐角频率随着破裂的传播而减小。每个子断层动力学拐角频率表达式为

$f_{c i j}(t)=4.9 \times 10^6 \beta\left(\Delta \sigma / M_{\text {0ave }}\right)^{1 / 3} \times N_R(t)^{-1 / 3} \times S $

(2)

其中，f_cij(t)为第ij个子断层的动力学拐角频率，t为第ij个子源被触发的时刻，N_R(t)为在时刻t已破裂子断层的累计数，M_0ave=M₀/N为子断层的平均地震矩，S为表示子断层辐射强度的一个常数。为了使断层总辐射能守恒，引入一个标度因子H。动力学拐角频率对子断层总辐射能的递减影响得到标度因子H的补偿，总辐射能保持不变，这样在保证辐射能守恒的前提下避免了子断层多次破裂这一非物理现象。至此，随机振动有限断层方法的理论框架得到最终完善。

1.2 场地条件对地震动模拟结果的影响

场地的浅层速度结构是引起地震场地效应的重要原因，直接影响地震灾害的严重程度。研究表明，基岩上的地震动参数值与考虑场地土层放大效应的参数值差别可达2.5倍(Borcherdt，1994)，因此在强地面运动模拟中需要考虑浅层速度结构的影响。精细的浅层速度结构数据往往依赖于钻孔数据，而由于造价成本等因素，钻孔勘探数量很少，无法满足大规模区域性的强地面运动模拟的需求。为解决这一问题，Wald等(2007)提出一种利用高精度地形数据计算地形坡度来求取浅层剪切波速V_s³⁰分布，进而得到近似场地校正放大因子的简单方法。这一方法已被成功应用于美国地质调查局(USGS)所属的地震烈度速报系统(ShakeMap)，并在实践中得到有效性的验证。V_s³⁰的全球数据可以在USGS官方网站下载，不同学者则提出不同的计算场地放大因子的方法。Borcherdt(1994)计算了不同场地类型条件下短周期(0.1~0.5s)及中周期(0.4~2.0s)4个基岩PGA和PGV输入档下的放大系数，这一结果无论在场地类型划分还是PGA或PGV分档方面均存在划分过粗的缺点。Boore等(2007)在建立地震动衰减关系NGA(Next Generation Attenuation)时认为，V_s³⁰对基岩PGA和PGV的放大作用可表达为

$ \left\{\begin{array}{l} \ln \left(\mathrm{PGA}^{\prime}\right)=\ln (\mathrm{PGA})-0.36 \ln \left(V_{\mathrm{s}}^{30} / 760\right) \\ \ln \left(\mathrm{PGV}^{\prime}\right)=\ln (\mathrm{PGV})-0.6 \ln \left(V_{\mathrm{s}}^{30} / 760\right) \end{array}\right. $

(3)

式中，PGA′和PGV′为考虑了场地放大效应后的强地面运动参数。式(3)避免了场地划分和地震动参数分档，适用于任何场地类型。通过式(3)可将V_s³⁰的影响加入到强地面运动模拟的结果中，使得模拟结果能够充分反映浅层速度结构的影响。图 1展示了本研究区内的地形及V_s³⁰分布情况，可以看到在平原、盆地等沉积层较厚的区域，V_s³⁰数值较低，而在丘陵地带、山区及基岩裸露地区，V_s³⁰数值较高。

2 模型参数 2.1 发震构造

发震断层的几何形态(包括走向、倾角、长度和宽度等)是随机振动有限断层模型的基本输入参数。灵山M6¾地震所在的华南地区由于历史强震较少，人类改造活动频繁且降水较多，地震遗址保存非常困难，因此很长一段时间内关于此次地震的发震断层存在争议。陈国达(1939)最早研究灵山M6¾地震时，认为灵山地震发震断层为NE走向的灵山—大容山断裂，而茹锦文等(1985)则认为发震断层应归属于NW向的昆仑关断裂南延部分。前人对该地震发震断层的研究大多从极震区震害和地裂缝的空间分布情况进行推测(陈国达，1939；茹锦文等，1985；李伟琦，1989；任镇寰等，1996)，缺乏直接的地质学证据。2015年研究人员首次在灵山断裂北段发现灵山地震2条地表破裂带遗迹(李细光等，2017a、2017b)，地表破裂带沿线开挖的探槽清晰揭露出断层对上覆土层的错动，这为灵山断裂北段晚更新世以来的活动提供了直接证据。探槽显示的最新一次地震事件距今约80年，与灵山地震发震时间吻合，证明了灵山断裂北段即为1936年灵山M6¾地震的发震断裂。另外，李伟琦(1989)和任镇寰等(1996)根据极震区的等震线形状和低烈度区的长轴方向，推测灵山地震可能为ENE向断裂与NNW向断裂共轭破裂的结果，其中ENE向构造为控震构造。李冰溯等(2018)综合极震区震害及地质学研究成果，认为灵山地震主发震构造为NE-ENE走向的灵山断裂，罗阳山南麓的泗州断裂同时发生共轭破裂，造成局部烈度增强。因此，本文在开展模拟时将发震断层设定为2条，即灵山断裂北段和泗州断裂(图 2)，2条断裂的几何参数见表 1。

① https://www.activefault-datacenter.cn/

2.2 断层滑动模型

断层滑动模型即断层面上的滑动量静态分布是进行强地面运动模拟的必要参数，一般通过震后地表形变数据或地震波数据反演获得(申文豪等，2019)。对于灵山M6¾地震这类没有观测数据的历史地震，已很难通过反演方式获得滑动模型。除了反演模型之外，地震学者还根据地震波观测总结出断层滑动模型的空间分布规律，提出了通过正演方式生成断层滑动分布的方法(王海云，2004)。本文采用应用广泛的k平方模型正演获得断层面滑动分布。Gallovič等(2004)提出的k平方模型二维波数谱为

$ D\left(k_x, k_y\right)=\frac{\bar{D} L W}{\sqrt{1+\left(\left(\frac{k_x}{k_{c x}}\right)^2+\left(\frac{k_y}{k_{c y}}\right)^2\right)^2}} \mathrm{e}^{i \mathit{\Phi}\left(k_x, k_y\right)} $

(4)

其中，L和W分别为断层破裂面的长度和宽度，D为平均滑移，Φ(k_x，k_y)为相位谱，k_cx和k_cy分别为走向方向滑移分布和下倾方向滑移分布的空间拐角波数。Causse等(2010)给出了震级M_W与空间拐角波数经验关系公式，即

$ \left\{\begin{array}{l} \lg \left(k_{c x}\right)=1.82-0.5 M_{\mathrm{W}} \\ \lg \left(k_{c y}\right)=1.93-0.5 M_{\mathrm{W}} \end{array}\right. $

(5)

本研究采取以下策略获取断层滑动模型：首先根据断层几何参数(表 1)，由震级与断层破裂长度经验关系(Thingbaijam等，2017)分别得到灵山断裂北段和泗州断裂的发震震级和平均滑动位移D，再由公式(5)得到2个断层滑动模型的空间拐角波数，最后根据公式(4)分别生成灵山断裂北段和泗州断裂的滑动模型(图 3)。

2.3 震源及路径传播参数

根据现有文献及地质资料基本可以确定1936年灵山M6¾地震发震断层的破裂长度、走向、倾角、埋深、震中位置等几何参数(表 1)。除此之外，为描述断层破裂和传播特征，还需要确定地震矩、应力降、kappa值、品质因子Q模型、几何衰减模型等参数。龙政强等(2015)选取广西数字地震台网记录到的2007—2013年广西主要断裂带M_L≥2.5地震事件，采用Atkinson方法反演得到介质品质因子Q=366.3f^0.47。根据Tsang等(2010)，Kappa值κ可以由经验关系式确定，即

$ \kappa=0.145-0.12 \ln \left(V_{u c}\right) $

(6)

其中，V_uc为上地壳4km深度内平均横波速度。根据广西中西部地区背景噪声成像结果，灵山地区V_uc约为3.25km/s(黄强强等，2022)，因此Kappa值计算结果为0.0036。杨晓瑜等(2021)利用中国国家地震台网336个固定地震台站记录的远震波形资料，通过P波接收函数分析估算了中国华南地区的地壳厚度，其结果显示灵山地震震源区地壳厚度约为28km，根据三段式几何衰减模型的定义(Atkinson等，1992)，该地区几何衰减模型可以表示为

$ Z(R)= \begin{cases}R^{-1} &, R \leqslant 42 \mathrm{~km} \\ 42^{-1} &, 42 \mathrm{~km}<R \leqslant 70 \mathrm{~km} \\ 42^{-1}(70 / R)^{0.5} &, R>70 \mathrm{~km}\end{cases} $

(7)

其中，R为震中距。龙政强等(2015)研究显示，广西地区10MPa及其以上的高应力降地震主要分布在桂东南，且历史上发生过5.0级以上地震的震源区及其附近仍具有高应力背景。2010年8月17日发生的灵山M_L3.6地震其应力降为11.54MPa，该地震位于1936年灵山M6¾地震和1958年9月25日灵山M5 ¾地震震中附近，这说明灵山地震区具有高应力背景，本研究将应力降的值取为8.0MPa。

3 结果分析与讨论

灵山M6¾地震虽然缺乏仪器记录，但具有丰富的烈度调查数据。本研究针对宏观区域(调查烈度覆盖区域及其扩展区域)和极震区(Ⅷ度和Ⅸ度区)分别划分不同的计算区域，并设置不同尺度的计算网格，利用随机振动有限断层模型计算每个网格点的PGA和PGV，通过经验转换关系(丁宝荣等，2017)分别转化为烈度，最后取两种烈度的加权平均为最终烈度值，并将其与野外调查结果进行比较。丁宝荣等(2017)给出的PGA、PGV与烈度MMI转换关系为

$ \left\{\begin{array}{l} \mathrm{MMI}=3.73 \lg (\mathrm{PGA})-1.23 \\ \mathrm{MMI}=3.61 \mathrm{lg}(\mathrm{PGV})+2.72 \end{array}\right. $

(8)

3.1 宏观烈度模拟

宏观烈度模拟区域划定在106.0°E~112.0°E、20.5°N~25.5°N范围内，区域网格大小为6.0°×5.0°，经度及纬度方向网格精度分别为0.06°和0.05°。为更好地展示强地面运动参数衰减趋势，绘制了lg(PGA)和lg(PGV)空间分布图，如图 4(a)、4(b)所示。从图中可以看到模拟的强地震动参数具有以下特征：①PGA最大值在900gal左右，PGV最大值达100cm/s，转化为烈度均已超过Ⅸ度；②PGA和PGV均随距断层距离增加而迅速衰减，且PGA的衰减趋势更快；③PGA和PGV均受到场地条件的影响，在平原等沉积层较厚区域存在不同程度的放大效应。

图 4(c)展示了模拟烈度与野外调查烈度(陈恩民等，1984)的对比，图中白色实线即为野外调查烈度等值线。可以看到，模拟获得的地震烈度在空间分布的大致走向、烈度等级范围与野外调查结果基本一致。模拟结果显示断层附近区域极震区的最大烈度已经超过Ⅸ度，其中灵山县城位于Ⅷ度区，南宁市区位于Ⅴ度区，玉林和钦州均位于Ⅵ度区，均与调查结果相同。另外需要注意的是，由于采样点的有限性，野外调查烈度图中的烈度等值线大多是经验性地圈划为圆形或者椭圆形，而本研究的数值模拟结果由于考虑了V_s³⁰的作用，烈度分布的形态明显受到场地条件的影响。其主要体现在一些沉积层较厚的冲积平原和盆地出现高烈度异常区，这些地区的V_s³⁰普遍较小，从而产生了局部场地放大效应，这可能是造成高烈度异常的主要原因。

3.2 极震区模拟

为更好地研究双断层共轭破裂造成的强地面运动特点，本研究针对极震区开展了网格尺度更精细的数值模拟，模拟区域划为109.1°E~109.7°E、22.2°N~22.8°N，区域网格大小为0.6°×0.6°，经度及纬度方向网格精度为0.006°。极震区地震动模拟结果如图 5所示，图中白色实线为野外调查Ⅷ度和Ⅸ度区等值线，黑色实线为模拟烈度Ⅷ度和Ⅸ度区等值线。国家地震局全国地震烈度区划编图组(1979)汇编的《中国地震等烈度线图集》中，灵山M6¾地震的震中烈度最大为Ⅸ度，并且形状为一长轴方向呈NE向的椭圆形。《中国地震目录》(顾功叙，1983)和《中国近代地震目录》(中国地震局震害防御司，1999)均采用这一结果。此后，陈恩民等(1984)、李伟琦(1989)以及任镇寰等(1996)对极震区烈度进行了重新调查和评定，结果均显示极震区的烈度分布有NEE和NWW两个优势方向，呈“T”型分布。由图 5可以看到，无论是PGA、PGV分布，还是烈度分布，本研究的模拟结果较好地反映了极震区烈度分布的这种“T”型特征，尤其是在Ⅸ度区。这一结果从模拟的角度验证了模型参数的合理性，也证实了灵山地震极震区分布是由灵山断裂及泗州断裂共轭破裂造成的。两断层发生共轭破裂的控制因素可能是罗阳山地区最大有效力矩准则约束下的以NWW-SEE向为主压应力方向的力学系统(张沛全等，2012)。图 5(c)还显示出与调查烈度不同的是，模拟烈度Ⅷ度区在北部及西南部有两处明显的烈度异常区，而这两处烈度异常区均为河流冲积平原，北部为郁江平原区，西南部为钦江平原区。这两处平原区沉积层较厚，V_s³⁰水平较低，从而导致明显的场地放大效应。

3.3 衰减关系对比

除了与调查烈度进行空间分布对比之外，本研究还将宏观区域和极震区模拟得到的地震动参数与衰减关系NGA(Boore等，2007)进行对比。图 6(a)、6(b)分别展示了宏观区域模拟PGA、PGV与NGA的衰减趋势，图 6(c)、6(d)分别展示了极震区内模拟PGA、PGV与NGA的衰减趋势，图中黑色实线为NGA衰减关系曲线，黑色虚线为95%置信区间曲线。从图 6可以看到，模拟结果与NGA结果的衰减趋势大体一致，在极震区内模拟结果与NGA曲线符合程度较高，但在震中距大于100km时，地震动模拟结果整体略高于NGA预测值。图 7展示了本研究模拟的烈度与王继等(2008)针对本地区烈度衰减公式的对比，其中蓝色实线为长轴衰减公式，蓝色虚线为短轴衰减公式。从图 7可以看到，模拟烈度与理论公式在不同震中距尺度上具有相似的衰减特征。在震中距小于5km的范围内，基于烈度衰减公式预测的最大烈度值为Ⅷ度强，显然低于模拟烈度和实际调查烈度。造成这一差异原因可能是王继等(2008)所使用的地震事件中超过6.0级的地震样本数较少，从而可能会低估大震级事件极震区内的烈度水平。在震中距大于200km处，模拟烈度则低于基于烈度衰减公式预测的烈度水平。考虑到在远场区域地震动主要以面波形式出现，而品质因子Q对其有着重要影响，在本研究过程中只采用了单一的Q衰减模型。事实上，许多研究已经表明Q值是依赖于震中距变化的(Tsang等，2010)，因此，在今后工作中需要建立随距离变化的Q值模型来改善模拟结果。

4 结论与展望

本研究收集整理了1936年灵山M6¾地震的地质资料、活动断层探测最新成果等，利用随机振动有限断层模型计算了宏观区域及极震区内网格点的峰值加速度及峰值速度等参数，并且加入了浅层横波速度结构对模拟结果的影响，最终得到此次地震的地震动分布，分析了地震动特征。本研究将模拟结果与历史地震烈度调查数据及地震动衰减关系进行对比，结果显示模拟结果与调查烈度值在整体特征、极震区的分布等方面均符合较好。

由于广西地区中强地震少发，缺乏足够的强地面运动记录回归拟合得到本地区地震动衰减关系，而模拟得到的强震数据可以补充强震数据的不足，为建立广西地区地震动衰减关系曲线提供有益补充。另外，由于随机振动有限断层模型未考虑地震波在介质中的传播过程，不需要像复合震源模型、有限差分法或者谱元法等方法耗费大量时间用于地质体网格剖分和格林函数计算，具备时效性强的突出特点，这一点正是地震烈度速报所需要的。因此，未来随机振动有限断层模型可以应用于本地区地震烈度速报，为震后灾情快速评估和应急救援提供有效帮助。