Spectral element-finite difference hybrid methods for seismic wave simulation
-
摘要:
为了充分发挥SEM和FDM各自在数值模拟复杂介质中地震波传播的优点，同时避免各自缺点，本文提出了模拟地震波传播的谱元-有限差分（SEM-FDM）混合方法．该混合方法在起伏地形附近采用谱元法计算地表附近的地震波传播；在远离起伏地形的区域之下，采用有限差分法计算地震波传播；通过构建数据交换区实现两种方法的耦合．结果表明，本文提出的谱元-有限差分混合方法能有效模拟任意非均匀介质中地震波传播，具有对复杂地形的刻画能力，能处理自由地表边界条件，且能适应模型内速度间断面．通过模型试算，并与谱元法进行对比，验证了该混合方法用于地震波模拟的有效性和高精度性．
Abstract:Efficient methods for seismic wave simulations play important roles in seismic waveform inversion. Currently, the spectral element method （SEM） and finite difference method （FDM） are mostly frequently used for simulating seismic waveforms in waveform inversion due to their efficiency. Compared with other methods, such as finite element method and pseudo spectral method, the computational costs of the SEM and FDM are relatively low, which is greatly important for seismic wave modeling in large-scale models and waveform inversion. Each method, the SEM or FDM, has its advantages and disadvantages when simulating seismic waveforms in complex medium. For example, the SEM can simulate seismic waves in arbitrary heterogeneous medium with velocity discontinuities due to its high-accuracy. In addition, the free-surface boundary can be automatically satisfied. However, designing a mesh with good quality is generally time-consuming. In some cases, it is difficult to construct an appropriate mesh for the SEM. The FDM generally uses regular meshes for seismic wave modeling. Therefore, it is extremely convenient to design a mesh for the FDM. However, free-surface boundary condition cannot be automatically satisfied and special treatments are required to incorporate the free-surface boundary condition. Although some special treatments are proposed in the past decades, algorithms for the approximation of free-surface boundary condition commonly have very low accuracy. Thus, it is difficult to accurately simulate surface waves using the FDM. To retain the advantages and avoid the disadvantages of the SEM and FDM, we develop a hybrid method, which is call SEM-FDM. For this method, the propagation of seismic waves near the free-surface boundary is simulated by the SEM, and in the area far away from the free-surface boundary is modeled by the FDM. A layer for data exchanging is used for coupling the two methods. The SEM-FDM hybrid method is efficient for the simulation of seismic waves in arbitrary heterogeneous media. The hybrid method has the properties for handling undulating topography, free-surface boundary condition and velocity discontinuities. Based on a plane wave analysis, we derive the stability conditions of the hybrid method, and present the stability conditions for spatial accuracy with different orders. In order to demonstrate the validity of the SEM-FDM hybrid method in seismic wave modeling, we construct four models for numerical tests. The first model is used to show the accuracy and efficiency of the SEM-FDM hybrid method. By a comparison with analytical solutions and results from the FDM and SEM, the SEM-FDM is indeed very suitable for seismic wave modeling. The second and third models are used to display the ability of the SEM-FDM hybrid method for modeling seismic wave in complex medium. Although relatively sample grids are adopted, the SEM-FDM can still accurately model seismic waves. The last is a really geological model, which is used to exhibit the usefulness of the proposed method in real cases. By a comparison with the results from the SEM, we demonstrate that the SEM-FDM hybrid method is a useful tool in real cases. Although the areas near the free-surface area are discretized by irregular elements, the process for generate irregular elements are not difficult. In our future work, we will use Fourier series to approximate the free surface and automatically generate a mesh based on the Fourier series and control parameters.
-
引言
印度板块与欧亚板块的碰撞导致青藏高原的快速隆升，而强烈的挤压作用对青藏高原周边地块的地质地貌及构造演化也造成了巨大影响。青藏高原东南缘由于地处高原的边界地区，地壳高度破碎，深大断裂发育，地震活动频繁且剧烈，在过去的30年里发生M6.0以上地震百余次，是现今地壳形变与地震活动最为强烈的地区之一。前人为了解释青藏高原的抬升过程与形变机制提出了刚性地块挤出模型（Tapponier，Molnar，1976；Tapponnier et al，1982；Avouac，Tapponnier，1993）和地壳流模型（Royden et al，1997；Clark，Royden，2000），这两种模型均认为青藏高原东南缘是高原地壳物质转移的重要场所与通道。因此，研究该地区的地壳结构，对进一步厘清地震孕育环境以及青藏高原构造动力学过程等均具有重要意义。
青藏高原东南缘地区因其独特的构造背景以及频繁的地震活动一直是研究的热点，近年来众多研究人员利用多种地球物理手段开展了大量研究，揭示了该区域地壳结构的基本特征，例如：体波层析成像的结果显示，研究区的中下地壳存在大范围的低速异常，多数强震发生在中上地壳内的高速异常区以及高、低速异常转换带上（王椿镛等，2002；韦伟等，2010）；面波与背景噪声层析成像的结果显示，研究区内壳内低速层的分布受断裂与构造边界的限制（Yao et al，2008），强震大多发生在中上地壳内的低速异常区以及高、低速异常分界面上（王琼，高原，2014；郑定昌等，2014；潘佳铁等，2015）；接收函数的研究结果则显示研究区内S波速度结构具有很强的横向不均匀性，低速层分布的深度和范围差别很大（吴建平等，2001；胡家富等，2003；李永华等，2009）。综合以上研究成果，在青藏高原东南缘的中下地壳存在低速层已得到公认，但是低速层的位置与分布形态仍然存在较大争议，地震分布与速度结构之间的关系也尚未明确，因此仍有必要继续开展该区域深部结构的研究工作。
为此，本文收集整理了云南及周边地区81个固定台站的观测数据，拟采用双差层析成像方法（Zhang，Thurber，2003，2006）进行地震重定位，获取青藏高原东南缘的三维地壳速度结构，分析地震分布与速度结构之间的关系，并试图结合地质构造背景来讨论速度结构特征及其包括的动力学意义，以期为该区的动力学模型提供更多地震学证据的支持。
1. 数据与方法
1.1 数据
本文采用的数据源于中国地震台网中心提供的2010—2016年期间云南及周边地区（图1）81个固定台站记录到的初至Pg波震相报告资料。为了保证反演结果的可靠性，我们对原始地震数据进行严格的筛选，仅选取观测报告中ML≥2.0且走时残差≤0.5 s的震相数据，并根据时距曲线拟合的方法剔除初始数据中存在明显异常的走时信息，即图2中两条绿线所括区间以外的走时。双差层析成像方法需要利用事件对之间的相对走时数据，而满足一定条件的两个事件才能够组成事件对。本文将事件对的参数选取标准设定为：事件对与台站的最大间距为800 km，事件对之间的最大间距为50 km，每个事件最多可以与50个事件组成事件对，每个事件对所需震相数的最小值与最大值分别为8和30。经过筛选，最终获得1万3 135个地震事件的9万7 908条绝对到时数据和209万1 189条相对到时数据参与反演，射线路径分布如图3所示。可见射线均匀地覆盖了整个研究区域，可以得到较好的反演结果。
![]()
图 1 研究区域构造背景及台站分布F1：怒江断裂；F2：澜沧江断裂；F3：南汀河断裂；F4：无量山断裂；F5：金沙江断裂；F6：红河断裂；F7：丽江断裂；F8：程海断裂；F9：绿汁江断裂；F10：安宁河断裂；F11：则木河断裂；F12：小江断裂Figure 1. Regional tectonic settings and distribution of seismic stations in the studied areaF1：Nujiang fault；F2：Lancangjiang fault；F3：Nantinghe fault；F4：Wuliangshan fault；F5：Jinshajiang fault；F6：Honghe fault；F7：Lijiang fault；F8：Chenghai fault；F9：Lüzhijiang fault；F10：Anninghe fault；F11：Zemuhe fault；F12：Xiaojiang fault1.2 双差层析成像方法
双差层析成像方法是Zhang和Thurber （2003，2006）在双差定位方法（Waldhauser，Ellsworth，2000）的基础上发展而来，该方法同时使用绝对走时数据和相对走时数据进行三维速度结构与震源位置联合反演，在反演过程中，先赋予绝对走时数据较高的权重以建立整个区域内的三维速度结构；几次迭代后赋予相对走时数据较高的权重来更好地约束速度结构与震源位置（肖卓，高原，2017）。由于考虑了介质速度结构的空间变化，该方法克服了双差定位对台站到事件对之间速度恒定的假设，所得定位结果精度更高（王长在等，2013）；而相对走时数据的加入，与传统的走时层析成像方法相比可以提高速度结构反演的精度，因此能够揭示更多的精细结构信息（Zhang，Thurber，2003）。
1.3 参数选择
反演所选用的初始模型参考了研究区域内人工地震测深的结果（Zhang，Wang，2009），并根据数据的残差分布适当调整而来（表1）。水平方向上采用0.5°×0.5°的网格间隔，竖直方向上的网格节点分别位于0，5，10，20，30，40 km。针对阻尼最小二乘问题，双差层析成像方法采用带阻尼的最小二乘正交分解（least squares QR factorization，简写为LSQR）算法，以总走时残差的2范数为目标函数进行迭代求解。由于平滑因子和阻尼参数的大小对反演结果的稳定性有较大影响，因此在反演之前需要对不同平滑因子和阻尼参数的组合进行权衡分析（Eberhart-Phillips，1986；王小娜等，2015），选取模型方差变化较小而数据方差显著降低时所对应的参数组合为最优值参与最终反演。本文利用L曲线法搜索最优参数值，设定平滑因子的搜索范围为1—600，阻尼参数的搜索范围为10—1 000，最终选取的最优平滑因子为40，阻尼参数为400 （图4）。
表 1 P波初始速度模型Table 1. Initial model of P wave velocity深度/km vP/（km·s−1） 深度/km vP/（km·s−1） 0 5.50 20 6.23 5 5.89 30 6.50 10 6.02 40 6.90 ![]()
图 4 利用L曲线法所选的最优平滑因子（a）和阻尼参数（b）Figure 4. The optimum smoothing parameter （a） and damping parameter （b） selected by L curve method1.4 分辨率测试
分辨率测试选用棋盘格测试方法（Spakman et al，1993）。该方法首先在初始模型的基础上添加正负相间的扰动合成理论棋盘模型，然后利用棋盘模型计算理论走时数据充当实际观测数据并反演模型参数，最后观察反演结果对理论棋盘模型的恢复程度。若在一些区域内恢复程度良好，则说明这些区域的结果是可信的。本文以0.5°×0.5°网格和±5%高低速异常的棋盘格模型作为理论模型正演计算理论走时数据，然后利用该数据和实际反演所用的初始速度模型进行反演。图5为棋盘格测试的恢复结果，图像显示在地壳5，10，20，30 km深度上，研究区内大部分区域均得到良好的恢复，证明反演结果相对可靠。
2. 主要结果
双差层析成像方法利用地震事件对的相对走时数据来提高震源间相对位置的精度，得到的结果更加准确，重定位前后的地震分布如图6所示。重定位前地震集中分布在5—10 km范围内，重定位后浅层地震数目增多，地震主要分布在20 km深度以上的中上地壳（图7）；而且重定位后所有地震事件的走时残差均方根显著减小，由0.380 s降至0.163 s。
![]()
图 6 重定位前（a）、后（b）的地震空间分布变化Figure 6. Spatial variation of seismic events before （a） and after （b） relocation![]()
图 7 重定位前（a）、后（b）的震源深度变化Figure 7. Focal depth variation before （a） and after （b） relocation图8给出了5，10，20，30 km深度处水平方向的P波速度扰动分布，可见在上地壳5—10 km的深度范围内，速度结构呈现明显的横向不均匀性，其分布特征与地表的地质构造相关。5 km深度上（图8a），西昌东部、大理地区、楚雄地区、景谷—思茅地区均表现为低速异常，这些地区处于不同的盆地内部，表明这些盆地内具有较厚的中生代沉积。在川西北次级地块内显示的低速特征可能与该区域的复理石杂岩沉积有关（Yin，Harrison，2000；潘佳铁等，2015）。熊绍柏等（1993）得到的丽江—者海深地震剖面显示丽江至永胜之间为低速异常，攀枝花以东至安宁河断裂为高速异常，安宁河断裂至小江断裂之间为低速异常，本文的成像结果与之基本吻合。小江断裂的东川地区所表现的低速异常可能是受到了深部热作用的影响。高速异常多集中在攀枝花地区、昭通地区以及南汀河断裂与澜沧江断裂的交会处。
![]()
图 8 不同深度h处P波速度扰动和地震分布图（粗实线AA′ ，BB′ ，CC′ 为地震剖面的位置）Figure 8. P wave velocity disturbance and earthquake distribution at different depth hThe thick lines AA′ ，BB′ ，and CC′ represent the position of the seismic sections（a） h＝5 km；（b） h＝10 km；（c） h＝20 km；（d） h＝30 km；10 km深度上（图8b），高低速异常分布主要受到断裂带的控制，例如：红河断裂、南汀河断裂、澜沧江断裂等整体处于低速异常；小江断裂两侧差异明显，其东侧表现为高速异常，西侧的昆明—玉溪表现为明显的低速异常，这与吴建平等（2013）的层析成像结果相近。腾冲地区的低速异常反映了新生代火山的活动特征，而攀枝花附近的区域仍然表现为高速异常。
在中下地壳20—30 km的速度扰动图上出现大范围的低速异常。20 km深度上（图8c），红河断裂的东西两侧速度结构差异明显。低速异常主要分布在红河断裂东侧的两条带状区域内，其中：一条低速带沿着安宁河断裂、小江断裂分布在川滇菱形地块东侧；另一条主要分布在川西北次级地块内，并穿过丽江断裂向南延伸。两条低速带之间的攀枝花地区表现为高速异常。红河断裂的西侧，除了临沧和腾冲地区表现为低速异常外，其余大部分地区均表现为高速异常。30 km深度上（图8d），高低速异常分布与20 km深度相近，低速带的分布范围较20 km深度的更广。前人的近震层析成像结果同样显示，在研究区的中下地壳内存在大范围的低速带（王椿镛等，2002；Huang et al，2002；韦伟等，2010），但是受观测资料的限制，不同研究显示的低速带位置及形态相差较大。
3. 讨论与结论
3.1 地震分布与速度结构
为了讨论研究区内的地震分布与速度结构之间的关系，我们将重定位后的震源位置分别投影至不同深度的速度扰动剖面上，即将震源深度处于0—7.5 km的地震投影至5 km深度的剖面上，7.5—15 km的地震投影至10 km深度的剖面上，15—25 km的地震投影至20 km深度的剖面上，大于25 km的地震投影至30 km深度的剖面上。从图8可以看出，区域内的地震活动与速度结构的变化明显相关，大多数地震发生在中上地壳的低速异常区内以及高、低速异常区域之间。钱晓东等（2011）统计分析了研究区域内的震源机制解参数认为，该地区主要受到3个方向的应力作用：印缅地块北东向的作用力、川滇菱形地块东南向的作用力以及华南地块北北西向的作用力。由于低速异常区域的强度相对较弱，在这种横向挤压的应力场作用下，这些低速异常区域更容易发生破裂导致地震。一些研究认为地震活动集中分布在中上地壳的低速区域内，可能与断层或微裂隙中存在流体有关（Zhao et al，1996；2000；齐诚等，2006；李永华等，2014；潘佳铁等，2015）。在研究区域内，地壳高度破碎，断层广泛发育，流体的存在会降低断层的强度从而造成孕震区的弱化，使地震更容易发生。另外，高、低速异常交界的区域往往是介质物性发生变化的区域，由于结构和应力场的变化，这些相对脆弱的区域更有利于能量的积累与释放进而引发地震。研究区的中下地壳存在大范围的低速异常区域，黄金莉等（2001）认为下部的低速区域是上部发生强烈地震的重要构造背景，王椿镛等（2002）认为中下地壳低速区域的存在有利于应力集中于上部地壳从而导致地震的发生。
将距离剖面两侧25 km范围内重定位后的地震震中投影至图8c所示的纵向速度剖面图中，以进一步分析地壳速度结构与地震分布的关系，结果如图9所示，可以看出，重定位后的地震沿着断裂带呈现明显的带状或簇状分布特征，地震大多发生于低速异常区域内以及高、低速异常区域之间。青藏高原东南缘地区的地壳结构复杂，地震活动频繁且剧烈，研究分析地震分布与速度结构的对应关系，对于进一步研究地震孕育环境等具有重要意义。
![]()
3.2 低速异常带与地壳流
Royden等（1997）提出的地壳流模型认为，在青藏高原东缘中下地壳存在一个黏性流动层，青藏高原的物质沿此向外流出。这一模型很好地解释了青藏高原东缘的地形变化以及上地壳无明显缩短的现象（Royden et al，1997；Clark，Royden，2000），获得了较多研究人员的支持，而寻找地壳流存在的证据则成为深部地球物理学必须面对的一个科学问题。
判断地壳流的主要依据是偏低的地震波速度，Yao等（2008）利用背景噪声成像的结果显示，研究区壳内低速层分布十分复杂，认为中下地壳物质不会发生大规模的流动，而是被断裂及构造边界限制在一定的区域内；韦伟等（2010）的近震P波层析成像结果显示，研究区中下地壳深度处的川滇菱形地块内低速异常广泛发育，认为该低速异常很可能就是青藏高原下地壳流的通道；Liu等（2014）利用背景噪声与接收函数联合反演的结果显示，松潘—甘孜地块内存在大范围的低波速、高泊松比区域，认为是高原内部物质流出所致。
近年来一些研究表明，青藏高原东南缘可能存在两条地壳流通道（Bai et al，2010；Bao et al，2015）。Bai等（2010）的大地电磁观测结果显示，青藏高原东部存在两条壳内高导带，一条环绕喜马拉雅东构造结向南转折通过滇西南地区，另一条沿鲜水河断裂、安宁河断裂、小江断裂向南延伸至印支地块，且这两条高导带被认为是两条独立的地壳流通道；Bao等（2015）关于面波频散与接收函数联合反演的结果表明，研究区壳内20—40 km深度存在两条S波低速带，且这两条低速带的分布较好地对应了上述高导带的位置。然而，另外一些研究人员采用面波频散与接收函数联合反演的方法得到的结果所显示的低速带位置和形态分布均存在差异，例如：郑晨等（2016）的结果显示两条低速带仅向南延伸至24°N左右，并未穿过红河断裂；而Li等（2016）的结果表明两条低速带穿过红河断裂，且最终交会于云南南部。这些结果的不一致一定程度上说明了研究区下方低速体的形态分布和位置较为复杂。
本文结果显示，在研究区的中下地壳同样存在着两条低速带，其中一条沿安宁河断裂、小江断裂分布在川滇菱形地块的东侧，另一条主要分布在川西北次级地块内，并穿过丽江断裂向南延伸，这两条低速带在前人的研究中均有所体现（Bao et al，2015；郑晨等，2016；Li et al，2016）。在低速带分布的位置同样存在高衰减（周龙泉等，2009；Zhao et al，2013）、高热流（Hu et al，2000）、高vP/vS （Sun et al，2014）等异常，这表明壳内物质存在部分熔融，其在重力的作用下可发生塑性流动。此外，王琼等（2015）基于背景噪声的瑞雷波方位各向异性的研究显示，在中下地壳深度范围内，川西北次级地块和小江断裂附近表现为强各向异性，本文中低速带的延展方向与其快波优势方向基本一致，进一步证明了处于低速带分布位置的壳内物质可以发生塑性流动，造成矿物的定向排列从而导致各向异性现象的产生。此外，这两条低速带并未完全连通，而是被攀枝花及周边区域的高速异常体分隔。地质研究表明，中晚二叠纪时期发生了一次巨大的玄武岩喷发事件，形成了著名的峨眉山大火成岩省（徐义刚，钟孙霖，2001；张招崇等，2006），这次巨大的玄武岩喷发事件很可能是由古地幔柱的活动作用所引起，攀枝花地区位于玄武岩喷发前地幔柱活动的核心部位（He et al，2003）。地幔柱活动导致的基性和超基性岩体残留在地壳内部，表现为高速异常（吴建平等，2013；徐涛等，2015），这些高速物质可能对青藏高原物质向南转移起到了一定的阻挡作用。
在图9所示的纵向速度剖面中均能观测到壳内低速带的存在，且低速带主要沿着研究区内绿汁江断裂、小江断裂等走滑断裂分布。Leloup等（1999）认为走滑断裂相对运动产生的热量能够降低中下地壳的黏度和地震波速度，从而形成低速层；而低黏度的低速层也易与上地壳发生相对运动。本文地震重定位的结果显示少量地震发生在低速带内部，绝大多数地震发生在低速带的边界区域，表明壳内低速带的存在可能会促使断层发生运动而引发地震。另外，沿小江断裂分布的低速带的延伸方向在26°N发生偏转，由与断裂平行的近南北向逐渐转变为西南向，与全球定位系统显示的速度场方向（相对于华南地块）一致，这表明低速带的分布可能对上地壳运动产生影响。张培震（2008）认为造成该地区地壳发生顺时针旋转的深部动力学机制可能是中下地壳软弱层物质的流动。结合前人的研究结果，本文认为青藏高原东南缘存在两条低速带，它们可能是青藏高原物质向南逃逸的两条通道。
3.3 结论
本文采用双差层析成像方法，利用2010—2016年中国地震台网中心提供的云南及周边地区Pg波走时资料，进行地震重定位并获得了青藏高原东南缘的三维地壳速度结构，主要结果如下：
1） 经重定位后，走时残差明显减小，地震主要分布在20 km深度以上的中上地壳；地震分布与速度结构存在一定关系，大多数地震发生在中上地壳的低速异常区内以及高、低速异常区域之间；
2） 层析成像结果显示，研究区上地壳速度结构存在明显的横向不均匀性，速度分布特征与地表的地质构造相关；
3） 研究区的中下地壳存在两条低速带，一条沿安宁河断裂、小江断裂分布在川滇菱形地块的东侧；另一条主要分布在川西北次级地块内，并穿过丽江断裂向南延伸。两条低速带被攀枝花及周边区域的高速异常体分隔，它们可能是青藏高原中下地壳物质向南逃逸的两条通道。
中国地震台网中心提供了震相数据，中国科学技术大学张海江教授提供了tomoDD-SE程序，两位评审专家在稿件修改过程中提出了宝贵意见，本文图件绘制采用GMT软件和Matlab软件，作者在此一并表示衷心的感谢。
-
![]()
图 1 SEM-FDM混合方法的网格模型示意图
Figure 1. Schematic of the grid model in the SEM-FDM hybrid method
![]()
图 2 数据交换区有限差分（a）和GLL （b）网格点上网格点上的波场数据插值示意图
（a） 数据交换区的位移插值；（b） 数据交换区底边界GLL点上的位移插值
Figure 2. Schematic of data interpolations for the wavefields on the FD （a） and GLL （b） grid nodes in the data exchange domain
![]()
图 3 SEM-FDM混合方法模拟地震波传播示意图
Figure 3. Schematic of wavefield simulation using the SEM-FDM hybrid method
![]()
图 4 用于检验SEM-FDM混合方法计算效率的均匀介质模型
模型0—1 km深度部分用SEM计算；1—2 km部分用FDM计算。黑色五角星代表震源；两个黑色三角形对应两个接收器
Figure 4. A homogeneous model used to verify the efficiency of the SEM-FDM hybrid method
The wavefields in upper （0−1 km depths） and lower （1−2 km depths） parts of the model are calculated by the SEM and FDM，respectively. A star represents a source；two triangles denote two receivers R1 and R2
![]()
图 5 均匀介质模型中t＝0.5 s时刻的波场水平分量（左）和垂直分量（右）快照
（a，b） SEM-FDM方法计算；（c，d） SEM方法计算；（e，f） FDM方法计算
Figure 5. Snapshots of displacement wavefield of horizontal （left） and vertical （right） components in a homogeneous model at t＝0.5 s generated by SEM-FDM （a，b），SEM （c，d） and FDM （e，f），respectively
![]()
图 6 由解析解得到的合成地震记录及其与三种数值方法得到的地震记录对比
（a，b） R1处位移水平和垂直分量；（c，d） R1处合成波形水平和垂直分量的误差；（e，f） R2处位移水平和垂直分量；（g，h） R2处合成波形水平和垂直分量的误差. 黑线、蓝线、红线和青色为解析解、SEM-FDM、SEM和FDM
Figure 6. Comparison of synthetic seismograms generated by analytical solution and three numerical methods
（a，b） The horizontal and vertical components of displacements recorded at the R1；（c，d） the errors of the horizontal and vertical components of the displacements at the R1；（e，f） the horizontal and vertical components of displacements at the R2，respectively；（g，h） the errors of the horizontal and vertical components of the displacements at the R2
![]()
图 7 用于地震波场模拟的层状模型
模型内部存在一倾斜界面（黑色实线）. 五角星代表震源；黑色三角形为两个接收器R1和R2
Figure 7. A layered model used for simulation of seismic wavefields
An inclined interface （solid black line） is in the model. The star located at （1000 m，100 m） represents a source；black triangles R1 and R2 denotes two receivers
![]()
图 8 层状模型中t＝0.5 s时刻的波场水平分量（左）和垂直分量（右）快照
（a，b） SEM-FDM计算所得；（c，d） SEM计算所得.
Figure 8. Snapshots of displacement wavefields of horizontal component （left） and horizontal component （right） in a layered model at t＝0.5 s generated by the SEM-FDM （a，b） and SEM （c，d）
The left plots show the displacement wavefields of horizontal component；left plots display the displacement wavefields of vertical component. The upper and lower plots are generated by the SEM-FDM and SEM，respectively
![]()
图 9 混合方法得到的合成地震记录与参考解对比
（a，b） R1处位移水平和垂直分量；（c，d） R1处合成波形水平和垂直分量的误差；（e，f） R2处位移水平和垂直分量；（g，h） R2处合成波形水平和垂直分量的误差
Figure 9. Comparison of synthetic seismograms generated by the SEM-FDM hybrid method and reference solution
（a，b） The horizontal and vertical components of displacements recorded at the R1，respectively；（c，d） Errors of the horizontal and vertical components of the displacements at the R1；（e，f） the horizontal and vertical components of displacements at the R2；（g，h） the errors of the horizontal and vertical components of the displacements at the R2
![]()
图 10 含起伏地表的层状模型中地震波场模拟
地表曲线代表起伏地形; 模型内部实线表示倾斜界面. 五角星表示震源，其位置为（1 000 m，100 m）. 三角形表示接收器，其平均水平间隔为16 m，为了显示需要，接收器每隔四个显示一个
Figure 10. Simulation of seismic wavefields in a layered model with topography
The topography is marked by the black line. The star located at （1 000 m，100 m） represents a source. The triangles denote receivers，the average horizontal distance of which is about 16 m. For visualization purposes，only one-fifth stations are displayed
![]()
图 11 SEM-FDM混合方法计算的单炮地震记录与SEM计算结果对比
图（a）和（b）分别为SEM-FDM和SEM方法计算的位移水平分量，图（c）为其相对差；图（d）和（e）分别为SEM-FDM和SEM方法计算的位移垂直分量，图（f）为其相对差
Figure 11. A comparison of single shot profiles generated by the SEM-FDM hybrid method and SEM
Plots （a） and （b） show displacement wavefields of horizontal component by the SEM-FDM and SEM，respectively，and plot （c） shows the relative differences between the wavefields shown in plots （a） and （b）. Plots （d） and （e） display displacement wavefields of vertical component by the SEM-FDM and SEM，respectively，and plot （f） displays the relative differences between the wavefields shown in plots （d） and （e）
![]()
图 12 龙门山断裂带地壳模型
F1：龙门山断裂带（LMSF）；F2：岷江断裂（MJF）；F3：虎牙断裂（HYF）；F4：龙泉山断裂（LQSF）；F5：安宁河断裂（ANHF）；SGB: 松潘-甘孜块体; SCB: 四川盆地. （a） 地质构造图. 红色线为AB垂直剖面所在位置，用于二维数值测试；（b） P波速度结构；（c） S波速度结构
Figure 12. The crustal model of the Longmenshan fault zone
F1：Longmenshan fault zone （LMSF）；F2：Minjiang fault （MJF）；F3：Huya fault （HYF）；F4：Longquanshan fault （LQSF）；F5：Anninghe fault （ANHF）. （a） The topography around the Longmenshan fault zone. The black lines denote faults. The red line AB represent a vertical profile，that is used for a 2D numerical test.（b） P-wave velocity model；（c） S-wave velocity model
![]()
图 13 地壳模型中t＝13 s时刻的波场快照. 左图为位移水平分量；右图为位移垂直分量
（a）和（b）为SEM-FDM混合方法计算得到；（c）和（d）为SEM计算得到
Figure 13. Snapshots of wavefields in the crustal model at t＝13 s. The left and right plots show horizontal and vertical components of displacement wavefields
Plots （a） and （b） are generated by the SEM-FDM hybrid method；plots （c） and （d） are computed by the SEM
![]()
图 14 地壳模型中合成地震记录
左图为位移水平分量；右图为位移垂直分量。图（a）和（b）为R1处合成地震记录；图（c）和（d）为R2处合成地震记录. 黑线为参考地震图; 蓝线为SEM-FDM混合方法计算得到
Figure 14. Synthetic seismograms in a crustal model
Left and right plots show synthetic seismograms of horizontal and vertical components，respectively；Plots （a） and （b） present seismograms recorded by the R1；plots （a） and （b） display seismograms recorded by the R2. The black lines are the reference seismograms. The blue lines are the synthetic seismograms computed by the SEM-FDM hybrid method
表 1 空间不同插值阶数SEM的稳定性参数a和b值
Table 1 The stability parameters a and b for the SEMs with different interpolation orders
空间插值阶数 a b 1 1 1 2 0 0.66 3 0 0.56 4 0 0.48 5 0 0.41 6 0 0.35 7 0 0.31 8 0 0.28 9 0 0.25 
下载: 导出CSV
表 2 空间不同差分阶数FDM的稳定性参数a和b值
Table 2 The stability parameters a and b for the different order FDMs
空间插值阶数 a b 2 1.08 1.00 4 1.07 0.80 6 0.93 0.73 8 0.73 0.69 10 0.70 0.67 12 0.55 0.65 14 0.36 0.63 16 0.31 0.62 18 0.21 0.61
下载: 导出CSV
表 3 三种地震波正演模拟方法运行内存和计算时间对比
Table 3 Computer memory and computation time of three forward modeling methods
正演模拟方法 运行内存/kb 计算时间/s SEM-FDM 68157 912.12 SEM 80740 961.31 FDM 28835 821.25
下载: 导出CSV
-
刘少林,杨顶辉,徐锡伟,李小凡,申文豪,刘有山. 2021. 模拟地震波传播的三维逐元并行谱元法[J]. 地球物理学报,64(3):993–1005. Liu S L,Yang D H,Xu X W,Li X F,Shen W H,Liu Y S. 2021. Three-dimensional element-by-element parallel spectral-element method for seismic wave modeling[J]. Chinese Journal of Geophysics,64(3):993–1005 (in Chinese).
刘少林,杨顶辉,孟雪莉,汪文帅,徐锡伟,李小凡. 2022. 模拟地震波传播的优化质量矩阵Legendre谱元法[J]. 地球物理学报,65(12):4802–4815. Liu S L,Yang D H,Meng X L,Wang W S,Xu X W,Li X F. 2022. A Legendre spectral element method with optimal mass matrix for seismic wave modeling[J]. Chinese Journal of Geophysics,65(12):4802–4815 (in Chinese).
刘有山,刘少林,张美根,马德堂. 2012. 一种改进的二阶弹性波动方程的最佳匹配层吸收边界条件[J]. 地球物理学进展,27(5):2113–2122. Liu Y S,Liu S L,Zhang M G,Ma D T. 2012. An improved perfectly matched layer absorbing boundary condition for second order elastic wave equation[J]. Progress in Geophysics,27(5):2113–2122 (in Chinese).
桑莹泉,刘有山,徐涛,白志明,解桐桐. 2021. 远震波场正演模拟方法及应用[J]. 地球与行星物理论评,52(6):569–586. Sang Y Q,Liu Y S,Xu T,Bai Z M,Xie T T. 2021. Forward modeling method and application of teleseismic wavefield[J]. Reviews of Geophysics and Planetary Physics,52(6):569–586(in Chinese).
苏波,李怀良,刘少林,杨顶辉. 2019. 修正辛格式有限元法的地震波场模拟[J]. 地球物理学报,62(4):1440–1452. Su B,Li H L,Liu S L,Yang D H. 2019. Modified symplectic scheme with finite element method for seismic wavefield modeling[J]. Chinese Journal of Geophysics,62(4):1440–1452 (in Chinese).
祝贺君,刘沁雅,杨继东. 2023. 地震学全波形反演进展[J]. 地球与行星物理论评(中英文),53(3):287–317. Zhu H J,Liu Q Y,Yang J D. 2023. Recent progress on full waveform inversion[J]. Reviews of Geophysics and Planetary Physics,54(3):287–317 (in Chinese).
Galis M,Moczo P,Kristek J. 2008. A 3-D hybrid finite-difference-finite-element viscoelastic modelling of seismic wave motion[J]. Geophys J Int,175(1):153–184. doi: 10.1111/j.1365-246X.2008.03866.x
Chen M,Niu F,Liu Q,Tromp J,Zheng X. 2015. Multiparameter adjoint tomography of the crust and upper mantle beneath East Asia:1. Model construction and comparisons[J]. J Geophys Res-Sol Ea,120:1762–1786. doi: 10.1002/2014JB011638
De Basabe J D,Sen M K,2007. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equation[J]. Geophysics,72(6):T81–T95.
De Hoop A T,1959. The surface line source problem[J]. Appl Sci Res,8(4):349–356.
Komatitsch D,Tromp J. 2002. Spectral-element simulations of global seismic wave propagation-Ⅰ. Validation[J]. Geophys J Int,149(2):390–412. doi: 10.1046/j.1365-246X.2002.01653.x
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]. B Seismol Soc Am,88(2):368–392. doi: 10.1785/BSSA0880020368
Lan H Q,Zhang Z J. 2011. Comparative study of the free-surface boundary condition in two-dimensional finite-difference elastic wave field simulation[J]. J Geophys Eng. 8:275–286.
Lei W,Ruan Y Y,Bozdağ E,Peter D,Lefebvre M,Komatitsch D,Tromp J,Hill J,Podhorszki N,Pugmire D. 2020. Global adjoint tomography-model GLAD-M25[J]. Geophys J Int,233(1):1–21.
Liu Y,2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling[J]. Geophys J Int,197(2):1033–1047.
Liu S L,Li X F,Wang W S,Xu L,Li B F. 2015a. A modified symplectic scheme for seismic wave modeling[J]. J Appl Geophys,116:110–120. doi: 10.1016/j.jappgeo.2015.03.007
Liu S L,Li X F,Wang W S,Zhu T. 2015b. Source wavefield reconstruction using a linear combination of the boundary wavefield in reverse time migration[J]. Geophysics,80(6):S203–S212. doi: 10.1190/geo2015-0109.1
Liu S L,Yang D H,Ma J. 2017. A modified symplectic PRK scheme for seismic wave modeling[J]. Comput Geosci-UK,99:28–36. doi: 10.1016/j.cageo.2016.11.001
Liu Y,Yu Z,Zhang Z Q,Yao H J,Wang W T,Zhang H J,Fang H J,Fang L H. 2023. The high-resolution community velocity model V2.0 of southwest China,constructed by joint body and surface wave tomography of data recorded at temporary dense arrays[J]. Sci China Earth Sci,66(10):2368–2385. doi: 10.1007/s11430-022-1161-7
Ma J,Yang D H,Tong P,Ma X. 2018. TSOS and TSOS-FK hybrid methods for modelling the propagation of seismic waves[J]. Geophys J Int,214(3):1665–1682. doi: 10.1093/gji/ggy215
Ma X,Yang D H,Song G J,Wang M X. 2014. A low-dispersion symplecitic partitioned Rung–Kutta method for solving seismic-wave equations:I. Scheme and theoretical analysis[J]. Bull Seismol Soc Am,104(5):2206–2225. doi: 10.1785/0120120210
Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[J]. Geophysics,49(5):533–549. doi: 10.1190/1.1441689
Moczo P,Robertsson J,Eisner L. 2007. The finite-difference time-domain method for modeling of seismic wave propagation[J]. Adv Geophys,48:421–516.
Rao Y,Wang Y. 2013. Seismic waveform simulation with pseudo-orthogonal grids for irregular topographic models[J]. Geophys J Int,194(3):1778–1788. doi: 10.1093/gji/ggt190
Shen W H,Yang D H,Xu X W,Yang S X,Liu S L. 2022. 3D simulation of ground motion for the 2015 MW7.8 Gorkha earthquake,Nepal,based on the spectral element method[J]. Nat Hazards,112:2853–2871. doi: 10.1007/s11069-022-05291-1
Tan P,Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation[J]. Geophys J Int,197:1250–1267. doi: 10.1093/gji/ggu077
Tong P,Zhao D,Yang D,Chen J,Liu Q. 2014. Wave-equation-based travel-time seismic tomography-Part 1:Method[J]. Solid Earth,5(2):1151–1168. doi: 10.5194/se-5-1151-2014
Tromp J,Tape C,Liu Q Y. 2005. Seismic tomography,adjoint methods,time reversal and banana-doughnut kernels[J]. Geophys J Int,160(1):195–216.
Zhang J H,Yao Z X. 2013. Optimized finite-difference operator for broadband seismic wave modeling[J]. Geophysics,78(1):A13–A18. doi: 10.1190/geo2012-0277.1
Zhang W Q,Zhang Z G,Fu H H,Li Z B,Chen X F. 2019. Importance of Spatial Resolution in Ground Motion Simulations With 3-D Basins:An Example Using the Tangshan Earthquake[J]. Geophys Res Lett,46(21):11915–11924. doi: 10.1029/2019GL084815
Zhou H,Liu Y,Wang J. 2021. Elastic wave modeling with high-order temporal and spatial accuracies by a selectively modified and linearly optimized staggered-grid finite-difference scheme[J]. IEEE T Geosci Remote,60(60):1–22.
-
期刊类型引用(24)
1. 郝梅煊，潘正洋，王武星，刘琦，赵国强. 青藏高原东南缘三维地壳形变特征与地震活动. 大地测量与地球动力学. 2025(03): 244-251 . 
百度学术
2. 茶文剑，王伟君，张圆缘，黑贺堂. 利用面波直接反演方法研究滇西北地壳三维S波速度结构. 地震工程学报. 2025(02): 458-467 . 百度学术
3. 武振波，邹昆，苏金蓉，滑玉琎，李萍萍，徐涛. 利用双差走时成像研究青藏高原东缘地壳速度结构. 地球物理学报. 2024(03): 871-888 . 百度学术
4. 茶文剑，金明培，王伟君，黑贺堂. 基于密集台阵的滇西北地区壳幔结构研究及孕震环境探讨. 大地测量与地球动力学. 2024(06): 611-617 . 百度学术
5. 杨建文，金明培，茶文剑，张天继，叶泵. 利用接收函数两步反演法研究小江断裂带及邻区地壳S波速度结构. 地震地质. 2023(01): 190-207 . 百度学术
6. 刘建明，吴传勇，余怀忠，王琼，刘代芹，赵彬彬，高荣，高歌，孔祥艳. 新疆天山中段地壳和上地幔顶部P波速度结构研究. 地球物理学报. 2023(04): 1348-1362 . 百度学术
7. 樊文杰，左可桢，赵翠萍，王光明. 2021年6月10日云南双柏地震序列发震构造分析. 地球物理学报. 2023(05): 1991-2006 . 百度学术
8. 宁铄现，陈永顺. 青藏高原东南缘存在连通的下地壳流吗?. 北京大学学报(自然科学版). 2023(03): 395-406 . 百度学术
9. 马永，张海江，高磊，宋程，毕金孟，高也. 滇西地区地壳三维精细结构成像与构造特征研究. 地球物理学报. 2023(09): 3674-3691 . 百度学术
10. 孙娅，邓士林，柳健新，Syed Muzyan SHAHZAD，陈波. 地震波成像揭示云南省的热液活动和成矿环境（英文）. Transactions of Nonferrous Metals Society of China. 2023(11): 3476-3486 . 百度学术
11. 高家乙，李永华，王志铄，张扬，贾漯昭，李大虎. 青藏高原东南缘地壳速度结构及云南漾濞M_S6.4地震深部孕震环境. 地球物理学报. 2022(02): 604-619 . 百度学术
12. 赵博，高原，马延路. 2021年5月21日云南漾濞M_S6.4地震序列重新定位、震源机制及应力场反演. 地球物理学报. 2022(03): 1006-1020 . 百度学术
13. 李洪丽，刘财，田有，闫冬，李志强. 地震走时层析成像方法在研究强震孕震、发生机制中的应用. 地球与行星物理论评. 2022(02): 148-158 . 百度学术
14. 周少贤，薛梅. 利用双差层析成像方法反演阿拉斯加地区岩石圈速度结构. 地震学报. 2022(03): 374-387 . 本站查看
15. 高天扬，丁志峰，王兴臣，姜磊. 利用接收函数、面波频散和ZH振幅比联合反演青藏高原东南缘地壳结构及其动力学意义. 地球物理学报. 2021(06): 1885-1906 . 百度学术
16. 蒋一然，宁杰远，李春来. 川滇地区地震分布及地壳结构的双差地震层析成像. 华北地震科学. 2021(02): 17-26 . 百度学术
17. 杜广宝，吴庆举，张雪梅. 云南漾濞M6.4地震震区三维速度结构. 地震学报. 2021(04): 397-409+533 . 本站查看
18. 刘森，边银菊，王婷婷，鲁志楠. 云南地区Lg波衰减成像研究. 地震学报. 2021(04): 410-426+533 . 本站查看
19. 杨峰. 滇西北地区主要断裂带及邻区三维P波速度结构的双差地震层析成像. 地震. 2021(03): 42-58 . 百度学术
20. 杜广宝，吴庆举，张雪梅，张瑞青. 四川威远及邻区中小地震活动特征及地壳精细结构研究. 地球物理学报. 2021(11): 3983-3996 . 百度学术
21. Xun Sun，Lianghui Guo. Crustal velocity, density structure, and seismogenic environment in the southern segment of the North-South Seismic Belt, China. Earthquake Science. 2021(06): 471-488 . 必应学术
22. 马永，张海江，高磊，毕金孟. 漾濞Ms6.4地震三维精细速度结构与地震活动研究（英文）. Applied Geophysics. 2021(04): 579-591+595 . 百度学术
23. 王月，孟令媛，韩颜颜，马亚伟，邓世广. 2018年云南通海M_S5.0震群序列重定位及震源区速度结构成像. 地震研究. 2020(02): 331-339 . 百度学术
24. 邓山泉，章文波，于湘伟，宋倩，王小娜. 利用区域双差层析成像方法研究川滇南部地壳结构特征. 地球物理学报. 2020(10): 3653-3668 . 百度学术
其他类型引用(9)
计量
- 文章访问数: 210
- HTML全文浏览量: 25
- PDF下载量: 69
- 被引次数: 33
