0.   引言

天水地区是我国典型的黄土丘陵区之一，黄土分布广，结构松散，且地势陡峻，不利于区内黄土斜坡体的稳定^[1~2]。该地区又位于青藏高原东北缘，区域内构造活动强烈(见图 1)，其南、北两侧分别发育礼县—罗家堡和西秦岭北缘全新世活动断裂，规模均可长达一、二百公里，历史上曾发生过8级地震或者具备发生8级地震的构造背景^[3~8]。斜坡地质条件和强震活动背景的双重作用下，奠定了天水地区地震黄土滑坡多发的基本态势。

图  1  天水盆地及周边地区地震构造背景

①—六盘山东麓断裂；②—清水河断裂；③—西秦岭北缘断裂；④—礼县-罗家堡断裂；⑤—通渭断裂  Ⅰ—龙门山地震带；Ⅱ—六盘山-祁连山地震带；Ⅲ—鄂尔多斯地震带；Ⅳ—长江中游地震带

Figure  1.  Seismotectonic background in Tianshui basin and its adjacent areas

下载: 全尺寸图片 幻灯片

第五代中国地震动参数区划图给出了天水盆地内基本地震动峰值加速度，并分别规定了罕遇地震动和极罕遇地震动峰值加速度对基本地震动峰值加速度的放大倍数，为单体滑坡机理研究、区域滑坡危险性定量评估提供了基础地震动参数。然而，区划图是综合所有潜在震源贡献的计算结果，与某个特定地震引起的地震动之间存在一定程度的差异。因此，运用数值模拟手段研究特定地震作用下的地震动场，对于单次地震引起的滑坡评估和灾害预防等更具有科学的指导意义。目前，基于大量强震动观测和地震动空间分布特性的研究，一般认为：低频地震动(＜1 Hz)是确定性的，可以基于震源模型和地下速度结构模型，采用确定性数值模拟方法(如有限差分法、有限元法等)得以重现；而高频地震动(＞1 Hz)由于受震源破裂细节和介质小尺度构造复杂性的影响，是随机的，可以采用随机振动合成方法模拟得到。鉴于近年来区域构造活动渐强的背景及天水地区历史大震发育特征及研究区设定地震模拟的空白，以礼县—罗家堡断裂为例，采用有限差分和随机振动合成相结合的复合方法，模拟该活动断裂上发生的Mw7.7级地震在天水盆地产生的地震动场，量化分析设定地震作用下盆地区域内的地震动参数特征，为区域地震滑坡研究提供合理的地震动参数。

1.   计算方法

采用有限差分和随机振动合成相结合的复合方法^[9~10], 模拟设定地震作用下天水盆地内的地表宽频带地震动。一方面采用有限差分法模拟1 Hz以下的长周期地震动时程，另一方面则采用随机振动合成方法模拟1 Hz以上的高频地震动时程，然后将两个地震动时程叠加，即可得到有效的宽频带地震动时程。

有限差分法是一种常用的确定性地震动数值模拟方法，通过联立求解一阶应力—速度方程组和本构关系方程组可得到任意地点的地震动时程。在求解过程中，通常以交错网格布局方式将各节点上的物理量进行定义，以空间四阶精度和时间二阶精度的差分格式对波动方程进行离散，通过应用水平地表应力镜像方案实现含起伏地形的地震波模拟，引入地球介质品质因子(Q值)描述地震波在介质中传播时的粘弹性衰减，采用吸收边界处理方法解决地震波到达计算边界处的反射和透射问题，并将震源地震矩转化成初始应力分量耦合到对应的网格点上^[11~15]。文中采用了在前述研究的基础上开发并逐步完善的曲线网格有限差分数值模拟方法^[16]；该方法对复杂起伏地形采用贴体网格进行剖分处理，对相对软弱介质采用局部空间网格加密的方法进行处理，在同类方法中综合性能较优，兼顾了计算精度与效率，并且采用高斯滤波克服了间断网格长时间模拟不稳定问题。

随机地震动合成高频地震动方法的原理是将平面断层划分为若干子源，每个子源考虑成一独立点源，每个点源对场点的贡献可根据谱的衰减关系计算得到，后经傅立叶逆变换获得该点源在场点产生的地震动时程。所有点源在场点产生的地震动时程进行叠加即为场点总的地震动时程^[17]。文中对该方法进行了改进，其一是突破了平面断层的限制，可以设置曲面形态的断层；其二是对描述震源破裂过程细节的几个参数考虑更充分，如可以设置复杂的子源破裂时间顺序、不同的子源上升时间、不同的子源滑动速率函数以及断层面上不均匀的滑动位错分布等，这使得模拟结果更加科学。

2.   天水盆地地下介质速度结构模型

天水盆地位于南北活动构造带中段北端，西秦岭断裂带和祁连山—六盘山断裂带的交汇处附近，也是青藏高原东北缘与黄土高原接壤地带，从地震区带划分上属于龙门山地震带北段(见图 1)。天水盆地为群山之间一条狭长的河谷盆地，具体范围为105.38°E~105.95°E，34.5°N~34.7°N，面积约为65 km×20 km(见图 1中黑虚线矩形框位置)。区内地势崎岖不平，落差较大，西侧海拔较东侧略高，山顶海拔最高为1900 m，盆地最低处海拔不足1100 m(见图 2)。

图  2  天水盆地区域地形图

Figure  2.  Topography of Tianshui basin

下载: 全尺寸图片 幻灯片

天水盆地第四系比较发育，主要为更新世和全新世堆积物。河谷区主要为砂砾石的冲洪积堆积，河谷两侧山体边坡主要为坡洪积物，沿山坡往上至山顶部以风积黄土为主。基于地脉动观测资料及钻孔资料进行反演，并利用克里金插值算法，得到了盆地内第四系覆盖层的厚度分布(见图 3)。区域内第四系覆盖层厚度介于5~90 m之间，但绝大部分区域厚度小于50 m。在中央几条冲沟的汇集处、西北侧地势低洼处、西南侧山间、以及东南侧较开阔的河谷等几个局部地区，覆盖层厚度较大，可超过50 m，而山体边坡及顶部的覆盖层厚度则相对较小。

图  3  天水盆地第四系覆盖层等厚线图

Figure  3.  Isopach of Quaternary covering layer in Tianshui basin

下载: 全尺寸图片 幻灯片

鉴于设定地震的发震断层规模较大，建立的有限差分计算模型尺度远大于盆地区域，计算模型范围为104°E~107°E，33°N~36°N(见图 1)。模型东西向和南北向均为300 km，竖向约40 km，共划分为5层不同的介质，包括第四系覆盖层、沉积层、上地壳、中地壳和下地壳，除了考虑地表第四系覆盖层的不均匀厚度分布及地形起伏外，其余下伏各层介质均假定为理想平层形态。此外，盆地区域外部因缺少覆盖层厚度的数据而假定为无覆盖层存在。深部介质的纵波速度以及几个平层分界面的埋深参考了甘肃东南部地壳一维速度结构模型^[18]。在设置纵波速度时，对参考的一维速度结构模型参数稍做修改，即未将单层介质设为均一波速而避免了整体上呈现出跳跃式突变的特点，而是考虑为由浅到深波速逐渐变化的趋势。横波速度根据甘肃东南部地壳介质的纵波速度与横波速度的比例关系(约1.70~1.75倍)计算给出^[19]，介质品质因子设定为横波速度的0.1倍，介质密度参考了邻区的研究结果^[20~21]。各层介质的计算参数见表 1。计算模型的网格划分方案为：横向上网格大小为100 m×100 m，以匹配100 m×100 m精度的起伏地形，同时适应第四系较低波速对计算稳定性的要求；纵向上为不均匀网格，近地表低速介质的网格尺寸较小，随着介质波速的增加，网格逐渐增大，这样可以降低计算量从而提高计算效率。

表  1  天水盆地地下介质的计算参数

Table  1.  Calculation parameters of underground media in Tianshui basin

介质层	第四系覆盖层	沉积层	上地壳	中地壳	下地壳
介质密度/g·cm-3	2.00	2.25~2.65	2.80	2.90	3.00
纵波速度/km·s-1	2.0	3.0~5.8	5.9~6.5	6.6~7.0	7.0~7.2
横波速度/km·s-1	1.0	1.5~3.3	3.4~3.8	3.8~4.1	4.1~4.2
介质品质因子	100	150~330	340~380	380~410	410~420
层介质埋深/km	0~0.09	0~6.5	6.5~24	24~38	38~40

下载: 导出CSV 

| 显示表格

随机合成方法模拟高频地震动时不要求精确的地下介质速度结构模型，只需输入地壳平均密度、平均剪切波速、粘弹性衰减因子三个计算参数。平均密度和平均剪切波速参考全球平均速度结构模型参数，计算中分别取2.8 g/cm³和3.7 km/s，粘弹性衰减因子(Q值)则利用甘肃东南地区近年来积累的大量M_L≥2.5数字地震观测记录反演得到^[22]，其随频率(f)变化的关系式为：

3.   震源破裂模型

距离天水盆地区最近且影响最大的主要有两条大型活动断裂(见图 1)，北西西向的西秦岭北缘断裂穿过盆地区域北侧，晚第四纪以来表现为强烈的左旋走滑和逆冲特征，易发生强度大的地震^[4~6]；北东东向的礼县—罗家堡断裂斜插入盆地区域的东南角，晚第四纪以来亦伴有强烈的地震活动发生^[7~8]。有史料记载以来，天水盆地及周边地区发生过天水西7级地震(734年)、天水南8级地震(1654年)、通渭7 1/2级地震等7级以上大震(1718年)，7级以下中强地震发生频繁。其中，天水南8级地震(1654年)发生于礼县—罗家堡断裂之上，天水西7级地震(734年)则发生于西秦岭北缘断裂之上^{[3, 7~8]}。

以礼县—罗家堡断裂作为发震断裂，探讨发生于该断裂上的大地震在天水盆地及周边地区产生的地震动场，可为该区黄土地震滑坡分析提供地震动参数。发震断裂的位置见图 1，设定地震的矩震级为Mw7.7级，标量地震矩为3.9811e+20牛顿·米。断层面长度为150 km，宽度根据矩震级与断层面面积的经验关系计算得出^[23]，约为30 km，断层走向约为NE60°，倾角84°，滑动角-160°。

震源模型的设置采用混合震源模型的方法，即在平均凹凸体模型的基础上叠加随机滑动分量以考虑滑动分布的不均匀性。关于平均凹凸体模型，Somerville在针对多个地震事件的震源破裂过程反演结果进行统计后，对断层面上滑动量分布特征(凹凸体分布)作了开创性的描述，他认为断层面上的平均凹凸体数目为2个，且凹凸体上的滑动量为背景滑动区的3倍^[24]；随后我国学者选择更多有代表性的地震样本进行统计分析，得到了断层面上凹凸体数量、大小、位置、以及初始破裂点位置的定量关系，对Somerville的研究结果做出了有价值的补充^[25~26]。然而，平均凹凸体模型的不足之处在于对断层面上滑动量分布设置得过于均匀，导致模拟的地震动仅包含某个频段的长周期成分，而添加随机滑动以后，则可以弥补这一缺憾，使模拟的结果更加合理。随机滑动分布遵循K方滑动模型，当波数大于空间拐角波数时，滑动波数谱以K^-2衰减，而当波数小于拐角波数时，滑动与波数无关^{[25, 27~30]}。

在上述研究成果的基础上，结合发震断裂中段为全新世活动而南北两端为晚更新世活动的性质(说明断裂中段可能是主要释放地震能量的地方)，设定了如图 4所示的断层面滑动分布模型。模型中共有2个凹凸体，主凹凸体位于断层面中段偏南，最大滑动量为6 m，次级凹凸体位于断层面中段偏北，最大滑动量为5 m左右，破裂的初始点位于主凹凸体上(见图 4中五角星)，以大约2.9 km/s的速度向四周扩展，直达边界处停止，断层面上端出露于地表。背景滑动区的上升时间(T_R)根据其与矩震级(M₀)的经验关系(式2)计算得到^[24]。从最近的研究结果来看，凹凸体的上升时间比背景滑动区要小，故本文取之为背景滑动区的0.6倍^[31]。

图  4  断层面上滑动量分布及初始破裂点位置

Figure  4.  Slip distribution on the fault plane and location of the initial rupture point

下载: 全尺寸图片 幻灯片

有限差分方法和随机振动合成方法在模拟地震动时，均采用上述震源模型，此外，随机振动合成方法尚需确定地震应力降参数。相关研究表明中国大陆中小地震释放的应力降主要为0.1~10 MPa，青藏块体东北缘及龙门山断裂带东北缘是高应力降地震集中的区域^[22]，并且全球M_b≥5.5级以上地震释放的应力降具有自相似性的特点，即应力降不随震级的增大而增大^[32]。基于此两方面的观点，本文在计算中将应力降设置为10 MPa。

4.   结果分析

以盆地区域内A、B、C、D四点(位置见图 2)作为示例，分别绘制了低频、高频和宽频带地震动加速度时程(见图 5)。A点位于区域最西侧山脉之中，距离断层最远，第四系覆盖层厚度小，因此覆盖层放大效应较弱。有限差分模拟的低频段PGA约为50 gal，随机法合成的高频段PGA约为210 gal，复合方法构建的宽频带PGA约为210 gal。B点距离断层较远，位于区域西侧的狭谷内，因此地形放大效应较弱。有限差分模拟的低频段PGA约为80 gal，随机法合成的高频段PGA约为250 gal，复合方法构建的宽频带PGA与随机法得到的PGA基本一致。C点位于区域中部地形起伏较缓的山脊之上，距离断层较近，第四系覆盖层厚度小，因此可忽略覆盖层对地震动的放大效应。有限差分模拟的低频段PGA大约为70 gal，随机法合成的高频段PGA约290 gal，复合方法构建的宽频带PGA与随机法得到的PGA基本相当。D点位于区域东南侧地形较为开阔的谷地中，距离断层最近，第四系覆盖层厚度较大，因此模拟结果理应反映出近断层地震动强的特点和覆盖层放大效应。有限差分模拟的低频段PGA为150 gal，随机法合成的高频段PGA为320 gal，复合方法构建的宽频带PGA约为370 gal。

图  5  A、B、C、D四点的地震动加速度时程

Figure  5.  Ground motion acceleration histories of station A, B, C and D

下载: 全尺寸图片 幻灯片

从A、B、C、D四点的地震动加速度傅立叶谱的计算结果可以看出，有限差分方法对1 Hz以上的高频地震动成分无法准确构建，随机合成方法则难以产生1 Hz以下的低频地震动，然而采用复合方法可以很好地弥补这两种方法各自的缺陷，模拟的宽频带地震动均包络了低频地震动和高频地震动，这说明复合方法在构建宽频带地震动时的科学有效性(见图 6)。其中，有限差分方法和随机合成方法模拟A点和B点的低频地震动相差较小，表明该两点处的覆盖层放大效应和地形放大效应不太显著，相比之下，两种方法模拟C点和D点的低频地震动傅立叶谱值相差较大，则分别说明了该处的地形效应和覆盖层效应比较突出，对1 Hz以下的低频地震动产生了强烈的放大作用。

图  6  A、B、C、D四点加速度傅立叶谱

Figure  6.  Fourier spectra of station A, B, C and D

下载: 全尺寸图片 幻灯片

盆地区域内东西向和南北向宽频带地震动PGA的分布情况(见图 7)表明：离断层较近的场点通常具有更大的地震动，随着断层距的增加，地震动PGA逐渐变小。如断层附近(区域东南角)两个水平方向上的PGA均接近900 gal，但在远离断层的区域西北侧，PGA减小至150 gal。此外，天水狭长型谷地的南侧，地震动PGA基本上大于400 gal，对应的地震烈度为Ⅸ~Ⅹ度，具备诱发滑坡的强大动力条件，而河谷北侧的PGA值小于400 gal，地震动强度相对较低。在第五代中国地震动参数区划图(GB 18306-2015)中，天水盆地的基本地震动PGA为300 gal，极罕遇地震动PGA为810~960 gal(基本地震动PGA的2.7~3.2倍)^[33]。与区划图的结果对比可知：盆地南部的模拟PGA均大于区划图中基本地震动PGA，而盆地北部则相反；离断层最近的盆地东南角的模拟PGA与区划图中极罕遇地震动PGA相当，而其它区域模拟结果均明显小于区划图结果。说明区划图是综合所有潜在震源贡献的计算结果，极罕遇地震动PGA代表了某个地区未来很长一段时间内的最大地震动强度，其值偏于保守，而模拟结果对于单次地震引起的滑坡评估和灾害预防等更具有科学的指导意义，两者可以形成有益的补充。

图  7  天水盆地内PGA分布

a-东西向  b-南北向

Figure  7.  Distribution of peak ground acceleration in Tianshui basin

下载: 全尺寸图片 幻灯片

PGV主要受1 Hz以下低频地震动成分的影响，文中的模拟结果能够体现出地形和覆盖层厚度对地震动的贡献。根据盆地区域内东西向、南北向及竖向宽频带地震动PGV的分布情况(见图 8)，可以总结出以下几点认识：(1)两水平方向上的PGV最大均为120 cm/s左右(地震烈度为Ⅹ度)，最小约为20 cm/s(地震烈度为Ⅷ度)；竖向PGV最大为40 cm/s(地震烈度为Ⅸ度)，最小约为10 cm/s(地震烈度为Ⅶ度)；(2)就水平向地震动而言，第四系覆盖层对地震动的放大效应和地形放大效应共同作用，使得盆地区域东侧和中部具有较大PGV，而西侧PGV相对较小；⑶与水平方向不同，竖向地震动在盆地区域东侧较弱，而在中部和西侧较强，特别是最西侧陡峭的山地上，PGV达到了最大值，显示受地形影响显著。此外，竖向地震动也体现了覆盖层厚度的影响，譬如在盆地区域南侧的中间部位，也具有较大PGV。

图  8  天水盆地区域内PGV分布

a-东西向  b-西向  c-竖向

Figure  8.  Distribution of peak ground velocity in Tianshui basin

下载: 全尺寸图片 幻灯片

5.   结语

采用复合方法模拟了当礼县—罗家堡断裂上发生矩震级Mw7.7级大地震时，在天水盆地区域产生的地震动场，分析了设定地震作用下盆地区域内的地震动参数PGA和PGV分布特征，量化了盆地区域内的地震动强度。主要结论有：(1)采用有限差分和随机振动合成相结合的复合方法，可以有效互补，得到盆地内地表宽频带地震动。这为该方法在其它地区的应用提供了经验支持；(2)地震在盆地区域产生了强烈地震动，PGA(峰值加速度)介于150~900 gal，离断层较近的区域东南角的PGA最大，随着断层距的增加，PGA逐渐变小。河谷南侧的PGA值基本上大于400 gal，对应的地震烈度为Ⅸ~Ⅹ度，具备诱发滑坡的强大动力条件，而河谷北侧的PGA值则小于400 gal，地震动强度相对较低；(3)盆地区域PGV(峰值速度)最大为120 cm/s。受第四系覆盖层放大效应和地形放大效应共同影响，水平向地震动在盆地区域东侧和中部具有较大PGV，而西侧PGV相对较小。竖向地震动在盆地区域东侧较弱，而在中部和西侧较强，特别是最西侧陡峭的山坡上，PGV达到了最大值。此外，竖向地震动明显受到覆盖层厚度的影响，譬如在盆地区域南侧的中间部位，也具有较大的PGV。

文中考虑了第四系覆盖层和地形因素对地震动的影响，为大比例尺地震滑坡危险性定量评估提供了依据，然而仍存在一些不足之处：其一是研究重点位于河谷两侧的盆地区域，对盆地外围斜坡地带的地震动响应没有进行研究；其二是位于斜坡地带的地脉动测点以及钻孔资料较少，因而反演得到的覆盖层厚度数据精度较低，这会给地震动参数的准确性带来一定影响，未来应该加强斜坡地带覆盖层厚度探测和地震动响应研究，以便提升区域地震动模拟的可靠性；此外，对于地震动参数与地形、覆盖层厚度等影响因素之间的定量关系，未做深入解释，下一阶段将增加这些方面的研究工作。