Comparison of seismic wavefield simulation between the curvilinear-grid finite difference method and ray method in undulated layered medium
-
摘要: 针对处理起伏地表(或含地下不规则波阻抗界面)条件下发展起来的地震波场数值模拟算法的模拟结果与解析解(大多数情形下无法得到)无法进行对比，且其有效性和正确性难以验证的情况，本文提出了一种可以相互验证波场数值模拟结果与射线追踪数值模拟结果的正确性和有效性的佐证方法，验证了参考射线追踪法.其中，波场数值模拟中采用曲线网格DRP/opt MacCormack有限差分法，射线追踪模拟则采用分区多步三角网格最短路径算法.通过系统对比上述两种方法得到的波场快照、单炮地震记录，以及合成理论地震图的结果显示，本方法相互作证了两种方法所得结果的正确性和有效性.双层和三层起伏层状模型的对比分析结果表明，这种方法不但可以加深理解地震波在复杂介质中的传播规律，同时射线法的引入为清晰识别和标定地震波场数值模拟中各种不同震相提供了一种便捷的途径.
-
关键词:
- 分区多步三角网格最短路径法 /
- 曲线网格有限差分法 /
- 地震波传播 /
- 数值模拟 /
- 起伏层状介质
Abstract: Commonly it is not possible to verify and test the accuracy and validity of the seismic wavefield modeling method, especially when the model surface (or subsurface) is undulated. For this reason, this paper tries an alternative way to mutually verify and test the computational accuracy and the solution correctness of both the ray theory (the multistage triangular shortest-path method) and the wave-equation simulation method (the curvilinear-grid finite difference method of the DRP/opt MacCormack format) in layered isotropic medium. Through the analysis and comparison of wave field snapshot, common source gather profile and synthetic seismogram, it is not only able to verify the accuracy and correctness of each of the methods at least for kinematic features, but also to thoroughly understand the kinematic and dynamic features of the wave propagation in isotropic media. The results show that the curvilinear-grid finite difference method is able to yield the correct results even for complex isotropic media; the multistage triangular shortest-path method is capable of predicting similar kinematic features as the wave-equation simulation method does, which can be used to mutually test each other for methodology accuracy and solution correctness. In addition, with the aid of the ray tracing results, it is easy to identify the multi-phases (or multiples) in the wave field snapshot, common source point gather seismic section and synthetic seismogram predicted by the wave-equation simulation method, which is a key issue for later seismic application. -
引言
考虑到大地震的记忆性，Utsu(1972)，Rikitake(1974)和Hagiwara(1974)基于Reid(1910)提出的弹性回跳理论，提出了一种更新模型，该模型假定断层上大地震的复发符合更新过程． 国外研究人员先后提出了多种概率分布来描述这种更新模型，具体的概率分布模型包括高斯分布(Rikitake，1974)、 韦布尔分布(Hagiwara，1974)、 对数正态分布(Nishenko，Buland，1987)、 布朗过程时间(Brownian passage time，简写为BPT)分布(Ellsworth et al，1999； Matthews et al，2002)等模型．
对于任何一种概率分布模型，若已知大震复发的概率密度函数f(T)和最近一次大震的离逝时间Te，则可利用下式计算断层源在未来一段时间ΔT内发生大震的条件概率(Wesnousky，1986)：
(1) 其中，概率密度函数f(T)的参数主要包括大震的平均复发间隔和反映复发间隔T的自然随机不确定性的参数，如标准差σT和变异系数α等．
在大震发生概率计算过程中存在大量的不确定性，例如： 由于活动断裂上的历史地震数据和古地震数据非常有限，由小样本数据所确定的模型参数存在很大的不确定性； 古地震数据的年代测定存在不确定性； 古地震的完整性存在不确定性； 大震发生概率模型本身也存在一定的不确定性． 然而，在实际大震发生概率计算过程中对不确定性的考虑仍不足，例如： 一般不考虑T，Te等参数的不确定性，而是较主观地采用一个固定值； 一般不考虑古地震的完整性，而是较武断地假定某时段内的古地震是完整的．
不确定性一般可分为认知不确定性和随机不确定性，本文在此基础上又对不确定性的定量化过程和处理过程分别进行了深入分析． 按照概率分布类型不同，将不确定性分布分为连续型分布和离散型分布;按照认知主体的不同，将不确定性的处理分为对个体自身认知不确定性的处理和对不同个体之间认知差异性的处理．最后进一步探讨了具体的不确定性处理方法和原则．
此外，大震发生概率计算过程中的不确定性构成往往很复杂，例如，不确定性因素较多，而不同不确定性因素之间又存在相关性，这使得利用数值方法实现起来比较困难． 针对该问题，本文首先提出利用蒙特卡罗模拟法对大震发生概率模型中的各种不确定性分布直接进行反复随机抽样； 然后，通过对大量的随机模拟结果进行统计分析，得到最终的大震发生概率计算结果； 最后，以东昆仑断裂带塔藏段的大震发生概率为计算实例，在充分考虑各种不确定性的基础上，建立了塔藏段的大震发生概率模型，并利用蒙特卡罗法计算塔藏段未来100年的大震发生概率．
1. 大震发生概率模型的不确定性分析
1.1 不确定性分类
在各种地震危险性模型中，对于不确定性的划分，一般可分为随机不确定性和认知不确定性(Bommer，2003)． 随机不确定性反映的是自然本身所固有的各种随机性； 认知不确定性反映的是基于不完备的信息和知识，人的主观认知和判断的不准确性．
大震危险性分析中的不确定性分析是对不确定性进行定量化处理，进而在大震危险性分析中考虑各种不确定性因素． 随机不确定性和认知不确定性有时可以分离，有时则很难分离． 因此，为了探讨不确定性的具体处理方法，本文在不确定性处理过程中，分别按照概率分布类型的不同和认知主体的不同进行划分，而未将随机不确定性和认知不确定性分离开进行讨论．
首先，按照概率分布类型不同，本文将不确定性分布分为连续型分布和离散型分布． 对于不确定性的定量化描述，根据具体情况的不同，可以采用连续型分布，也可以采用离散型分布． 其中，连续型分布更为客观，但建立相对可靠的分布模型需要足够的样本支持，同时连续型的概率分布不方便计算； 离散型分布的确定相对主观，但却便于计算． 其次，按照认知主体的不同，本文将不确定性的处理分为对个体自身认知不确定性的处理和对不同个体之间认知差异性的处理． 在不确定性的实际处理中，不仅要考虑研究个体，同时还需考虑不同研究人员之间所得出的不确定性定量化结果的差异性(Senior Seismic Hazard Analysis Committee，1997)． 这种差异性反映的是人与人之间存在的认知分歧，为了避免计算结果过于主观和片面，需要在大震危险性分析中考虑这种分歧．
1.2 不确定性的处理方法
为便于讨论不确定性的具体处理方法，根据本文所采用的分类方式，并参考前人在概率地震危险性分析中不确定性处理方面的研究成果(Senior Seismic Hazard Analysis Committee，1997)，本文对不确定性的处理过程参见图 1.
1.2.1 采用连续型分布来定量化不确定性
连续型分布较离散型分布更为精确，但其计算量较大，因此如果计算条件允许，可以采用连续型分布来定量化不确定性． 连续型分布的具体分布形态又分为均匀分布和非均匀分布．
1) 均匀分布. 客观信息量越少，不确定性越大，这是个不变的原则，必须体现在大震危险性分析中的每个具体环节中． 对于基本不可知因素或完全不可知因素的不确定性评价，一般采用均匀加权法，其本质上属于拉普拉斯决策法，即等可能决策法． 均匀分布在实际地震危险性分析中应用广泛，例如： 在潜在震源区内任何地方发生地震的可能性是相同的(Cornell，1968)； 古地震发生在其年代测定不确定性的上下限内任意一年的可能性也是相同的．
2)非均匀分布. 若要用非均匀分布来表示不确定性，首先必须假设一个非均匀分布模型，然后通过参数估计法确定该分布模型的参数． 具体的参数估计法包括矩估计、 最大似然估计、 贝叶斯估计等． 其中矩估计和最大似然估计多用于计算某些自然随机分布的参数，而贝叶斯估计则多用于计算某些模型参数的认知不确定性分布． 如果有连续性的证据或数据，则可以采用非均匀分布来表示不确定性，但必须先定义一个分布模型，甚至是出于个体的主观认识．
1.2.2 采用离散型分布来定量化不确定性
逻辑树是一种常用的离散化不确定性处理方法． 该方法一般分为两种： 一种是加权逻辑树法，一种是不加权逻辑树法．
1) 不加权逻辑树． 对于个体的认知不确定性处理，类似于连续型分布中的均匀分布，不加权逻辑树也是一种等可能加权． 对于不同个体之间的认知不确定性的差异性处理，可以采用不加权逻辑树，这意味着在决策过程中不同个体的认知被赋予相同的权重，即不同个体之间是平等的．
2)加权逻辑树． 利用加权逻辑树对个体认知不确定性的处理，是单独地研究个体基于有限的地震地质资料和信息，对大震危险性评估过程中某个环节的认知不确定性分布所作的离散处理，其中每个逻辑树分支均被主观地赋予一定的权重． 根据有限的地震地质资料和信息，对两种或两种以上可能(或离散型的参数值)的权重进行主观评定． 最终评定结果的确定须参考信息的完备程度，信息量越少，未知因素越多，不确定性也就越大． 例如，某个不确定性环节包括A和B两种可能，即使所有的已知信息均支持A，但如果信息量非常有限，也不可以给出确定的结论，而应参考信息量的完备程度，在充分考虑不确定性的基础上进行权重评估．
对于个体之间的认知不确定性的差异性处理，也可以采用加权逻辑树法，其综合了不同个体认知的权重不平等性和单独个体的认知不确定性． 此外，在不确定性处理过程中还需注意以下两点：
1) 必须保证权重赋予和不确定性分布范围确定过程的一致性原则． 一个独立的人或者独立的小组，在大震危险性评估过程中每个具体的不确定性环节，必须遵循相同的不确定分级原则或权重分配原则．
2) 在评估不确定性分布范围的过程中，需要保守考虑，但又不可以无限地保守． 一些概率极低的极端事件对计算结果的贡献很小，也很难进行概率评估和计算，所以，在实际大震危险性评估过程中可以对不确定性分布作适当的截断处理．
2. 计算实例
本文以东昆仑断裂带东段塔藏段的大震发生概率为计算实例，具体研究在大震发生概率计算过程中如何考虑不确定性．
2.1 东昆仑断裂带塔藏段的古地震
东昆仑断裂带东段位于巴颜喀拉地块与西秦岭地块边界断裂的东段，全长约330 km，自西向东可划分为3段： 玛沁段、 玛曲段和塔藏段(青海省地震局，中国地震局地壳应力研究所，1999)． 国内外研究人员在东昆仑断裂带东段上获得了许多古地震数据，为评估该段的大地震发生概率提供了重要依据． 其中，玛沁段和玛曲段上的古地震数据比较丰富，而塔藏段上的古地震数据量则比较少，完整性较差，最近一次大震的离逝时间也存在很大的不确定性．
根据李正芳等(2012)对前人古地震研究结果的分析，本文给出了塔藏段古地震的发生时间分别为距今(4693±151)年，(7304±500)年和(9136±131)年，其中存在4700年左右的大地震空白期，很可能存在历史地震漏记的现象． 另根据李正芳等(2012)对四川松潘(位于塔藏段东南方向约50 km)历史文献记载的分析，推断塔藏段最近一次地震的发生时间应在公元638年之前，即最近一次地震的离逝时间应介于1377—(4693±151)年之间．
2.2 不确定性分析
由于塔藏段古地震数据量比较少，完整性较差，不管是大震复发间隔还是大震离逝时间，均存在着非常大的不确定性． 具体的不确定性因素包括以下3个方面：
1) 在离逝时间Te的确定过程中，需要考虑古地震记录和历史地震记载的不完整性所造成的最近一次大地震发生年代的不确定性
首先，古地震的完整性研究一直是古地震研究的难题，不同断裂上的古地震研究深度不同，地层的完整性情况也不同，造成不同断裂和不同时段的古地震记录的完整性情况相差较大(冉勇康，邓起东，1999)； 其次，我国历史地震资料保存得相对较多，但就各个历史时期和不同区域来说，历史地震记载的完整性程度相差较大(陈春梅，任雪梅，2014)． 因此，无法通过客观的统计分析方法得到最近一次大地震发生时间的不确定性分布，也无法直接判定某时间段内的大震记录是否完整或缺失的数量，只能主观评估已知信息量能否足够(或充分)判定给定时间段的大地震记载是完整的．
本文提出采用逻辑树法来分析由古地震记录和历史地震记载的不完整性所造成离逝时间Te的认知不确定性． 逻辑树在此处共有两个分支： ① 根据已知的区域古地震研究成果和区域历史地震记载情况，无法确定给定时间段的大地震记载是完整的，即信息量有限，无法得出确定性的结论； ② 根据已知的区域古地震研究成果和区域历史地震记载情况，可以确定给定时间段的大地震记载是完整的，即信息量是足够充分的，可以得出确定性的结论． 根据已知的数据和资料，并征求相关方面专家的意见，对两个逻辑树分支的权重进行主观赋值．
2) 在平均复发间隔T的确定过程中，需要考虑古地震年代测定的不确定性. 由于每次古地震发生的年代均存在不确定性，而对于这种不确定性，根据多探槽古地震分析的逐次限定法可以确定古地震发生年代的上下限(毛凤英，张培震，1995)，但很难确定其具体的分布形态，因此可以采用均匀分布来处理此种不确定性．
3) 在平均复发间隔T的确定过程中需要考虑利用有限的复发间隔数据样本来确定平均复发间隔参数的客观认知不确定性．
在利用塔藏段有限的两个复发间隔数据估计T的过程中，需要同时计算基于小样本数据估计参数T的认知不确定性． 本文采用贝叶斯估计法进行评估，具体做法可以参考郭星和潘华(2016)提出的平均复发间隔T的认知不确定性评估方法．
2.3 考虑不确定性的塔藏段大震发生概率计算模型
在Reid(1910)弹性回跳理论的基础上，国内外研究人员先后提出了多种符合更新过程的概率模型，其中比较常用的概率分布模型有对数正态分布模型(Nishenko，Buland，1987)、 BPT(Ellsworth et al，1999； Matthews et al，2002)模型等．
本文未考虑模型本身的不确定性，直接采用BPT模型计算大震发生概率. BPT分布也被称作逆高斯分布(Seshadri，1983)，其概率密度函数为
(2) 式中： T为断裂段大震的平均复发间隔； α为大震复发间隔的变异系数，反映的是大震复发间隔的随机不确定性，可以通过对大量归一化复发间隔数据的统计分析得到． 在广泛搜集的39个原地复发的地震序列基础上，郭星和潘华(2015)统计分析得到了一个通用的变异系数α＝0.34，本文在发震概率计算过程中直接选取α=0.34．
塔藏段3次古地震的发生时间分别为距今(4693±151)年，(7304±500)年和(9136±131)年，本文利用这3个古地震数据来评估BPT模型的参数T，在参数确定过程中的各种不确定性则采用蒙特卡罗法来处理，最后利用式(1)计算塔藏段未来100年的大震发生概率．
参考对不确定性处理方法的分类讨论，本文采用加权逻辑树法来考虑古地震的不完整性． 考虑到每个逻辑树分支所对应的计算过程均存在一定差异，本文分成A和B两个相互独立的蒙特卡罗模拟过程分别进行计算． 其中分支A假定根据有限的信息，无法确定给定时间段的大地震记载是完整的； 分支B假定信息量足够充分，可以确定给定时间段的大地震记载是完整的． 根据已知的数据和资料，并征求相关方面专家的意见，对两个逻辑树分支的权重进行主观赋值．
分支A的具体实现步骤如图 2所示：
![]()
图 2 塔藏段大震发生概率计算的分支A流程图Figure 2. The process chart of branch A for calculating the occurrence probability of large earthquakes on the Tazang fault1) 考虑到古地震年代测定的不确定性所导致的平均复发间隔参数T的不确定性，在3次古地震发生年代的不确定性范围内生成均匀分布的随机数，这样就可以得到两个随机复发间隔数据. 同时还需要得出最近一次已知古地震的离逝时间Te1，这是考虑到最近一次古地震事件年代测定的不确定性，平均复发间隔参数与大震离逝时间参数之间存在一定的相关性，不可以分别得出.
2) 对于平均复发间隔的参数认知不确定性分布的定量化处理，可以参考郭星和潘华(2016)提出的基于两个复发间隔数据所确定的平均复发间隔T的参数认知不确定性研究，具体方法如下：
若已知两个大震复发间隔的实际数据T21和T22，对于任意平均复发间隔T，根据式(2)得到的复发间隔为T21的概率密度函数为
(3) 而复发间隔为T22的概率密度函数为
(4) 再根据贝叶斯估计法，即可得到T的认知不确定性分布的概率密度函数f(T)
(5) 3) 对于任意平均复发间隔T，其所对应的BPT模型的概率密度函数为f(T|T)，因此，根据式(1)其对应的发震概率计算公式需修改为
(6) 根据f(T)，可以随机生成一个平均复发间隔T，也就同时确定了一个f(T|T)．
4) 根据步骤3中确定的概率密度函数f(T|T)，随机生成下一次地震发生的时间Ts． 若要确定当前的大震离逝时间Te，则需进行如下判断： 若Ts＞Te1，当前的大震离逝时间Te=Te1; 若Ts＜Te1-1377，当前的大震离逝时间Te=Ts; 其它情况则需重新模拟，再随机生成下一次地震发生的时间Ts，直到Ts满足上述两种情况为止．
从保守和计算效率的角度考虑，判定过程中并未考虑记载空段内漏记两次或两次以上古地震的小概率特殊事件的发生．
5) 已知f(T|T)和大震离逝时间Te，可以根据式(6)计算得到未来100年塔藏段大震的发生概率Pn．
6)反复模拟10万次，计算出10万次模拟结果的平均值P．
分支B的具体实现步骤与分支A类似(图 3)，只是第4步有所差异，考虑到分支B假定不可能出现古地震漏记，则直接取大震离逝时间Te=Te1．
![]()
图 3 塔藏段大震发生概率计算的分支B流程图Figure 3. The process chart of branch B for calculating the occurrence probability of large earthquakes on the Tazang fault2.4 大震发生概率计算结果
分支A的大震发生概率计算结果为0.11； 分支B的大震发生概率计算结果为0.16． 而对于A和B两个逻辑分支的权重赋值则需要综合考虑以下几个方面： ① 研究区域的人类活动历史和历史地震记载的完整性情况； ② 研究断裂上的地层完整性情况； ③ 研究断裂上的古地震研究深度．
基于上述3个方面的信息，并征求相关方面专家的意见，本文给出的A和B两个逻辑分支的权重赋值WA和WB分别为0.8和0.2，最后计算出的塔藏段未来100年的大震发生概率为
(8) 3. 讨论与结论
本文对大震发生概率评估过程中的各种不确定性进行了系统分析，在认知不确定性和随机不确定性的基础上对不确定性的定量化处理过程进行了分类研究，并对不确定性的具体定量化方法进行了梳理．
考虑到大震复发概率计算过程中存在大量的不确定性，且不确定性的构成又较为复杂，数值方法实现起来比较困难，故此类不确定性往往被忽略掉． 本文首先明确了大震复发概率模型中的各种不确定性，特别是常被忽略的不确定性，然后利用蒙特卡罗模拟法对大震发生概率模型中的各种定量化的不确定性直接进行随机抽样，最后通过对大量随机模拟结果的统计分析得到了最终的大震发生概率计算结果． 虽然计算过程相对复杂，但可以尽可能地减少一些主观因素，使得计算结果更为客观．
本研究的重点在于探讨如何利用蒙特卡罗法考虑发震概率计算过程中的各种不确定性，对于各个不确定性分支的具体赋值过程则未进行更深入的研究． 而不确定性分支的权重将直接影响最终的计算结果，需要作大量的地震地质调查和统计分析工作，也需要采取科学有效的方法对多位专家的意见进行综合． 在实际大震危险性评估中，怎样考虑和处理这种人与人之间的认知分歧，也需进一步研究．
-
![]()
图 1 含起伏界面的层状模型示意图
(a)计算反射波、透射波、转换波的模型示意图；(b)模型参数化时所采用的三角单元十字表示震源；圆圈表示主节点；圆点表示次级节点
Figure 1. Schematic diagram of the undulating layered medium model
(a) A schematic explanation for computing the transmitted, reflected and converted arrivals; (b) A triangular cell used in the model parameterization The cross indicates earthquake source, circles indicate primary nodes, and dots indicate secondary nodes
![]()
图 2 起伏层状均匀介质中多震相射线路径
(a)一次纯P波反射波；(b)一次纯S波反射波；(c)一次P-S波反射转换波；(d)一次S-P波反射转换波
Figure 2. The multi-phase ray paths of undulating layered medium
(a) Purely primary reflected P-wave; (b) Purely primary reflected S-wave; (c) Primary reflected and converted from P to S wave; (d) Primary reflected and converted from S to P wave
![]()
图 3 曲线网格和计算网格的网格变换关系
Figure 3. Mapping between the boundary-conforming grid and uniform grid
![]()
图 5 双层介质模型中t=0.6 s时x分量波场快照.虚线表示模型速度界面
(a)曲线网格DRP/opt MacCormack差分法模拟结果；(b) MTSPM射线追踪法结果
Figure 5. Snapshots of x component wave field at the moment of 0.6 s for the two-layer model where the dashed lines indicate the model subsurface
(a) Results by the curvilinear-grid finite difference method of the DRP/opt MacCormack; (b) Results by MTSPM ray tracking algorithm ① P1; ② S1; ③ P1P1; ④ P1S1; ⑤ S1P1; ⑥ S1S2; ⑦ P1S2; ⑧ S1P2; ⑨ P1P2; ⑩ S1S1; ⑪ P1S1; ⑫ P1P1
![]()
图 6 双层模型中x分量单炮地震记录剖面
(a)曲线网格DRP/opt MacCormack差分法模拟结果; (b) MTSPM射线追踪法结果
Figure 6. Common source point gather of x component for two-layer model
(a) Results by the curvilinear-grid finite difference method of the DRP/opt MacCormack; (b) Results by MTSPM ray tracking algorithm ① P1; ② S1; ③ P1P1; ④ S1P1; ⑤ P1S1; ⑥ P1P1; ⑦ S1S1; ⑧ P1S1; ⑨ S1P1; ⑩ S1S1
![]()
图 7 基于曲线网格DRP/opt MacCormack差分法得到的双层介质模型中x分量(a)和z分量(b)的合成理论地震图
垂线为MTSPM算法所预测的各震相理论初至时间
Figure 7. Synthetic seismogram for two-layer model with the curvilinear-grid finite difference method of the the DRP/opt MacCormack for x component (a) and z component (b)
Vertical lines indicate the onset times for possible phases by MTSPM algorithm. ① P1; ② S1; ③ P1P1; ④ S1P1; ⑤ P1P1; ⑥ S1S1; ⑦ P1S1; ⑧ S1P1; ⑨ P1P1P1; ⑩ S1S1; ⑪ P1P1P1
![]()
图 9 三层均匀介质模型中t=0.6 s时x分量波场快照.虚线表示模型速度界面
(a)曲线网格DRP/opt MacCormack差分法模拟结果; (b) MTSPM射线追踪法结果
Figure 9. Snapshots of x component wavefield at the moment of 0.6 s for the three-layer model where the dashed lines indicate the model subsurface
(a) Results by the curvilinear-grid finite difference method of the DRP/opt MacCormack; (b) Results by MTSPM ray tracking algorithm ① P1; ② P1P1; ③ P1S1; ④ S1; ⑤ S1P1; ⑥ S1S1; ⑦ S1S2; ⑧ S1P2; ⑨ P1S2; ⑩ P2P2S3; ⑪ P1P2P3; ⑫P1P2P2; ⑬ P1P2S2; ⑭ P1P1; ⑮ P1S1
![]()
图 10 三层模型中x分量单炮地震记录剖面
(a)曲线网格DRP/opt MacCormack差分法模拟结果；(b) MTSPM射线追踪法结果
Figure 10. Common source point gather of x component for three-layer model
(a) Results by the curvilinear-grid finite difference method of the DRP/opt MacCormack; (b) Results for MTSPM ray tracking algorithm ① P1; ② S1; ③ P1P1; ④ S1P1; ⑤ P1S1; ⑥ P1P1; ⑦ S1S1; ⑧ P1P2P2P1; ⑨ P1S1; ⑩ P1S2P2P1; ⑪ S1P1; ⑫ P1P2P2S1
![]()
图 11 基于曲线网格DRP/opt MacCormack差分法得到的三层介质模型中x分量(a)和z分量(b)的合成理论地震图
垂线为MTSPM算法各震相的理论初至时间
Figure 11. Synthetic seismogram for three-layer model with the curvilinear-grid finite difference method of the DRP/opt MacCormack for x component (a) and z component (b)
Vertical lines indicate the onset times for possible phases by MTSPM algorithm. ① P1; ② S1; ③ P1P1; ④ S1P1; ⑤ P1S1; ⑥ P1P1; ⑦ S1S1; ⑧ P1P2P2P1; ⑨ P1S1; ⑩ P1P2S2P1; ⑪ S1P1; ⑫ P1P2P2S1; ⑬ P1P2S2S1; ⑭ P1S2S2S1; ⑮ S1S2S2S1
-
黄国娇. 2014.各向同性/各向异性介质中多震相走时同时反演方法技术研究[D].西安:长安大学地质工程与测绘学院: 16-19. http://cdmd.cnki.com.cn/Article/CDMD-10710-1014070093.htm Huang G J. 2014. Simultaneous Traveltime Inversion in Complex Isotropic/Anisotropic Medium Using Multi-Phase Arrivals[D]. Xi'an: College of Geological Engineering and Surveying, Chang'an University: 16-19 (in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10710-1014070093.htm
黄国娇, 白超英, 钱卫. 2015.球坐标系下多震相走时三参数同时反演成像[J].地球物理学报, 58 (10): 3627-3638. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201510016.htm Huang G J, Bai C Y, Qian W. 2015. Simultaneous inversion of three model parameters with multiple phases of arrival times in spherical coordinates[J]. Chinese Journal of Geophysics, 58 (10): 3627-3638 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201510016.htm
蒋丽丽, 孙建国. 2008.基于Poisson方程的曲网格生成技术[J].世界地质, 27 (3): 298-305. http://www.cnki.com.cn/Article/CJFDTOTAL-SJDZ200803012.htm Jiang L L, Sun J G. 2008. Source terms of elliptic system in grid generation[J]. Global Geology, 27 (3): 298-305 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-SJDZ200803012.htm
李晓玲, 白超英, 胡光义. 2013.起伏层状TI介质中多次波射线追踪[J].石油地球物理勘探, 48 (6): 924-931. http://www.cnki.com.cn/Article/CJFDTOTAL-SYDQ201306011.htm Li X L, Bai C Y, Hu G Y. 2013. Multiple ray tracing in an undulated layered TI media[J]. Oil Geophysical Prospect-ing, 48 (6): 924-931 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-SYDQ201306011.htm
王涛, 张延腾, 白超英. 2014. 2D各向异性TTI介质中地震波场的射线法和波动方程数值模拟对比研究[J].地球物理学进展, 29 (2): 580-587. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201402014.htm Wang T, Zhang Y T, Bai C Y. 2014. Seismic wavefield propagation on 2D anisotropic TTI media through ray theory and wave-equation simulation and comparison[J]. Progress in Geophysics, 29 (2): 580-587 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201402014.htm
王涛, 白超英, 杨尚倍. 2016.各向同性介质中射线理论与有限差分地震波场模拟方法比较研究[J].地球物理学进展, 31 (2): 606-613. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201602013.htm Wang T, Bai C Y, Yang S B. 2016. Comparison of ray theory and FDM for simulating seismic wavefield in isotropic media[J]. Progress in Geophysics, 31 (2): 606-613 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201602013.htm
王童奎, 张美根, 李小凡, 李瑞华, 龙桂华. 2007. PS转换波界面二次源法射线追踪[J].地球物理学进展, 22 (1): 165-170. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ200701022.htm Wang T K, Zhang M G, Li X F, Li R H, Long G H. 2007. Interface points as secondary sources method raytracing PS converted waves[J]. Progress in Geophysics, 22 (1): 165-170 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ200701022.htm
赵爱华, 张中杰. 2004.三维复杂介质中转换波走时快速计算[J].地球物理学报, 47 (4): 702-707. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200404023.htm Zhao A H, Zhang Z J. 2004. Fast calculation of converted wave traveltime in 3-D complex media[J]. Chinese Journal of Geophysics, 47 (4): 702-707 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200404023.htm
赵爱华, 丁志峰. 2005.宽角反射地震波走时模拟的双重网格法[J].地球物理学报, 48 (5): 1141-1147. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200505022.htm Zhao A H, Ding Z F. 2005. A double-grid algorithm for calculating traveltimes of wide-angle reflection waves[J]. Chinese Journal of Geophysics, 48 (5): 1141-1147 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200505022.htm
朱多林, 白超英. 2011.基于波动方程理论的地震波场数值模拟方法综述[J].地球物理学进展, 26 (5): 1588-1599. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200505022.htm Zhu D L, Bai C Y. 2011. Review on the seismic wavefield forward modelling[J]. Progress in Geophysics, 26 (5): 1588-1599 (in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200505022.htm
Bai C Y, Li X L, Wang Q L, Peng J B. 2012. Multiple arrival tracking within irregular triangular or tetrahedral cell model[J]. J Geophys Eng, 9 (1): 29-38. doi: 10.1088/1742-2132/9/1/004
Carcione J M, Herman G C, Kroode A P E. 2002. Seismic modeling[J]. Geophysics, 67 (4): 1304-1325. doi: 10.1190/1.1500393
Červený V. 2001. Seismic Ray Theory[M]. Cambridge: Cambridge University Press: 117-118.
Crampin S. 1984. Effective anisotropic elastic constants for wave propagation through cracked solids[J]. Geophys J Int, 76 (1): 135-145. doi: 10.1111/j.1365-246X.1984.tb05029.x
Hixon R. 1998. Evaluation of a high-accuracy MacCormack-type scheme using benchmark problems[J]. J Comp Acoust, 6 (3): 291-305. doi: 10.1142/S0218396X9800020X
Hixon R, Turkel E. 2000. Compact implicit MacCormack-type schemes with high accuracy[J]. J Comp Phys, 158 (1): 51-70. doi: 10.1006/jcph.1999.6406
Igel H, Mora P, Riollet B. 1995. Anisotropic wave propagation through finite-difference grids[J]. Geophysics, 60 (4): 1203-1216. doi: 10.1190/1.1443849
Komatitsch D, Barnes C, Tromp J. 2000. Simulation of anisotropic wave propagation based upon a spectral element method[J]. Geophysics, 65 (4): 1251-1260. doi: 10.1190/1.1444816
Moser T J. 1991. Shortest path calculation of seismic rays[J]. Geophysics, 56 (1): 59-67. doi: 10.1190/1.1442958
Tam C K W, Webb J C. 1993. Dispersion-relation-preserving finite difference schemes for computational acoustics[J]. J Comp Phys, 107 (2): 262-281. doi: 10.1006/jcph.1993.1142
Zhang W, Chen X F. 2006. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J]. Geophys J Int, 167 (1): 337-353. doi: 10.1111/gji.2006.167.issue-1
Zhang W, Zhang Z G, Chen X F. 2012. Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids[J]. Geophys J Int, 190 (1): 358-378. doi: 10.1111/gji.2012.190.issue-1
Zhou H, Chen X F. 2008. Localized boundary integral equation-discrete wave number method for simulating wave propagation in irregular multiple layers, part I: Theory[J]. Bull Seismol Soc Am, 99 (3): 1984-1994.
Zhu J L, Dorman J. 2000. Two-dimensional, three-component wave propagation in a transversely isotropic medium with arbitrary orientation-finite-element modeling[J]. Geophysics, 65 (3): 934-942. doi: 10.1190/1.1444789
-
期刊类型引用(3)
1. 高曹珀，王自法，周阳，WANG Jianming. 基于大样本模拟的巨灾损失计算方法. 保险研究. 2023(03): 105-116 . 
百度学术
2. 邵霄怡，王晓青，窦爱霞，袁小祥. 基于蒙特卡洛方法的地震目录模拟及相符性检验——以汾渭地震带为例. 中国地震. 2018(01): 35-47 . 百度学术
3. Shao Xiaoyi，Wang Xiaoqing，Dou Aixia，Yuan Xiaoxiang. Seismic Inventory Simulation and Conformance Test Based on the Monte Carlo Method——A Case Study of the Fenhe-Weihe Belt. Earthquake Research in China. 2018(04): 533-548 . 必应学术
其他类型引用(0)
下载:
计量
- 文章访问数: 675
- HTML全文浏览量: 273
- PDF下载量: 21
- 被引次数: 3
