引言

成都位于四川省中部,四川盆地西部,介于102°54′E—104°53′E和30°05′N—31°26′N之间,全市东西长192km,南北宽166km,总面积14605km2。其地理位置非常重要,是我国西南政治、经济、文化、旅游和交通的中心。成都市人口众多,各种基础建设较为丰富,大地震的发生必将给该地区经济、人文带来巨大的损失。成都位于南北地震带上,2017年开展的“成都市地震活断层普查详勘”,完成了深部地震构造环境探测、地震活断层详细探测与综合制图等地震探测工作,确定了双石-大川断裂发生最高7.6级地震的可能性。依据这一结论,并结合研究区域地壳三维速度结构等资料,开展了双石-大川7.6级地震的强地面运动模拟,给出了潜在地震在本地区产生的强地面运动结果。

强地面运动的短周期地震动利用随机有限断层方法模拟,该方法从点源随机法(SMSIM)发展起来(Boore, 1983, 2003),已用于多地的地震动预测中(Zonno等,2005, 2006Satyam等,2009),为了提高其模拟效果,学者们对其进行了多次修正(Motazedian等,2005Atkinson等,2009Boore等,2015)。强地面运动的长周期地震动利用谱元方法计算。谱元法(SEM)是由Patera(1984)引入到流体计算中,并被Seriani等(1992)首次用于地震波传播模拟运算当中,目前已被成功地应用于2D和3D弹性、非弹性介质地震波传播、激发的问题中(Priolo等,1994Faccioli等,1997Komatitsch等,1998, 1999)。最终,采用2种结果合成统一的强地面运动(Graves, 2004, 2010Frankel,2009Pacor等,2005Motazedian等,2006Atkinson等,2011Shahjouei等,2015Zhou等,2015),其中长周期利用确定性谱元方法结果获得,而短周期利用随机方法模拟结果产生。

1 地震动预测模拟方法介绍

从理论基础划分,近场强地震动模拟预测方法主要可以归纳为3种:依据波动方程理论的确定性方法,主要适合于模拟长周期地震动;依据随机理论的随机有限断层法,主要适合于模拟短周期地震动;两者结合模拟宽频地震动的混合方法。

1.1 谱元法

对于长周期地震动预测,主要采用谱元方法计算(Komatitsch等,1998)。谱元方法即勒让德(Legendre)谱元方法,由Patera(1984)引入到流体计算中,并被Seriani等(1992)首次用于地震波传播模拟运算当中。目前,谱元法广泛应用于地震模拟领域。

地震产生的位移场u可由波动方程表示:

$ \left\{ \begin{matrix} \rho \bf{\ddot{u}}=\nabla \cdot \bf{ \pmb{\mathsf{ σ}} }+\bf{f} \\ \bf{ \pmb{\mathsf{ σ}} }\text{=}\bf{C}:\bf{ \pmb{\mathsf{ ε}} } \\ \bf{ \pmb{\mathsf{ ε}} }\text{=}\frac{1}{2}[\nabla \bf{u}+{{\left(\nabla \bf{u} \right)}^{\text{T}}}] \\ \end{matrix} \right. $ (1)

其中,u为位移向量,σ为应力张量,C为刚度张量,ρ是密度,f是外力。

波动方程(1)可运用积分形式求解(Komatitsch等,1999)。对于1个有限求解区域,在计算区域内将波动方程(1)乘以1个任意测试向量w,再对整个求解区域Ω进行分步积分,得到:

$ \int_\Omega {\rho {\bf{w}} \cdot {\bf{\ddot u}}{\rm{d\Omega }}} + \int_\Omega {\nabla {\bf{w}}:{\bf{C}}:\nabla {\bf{u}}{\rm{d\Omega }}} {\rm{ = }}\int_\Omega {{\bf{w}} \cdot {\bf{f}}{\rm{d\Omega }}} + \int_\Gamma {{\bf{w}} \cdot {\bf{T}}{\rm{d}}\Gamma } $ (2)

式中,$ \Gamma $是人工吸收边界,由于地表应力为零,因此满足地表条件。谱元法将区域Ω分成不重叠的六面体空间子单元Ωe(三维情况),在单元里取高斯罗巴托勒让德(Gauss-Lobatto- Legendre,GLL)点剖分,并以GLL节点做多项式插值,结合GLL积分规则,用U表示介质位移未知量,可以得到方程:

$ {\bf{M\ddot U}} + {\bf{KU}} = {\bf{F}} $ (3)

其中,M是对角线质量矩阵,K是带宽刚度矩阵,F是源项。M的对角性使得谱元法具有高效率、易于并行运算的计算优势。

1.2 随机有限断层方法

随机有限断层模拟方法的基本步骤是将发震断层划分成一系列小断层子源,由于子源尺度较小,每个子源可以当作点源,利用随机方法(Boore,1983)计算每个子源的地震动,再基于小震合成大震的原则得到大断层的地震动,即:

$ a\left(t \right) = \sum\limits_{i = 1}^{{N_{\rm{L}}}} {\sum\limits_{j = 1}^{{N_{\rm{W}}}} {{a_{ij}}\left({t + \Delta {t_{ij}}} \right)} } $ (4)

其中,NLNW分别是沿着断层走向和倾向的子断层数,aij表示第(ij)个子源引起的场地地震动,$ \Delta {t_{ij}} $表示破裂传播到第(ij)个子源引起的时间滞后和从第(ij)个子源到场地间由于传播距离的不同引起的时间滞后。Motazedian等(2005)提出了“动力学拐角频率”的新概念,使得剖分的子断层大小不再影响模拟结果。

动拐角频率与破裂面积成反比,其计算公式为:

$ {{f}_{{{\text{c}}_{ij}}}}\left(t \right)=\text{4}\text{.9}\times {{10}^{6}}\beta {{\left(\Delta \sigma /{{M}_{\text{0ave}}} \right)}^{{1}/{3}\;}}\cdot {{N}_{\text{R}}}{{\left(t \right)}^{{\text{ }-1}/{3}\;}}\cdot S $ (5)

其中,$ {{f}_{{{\text{c}}_{ij}}}}(t) $是第(ij)个子断块的动力学拐角频率,t是第(ij)个子源被触发的时刻,NRt)是在t时刻已经破裂的子断块数目,M0ave是子断层的平均地震矩,S是表示子断层辐射强度的常数。随着破裂的发展,破裂子断层数增加,$ {{f}_{{{\text{c}}_{ij}}}}(t) $发生变化。

动拐角频率法消除了断层子断块划分对模拟结果的影响,可以更好地反映断层面上滑动分布的不均匀性。

1.3 混合计算方法

谱元法适合计算长周期地震动,随机有限断层法适合计算短周期地震动。为最大限度地发挥2种方法的长处,分别对谱元法计算的长周期地震波和随机有限断层方法计算的短周期地震波在频带之间进行滤波,然后在频域内将滤波后的长周期地震波和短周期地震波进行分段组合,即合成记录中0—1Hz的傅立叶谱由谱元计算结果提供,大于1Hz的傅立叶谱由随机有限断层的计算结果提供,将该傅立叶谱经傅立叶逆变换得到宽频带的最终地震动时程。

2 三维地下介质计算模型

根据成都及邻区的三维速度探测结果1,及所收集的第四系等厚线图、DEM数据、芦山科考地层探测结果,建立了成都市及邻区深部地壳结构三维地下介质计算模型(图 1,红色线为双石-大川断层延地表的展布)。图中显示计算区域内地下P波纵向、横向都呈现非均匀性的分布特征,向下随着深度的增加,P波速度逐步增加,30km深度之下,达到7km/s;而四川盆地内浅层速度较小,在1—3km/s之间变化,具有沉积体地震波低速的特性;地表观测计算区域呈现西高东低的态势,成都地区位于地势较低的四川盆地之中,地下速度构造较为复杂。计算模型范围为29°—31.5°N、102°—105°E,深度为50km至地表。图 1主要显示了深度30km内的介质构造。

1  四川赛思特科技有限责任公司,2017.成都市主要活动断裂地震危险性评价报告.成都市地震活断层普查详勘项目报告.


图 1 三维P波速度结构 Fig. 1 3D structure of P wave velocity
3 震源模型

确定断层的震源机制,建立震源模型,是近断层强地面运动场数值模拟的另一前提。在目标区内,四川盆地西部边界的前山断裂主要分布有双石-大川断裂和彭灌断裂。在双石-大川断裂上,1970年发生过大邑6.2级地震,彭灌断裂则在2008年汶川8.0级地震时发生同震地表破裂,2条断裂的地震活动强度有所不同,双石-大川断裂和彭灌断裂存在级联破裂的可能性较小,因此仅考虑双石-大川断裂发生破裂的情况。依据地震活动性的多种分析方法综合评估,最终确定双石-大川断裂可能发生自南向北的全段破裂,震级为MS 7.6。双石-大川断裂全长约135km,根据断层的空间位置1,在模拟中按走向变化将其分成6段(表 1)。图 1中红色曲线显示了6段弯曲,表 1给出了6段对应的长度、起始端点的经、纬度坐标及产状。双石-大川断裂6段的走向变化范围较大,最小208°,最大230°,有20多度的变化幅度,而6段的倾向变化不大,倾角围绕44°做小幅变化。

1  四川赛思特科技有限责任公司,2017.成都市主要活动断裂地震危险性评价报告.成都市地震活断层普查详勘项目报告.

表 1 双石-大川断裂各段参数 Table 1 The fault parameters of Shuangshi-Dachuan fault

因双石-大川断裂自南向北破裂,依据Somerville等(1999)的研究结论,设定震源的凹凸体位置、大小等计算参数(表 2)。依据表中参数,将破裂地震矩、破裂发生时间过程投影到地表上断层面投影面内(图 2),图 2中最强的位错发生在断层投影面上的浅蓝区域,对应凹凸体的分布区域。图 3显示断层不同地点的破裂时间,破裂由南向北传播,北端大概在地震发生40s后发生破裂。利用谱元法模拟长周期地震动和随机有限断层方法模拟短周期地震动中所采用的设定震源,具有相同的位错分布、破裂速度、破裂形式、子断层分布和凹凸体分布,谱元法采用的震源破裂时间函数为上升时间3.01s的Ramp函数。

表 2 双石-大川断裂计算参数 Table 2 The computation parameters of Shuangshi-Dachuan fault

图 2 凹凸体空间分布 Fig. 2 The distribution of asperities

图 3 双石-大川破裂起始时间 Fig. 3 The rupture initial time of Shuangshi-Dachuan fault
4 计算结果分析
4.1 地震动时程特征

整个计算体(图 1)南北向273km,东西向288km,基于地下结构参数(表 1)、震源参数(表 2)和震源破裂方式(图 23),模拟计算了计算区域内的强地面运动。完成长周期地震波模拟和短周期地震动计算后,利用混合法合成了宽频带模拟地震动。图 4显示了成都市中心的地震动时程,由图可明显看出,合成的位移包含了更多的长周期能量,而合成之后的速度与加速度波形受长周期影响可以忽略,主要由短周期控制。由合成的成都加速度时程上可以看到成都峰值接近300gal。


图 4 合成的成都市宽频带地震动时程 Fig. 4 The simulated ground motion of Chengdu city
4.2 地震动峰值空间分布特征

依据地震动合成记录计算了研究区域的地震动峰值,绘制成PGA、PGV和PGD分布图,见图 57。在PGA分布图(图 5)中,蓝色不规则图框为断层面在地表的投影区域,棕色不规则框为成都地区的行政区域范围,断层投影面里PGA最大值超过1000gal,其位置处于凹凸体之上的地表位置;围绕断层区域,PGA呈现出沿断层走向长圆形空间分布形态,断层地表投影面内PGA最大,从断层两侧远离断层,PGA逐步减小,到达成都中心位置PGA降到280gal左右,向东再逐渐减小,成都东边界处的PGA也达到了200gal,在设定的震源影响下,成都大部分地区的加速度峰值为200—600gal。断层破裂长度长、破裂面积大,断层投影也宽大,PGA的分布主要与断层破裂方式、凹凸体分布相关,断层距离成都中心已有一段距离,但是对成都地区的西部山区影响较大,PGA较大可以引起山区滑坡等地震此生灾害。双石-大川断裂向西倾斜,上盘在断层的西部,断层的向西产状特性减小了对于成都(位于断层下盘)的危害。图 67分别为计算地区的速度峰值PGV和位移峰值PGD空间分布,两者基本形态与PGA相似,都是呈沿断层走向呈长圆形分布。成都市区的PGV为12cm/s,向西增大至50cm/s,向东减小至10cm/s;研究区域内PGD变化范围1—12cm,成都中心的位移峰值约为1.5cm,显示双石-大川断裂带列发生的7.6级地震在成都市区产生的位移不大。


图 5 计算区域加速度峰值(PGA)分布 Fig. 5 The distribution of the peak ground acceleration (PGA) in the computation zone

图 6 计算区域速度峰值(PGV)分布 Fig. 6 The distribution of the peak ground velocity (PGV) in the computation zone

图 7 计算区域位移峰值(PGD)分布 Fig. 7 The distribution of the peak ground displacement (PGD) in the computation zone
4.3 反应谱空间分布分析

图 8为加速度0.3s及1.0s反应谱(PSA)的空间分布。0.3s反应谱显示成都PSA超过1000gal,全境处于800—3000gal之间,在断层投影面内部,PSA超过2000gal,最大超过5000gal。与PGA、PGD和PGV分布相比较,0.3s反应谱分布复杂,与断层走向的相关性减弱,说明混合法模拟地震动较高的频率成分随机性增大,受断层破裂的影响或控制作用变弱,特别是与断层相隔一定距离后(如30km外),断层破裂特性影响性明显降低。1.0s反应谱由于频率较低,断层破裂方式对其控制作用很强,断层面地表投影面内PSA值较大,600—1600gal等值线呈椭圆形沿断层走向展布,最强值出现在凹凸体对应的地表之上,与地震动峰值相比断层的控制作用也有所减弱,不再是完全围绕断层地表投影面呈闭合分布。成都中心的1.0s反应谱值为300—400gal,处于等值线变化减弱的区域,显示成都市区内1.0s反应谱值变化不显著。


图 8 计算区加速度反应谱(PSA)分布 Fig. 8 The distribution of pseudo-spectral acceleration (PSA) in the computation zone
5 结论

本文依据成都活断层探测的资料建立了研究区域地下三维构造模型,构建了断层几何参数、地质参数、破裂参数、震级等震源模型。针对该震源模型和介质模型,利用谱元方法和随机有限断层方法分别模拟了长周期地震动和短周期地震动,两者在频率域结合成宽频地震动,最后依据该宽频地震动结果,研究了PGA、PGV、PGD及反应谱分布特征。

双石-大川断裂全长135km,宽18km,是1个狭长破裂断层,断层面地表投影面狭长,模拟地震动以断层地表投影面为中心呈狭长带状分布,说明断层破裂对于地震动具有强大的控制作用。因断层倾角约44°,地震动分布表现出了上下盘效应,上盘地震动强,下盘远离断层30km,地震动能量较快速度衰减。成都位于断层下盘,市区PGA接近300gal,PGV达到10cm/s,PGD约1.5cm。

依据地震动结果研究了0.3s、1.0s反应谱分布,1.0s反应谱分布特征主要沿断层走向呈现椭圆分布,而0.3s反应谱沿断层走向分布的特征减弱,随机性增强,这说明断层作用主要控制地震动长周期成分。

致谢: 感谢国家超级计算天津中心提供超级计算平台。
参考文献
Atkinson G. M., Assatourians K., Boore D. M., et al, 2009. A guide to differences between stochastic point-source and stochastic finite-fault simulations[J]. Bulletin of the Seismological Society of America, 99(6): 3192-3201. DOI:10.1785/0120090058
Atkinson G. M., Goda K., Assatourians K., 2011. Comparison of nonlinear structural responses for accelerograms simulated from the stochastic finite-fault approach versus the hybrid broadband approach[J]. Bulletin of the Seismological Society of America, 101(6): 2967-2980. DOI:10.1785/0120100308
Boore D. M., 1983. Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra[J]. Bulletin of the Seismological Society of America, 73(6): 1865-1894.
Boore D. M., 2003. Simulation of ground motion using the stochastic method[J]. Pure and Applied Geophysics, 160(4): 635-676.
Boore D. M., Thompson E. M., 2015. Revisions to some parameters used in stochastic-method simulations of ground motion[J]. Bulletin of the Seismological Society of America, 105(2A): 1029-1041. DOI:10.1785/0120140281
Faccioli E., Maggio F., Paolucci R., et al, 1997. 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method[J]. Journal of Seismology, 1(2): 237-251.
Frankel A., 2009. A constant stress-drop model for producing broadband synthetic seismograms:comparison with the Next Generation Attenuation relations[J]. Bulletin of the Seismological Society of America, 99(2A): 664-680. DOI:10.1785/0120080079
Graves R. W., Pitarka A., 2004. Broadband time history simulation using a hybrid approach. In: Proceedings of the 13th World Conference Earthquake Engineering. Vancouver, Canada, 1-6 August.
Graves R. W., Pitarka A., 2010. Broadband ground-motion simulation using a hybrid approach[J]. Bulletin of the Seismological Society of America, 100(5A): 2095-2123. DOI:10.1785/0120100057
Komatitsch D., Vilotte J. P., 1998. The spectral element method:an efficient tool to simulate the seismic response of 2D and 3D geological structures[J]. Bulletin of the Seismological Society of America, 88(2): 368-392.
Komatitsch D., Tromp J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation[J]. Geophysical Journal International, 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x
Motazedian D., Atkinson G. M., 2005. Stochastic finite-fault modeling based on a dynamic corner frequency[J]. Bulletin of the Seismological Society of America, 95(3): 995-1010. DOI:10.1785/0120030207
Motazedian D., Moinfar A., 2006. Hybrid stochastic finite fault modeling of 2003, M 6.5, Bam earthquake (Iran)[J]. Journal of Seismology, 10(1): 91-103. DOI:10.1007/s10950-005-9003-x
Pacor F., Cultrera G., Mendez A., et al, 2005. Finite fault modeling of strong ground motions using a hybrid deterministic-stochastic approach[J]. Bulletin of the Seismological Society of America, 95(1): 225-240. DOI:10.1785/0120030163
Patera A. T., 1984. A spectral element method for fluid dynamics:laminar flow in a channel expansion[J]. Journal of Computational Physics, 54(3): 468-488. DOI:10.1016/0021-9991(84)90128-1
Priolo E., Carcione J. M., Seriani G., 1994. Numerical simulation of interface waves by high-order spectral modeling techniques[J]. The Journal of the Acoustical Society of America, 95(2): 681-693. DOI:10.1121/1.408428
Satyam N., Rao K. S., 2009. Estimation of peak ground acceleration for Delhi NCR using FINSIM, a finite fault simulation technique[J]. International Journal of Geotechnics and Environment, 1(2): 147-159.
Seriani G., Priolo E., Carcione J. M., et al, 1992. High-order spectral element method for elastic wave modeling. In:Proceedings of the 62th SEG Annual Meeting[M]. New Orleans, Louisiana: Society of Expanded Abstracts, 1285-1288.
Shahjouei A., Pezeshk S., 2015. Synthetic seismograms using a hybrid broadband ground-motion simulation approach:application to central and eastern united states[J]. Bulletin of the Seismological Society of America, 105(2A): 686-705. DOI:10.1785/0120140219
Somerville P., Irikura K., Graves R., et al, 1999. Characterizing crustal earthquake slip models for the prediction of strong ground motion[J]. Seismological Research Leteers, 70(1): 59-80. DOI:10.1785/gssrl.70.1.59
Zhou H., Jiang H., 2015. A new time-marching scheme that suppresses spurious oscillations in the dynamic rupture problem of the spectral element method:the weighted velocity Newmark scheme[J]. Geophysical Journal International, 203(2): 927-942. DOI:10.1093/gji/ggv341
Zonno G., Carvalho A., Franceschina G., et al., 2005. Simulating earthquake scenarios using finite-fault model for the Metropolitan Area of Lisbon (MAL). In: Proceedings of the 250th Anniversary of the 1755 Lisbon Earthquake. Portugal, 1-4 November.
Zonno G., Carvalho A., 2006. Modeling the 1980 Irpinia earthquake by stochastic simulation. Comparison of seismic scenarios using finite-fault approaches. In: Proceedings of the 1st European Conference on Earthquake Engineering and Seismology. Geneva, Switzerland, 3-8 September.