2. 长安大学理学院, 西安 710064
2. School of Science, Chang'an University, Xi'an 710064, China
0 引言
瑞利面波凭借其信噪比高、抗干扰能力强以及在层状介质中的频散特性在地震勘探中应用广泛。现有的常规多道面波分析技术(MASW)基于水平层状假设，已经不能满足横向不均匀复杂介质的要求，而全波形反演技术为其提供了解决方案：不需要根据频散方程计算频散曲线，理论上避免了常规多道面波技术横向分辨率不足的问题[1]。全波形反演技术在实现中需要大量的正演计算：一方面是当前模型波场的正演计算，另一方面是野外数据与模拟数据残差作为震源的波场逆推计算，然后将这两种波场结合求弗雷歇(Frechet)导数进行模型更新[2]。全波形反演利用了所有波至的振幅和相位信息，因此反演过程中必须进行波场正演计算才能得到相应波形的振幅和相位信息。全波形反演的难点之一就是海量计算，主要由多震源的正演模拟导致。由于瑞利波速度较小，模拟瑞利面波在地表附近传播时需要有较小的网格间距才能避免数值频散，这无疑大大增加了正演模拟的计算量。而全波形反演的计算量及计算效率主要集中在反复迭代的正演算法上，因而快速高效的瑞利波并行正演方案设计与实现具有重要的作用和现实意义。
MASW是通过频散函数来进行频散曲线正演计算的，如Thromson-Haskell法[3]、Schwab-Knopoff法等[4]。而波形反演方法则需要从波动理论出发进行波场的正演计算，常见的方法主要有直接法、积分方程法和近似法[5]。直接法包括有限差分法、伪谱法和有限元法等；积分方程法基于惠更斯原理，主要有时间域的积分方程法和边界积分方程法；近似法主要指射线追踪法，基于程函方程，计算效率高。在以上方法中，有限差分法尤其是交错网格有限差分法因其模拟精度高、内存占用小、计算效率高、易于并行实现[6]得到了广泛应用。瑞利面波波动方程与弹性波波动方程一致，不同之处在于地表自由边界条件的处理。Virieux[7]提出弹性波的交错网格有限差分模拟方案，获得了良好的模拟结果。Levander[8]采用空间四阶交错网格法，并在地表应用自由边界条件，为高阶交错网格提供了基础。Berenger[9]提出的完全匹配吸收层(PML)及Francis等[10]提出的卷积完全匹配吸收层(CPML)方法，减少了人工边界的反射。Robersson等[11]对黏弹地震波模拟进行了研究，探讨了广义标准线弹性固体(SLS)模型的计算效果。Zahradnık等[12-13]用真空法处理地表自由边界条件，模拟了瑞利波波场。Robersson[14]将真空法应用到起伏地表黏弹介质的面波模拟中。Rune[15]采用改进的真空法，对地表介质参数进行修正来模拟面波。周竹生等[16]讨论了弹性介质不同模型对应面波记录的频散情况，并讨论了边界和角点的处理。上述学者在波场模拟时采用常规串行算法进行计算，重点在于提高地震波模拟的精度，较少涉及并行算法。
并行方案的发展依赖于计算机技术的提高，较为实用的并行方案出现在20世纪末，其主要标志为消息传递接口(message passing interface, MPI)并行语言的广泛应用。Thomas等[17]提出了基于MPI的三维弹性和黏弹性波动方程模拟方案。Felix等[18]提出了基于多GPU的广义交错网格弹性波模拟法。Liu等[19]提出了基于图形处理器(GPU)并行的不规则自由地表下二维弹性波模拟方案。周洲[20]对比了CPU和异构CPU&GPU结构下三维弹性波波场模拟计算效率。张明财等[21]进行了三维面波模拟讨论，在弹性介质中实现了基于MPI并行的三维瑞利面波模拟，但没有涉及黏弹介质。以上研究主要集中在反射波场的并行计算或三维面波的并行计算，各位学者的研究对波场模拟的快速计算起到了很大的推动作用。从目前瑞利波波形反演技术的发展和应用来看，仍然以二维波形反演为主[1]，因此研究二维瑞利波正演方案具有更大的现实意义。现有的并行计算技术应用广泛，MPI标准以其强大的可移植性、底层封装和不易出错的特点备受青睐。MPI标准有各种不同的版本，如MPICH和LAMMMPI等，本文选择MPICH3.2.1版本实现瑞利波场并行计算。
在面波勘探中，通常以地下介质为完全弹性的假设为前提进行数据处理和解释。而在实际中，地表土壤等松散介质更接于黏弹性介质模型[22]。因此，在正演计算时应当考虑品质因子Q对波场的影响。地震波场模拟中，黏弹性介质模型主要有两种类型：Kelvin模型[22]和SLS模型[14]。本文拟以MPI并行标准为基础，根据Foster[23]的并行方案进行面波并行化设计，合理安排CPU计算和通信之间的关系，对完全弹性模型、Kelvin黏弹性模型和SLS黏弹性模型3种情况进行模拟计算，详细叙述并行中的区域分解、坐标转换、区域通信、波场合并等方案，从而进一步验证并行算法的可行性和有效性。
1 方法原理 1.1 波动方程有限差分算法弹性波波动方程[24]为
(1) 式中：vx, vz分别为质点在x和z方向的速度；σxx, σzz为正应力；σxz为切应力；λ, μ为拉梅系数；ρ为介质密度；ρx, ρz分别为x和z方向的平均密度；μ为剪切模量平均值；t为时间。
Kelvin黏弹性模型波动方程的速度更新公式与弹性波动方程相同，应力计算公式如下：
(2) 其中：
式中：vP、vS分别为纵、横波速度；QP、QS分别为纵、横波品质因子；ω为角频率。
单个SLS黏弹性模型波动方程的速度更新公式与弹性波动方程相同，应力计算公式如下：
(3) 式中：τσ为应力松弛时间；τεP、τεS分别为P波和S波的应变松弛时间；rxx, rxz, rzz为中间变量。rxx, rxz, rzz计算公式如下：
(4) 采用交错网格有限差分模拟，网格形式与Virieux[7]相同，交错网格波场及参数变量分布示意图如图 1所示。速度和应力分量交错分布在网格上，各相邻分量距离为半个网格长度(假定x和z方向网格间距相同)；同时在时间上各分量也是交错的，速度和应力分量差半个时间步长。对方程离散得到对应的离散表达式，差分系数采用Taylor展开得到。
|
| 图 1 交错网格参数变量分布示意图 Fig. 1 Sketch map of parameter variable distribution in a staggered grid |
|
|
对弹性波方程中水平方向的速度进行离散化表示，得到：
(5) 式中：CnN是与差分阶数有关的系数[25]；2N为差分精度；Δt为时间步长；k为离散时间点，t=kΔt；Δx、Δz分别为水平和垂直方向的网格点间距；i、j为离散空间点，x=iΔx，z=jΔz。其他速度和应力分量的离散表达式可用类似方法得到。
对地表采用镜像法进行自由边界处理，令自由地表处垂向应力和切向应力为零，可得到：
(6) 在SLS黏弹性模型中，自由地表的处理还需要rxx、rzz为零。
边界采用CPML[10]。该方法与PML相比不需要将波场分解，计算效率高，且对于掠入射效果较好。CPML的加载如下：
(7) 
(8) 式中：∂x为应力和速度在原始区域内的导数；
从式(7)、(8)可知，在CPML中需要使用一个中间变量φx：先计算x方向不加边界处理的导数，通过式(8)计算得到中间变量，再通过式(7)计算得到经过边界处理的x方向导数。
1.2 基于MPI的并行有限差分计算策略并行算法的核心在于串行算法的并行化，即将可以同时计算的代码分配到不同的进程或CPU中，每个进程同时计算，以此来实现计算效率的加速。由式(5)可知，当空间精度为2阶(N=1)时，点(i, j)处的波场vx只与该点附近4个点的值有关, 即
根据Foster[23]的并行算法设计流程：划分、通信、聚集和映射对波场模拟过程进行讨论。第一步：划分。区域分解方式一般有行分解、列分解和矩形分解，不同的分解方式对应不同的算法流程；常用的是矩形区域分解的方案。第二步：通信。在模拟过程中合理地采用全局通信和局部通信来减少内存和计算开销。第三步：聚集。将相似的工作任务聚集安排给一个处理器，减少处理器之间的通信开销。第四步：映射。将任务和处理器建立连接，保持通信最小化和处理器利用最大化之间的平衡。
然后，将按照Foster的思路对有限差分计算过程进行并行化设计，实现快速准确的波场并行计算。图 2为MPI并行实现的算法流程图。
|
| 图 2 面波MPI并行实现流程图 Fig. 2 Flow chart of MPI parallel implementation of surface wave |
|
|
1) 区域分解
在MPI并行计算中，使用MPI_Comm_rank和MPI_Comm_size两个函数就可确定进程的数目和序号，重要的是将进程与计算区域对应。用两个数组POS(2)和INDEX(4)来表征进程之间的位置关系：POS数组的两个值表示进程对应区域的行和列；INDEX数组的4个值表示该区域上下左右4个区域对应的进程。如图 3所示，采用6个进程进行区域分解，x方向有3个进程，z方向有2个进程，记为Npx×Npz=3×2。
|
| 图 3 区域分解示意图 Fig. 3 Sketch map of region decomposition |
|
|
2) 坐标位置转换
检波器和震源位置，速度和拉梅常数等都是相对于整个模型空间的位置变量，在并行中需要进行坐标转换，将全局坐标转换为对应进程的局部坐标。如模型大小为300×300，网格尺寸为1 m，有两个炮点(5 m, 1 m)和(100 m, 1 m)，炮点纵坐标为zs、横坐标为xs，区域分解形式为Npx×Npz=3×2，则每个矩形区域对应的大小为Ix×Iz=100 m×150 m。进程0对应区域为x=0~99 m、z=0~149 m，进程1对应区域为x=100~199 m、z=0~149 m，以此类推。此时两个炮点的坐标位置转化为(5 m, 1 m)和(1 m, 1 m)，分别表示为xloc和zloc。计算过程如下：
(9) 式中，MOD表示取余数。
介质参数也需要进行类似的坐标转换，得到每个进程对应区域的参数信息。
3) 计算波场
在波场计算时，需要进程间的通信。进程需要将对应区域边界附近的波场值传递给邻近的区域，而传递波场值的数量与差分阶数以及交错网格的定义有关。我们将波场位置按照所在区域划分为3类：①区域内部的点，不需要波场值传递；②区域边界上的点，不包括角点，需要将波场传递给水平或者垂直方向的临近区域；③角点，需要同时传递给水平和垂直方向的临近区域。
现以差分公式空间精度等于4阶为例来说明区域间的波场值传递过程。图 4为进程0区域的右下角：D区域蓝色点代表内部点，不需要传递；A区域红色点代表区域边界上的点，需要传递给右侧区域(进程1)；C区域红色点则需要传递给下侧区域(进程3)；B区域黄色点代表角点，传递给下侧和右侧区域。
|
| 图 4 区域点分类示意图 Fig. 4 Sketch map of regional point classification |
|
|
在计算应力时需要先计算速度vx、vz的垂直方向和水平方向差分。由图 1交错网格的定义可知，在计算正应力时，∂vx/∂x, ∂vz/∂z需要后向差分，而在计算切应力时，∂vx/∂z, ∂vz/∂x需要前向差分。在B和C区域中：vx在z方向是前向差分，所以只需要将边界上一层L1上的vx点值传递给下侧区域，相对的下侧区域要将边界上两层的vx点值传递给进程0区域；vz在z方向是后向差分，所以需要将L1和L2两层的vz点值传递给下侧区域，下侧区域则要将边界的上一层vz点值传递给进程0区域。同理，在B和A区域，vx在x方向是后向差分，需要将边界上W1和W2两层vx的点值传递给右侧区域，相对的右侧区域需要将边界上一层vx点值传递给进程0区域；vz在x方向是前向差分，需要将边界上W1一层的vz点传递给右侧区域，而右侧区域需要将边界上两层的vz点值传递给进程0区域。区域上其他位置的速度场值传递情况可类似得到。
在计算速度时需要先计算∂σxx/∂x、∂σxz/∂x、∂σzz/∂z、∂σxz/∂z。由图 1的定义可知，正应力的差分是前向差分，而切应力的差分是后向差分。如图 4所示，在C和B区域：波场值σzz只需要将边界上一层L1上的点传递给下侧区域，对应的下侧区域需要将边界上两层σzz传递给进程0区域；σxz则需要将边界上L1和L2两层的点传递给下侧区域，对应下侧区域将边界上一层σxz值传递给进程0区域。在A和B区域：波场值σxx只需要将边界上一层W1上的点传递给右侧区域，对应的右侧区域需要将边界上两层σxx传递给进程0区域；σxz则需要将边界上W1和W2两层的点传递给右侧区域，对应右侧区域将边界上一层σxz值传递给进程0区域。区域其他位置处应力传递的情况可类比得到。
根据区域所处的位置判断区域边界是否需要通信。如图 3：进程处于0行，上边界不需要通信，处于0列左边界不需要通信，如进程0；进程处于Npx-1列，右边界不需要通信，处于Npz-1行下边界不需要通信，如进程5。这可以通过数组POS(2)来进行判断。在实际的通信过程中，区域边界需要“加边”来存储交换的波场信息，称为缓冲层。例如模型大小为400×400，区域分解形式为Npx×Npz=2×2，则每个区域的计算量为200×200，但在编程实现过程中每个区域的计算量为[200+2×(N+1)]2。如空间精度为4阶，每个区域的计算量为206×206，区域通信示意图如图 5所示。图 5a中左侧区域将边界上的波场值传递给右侧区域的缓冲层，而右侧区域将边界上的波场信息传递给左侧的缓冲层。图 5b为上下区域之间波场传递示意图。
|
| a.左右区域；b.上下区域。 图 5 区域通信示意图 Fig. 5 Sketch map of regional communication |
|
|
4) 波场输出
波场快照的输出较为简单，在需要输出的时刻将本区域计算的波场直接输出为文件，文件名中包含时刻、炮序和进程位置信息。如第一炮，在第1 000个时刻，进程0输出的vx波场快照文件名可以是：snap_1_1000_vx.0.0。若区域分解形式为Npx×Npz=2×2，则vx在1 000时刻有4个文件，分别是：snap_1_1000_vx.0.0、snap_1_1000_vx.0.1、snap_1_1000_vx.1.0以及snap_1_1000_vx.1.1。当该炮计算完成后，根据区域分解的形式，将4个文件合成为一个完整的文件，就可以得到正确的波场快照。
单炮记录的输出思路如下：设置一个与检波器数目相同的整型数组switch(ntr)，其中ntr表示检波器数目。在坐标转换时，若第i个检波器属于某个区域，则switch(i)等于1；否则switch(i)等于0。如此，对于顶层的区域该数据就成了由0和1组成。假定有5个检波器，位置分别为(5 m, 2 m)、(80 m, 2 m)、(102 m, 2 m)、(164 m, 2 m)、(187 m, 2 m)，网格大小为200×200，区域分解形式为Npx×Npz=2×2，网格尺寸为1 m。经过坐标转换：前两个检波器在进程0区域，且变成了(5 m, 2 m)和(80 m, 2 m)，数组switch为(1, 1, 0, 0, 0)；后3个检波器在进程1区域，且变成了(2 m, 2 m)，(64 m, 2 m)和(87 m, 2 m)，数组switch为(0, 0, 1, 1, 1)。按照一般的方式存储检波器位置处的波场，则会得到两个数组，数据大小分别为2nt和3nt(nt为模拟计算的采样点数)。按照数组switch可以将数据保存在临时数据temp里面，再用MPI_Allreduce函数对该临时数组进行求和归约就可以得到正确的单炮记录，最后按照sgy或者su等格式输出即可。
2 数值算例建立一个两层介质模型，模型参数如下：第一层为0~150 m，vP=2 800 m/s，vS=1 400 m/s，ρ=1.95 g/cm3；第二层为151~299 m，vP=3 000 m/s，vS=1 500 m/s，ρ=2.40 g/cm3。在进行黏弹性介质计算时QP=80，QS=50。
波场模拟时选择的相关参数如下：模型大小为300×300，Δx=Δz=1 m，Δt=0.1 ms；接收长度为0.4 s、震源为30 Hz雷克子波；将炮点放置在(20 m, 2 m)和(280 m, 2 m)两个位置进行计算对比。图 6a、b分别是弹性两层介质模型对应两个单炮记录的垂直分量。由于在地表采用自由边界处理，面波部分的能量较强，而反射部分的能量很弱；这与实际情况相符。黏弹性模型面波单炮记录与之类似，此处不再一一给出。
|
| a.第一炮；b.第二炮。 图 6 弹性介质单炮记录垂直分量 Fig. 6 Vertical component of single shot for elastic medium |
|
|
并行计算时，采用不同的CPU数目和区域分解形式进行算法测试，运行时间如表 1和图 7所示。
| 编号 | CPU数目 | 区域分解形式 | 计算时间/s | ||
| 弹性模型 | Kelvin黏弹性模型 | SLS黏弹性模型 | |||
| 1 | 1 | 1×1 | 253.57 | 484.51 | 502.28 |
| 2 | 2 | 1×2 | 135.38 | 220.32 | 227.66 |
| 3 | 2 | 2×1 | 131.35 | 240.25 | 237.81 |
| 4 | 3 | 1×3 | 98.10 | 147.52 | 156.43 |
| 5 | 3 | 3×1 | 93.95 | 176.43 | 174.02 |
| 6 | 4 | 2×2 | 67.68 | 111.25 | 103.56 |
| 7 | 5 | 1×5 | 59.17 | 89.85 | 84.12 |
| 8 | 5 | 5×1 | 52.73 | 107.30 | 96.57 |
| 9 | 6 | 3×2 | 162.27 | 269.12 | 216.41 |
| 10 | 8 | 4×2 | 212.35 | 325.46 | 261.41 |
| 11 | 10 | 5×2 | 254.20 | 372.90 | 287.01 |
| 12 | 12 | 6×2 | 349.07 | 535.70 | 421.54 |
|
| 图 7 不同并行方案计算时间对比 Fig. 7 Elapsed time comparison of different parallel schemes |
|
|
以弹性模型计算时间为例，可以看到两炮模拟计算的时间随着处理器数目的增加，先逐渐减少，再逐渐增大。并行程序主要分为两部分：必须串行执行的部分以及可以并行执行的部分。增加处理器数目，串行部分的执行时间几乎不变，而并行部分的计算时间会变化。并行部分的计算时间主要由两部分构成：处理器的纯计算时间和处理器之间的通信时间。采用一个处理器模拟计算时，处理器通信时间最少(几乎不需要通信)，计算时间占大多数。当处理器数目增加到通信时间大于计算时间后，继续增加处理器数目，总的计算时间会增加。在图 7中, 对于弹性模型，处理器数目为6之后即为通信时间大于计算时间的情况。当采用5个处理器且区域分解形式为1×5时，计算时间为59.17 s，约为采用一个处理器计算时间253.57 s的23%，约等于处理器数目之比。在表 1中，编号2与3、编号4与5、编号7与8的处理器数目相同，但区域分解形式不同，计算的时间略有不同；这说明区域分解形式对并行的效率有影响。
两种黏弹性介质模型在不同处理器数目情况下，计算时间都大于弹性介质模型的计算时间，这与波动方程的复杂性一致。在处理器数目小于等于5时，两种黏弹性介质模型的计算时间几乎相等；当继续增大处理器数目时，Kelvin黏弹性模型计算时间明显增大，这说明Kelvin的计算效率低于SLS黏弹性模型(图 7)。
为对比不同模型的模拟情况，现对单道记录进行对比分析。为了说明情况，选取第一炮z分量不同偏移距的3道进行对比。图 8a、b、c分别是第40、150、280道3种模型的单道对比图，分别表示了近偏移距、中偏移距和远偏移距3种情况。从图 8中可以看到：在近偏移距时，3种模型记录很贴近，几乎相同；在中偏移距时，弹性模型振幅较大，两种黏弹性模型振幅较小；在远偏移距时，弹性模型振幅最大，SLS次之，Kelvin模型最小，且两种黏弹性模型记录相位略有不同。这与实际情况相同，弹性模型没有考虑吸收衰减和偏移距的影响，因此振幅不随之变化；黏弹性模型考虑了吸收衰减的影响，偏移距越大时振幅越小。
|
| a.第40道；b.第150道；c.第280道。 图 8 垂直分量单道记录对比 Fig. 8 Comparison of vertical component for a single trace |
|
|
1) 利用MPI标准实现了弹性和黏弹面波的并行模拟，采用并行策略可以有效减少计算时间，提高计算效率，为进一步进行瑞利波波形反演提供了高效的计算方案。
2) 与单处理器计算时间相比，增加处理器数目可以明显减少计算时间。但随着处理器数目的增加，不同处理器之间的通信时间也增大。当处理器间用于通信时间大于并行计算的时间后，增加处理器的数目，总的计算时间会逐渐增加。这要求在并行时需要选择合适的处理器数目，盲目地选择较多数目的处理器有时反而会降低计算效率。
3) 实际介质中存在着吸收衰减，黏弹性介质模型更符合实际情况。通过对比两种黏弹性模型的计算时间，SLS黏弹性模型的计算效率优于Kelvin黏弹性模型。
| [1] |
吴华, 李庆春, 邵广周. 瑞利波波形反演的发展现状及展望[J]. 物探与化探, 2018, 42(6): 1103-1111. Wu Hua, Li Qingchun, Shao Guangzhou. Development Status and Prospect of Rayleigh Waveform Inversion[J]. Geophysical and Geochemical Exploration, 2018, 42(6): 1103-1111. |
| [2] |
Tarantola A. Inversion of Seismic Reflection Data in the Acoustic Approximation[J]. Geophysics, 1984(49): 1259-1266. |
| [3] |
Haskell N A. The Dispersion of Surface Waves on Multilayered Media[J]. Bulletin of the Seismological Society of America, 1953, 43(1): 17-34. |
| [4] |
Schwab F. Surface Wave Dispersion Computation:Knopoff's Method[J]. Bulletin of the Seismological Society of America, 1970, 60(5): 1491-1520. |
| [5] |
Carcione J M, Herman J C, Kroode A P E. Seismic Modeling[J]. Geophysics, 2002, 67(4): 1304-1325. DOI:10.1190/1.1500393 |
| [6] |
裴正林, 牟永光. 地震波传播数值模拟[J]. 地球物理学进展, 2004, 19(4): 933-941. Pei Zhenglin, Mou Yongguang. Numerical Simulation of Seismic Wave Propagation[J]. Progress in Geophysics, 2004, 19(4): 933-941. DOI:10.3969/j.issn.1004-2903.2004.04.038 |
| [7] |
Virieux J. P-SV Wave Propagation in Heterogeneous Media:Velocity-Stress Finite-Difference Method[J]. Geophysics, 1986, 51(4): 889-901. DOI:10.1190/1.1442147 |
| [8] |
Levander A R. Fourth-Order Finite-Difference P-SV Seismograms[J]. Geophysics, 1988, 53(11): 1425-1436. DOI:10.1190/1.1442422 |
| [9] |
Berenger J P. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Computational Physics, 1994, 114(2): 185-200. DOI:10.1006/jcph.1994.1159 |
| [10] |
Collino F, Tsogka C. Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media[J]. Geophysics, 2001, 66(1): 294-307. DOI:10.1190/1.1444908 |
| [11] |
Johan O A R, Joakim O B, William W S. Viscoelastic Finite-Difference Modeling[J]. Geophysics, 1994, 59(9): 1444-1456. DOI:10.1190/1.1443701 |
| [12] |
Zahradnık J, Priolo E. Heterogeneous Formulations of Elastodynamic Equations and Finite Difference Schemes[J]. Geophys J Int, 1995, 120(3): 663-676. DOI:10.1111/j.1365-246X.1995.tb01844.x |
| [13] |
Bohlen T, Saenger E H. Accuracy of Heterogeneous Staggered-Grid Finite-Difference Modeling of Rayleigh Waves[J]. Geophysics, 2006, 71(4): T109-T115. DOI:10.1190/1.2213051 |
| [14] |
Robertsson O A J. A Numerical Free-Surface Condition for Elastic/Viscoelastic Finite-Difference Modeling in the Presence of Topography[J]. Geophysics, 1996, 61(6): 1921-1934. DOI:10.1190/1.1444107 |
| [15] |
Rune M. Free-Surface Boundary Conditions for Elastic Staggered-Grid Modeling Schemes[J]. Geophysics, 2002, 67(5): 1616-1623. DOI:10.1190/1.1512752 |
| [16] |
周竹生, 刘喜亮, 熊孝雨. 弹性介质中瑞雷面波有限差分法正演模拟[J]. 地球物理学报, 2007, 50(2): 567-573. Zhou Zhusheng, Liu Xiliang, Xiong Xiaoyu. Finite Difference Modeling of Rayleigh Surface Wave in Elastic Media[J]. Chinese Journal of Geophysics, 2007, 50(2): 567-573. DOI:10.3321/j.issn:0001-5733.2007.02.030 |
| [17] |
Thomas B, Tobias M M, Bernd M. Parallel Finite-Difference Modeling of Seismic Wave Scattering in 3-D Elastic Random Media[C]//SEG Technical Program Expanded Abstracts 2001. Tulsa: SEG, 2001: 1147-1150.
|
| [18] |
Felix R, Mauricio H, Albert F, et al. Generalized Elastic Staggered Grids on Multi-GPU Platforms[C]//SEG Technical Program Expanded Abstracts 2012. Tulsa: SEG, 2012: 1-5.
|
| [19] |
Liu Xiaobo, Chen Jingyi, Lan Haiqiang, et al. Numerical Modeling of Wave Propagation with an Irregular Free Surface and Graphic Processing Unit (GPU) Implementation[C]//SEG Technical Program Expanded Abstracts 2015. Tulsa: SEG, 2015: 3704-3709.
|
| [20] |
周洲. 基于多线程并行计算的有限差分法弹性波数值模拟[J]. 硅谷, 2015, 4: 76. Zhou Zhou. Numerical Simulation of Elastic Waves Using Finite Difference Method Based on Multithread Parallel Computation[J]. Silicon Valley, 2015, 4: 76. |
| [21] |
张明财, 熊章强, 张大洲. 基于MPI的三维瑞雷面波有限差分并行模拟[J]. 石油物探, 2013, 52(4): 354-362. Zhang Mingcai, Xiong Zhangqiang, Zhang Dazhou. Finite Difference Parallel Simulation of 3D Rayleigh Surface Wave Based on MPI[J]. Geophysical Prospecting for Petroleum, 2013, 52(4): 354-362. DOI:10.3969/j.issn.1000-1441.2013.04.004 |
| [22] |
钟飞, 张伟, 焦标强, 等. 可控震源粘弹性波动方程有限差分模拟[J]. 煤田地质与勘探, 2011, 39(2): 57-65. Zhong Fei, Zhang Wei, Jiao Biaoqiang, et al. Finite Difference Simulation of Viscoelastic Wave Equation in Vibroseis[J]. Coal Geology & Exploration, 2011, 39(2): 57-65. DOI:10.3969/j.issn.1001-1986.2011.02.013 |
| [23] |
Foster I. Designing and Building Parallel Programs:Concepts and Tools for Parallel Software Engineering[M]. New Jersey: Addison-Wesley, 1995.
|
| [24] |
杨庆节, 刘财, 耿美霞, 等. 交错网格任意阶导数有限差分格式及差分系数推导[J]. 吉林大学学报(地球科学版), 2014, 44(1): 375-385. Yang Qingjie, Liu Cai, Geng Meixia, et al. Staggered Grid Finite Difference Scheme and Coefficients Deduction of Any Number of Derivatives[J]. Journal of Jilin University (Earth Science Edition), 2014, 44(1): 375-385. |
| [25] |
崔永福, 李国发, 吴国忱, 等. 基于面波模拟和曲波变换的去噪技术[J]. 吉林大学学报(地球科学版), 2016, 46(3): 911-919. Cui Yongfu, Li Guofa, Wu Guochen, et al. Seismic Denoising Technique Based on Surface Wave Modeling and Curvelet Transform[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(3): 911-919. |
| [26] |
Roden J A, Gedney S D. Convolutional PML(CPML):An Efficient FDTD Implementation of the CFS_PML for Arbitrary Media[J]. Microwave and Optical Technology Letters, 2000, 27(5): 334-339. DOI:10.1002/1098-2760(20001205)27:5<334::AID-MOP14>3.0.CO;2-A |
| [27] |
Damir P, Ray M. Convolutional Perfectly Matched Layer for Isotropic and Anisotropic Acoustic Wave Equations[C]//SEG Technical Program Expanded Abstracts 2010. Tulsa: SEG, 2010: 2925-2929.
|