Accurate Mesh Generation of Undulating Interface Based on Truegrid and the Ground Motion Simulation
-
摘要: 利用谱元法的规则六面体单元进行网格剖分时,界面起伏较大处会出现阶梯状网格而导致模拟时产生数值散射。为消除阶梯状网格对起伏界面地震动模拟计算的影响,本文基于TrueGrid软件编写了应用程序,提出了起伏界面处六面体单元网格剖分方式,通过该程序可快速建立起伏界面处均匀的六面体网格模型。本文采取了删除四纵列拐角单元、删除一纵列拐角单元以及构造单元过渡环3种剖分方式,解决两个方向上(x-z与y-z方向)单元二合一过渡后拐角处产生扭曲单元的问题。将阶梯状网格经二合一处理后变为斜面网格,并投影至起伏界面,使得网格完全贴合起伏界面,改善了用台阶状网格近似描绘起伏界面的问题。将3种模型通过谱元法进行数值模拟计算验证了该剖分方式的正确性,对比结果发现删除拐角处一纵列单元方式与设置过渡环的方式均可使用,删除四纵列拐角单元方式不推荐使用,本文提出的方案有助于提高谱元法处理起伏界面问题的灵活性。Abstract: When the regular hexahedral element of the spectral element method is used for meshing, a stepped grid will appear where the interface fluctuates greatly, which leads to numerical scattering during simulation. In order to eliminate the influence of the stepped grid on the ground motion simulation calculation of the undulating interface, an application program is written based on the TrueGrid software, and a method of meshing hexahedral element at the undulating interface is proposed. By using this program, a uniform hexahedral mesh model at the undulating interface can be quickly established. This paper adopts three meshing methods: deleting four-column corner elements, deleting one-column corner elements, and constructing element transition rings to solve the problem of twisting elements at the corners after the two-in-one transition of the elements in the two directions (x-z and y-z directions). The stepped grid is transformed into an inclined grid after two-in-one processing, and projected to the undulating interface, so that the grid fits the undulating interface completely, and the problem of using a stepped grid to approximate the undulating interface is improved. Numerical simulations of the three models verify the correctness of the subdivision method. Comparing the results, it is found that both the method of deleting one column of elements at the corner and the method of setting a transition ring are effective,but the method of deleting four columns of corner elements is not recommended. The schemes proposed in this paper can help to improve the flexibility of the spectral element method to deal with undulating interface problems.
-
Key words:
- Mesh generation /
- Hexahedral element /
- Spectral element method /
- Numerical simulation
-
引言
据河北省地震台网测定,2018年2月12日18时31分在河北省廊坊市永清县(39.37°N,116.67°E)发生4.3级地震,震源深度20km。此次地震未造成人员伤亡,在震中区域造成个别房屋墙皮脱落,河北省廊坊市、固安市震感较强,北京市、天津市普遍有感。地震发生后,中国地震局工程力学研究所对永清4.3级地震强震动记录进行了远程和现场记录回收。强震动记录是地震能量在地表或建筑结构上的真实反映,是工程地震、地震学研究的基础数据资料(朱升初等,2018),本文对此次地震构造背景及强震动记录进行了分析,对永清4.3级地震烈度采用建立虚拟台站方法进行了校正,为河北地区烈度调查、地震区划等提供基础资料。
1. 发震构造背景
永清4.3级地震区域构造主要以NE、NW向断裂为主,地震发生在廊固凹陷内(图 1),廊固凹陷是冀中坳陷北部的一个古近系湖盆,整个凹陷大体呈NE走向,南北长约90km,东西宽20—40km,面积约2600km2。该凹陷是受区域伸展并伴有走滑作用形成的典型箕状凹陷,整个凹陷表现为北断南超、西断东超的构造格局。由于受NE向断裂构造控制,表现出整体不均匀下降的箕状断陷特征。该凹陷内曾于1621年发生廊坊永清$ 5\frac{1}{2} $级地震(国家地震局震害防御司,1995)。
图 1 永清4.3级地震构造背景Figure 1. The seismotectonic background map of the Yongqing MS4.3 earthquake(1)施庄断裂(2)南口山前断裂(3)南口-孙河断裂(4)小汤山-东北旺断裂(5)黄庄-高丽营断裂 (6)顺义-前门-良乡断裂(7)南苑-通县断裂(8)新夏垫断裂(9)夏垫断裂(10)大兴凸起东缘断裂 (11)蓟县山前断裂(12)宝坻断裂(13)桐柏断裂(14)河西务断裂(15)天津断裂(16)大城东断裂 (17)沧东断裂(18)海河断裂(19)牛驼镇凸起东缘断裂(20)徐水断裂(21)东垒子-涞水断裂 (22)永定河断裂(23)徐水南断裂本次地震发生在廊固凹陷东侧控制断裂河西务断裂带附近,该断裂北起河西务西北,向西南经大王务,延伸至永清东码头镇,全长近50km,走向NE,倾向SE,为上陡下缓的铲形正断层,断层产状与震源机制参数基本一致(王晓山等,2018),所以永清4.3级地震发震构造可能为河西务断裂。周月玲等(2018)通过浅层地震勘探和钻孔联合剖面探测,结合年代样品测试,对河西务断裂活动性进行了综合研究,揭示断裂上断点埋深约150m或以浅,第四系底界面垂直断错为20—45m,断裂最新活动时代为晚更新世早期。
2. 观测台站及记录获取
2.1 强震动台站分布
本次地震强震动数据主要由中国地震局工程力学研究所提供,共获取3个分量自由场记录74组,常规处理后得到222条3分向加速度记录,获取记录的强震动台站均为土层台,记录台站分布如图 2所示,由图 2可知,触发台站主要分布在震中北部,南部无台站触发,这可能是因为首都圈台站较密集,而河北省境内强震动台站稀疏,或与仪器型号不同且很多台站老旧未改造有关。触发台站的震中距为29.4—185.5km,其中震中距<50km的记录有5组,震中距为50—100km的有29组,震中距为100—150km的有37组,震中距为150—200km的有3组,距震中最近的台站为长子营台(11ZZY),震中距为29.4km。虽然此次地震获取的强震动记录较全面,但由于震级较小,未造成大的灾害损失,且震中附近台站相对稀疏,能供工程使用的强震动记录较少。
2.2 强震动记录获取
在强震动记录使用分析前,首先对波形数据进行消除基线漂移,即原始强震动记录通过减去事前记录平均值的方法进行零线调整(吕俊强等,2015),通过分析计算得到震中距R<50km的5个台站强震动记录基本信息及相关参数,见表 1。其中距震中最近的长子营台记录到的最大峰值加速度为7.9cm/s2,觅子店台记录到的东西向最大峰值加速度为36.1cm/s2、南北向最大峰值加速度为40.4cm/s2,红星中学台记录到的垂直方向最大峰值加速度为10.8cm/s2,5个台站加速度时程曲线如图 3。由表 1可知,离震中较近的长子营台(震中距29.4km)较离震中较远的觅子店台(震中距38.4km)和于家务台(震中距36.8km)各分量加速度小,影响地震动强度的因素除震中距外,还包括仪器型号、场地特征条件和地质构造、场地土层结构、覆盖层厚度、地形及断裂分布等(王海云等,2010;温瑞智等,2013;郭秋娜等,2018)。本文根据温瑞智等(2015)提出的适用于我国强震动台站场地划分经验曲线,采用H/V谱比法对本次地震触发台站进行分类,分类结果见图 4。由图 4可知,记录到的本次地震台站大部分为Ⅰ类场地(51.35%),其次为Ⅱ类场地(35.14%)和Ⅲ类场地(13.51%),其中,觅子店台和红星中学台为Ⅰ类场地,于家务台和张家湾中学台为Ⅱ类场地,长子营台为Ⅲ类场地。长子营台与觅子店台和于家务台记录差异较大,这可能与场地类型有关。
表 1 强震动记录(震中距<50km)及相关参数Table 1. Strong motion records(epicenter distance < 50km) and related parameters台站名称 台站代码 仪器型号 震中距/km PGA/cm·s-2 EW SN UD 觅子店 11MZD REFTEK/SLJ-100 38.4 36.1 40.4 8.5 红星中学 11YHZ REFTEK/SLJ-100 43.4 -5.4 -4.6 10.8 于家务 11YJW ETNA/ES-T 36.8 -21.7 -18.5 -8.1 张家湾中学 11ZJW REFTEK/SLJ-100 47.9 -5.8 4.2 -4.4 长子营 11ZZY REFTEK/SLJ-100 29.4 5.0 -7.9 -3.5 3. 谱分析
加速度反应谱是目前国内外工程抗震设计的重要依据,反应谱形状、幅值及卓越周期的差异体现了场地条件、震级对地面运动反应谱的影响(温瑞智等,2013;徐钦等,2017)。为避免地震时建筑结构发生共振,确保结构安全,有必要对加速度反应谱成分进行分析。本文选取震中距<50km的5个典型强震动台站进行频谱分析,阻尼比取0.05,得到本次地震EW、SN、UD向加速度反应谱(图 5)。
由图 5可知,觅子店台记录到了本次地震EW向和SN向加速度反应谱最大值,红星中学台记录到了本次地震SN向加速度反应谱最大值。觅子店台与长子营台在EW和SN向差别最大,觅子店台与张家湾中学台在SN向差别最大。5个典型台站记录到的加速度反应谱高频成分丰富,谱卓越周期主要为0.1—0.3s,0.3s后迅速下降。根据现场调查,震区主要房屋类型为土木结构、木房顶砌体结构、预制板砌体结构、钢结构简易房屋4种,本次地震除少数土木结构房屋原有裂缝扩大、墙皮开裂外,其他类型房屋未发现新开裂、屋顶掉瓦、墙皮掉灰现象,中小城市主要建筑物结构自振周期为0.3—1s(王文才等,2018),因此本次地震对该震区自振周期范围内的建筑物破坏影响较小。通过对比不同震中距台站各方向反应谱特征周期发现,震中距对特征周期的影响不明显。
傅氏谱能反映地震波在频域上的幅值特性,地震动记录中包含多种振动频率,不同地区的主要振动频率不同,结构自身具有不同频率、幅频响应,因而在不同频率振动作用下,结构振动不同。为避免地震作用下结构发生共振现象,确保结构安全,有必要对地震频率成分进行分析(冉志杰等,2012)。
本文计算了永清4.3级地震震中距50km以内的5个典型台站强震动记录傅氏谱(图 6),由图 6可知典型台站傅氏谱谱型主要为单峰值型,强震动记录EW分量傅氏谱最大峰值为3.5—8.5Hz,SN分量傅氏谱最大峰值为2.9—6.2Hz,UD分量傅氏谱最大峰值为2.1—9.8Hz。总体来看,UD分量频率范围分布较宽,EW分量傅氏谱最大峰值大于SN分量和UD分量傅氏谱峰值。一般而言,随着震中距的增大,地震动记录长周期(低频)越显著,但本次地震5个典型台站特征不明显,可能与仪器型号及场地特征条件有关。另外,EW分量傅氏谱最大峰值随着震中距的增大呈先减小后增大的趋势,而SN分量和UD分量无该特征。
本文以东西向为正方向,统计峰值加速度反应谱与永清4.3级地震触发台站同震中所处方位角和震中距的关系(图 7)。由图 7可知,峰值加速度反应谱在同一震中距下整体随着方位角的增大而增大,达70°左右后下降,这可能与附近河西务断裂构造走向有关;峰值加速度反应谱值随着震中距的增加而减小。
4. 地震烈度分析
地震烈度是衡量地震破坏和地震动强弱程度的重要指标。Wald等(1999)认为当烈度较小时,用峰值加速度计算烈度较合适,因为人对加速度较敏感。随着强震动观测技术的发展,由于以各观测点峰值加速度作为依据可真实反映地震发生时各点地面运动强度(田秀丰等,2013),基于地震动记录快速确定地震仪器烈度,越来越受到研究者的重视(崔建文等,2016)。目前国内外基于强震动记录资料的地震烈度计算方法较多,不同方法利用不同的地震动参数,包括时间、频度和强度等,建立地震烈度与地震动参数之间的关系。本文采用美国ShakeMap烈度计算方法(李俊等,2010),计算永清4.3级地震事件记录台站烈度,并采用克里格插值方法进行网格化插值,绘制等值线图(图 8)。Wald等(1999)的研究认为,当烈度低于Ⅴ度时峰值加速度PGA与烈度相关性明显高于峰值速度PGV,而当烈度为Ⅴ—Ⅶ度时,PGA和PGV与烈度相关程度相当,考虑本次地震烈度较低,故采用峰值加速度计算烈度。由图 9可知,本次地震仪器烈度等震线长轴沿NW走向,且极震区位于震中东北部,而不是震中附近,但由现场宏观烈度调查(图 8)与震源机制解(王晓山等,2018)可知该地震等震线长轴沿NE向展布,这种差异现象的出现可能与台站分布不均匀有关,震中南部无强震动记录,所以图 8中震中北部附近出现数据空白现象,不利于获取实际烈度长轴方向真实性。
可靠快速生成地震烈度可为人员伤亡估计、应急救援决策和工程抢险修复决策提供重要依据,为解决由于缺少地震动观测值出现空白区域现象及烈度长轴方向和实际方向不一致的问题,本文首先利用永清4.3级地震中获取的所有台站加速度峰值,按地震动峰值衰减模型(王国新,2001)确定该地区地震动衰减关系,然后在整个地震区内建立空间随机假设台站进行补点插值,得到新的地震烈度分布图。地震动峰值衰减模型(王国新,2001)如下:
$$ \lg Y = {{\rm{C}}_1} + {{\rm{C}}_2}M + {{\rm{C}}_3}M\lg \left({R + {R_0}} \right) + {{\rm{C}}_4}\lg \left({R + {R_0}} \right) \pm \varepsilon $$ (1) 式中,Y为PGA(峰值加速度),M为震级,R为震中距,C1、C2、C3、C4为待求常数,R0为近场饱和因子,ε为统计分析误差项。采用通用全局优化法拟合得到地震动峰值加速度随震级和震中距的变化规律,考虑本次地震震级较小,故选取震中距<100km的Ⅰ类场地强震动数据进行拟合,其中水平向按EW、NS向加速度峰值矢量拟合回归分析,拟合曲线如图 9,计算得到的衰减关系式为:
$$ \lg Y = 1.6882 + 0.3926M - 0.1546M\lg \left({R + 5} \right) + 0.7000\lg \left({R + 5} \right) \pm 0.2287 $$ (2) 峰值加速度随震中距的变化如图 10所示,由图 10可知峰值加速度随着震中距的增加而减小。同一震中距不同台站测得的峰值加速度离散性较大。通过拟合各台站实际PGA值,与俞言祥等(2013)地震动衰减关系进行对比,可知20km范围内2条曲线吻合较好,而20km外拟合曲线呈衰减相对缓慢状态,这可能与本次地震震级较小、记录台站均为土层台、远场有一定放大作用有关(王海云,2011)。由观测值得到的预测方程虽存在一定离散性,但整体趋势与实际结果保持一致。
在震区空白区域随机产生N个空间坐标的假设台站,按式(2)计算水平峰值加速度,按美国ShakeMap烈度计算方法转换成仪器烈度,绘制增加假设台站的烈度分布图(图 11),由图 11可知,增加假设台站的仪器烈度分布图与现场烈度调查分布图极震区方向基本一致,进一步验证了河西务断裂为本次地震发震构造的可能,解决了台站分布不均匀导致地震烈度影响场计算出现缺值现象的问题,为以后缺少台站记录震区提供准确快速制作烈度分布图的思路,为震害调查和地震应急救援提供重要依据。
5. 结论
本文介绍了河北永清4.3级地震发震构造背景,并对国家强震动台网中心在本次地震中获取的加速度记录特征进行了分析,为河北省地震动特征和震害调查提供了一定参考,结论如下:
(1) 河北永清4.3级地震在廊固凹陷东侧控制断裂河西务断裂带附近,地震烈度极震区走向与该断裂展布方向一致,结合震源机制解初步判定本次地震发震构造为河西务断裂。
(2) 74个强震动台站获取到强震动加速度记录,震中距位于29.4—185.5km,其中觅子店台记录到的本次地震东西向最大峰值加速度为36.1cm/s2,垂直方向最大加速度为40.4cm/s2,于家务台记录到的南北向最大峰值加速度为18.5cm/s2。采用H/V谱比法对触发台站进行分类,可知Ⅰ类场地居多。
(3) 通过对5个典型台站反应谱分析发现永清4.3级地震加速度反应谱高频成分丰富,谱卓越周期主要为0.1—0.3s。傅氏谱谱型主要为单峰值型,UD分量频率范围分布较宽,EW分量傅氏谱最大峰值大于SN分量和UD分量。
(4) 通过分析峰值加速度反应谱与震中距、方位角的关系,发现峰值加速度反应谱整体随着方位角的增大而增大,达70°左右后下降,这可能与附近河西务断裂构造走向有关,有待进一步研究。其次,峰值加速度反应谱值整体随着震中距的增加而减小。
(5) 由本次地震仪器烈度分布图可知,台站分布不均匀性直接影响烈度分布的真实性,为避免此类问题,拟合出该地区地震动衰减关系,随机产生空间假设台站进行弥补,为以后缺少台站记录震区提供准确快速制作烈度分布图的思路,为震害调查和地震应急救援提供重要依据。由于收集到的强震动数据有限,文中仅采用本次地震获取的峰值加速度拟合当地地震动衰减关系,拟合结果可能存在误差,今后将收集更多的强震动数据进行拟合,以校正该地区地震动衰减关系。
-
表 1 模型基本参数
Table 1. Basic parameters of the model
模型尺寸/m 总步长/步 格点间距/m 时间步长/s 震源时间函数主频/Hz 震源时间函数 6 000×5 000×1 100(长×宽×深) 4 000 100 0.001 4 主频为4 Hz的Ricker子波 -
[1] 董良国, 马在田, 曹景忠等, 2000. 一阶弹性波方程交错网格高阶差分解法. 地球物理学报, 43(3): 411—419 doi: 10.3321/j.issn:0001-5733.2000.03.015Dong L. G. , Ma Z. T. , Cao J. Z. , et al. , 2000. A staggered-grid high-order difference method of one-order elastic wave equation. Chinese Journal of Geophysics, 43(3): 411—419. (in Chinese) doi: 10.3321/j.issn:0001-5733.2000.03.015 [2] 董良国, 李培明, 2004. 地震波传播数值模拟中的频散问题. 天然气工业, 24(6): 53—56 doi: 10.3321/j.issn:1000-0976.2004.06.016Dong L. G. , Li P. M. , 2004. Dispersive problem in seismic wave propagation numerical modeling. Natural Gas Industry, 24(6): 53—56. (in Chinese) doi: 10.3321/j.issn:1000-0976.2004.06.016 [3] 付长华, 高孟潭, 陈鲲, 2012. 北京盆地结构对长周期地震动反应谱的影响. 地震学报, 34(3): 374—382 doi: 10.3969/j.issn.0253-3782.2012.03.009Fu C. H. , Gao M. T. , Chen K. , 2012. A study on long-period response spectrum of ground motion affected by basin structure of Beijing. Acta Seismologica Sinica, 34(3): 374—382. (in Chinese) doi: 10.3969/j.issn.0253-3782.2012.03.009 [4] 高孟潭, 俞言祥, 张晓梅等, 2002. 北京地区地震动的三维有限差分模拟. 中国地震, 18(4): 356—364 doi: 10.3969/j.issn.1001-4683.2002.04.005Gao M. T. , Yu Y. X. , Zhang X. M. , et al. , 2002. Three-dimensional finite-difference simulations of ground motions in the Beijing area. Earthquake Research in China, 18(4): 356—364. (in Chinese) doi: 10.3969/j.issn.1001-4683.2002.04.005 [5] 国家地震局, 1986. 中国地壳上地幔地球物理探测成果. 北京: 地震出版社. [6] 金安蜀, 刘福田, 孙永智, 1980. 北京地区地壳和上地幔的三维P波速度结构. 地球物理学报, 23(2): 172—182 doi: 10.3321/j.issn:0001-5733.1980.02.006Jin A. S. , Liu F. T. , Sun Y. Z. , 1980. Three-dimensional P velocity structure of the crust and upper mantle under Beijing region. Chinese Journal of Geophysics, 23(2): 172—182. (in Chinese) doi: 10.3321/j.issn:0001-5733.1980.02.006 [7] 罗宏伟, 郭鹏, 王树山等, 2013. 基于TrueGrid软件的射孔弹有限元仿真模型参数化建立方法. 测井技术, 37(2): 219—222 doi: 10.3969/j.issn.1004-1338.2013.02.022Luo H. W. , Guo P. , Wang S. S. , et al. , 2013. Parametric modeling method for finite element simulation of perforating charge based on TrueGrid software. Well Logging Technology, 37(2): 219—222. (in Chinese) doi: 10.3969/j.issn.1004-1338.2013.02.022 [8] 蒙茂洲, 2010. 功能强大的网格生成软件——TrueGrid. CAD/CAM与制造业信息化, (1): 57—60. [9] 孙建国, 2007. 复杂地表条件下地球物理场数值模拟方法评述. 世界地质, 26(3): 345—362 doi: 10.3969/j.issn.1004-5589.2007.03.015Sun J. G. , 2007. Methods for numerical modeling of geophysical fields under complex topographical conditions: a critical review. Global Geology, 26(3): 345—362. (in Chinese) doi: 10.3969/j.issn.1004-5589.2007.03.015 [10] 孙若昧, 赵燕来, 吴丹, 1996. 京津唐地区地壳结构与强震的发生──Ⅱ. S波速度结构. 地球物理学报, 39(3): 347—355Sun R. M. , Zhao Y. L. , Wu D. , 1996. Crust structures and strong earthquakes in the Beijing, Tianjin, Tangshan area: Ⅱ. S wave velocity structure. Acta Geophysica Sinica, 39(3): 347—355. (in Chinese) [11] 王夫运, 张先康, 陈棋福等, 2005. 北京地区上地壳三维细结构层析成像. 地球物理学报, 48(2): 359—366 doi: 10.3321/j.issn:0001-5733.2005.02.018Wang F. Y. , Zhang X. K. , Chen Q. F. , et al. , 2005. Fine tomographic inversion of the upper crust 3-D structure around Beijing. Chinese Journal of Geophysics, 48(2): 359—366. doi: 10.3321/j.issn:0001-5733.2005.02.018 [12] 王瑞, 高曙明, 吴海燕, 2020. 六面体网格生成和优化研究进展. 计算机辅助设计与图形学学报, 32(5): 693—708Wang R. , Gao S. M. , Wu H. Y. , 2020. Progress in hexahedral mesh generation and optimization. Journal of Computer-Aided Design & Computer Graphics, 32(5): 693—708. [13] 赵金仁, 张先康, 张成科等, 1999. 香河-北京-涿鹿及其相邻地区壳幔构造与速度结构特征. 地震地质, 21(1): 29—36 doi: 10.3969/j.issn.0253-4967.1999.01.004Zhao J. R. , Zhang X. K. , Zhang C. K. , et al. , 1999. The crust-mantle tectonics and velocity structure characteristics in Xianghe-Beijing-Zhulu and its adjacent areas. Seismology and Geology, 21(1): 29—36. doi: 10.3969/j.issn.0253-4967.1999.01.004 [14] Alford R. M. , Kelly K. R. , Boore D. M. , 1974. Accuracy of finite-difference modeling of the acoustic wave equation. Geophysics, 39(6): 834—842. doi: 10.1190/1.1440470 [15] Alterman Z. , Karal Jr F. C. , 1968. Propagation of elastic waves in layered media by finite difference methods. Bulletin of the Seismological Society of America, 58(1): 367—398. [16] Crase E. , Pica A. , Noble M. , et al. , 1990. Robust elastic nonlinear waveform inversion; application to real data. Geophysics, 55(5): 527—538. doi: 10.1190/1.1442864 [17] Dablain M. A. , 1986. The application of high-order differencing to the scalar wave equation. Geophysics, 51(1): 54—66. doi: 10.1190/1.1442040 [18] Furumura T. , Kennett B. L. N. , Furumura, M. , 1998. Seismic wavefield calculation for laterally heterogeneous whole earth models using the pseudospectral method. Geophysical Journal International, 135(3): 845—860. doi: 10.1046/j.1365-246X.1998.00682.x [19] Graves R. W. , 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bulletin of the Seismological Society of America, 86(4): 1091—1106. [20] Hestholm S. , Ruud B. , 1994. 2D finite-difference elastic wave modelling including surface topography. Geophysical Prospecting, 42(5): 371—390. doi: 10.1111/j.1365-2478.1994.tb00216.x [21] Hestholm S. , Ruud B. , 2002. 3D free-boundary conditions for coordinate-transform finite-difference seismic modelling. Geophysical Prospecting, 50(5): 463—474. doi: 10.1046/j.1365-2478.2002.00327.x [22] Jastram C. , Tessmer E. , 1994. Elastic modelling on a grid with vertically varying spacing. Geophysical Prospecting, 42(4): 357—370. doi: 10.1111/j.1365-2478.1994.tb00215.x [23] Koketsu K. , Fujiwara H. , Ikegami Y. , 2004. Finite-element simulation of seismic ground motion with a voxel mesh. Pure and Applied Geophysics, 161(11—12): 2183—2198. [24] Komatitsch D. , Vilotte J. P. , 1998. The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures. Bulletin of the Seismological Society of America, 88(2): 368—392. [25] Komatitsch D. , Tromp J. , 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophysical Journal International, 139(3): 806—822. doi: 10.1046/j.1365-246x.1999.00967.x [26] Komatitsch D. , Barnes C. , Tromp J. , 2000. Simulation of anisotropic wave propagation based upon a spectral element method. Geophysics, 65(4): 1251—1260. doi: 10.1190/1.1444816 [27] Komatitsch D. , Tromp J. , 2002. Spectral-element simulations of global seismic wave propagation—I. Validation. Geophysical Journal International, 149(2): 390—412. doi: 10.1046/j.1365-246X.2002.01653.x [28] Komatitsch D., Tsuboi S., Tromp J., 2005. The spectral-element method in seismology. In: Levander A, Nolet G, eds., Seismic Earth: Array Analysis of Broadband Seismograms. Washington: American Geophysical Union, 157. [29] Lan, H. Q. , Zhang, Z. J. , 2011. Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface. Bulletin of the Seismological Society of America, 101(3): 1354—1370. doi: 10.1785/0120100194 [30] Levander A. R. , 1988. Fourth-order finite-difference P-SV seismograms. Geophysics, 53(11): 1425—1436. doi: 10.1190/1.1442422 [31] Madariaga R. , 1976. Dynamics of an expanding circular fault. Bulletin of the Seismological Society of America, 66(3): 639—666. doi: 10.1785/BSSA0660030639 [32] Moczo P. , 1989. Finite-difference technique for SH-waves in 2-D media using irregular grids-application to the seismic response problem. Geophysical Journal International, 99(2): 321—329. doi: 10.1111/j.1365-246X.1989.tb01691.x [33] Moczo P. , Kristek J. , Bystrický E. , 2001. Efficiency and optimization of the 3-D finite-difference modeling of seismic ground motion. Journal of Computational Acoustics, 9(2): 593—609. doi: 10.1142/S0218396X01000681 [34] Olsen K. B. , Archuleta R. J. , Matarese J. R. , 1995. Three-dimensional simulation of a magnitude 7.75 earthquake on the San Andreas fault. Science, 270(5242): 1628—1632. doi: 10.1126/science.270.5242.1628 [35] Olsen K. B. , Madariaga R. , Archuleta R. J. , 1997. Three-dimensional dynamic simulation of the 1992 landers earthquake. Science, 278(5339): 834—838. doi: 10.1126/science.278.5339.834 [36] Olsen K. B. , 2000. Site amplification in the Los Angeles basin from three-dimensional modeling of ground motion. Bulletin of the Seismological Society of America, 90(6B): S77—S94. doi: 10.1785/0120000506 [37] Patera A. T. , 1984. A spectral element method for fluid dynamics: laminar flow in a channel expansion. Journal of Computational Physics, 54(3): 468—488. doi: 10.1016/0021-9991(84)90128-1 [38] Pitarka A. , 1999. 3D elastic finite-difference modeling of seismic motion using staggered grids with nonuniform spacing. Bulletin of the Seismological Society of America, 89(1): 54—68. doi: 10.1785/BSSA0890010054 [39] Priolo E. , Carcione J. M. , Seriani G. , 1994. Numerical simulation of interface waves by high-order spectral modeling techniques. The Journal of the Acoustical Society of America, 95(2): 681—693. doi: 10.1121/1.408428 [40] Robertsson J. O. A, 1996. A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography. Geophysics, 61(6): 1921—1934. doi: 10.1190/1.1444107 [41] Seriani G. , 1998. 3-D large-scale wave propagation modeling by spectral element method on Cray T3E multiprocessor. Computer Methods in Applied Mechanics and Engineering, 164(1—2): 235—247. doi: 10.1016/S0045-7825(98)00057-7 [42] Smith W. D. , 1975. The application of finite element analysis to body wave propagation problems. Geophysical Journal International, 42(2): 747—768. [43] Tarrass I. , Giraud L. , Thore P. , 2011. New curvilinear scheme for elastic wave propagation in presence of curved topograph. Geophysical Prospecting, 59(5): 889—906. doi: 10.1111/j.1365-2478.2011.00972.x [44] Tessmer E. , Kosloff D. , Behle A. , 1992. Elastic wave propagation simulation in the presence of surface topography. Geophysical Journal International, 108(2): 621—632. doi: 10.1111/j.1365-246X.1992.tb04641.x [45] Xue Z. , Dong L. G. , Li X. B. , et al. , 2014. Discontinuous galerkin finite-element method for elastic wave modeling including surface topography. Chinese Journal of Geophysics, 57(4): 1209—1223. [46] Zahradník J. , Moczo P. , Hron F. , 1993. Testing four elastic finite-difference schemes for behavior at discontinuities. Bulletin of the Seismological Society of America, 83(1): 107—129. [47] Zhang W. , Chen X. F. , 2006. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation. Geophysical Journal International, 167(1): 337—353. doi: 10.1111/j.1365-246X.2006.03113.x [48] Zhang W. Q. , Zhang Z. G. , Fu H. H. , et al. , 2019. Importance of spatial resolution in ground motion simulations with 3-D basins: an example using the Tangshan earthquake. Geophysical Research Letters, 46(21): 11915—11924. doi: 10.1029/2019GL084815 [49] Zhou H. , Chen X. F. , 2010. A new technique to synthesize seismography with more flexibility: the Legendre spectral element method with overlapped elements. Pure and Applied Geophysics, 167(11): 1365—1376. doi: 10.1007/s00024-010-0106-0 -