地球物理学报  2012, Vol. 55 Issue (2): 547-559   PDF    
求解3D弹性波方程的并行WNAD方法及其TI介质中的波场模拟
宋国杰1 , 杨顶辉1 , 童平1 , 练永生2     
1. 清华大学数学科学系, 北京 100084;
2. Department of Mechanical Engineer, University of Louisville, Louisville, KY 610500, US;
3. 清华大学计算机科学与技术系, 北京 100084
摘要: 准确模拟TTI介质中弹性波的传播是研究地震各向异性、AVO反演的基础. 在二维加权近似解析离散化(WNAD)算法的基础上, 本文发展的并行WNAD算法是一种研究三维横向各向同性(TI)介质中弹性波传播的、快速高效的数值模拟方法. 我们首先介绍三维WNAD方法的构造过程, 然后与经典的差分格式——交错网格(SG)算法进行了比较. 理论分析和数值算例表明, WNAD算法比交错网格算法更适合在高性能计算机上进行大规模弹性波场模拟. 同时, 本文利用并行的WNAD方法研究了弹性波在TTI介质中的传播规律, 观测了TI介质中弹性波传播的重要特征:横波分离、体波耦合和速度各向异性等. 在TTI介质分界面处, 弹性波产生更加复杂的折射、反射和波型转化, 使得波场非常复杂, 研究和辨别不同类型的波能够加深我们对由裂隙诱导的各向异性介质的认识.
关键词: 波场模拟      加权的近似解析离散化方法      TTI介质      数值频散      并行计算     
Parallel WNAD algorithm for solving 3D elastic equation and its wavefield simulations in TI media
SONG Guo-Jie1, YANG Ding-Hui1, TONG Ping1, LIAN Yong-Sheng2     
1. Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China;
2. Department of Mechanical Engineering, University of Louisville, Louisville, KY 610500, US;
3. Department of Computer Science and Technology, Tsinghua University, Beijing 100084, China
Abstract: The accurate modeling of elastic wave propagation in TTI media is significantly important for the study of the anisotropy, amplitude versus offset (AVO) analysis, and other inverse problems. A new numerical algorithm named parallel weighted nearly-analytical discrete (WNAD) method is presented in this paper. As an extension of 2D WNAD, the parallel WNAD method is a fast and efficient algorithm for simulating the elastic wave propagation in TTI media. We first show how to formulate the 3D WNAD method, and then compare the wavefield generated by this method with the wavefield computed by a conventional staggered-grid method. The numerical results show that the WNAD algorithm is more suitable for the simulations of the large-scale seismic wavefield by using the high-performance computers. Using the parallel WNAD algorithm, we study the elastic wave propagation in the TTI media and observe the important feature of TI media: shear-wave splitting, the quasi body wave coupling and velocity anisotropy. The reflected, refracted, and converted waves are generated at the interfaces. This makes the wavefield complex. To better understand the anisotropic media induced by fracture, it's useful to study and identify those waves..
Key words: Wave-field simulation      Weighted nearly-analytic discrete method      TTI media      Numerical dispersion      Parallel computing     
1 引言

地震观测表明,由裂隙诱导的各向异性是广泛存在的.如在页岩、周期性薄互层等所谓的横向各向同性(TI)介质中,垂直于裂隙方向传播的地震波波速相对较小,而沿着平行于裂隙方向传播的地震波波速相对较大.主要由沉积作用形成的TI介质的对称轴垂直于水平地表,此时我们称其为具有垂直对称轴的横向各向同性(VTI)介质.在漫长的地质变化过程中,由于褶皱、上冲等地质作用的影响,TI介质的对称轴方向会发生变化,可能不再垂直于水平面,而是与铅垂线方向存在一定的夹角,形成所谓的具有倾斜对称轴的横向各向同性(TTI)介质.目前,国内外许多地学工作者都对TTI介质进行了研究.1986年Thomsen[1]系统地分析了TI介质中的群速度和相速度之间的关系,并给出了Thomsen系数来表示介质的各向异性特征.牛滨华等[2]使用方位矢量波动方程求解六方各向异性介质相速度.吴国忱等[3]研究了TTI介质中弹性波的相速度及偏振特征.郝重涛等[4]研究了任意空间取向TI介质中体波速度的角散和方位变化,杜向东等[5]研究了TI介质偏移的速度建模.

在数值模拟方面,包括Crampin[6],Carcione[7],Igel[8],Juhlin[9]等在内的国际知名地球物理学家都使用数值模拟算法研究了TI介质,尤其是VTI介质中弹性波的传播.在国内,侯安宁等[10]研究了各向异性介质中弹性波高阶差分方法及其稳定性问题,吴国忱[11]研究了三维TTI介质的相速度和群速度,牟永光等[12],郝奇[13]用交错网格方法进行了二维三分量地震波场模拟.

因为有限差分方法编程简单、存储量小、计算效率高,所以有限差分算法是目前研究TI介质使用最多的数值算法.数值频散是有限差分方法的固有现象,能否有效压制数值频散成为判定一种数值模拟方法是否实用的重要判定标准之一.由于数值频散产生的原因在于差分方程和微分方程之间所存在的近似误差,所以采用精细网格进行计算和使用高精度的差分格式是消除数值频散两种基本策略,但是对于目前的软硬件条件来说,这两种策略都存在一定的不足.采用精细网格会使得单位波长内的采样点数大大增加,虽然可以消除数值频散,但是由于目前计算机的存储能力和计算能力限制,这也使得细网格条件下能够计算的区域相对较小,难以应用于大规模地震波场模拟;经典的高阶有限差分格式根据其网格上应力以及位移(或速度)分布位置的不同,又分为规则网格高阶差分方法(如Lax-WendroffCorrection (LWC))[14]和交错网格高阶差分方法(StaggeredGrid (SG))等方法[15].其中规则网格上的高阶差分方法应力位移均分布在整网格点上,每个方向上至少需要2n+1个网格点才能保证差分格式达到2n阶精度;而交错网格方法应力位于半网格点上、位移/速度分布在整网格上,只需要使用2n个网格点即可达到2n阶精度,并且交错网格方法的收敛性和局部性都优于规则网格上的同阶差分算法[12].这也是目前地学领域中广泛使用交错网格方法进行地震波研究的原因之一.

总体来说,这两类高阶差分方法都是通过在空间上每个方向使用较多的网格点来获得高精度,以达到消除数值频散的目的.每个方向上使用较多的网格点增加了边界条件的处理难度,容易降低边界点的计算精度和算法的稳定性,而边界处的低精度和较差的稳定性又会反过来影响全局的计算精度和稳定性,从而导致整个计算结果不能尽如人意.另外在利用并行计算机计算时,较长的空间算子破坏了差分格式良好的局部性,使得不同进程之间的通信量大大增加,严重降低了并行计算效率.

杨顶辉等[16-19]提出的近似解析离散化(NAD)及其改进方法是一类新的扰动差分计算方法.同传统高阶差分方法不同的是,这类方法在近似波位移的高阶空间偏导数时候,不仅使用位移,同时也使用位移的梯度场共同重构波位移的高阶空间偏导数.我们知道,梯度能够反映函数变化趋势,这一重要的波场信息保留在NAD 类算法的空间高阶偏导数逼近表达式中,这使得NAD 算法即使是在粗网格条件下也能有效压制数值频散,极大地提高了波场模拟的计算效率、节省了存储空间.

模拟某一具体问题,我们希望能够在尽可能短的时间内得到模拟结果.要想尽可能地缩短正演模拟时间,一方面需要采用高效的模拟算法,另一方面也需要紧密结合计算机硬件的发展,引入先进的计算工具加快计算速度.目前,我国在高性能计算机硬件发展方面已经取得了举世瞩目的成就,由国防科技大学研制的“天河一号A"已经成为当前世界上计算速度最快的高性能计算机[20].在这样的背景下,掌握并利用大规模计算机进行地球科学研究已经成为新一代地学工作者必备的基础知识之一.目前,国内外地学领域有许多优秀的成果都是借助于高性能计算机的卓越计算性能完成的.早在1995年,Olsen等[21]使用了512 个CPU 进行了网格规模为576×352×116 的三维交错网格并行计算;Bohlen[22]结合MPI实现了三维黏弹性交错网格(精度为空间4阶,时间2阶)并行计算,并给出了并行加速比分析.Komatitsh等[23]使用贝奥武夫机群结合谱元法研究了三维地球模型.高红伟等[24]使用并行化的格子法模拟了三维各向异性介质中弹性波传播,Tape等[25]使用基于伴随方法的地震层析成像,耗时超过80万CPU 机时研究了美国加利福尼亚州的地壳情况.

在本文中,我们首先给出并行化的三维近似解析离散化方法,然后分析它的计算效率和加速比,并使用并行WNAD 方法[1719]来研究弹性波在单层、双层横向各向同性(TI)介质中的传播规律.

2 TI介质中的弹性波方程和WMAD算法 2.1 TI介质中的弹性波方程

三维弹性介质中传播的弹性波方程可以写为

(1)

其中C= {cijkl(xyz)}表示四阶弹性张量,ρ 表示介质密度,ui为第i个方向上的波位移分量,fi为第i个方向上的震源分量.变量上的“·"表示对时间求一阶偏导数,下标中记号“,j"表示对xj求偏导.

对于一般的完全各向异性介质中的弹性波传播方程,四阶弹性张量C含有21个独立参数[26],非常复杂,我们很难分析弹性波在这种介质中的传播规律.同时,地球介质含有较多的定向排列的裂隙和薄层结构,使得在这种介质中传播的地震波具有明显的横向各向同性特征,所以TI介质是目前地球物理领域内研究最多的弹性介质之一.

根据对称轴的不同,TI介质可以分为三类:(1)具有垂直对称轴的横向各向同性(VTI)介质,其对称轴垂直于水平地表;(2)具有水平对称轴的横向各向同性(HTI)介质,其对称轴平行于水平地表;(3)具有倾斜对称轴的横向各向同性(TTI)介质,其对称轴既不垂直于水平地表,也不平行于水平地表,而是与水平地表之间存在一定大小的夹角.

对于TTI介质,如果我们以其对称轴为z1 轴,建立笛卡尔坐标系x1y1z1 (称这个坐标系为自然坐标系,以区别普通的、以铅垂线为z轴建立的观测坐标系xyz),则在自然坐标系下TTI介质变为VTI介质.反之,任何一个TTI介质可由自然坐标系x1y1z1 下的VTI介质经过旋转得到.根据对称性,VTI介质的弹性张量可以写为

其中).

不妨设TTI介质对称轴z1 的走向角(Strikeangle)为φ,倾向角(Dipangle)为θ,则此TTI介质可以通过两次旋转(首先沿y轴旋转θ 角,然后绕z轴旋转$\frac{\pi }{2}-\phi $角)变成VTI介质,对应的弹性张量由下式给出[26]:

其中Cmnpq表示自然坐标系下VTI介质的四阶弹性张量,Cijkl*表示观测坐标系下TTI介质的弹性张量,Ωim矩阵

的对应元素.

特别地,当θ =0,φ = π/2时,TTI介质变化为VTI介质;当θ = π/2,φ = π/2时,TTI介质变化为HTI介质,其弹性张量为

2.2 WNAD算法

为了叙述方便,我们引入记号U= {u1u2u3u1,1u1,2u1,3u2,1u2,2u2,3u3,1u3,2u3,3}T 表示波位移及其梯度分量组成的向量,变量$W=\dot{U}=\partial U/\partial t$称为粒子“速度",变量$P=\dot{W}=\ddot{U}={{\partial }^{2}}U/\partial {{t}^{2}}$称为粒子“加速度".

与Yang等[16]的思路一样,利用4 次截断的Taylor展开式,可以得到由第n时间层上的位移、粒子速度及其梯度表示的第n+1 时间层的位移和粒子速度表达式分别为

(2)

(3)

对于公式(3),考虑到已由公式(2)计算得到了最新的粒子速度Wijkn+1,因此,我们可借鉴求解线性代数方程组Gauss-Seidel方法的类似思想,采用Wijkn+1Wijkn的一个加权组合来替代公式(3)中的Wijkn,从而得到新的位移更新公式[1719]:

(4)

其中ω 为一加权常数,满足0≤ω ≤1,在本文的计算中,我们取ω =0.4.特别地,当ω =1时,即为精化的近似解析离散化方法[18].

要联合使用公式(2)和(4)实现时间方向上的推进计算,显然必须首先计算UWP对时间的高阶偏导数.如果直接离散公式(2)和公式(4)中的高阶时间偏导数,则需要存储多个时间层上的波位移、粒子速度及其梯度值,内存需求量庞大.注意到在公式(2)和(4)中,只含有波位移和粒子速度关于时间的偶数阶偏导数,所以我们可采取Dablain的思想,利用方程(1),首先将时间的高阶偏导数转化为空间的高阶偏导数[1416],然后再利用Yang等给出的高阶空间偏导数的逼近公式[27]来计算公式(2)和(4)中所涉及到的这些高阶空间偏导数.

3 二维分解求解波动方程的差分方法之并行加速比分析

加速比是并行计算中的一个重要概念,它反映了并行程序的执行速度相对于串行程序快了多少倍.本小节将定量地分析并行差分格式的加速比,根据Crama[28]给出的加速比定义:

(5)

其中TS 表示使用最好的串行算法完成模拟所需要的计算时间,TP 表示并行计算时间.在实际的分析中,公式(5)中的TS 很难精确找到能够得到公认的“最好的算法",同时为了方便分析,我们将上式改写为

(6)

其中T1 表示对于某一个具体的算法,仅使用一个进程进行计算所需要的计算时间.根据TST1 的定义知,通常TST1.

下面我们分析公式(6)中的两个因子.对于第一个因子TS/T1,由于TS 表示使用最好的串行算法完成模拟所需要的计算时间,所以TS 对于具体问题是固定不变的,而T1 表示使用1个进程进行计算所需要的时间,不同的计算方法的T1 是不同的.因为单机上面最好的串行程序并不一定适合并行化后在并行机上进行计算,所以我们需要找到合适的串行程序进行并行化.下面我们来分析公式(6)中的第二个因子.

并行计算的墙上时间可以分为两部分:真正用于计算的时间和不同进程之间的通讯时间(进程为了完成计算需要用到分布在其他进程上的数据,从其他进程发送相应的数据并接受到本进程上,此过程在高性能计算中被称为通信,其消耗的时间称为通信时间).首先,我们分析计算时间.假定整体的网格计算规模为N1×N2×N3,划分采用沿x方向和y方向进行二维分解,并假设有p=p1 ×p2 个进程执行任务,则每个进程上的需计算的网格点数为$\frac{{{N}_{1}}{{N}_{2}}{{N}_{3}}}{{{p}_{1}}{{p}_{2}}}$.若每个网格点需要的计算时间为tgrid, 则在一个时间迭代步内,每个进程需要的CPU 计算时间为$\frac{{{N}_{1}}{{N}_{2}}{{N}_{3}}}{{{p}_{1}}{{p}_{2}}}{{t}_{\text{grid}}}$;接下来,我们分析单个时间层内的通信时间.它包括通信建立时间和数据传输时间.若差分格式算子长度为m+1,则x方向需通信m次.每次通信需要建立一对通信和传输计算子区域yz侧面上的一层数据(共计$\frac{{{N}_{2}}}{{{p}_{2}}}{{N}_{3}}$个双精度(单精度)数据类型).假设建立一对通信需要时间tsettle, 每传递一个双精度(或单精度)数据需要时间tcom, 则x方向需要的通信时间为$m\left( {{t}_{\text{settle}}}+\frac{{{N}_{2}}}{{{p}_{2}}}{{N}_{3}}{{t}_{\text{com}}} \right)$.类似地分析可知,y方向需要的通信时间为$m\left( {{t}_{\text{settle}}}+\frac{{{N}_{1}}}{{{p}_{1}}}{{N}_{3}}{{t}_{\text{com}}} \right)$.所以每个进程在每次时间迭代中需要的通信时间(含通信建立时间)为

(7)

由于通信可以并发的进行,所以在不考虑通信阻塞的前提下,全部进程完成通信所需要的时间仍然为Tcom.通信时间公式(7)表明,对于二维分解,随着x方向和y方向划分的块数的增加,每个进程计算子区域的表面积减小,通信量随之下降,通信次数为4m次(注:为分析简便起见,这里不考虑将每个层上所有的数据打包后统一传输,也未考虑通信优化措施).

所以,单个时间层上计算需要的墙上时间=计算时间+通信时间,写成数学表达式即为

对于串行计算,通信时间为0,单个时间层内计算所需要的墙上时间T1 = N1N2N3tgrid.所以

(8)

公式(8)表明,影响加速比因子T1/Tp 的主要因素包括:

(1) 并行计算平台的性能(tsettle, tcom).并行计算平台的性能越好,建立通信、传输数据所需要的时间越短,并行加速比越快.理论上,当tsettle/tgrid =0,tcom/tgrid =0,即建立通信、传输数据不花费时间,此时T1/Tp 达到理想值p1p2

(2) 问题的规模(Ni).在进程数固定的前提下,T1/Tp 随着计算规模的增大单调上升,且渐渐趋于p1p2.换言之,大规模并行计算机适合计算“足够大"规模的问题;

(3) 计算区域的划分.当总的进程数p1p2 固定时,若p1p2 满足p1/p2 = N1/N2(即p1/N1 =p2/N2,在进行区域分解的时候,每个进程计算的三维子区域在xoy面的投影尽可能是正方形),T1/Tp 达到最大;

(4) 使用的算法(算子长度m+1).一个具有较高效率的串行算法并行化后未必能得到一个很好的并行算法,但是一个低效率串行算法并行后一定不会是一个优秀的并行算法.如果一个算法的算子长度很长(每个方向上使用较多的网格点),即使它在串行条件下有很好的效果,也很难保证它的并行化算法仍然是一个优秀的并行算法.因为算子长度越长,其局部性越差,不同进程之间的通信量越大,通信次数也越多,这些都会降低并行加速比.

可见,影响加速比的几个重要因素依次为:计算机性能,问题规模,并行划分,算法的串行计算效率和算法的算子长度.其中计算机的性能是无法改变的,除非将计算机硬件更新换代;问题的规模取决于问题本身,一般也不是我们所能控制的,尤其是目前的地球物理问题对于目前的计算机硬件而言,其规模总是“足够大"的;最佳的剖分要求每个子区域的横截面是正方形的,在实际中未必一定能做到每个子区域的横截面都是正方形的,这时要尽量做到近似正方形,避免出现长、宽比值过大的情况,这样可以提高并行计算速度.所以提高加速比的主要方法是采用好的数值算法和对计算区域进行合理的划分;串行WNAD 方法的计算效率高于一般差分算法[171827];同时,差分算子长度越短,则算法的并行加速比就会越好,在这方面,WNAD 方法在每个方向只使用3个网格点,具有非常好的局部性,这些都说明WNAD 方法具有自然的并行优越性.

下面的小节中,本文将通过比较WNAD 方法和交错网格算法的计算效率、并行加速比,进一步阐明并行WNAD 算法在进行弹性波动方程正演方面的优势.本文的所有计算均在清华大学数学科学系的PC-Cluster上进行.该机群有32 个计算节点,每个节点上拥有8个主频为2.34GHz的核,每个节点上内存为12GB.在我们的计算中,计算区域的分解方式为按x方向和y方向进行二维区域分解,使用Fortran语言编写并行代码.

4 数值算例

在本节中,我们分别设计了5 个计算模型:均匀各向同性模型、均匀的VTI介质模型,双层的TI介质模型和方位角连续变化的非均匀介质模型来检验并行WNAD 算法的有效性.同时,我们将通过研究弹性波在平行于坐标面的三个截面上的波场快照来考察弹性波在TI介质中的传播规律.如无特别说明,xoy截面指的是在深度为介质深度1/2处,平行于xoy平面的截面;xoz截面指的是y分量为介质宽度1/2处,平行于xoz平面的截面;yoz截面指的是x分量为介质长度1/2处,平行于yoz平面的截面.

4.1 模型1 三维均匀各向同性模型

本文第一个模型为三维均匀各向同性模型.我们将分别使用4阶交错网格方法和WNAD 方法计算,通过对比两者的计算结果和计算代价来说明并行WNAD 方法的计算效率.

整个模型的计算区域为0≤xyz≤7.17km, P波和S波速度分别为vp=5.5km/s, vs=3.0km/s, 介质密度为ρ=2.63g/cm3.震源位于计算区域的中心,坐标为(3.57km, 3.57km, 3.57km).震源是一个主频率为f0=30 Hz的Rick 子波,其随时间变化的震源函数为

图 1(a, b, c)给出了在粗网格条件下(Δx= Δyz=30m)使用WNAD 方法计算得到的u3 分量在xoyxozyoz截面上T=0.66s时刻的波场快照.波场快照显示了明显的各向同性特征,P 波和S波的波前以震源为中心形成环形,整个波场非常清晰,无任何可见数值频散.

图 1 u3分量在xoy,xozyyoz截面上T=0.66 s时刻的波场快照 (a, b, c)粗网格条件下WNAD算法得到的u3分量;(d, e, f)粗网格条件下交错网格方法得到的吨分量;(g, h, i)精细网格条件下交错网格方法得到的吨分量. Fig. 1 Thewave-field snapshots of displacement component u3 cm the xoy,xozyand yoz-plane at T= 0.66 s (a, b, c) The wavefield snapshots of u3 with coarse grid by using the WNAD algorithm; (d, e, f) The wavefield snapshots of u3 with coarse grid by using the staggered grid method; (g, h, i) The wavefield snapshots of u3 with refined grid by using the staggered grid method.

图 1(d, e, f)给出了在相同计算条件下,使用交错网格方法计算得到的同时刻u3 分量在三个截面上的波场快照.容易看出,分别由WNAD 算法和交错网格方法计算得到的波场快照基本一致,但是WNAD 格式得到的波场快照非常清晰,而由交错网格方法得到的波场快照中出现了非常明显的数值频散现象.

为了消除交错网格方法引起的数值频散,我们逐渐缩小交错网格的空间步长和时间步长,当空间网格步长减小到Δxyz=10m 时,交错网格方法可以有效消除数值频散,得到清晰的波场快照图 1(g, h, i).可以看出,此时的快照和粗网格条件下由WNAD 方法计算得到的快照图 1(d, e, f)几乎没有差别.说明交错网格方法需要使用更多的网格点数才能获得WNAD 算法在粗网格下一样的结果.

由于WNAD 方法可以在粗网格条件下消除数值频散,所以WNAD 方法可以快速高效地进行地震波场模拟.事实上,得到如图 1(a, b, c)中所示的无明显可见数值频散的波场快照,WNAD 方法只需要5383s, 而在同样的计算环境下,使用加细的交错网格方法计算则需要34347 s, 这说明WNAD 方法的计算效率是交错网格方法的6.38倍;由于网格点数的减少,WNAD 方法所需要的存储量和不同进程之间的通信量也大大减少.本算例中,WNAD方法的存储量只需要加细交错网格方法的11.85%;并行计算时,每步时间层更新时,不同进程之间所需通信量只需要加细交错网格方法的14.81%.这说明,WNAD 方法比交错网格方法更适合在并行计算机上进行大规模地震波场模拟.

下面通过数值算例来比较WNAD 方法和交错网格方法的加速比.为此,我们模拟了区域为0≤xyz≤7.17km 内的弹性波传播.使用不同的进程数计算此问题,得到WNAD 方法和交错网格方法的并行加速比(S)-进程数(p)如图 2.图 2中横轴表示计算使用的进程数,纵轴表示加速比.图 2中的虚线表示理想加速比,实线表示WNAD 方法的加速比,而点划线表示交错网格方法的加速比.一般来说,算法的加速比总是小于理想加速比,越接近理想加速比越说明算法在并行计算时候具有优势.由图 2 可知,随着计算所用的进程数的增加,WNAD 方法和交错网格方法的加速比都在增大,说明对于这两种算法而言,随着参与计算的进程数增大,都会减少完成计算所需要的计算时间.但是WNAD 方法的加速比近似地正比于进程数,比率大约为0.5,这意味着使用WNAD 方法进行波场模拟时,使用n个进程进行并行计算,其计算速度将是同规模问题使用串行WNAD 算法的n/2 倍;而交错网格方法随着进程数的增加,并行加速比增长相对缓慢.这样要提高同样倍数的加速比,交错网格方法需要使用更多的进程参与计算,这说明WNAD 算法可以使用较少的计算硬件实现较大的计算加速,更加适合并行计算.(注:因为交错网格方法需要的内存较大,所以当使用较少的进程计算时,会因为计算机内存不够大而不能开展运算,所以当p较小时,无法得到计算时间,从而无法给出对应的加速比.在图 2中表现为交错网格的加速比曲线从中间(p=16)开始,而不是从p=1开始.)

图 2 WNAD算法(实线)和交错网格算法(点划线)的并行加速比 Fig. 2 The speedup of WNAD algorithm (olid line) and staggered grid (dot dash line) method
4.2 模型2 VTI模型

横波分裂现象是裂缝介质的重要特征之一,本文的第二个数值算例研究弹性波在三维VTI介质中的传播,并观察VTI介质中的横波分离、速度各向异性等.

模型2的计算区域大小为0≤xyz≤7.98km, 介质的弹性参数为:C11=25.2GPa, C33=15.0GPa, C13=6.11GPa, C44=4.38GPa, C66=6.57GPa;密度ρ=2.1g/cm3.震源位于计算区域的中心,随时间变化的震源函数同数值算例1中的震源函数一致.计算参数为Δxyz=20m, Δt=1ms.

图 3给出了T=1.45s时刻三个坐标剖面上的波场快照.从快照图 3可以看到,在xoy截面,三个位移分量表现出明显的各向同性特征:P 波和S波的波前均以震源为中心形成一个圆形.这是因为对于VTI介质而言,xoy截面是一个各向同性面.而在xoz截面和yoz截面,波场表现出明显的各向异性特征:拟P波(qP)和拟S波(qS)沿着不同的方向速度不同,使得波前不再是一个圆形,且波场中出现四个因体波耦合而产生的三叉区(Cusp).

图 3 VTI介质中的弹性波场快照 (a, b, c) ulu2 u3 在:xoy 截面上;(d, e, f) ulu2 u3 在:xoz 截面上;(g, h, i) ulu2 u3 在:yoz 截面上. Fig. 3 The elastic wavefield snapshots in VTI meida (a, b, c) Displacement component ul,u2 and cin xoy-plane; (d, e, f) Displacement component ul,u2 and u3 in xoz-plane; (g, h, i) Displacement component ul,u2 and u3 in yoz-plane.

同时我们注意到:u1 分量在xoz截面、u2 分量在yoz截面内的波场快照出现了横波分离现象,为了更一步验证计算的正确性,我们给出了图 4.在图 4中所示的虚线为u2 分量在yoz截面内的、解析的相速度曲线.可以看出,解析的速度曲线和我们数值模拟得到的波前曲线吻合得非常好,再次说明了WNAD 算法的正确性.

图 4 u2分量在yoz似截面上的波场快照与解析相速度曲线对比图 Fig. 4 Comparison between numerical wavefield and analyticalphase velocity of u2 component in yoz-plane
4.3 模型3 双层模型

双层介质模型是研究多层地球介质问题的基础.CarcioneJ[7]、杨顶辉[18]等均在他们的论文中研究了具有明显速度反差的断层模型.在他们的工作中,由于断层两层由不同的介质构成,所以在断层两侧具有明显的速度反差.本文考虑另外一种双层结构,断层两侧由相同的介质构成,但是介质的对称轴发生改变.在这种情况下,弹性波在断层两侧传播的最大速率和最小速率都没有发生改变,但是最大速度和最小速度的传播方向却发生了改变.接下来的数值算例即是用来研究弹性波在这种仅有裂隙方位角发生变化的双层断层处的反射和折射规律.

图 5给出了模型3的示意图.模型的大小为0≤xyz≤11.98km, 模型的上半部分为VTI介质(z≤6.58km),下半部分为TTI介质(z>6.58km).在自然坐标系下,两层介质具有相同的弹性参数:C11=25.2GPa, C33=15.0GPa, C13=6.11GPa, C44=4.38GPa, C66=6.57GPa;密度ρ=2.1g/cm3.但是,下层介质的对称轴发生倾斜,其θ = π/4,φ =π/6.震源位于(5.98km, 5.98km, 5.78km),震源函数与模型1中使用的震源一致.

图 5 双层模型(模型3),上层为VTI介质,下层为TTI介质.ABCD和EFGH分别表示上下层介质的各向同性面,n1n2分别表示上下层介质的对称轴 Fig. 5 Two-layer model (model 3),the upper layer is VTI and the lower is TTI media.ABCDand EFGH denote the isotropic plane, n1 and n2 denote the symmetry axis of the upper and lowrr media respectively

图 6给出了T=1.62s时弹性波在双层介质中传播的波场快照.从图 6 可以看出,在xoy截面,对比单层的VTI介质中的波场快照图 3图 6xoy截面的快照中直达P波和直达S波基本呈现出各向同性特征,由于P波和S波在层界面处发生反射、折射和转换,所以在快照中可以看到直达的P波和S波后面有能量相对较弱的反射P 波、反射S波以及转换P波和转换S波,而xoy截面恰好是上层VTI介质的各向同性面,所以这些直达波、反射波和转换波均呈现出各向同性特征,整个波场由一系列圆形波阵面构成;在xoz截面和yoz截面,由于上层介质为VTI介质,所以上半部分的速度的最大方向在水平方向;而下层为TTI介质,所以下半部分波速度的最大方向与水平方向有一定的夹角,从而使得波场的外部轮廓(qP波及其折射波的波前面)变得非常不规则.直达拟SH 波(qSH)在介质突变面发生折射、反射,直达qP 波和拟SV 波(qSV)在介质突变面发生反射、折射和转换,使得多种波形混合在一起,能量发生积聚,在qS 波附近的能量较大,同时也导致了我们很难一一辨认波场中的各种波.与速度反差较大的断层面处的传播不同的是,在图 6中,首波成分相对较少且不明显,直达的qP波与它的折射波直接相连,直达的qS波和其折射波也直接相连.

图 6 双层介质模型3中u1u2u3的弹性波场快照(T=1.62 s) (a, b, c)在:xoy截面上;(d, e, f)在:xoz截面上;(g, h, f)在yoz截面上. Fig. 6 The wavefield snapshots of displacement component u1,u2 and u3 in two-layer model 3 (T=1.62 s) (a, b, c) In xoy-plane; (d, e, f) In xoz-plane; (g, h, f) In yoz-plane.
4.4 模型4 对称轴连续变化模型

为了能够更进一步考察弹性波在地下介质对称轴连续变化情形下地震波的传播特征,我们设计了模型4.具体的说,在三维模型4中,任何一个子区域中的弹性介质在自然坐标下都是VTI介质,且弹性参数均为:C11 =25.2 GPa, C33 =15.0 GPa, C13=6.11GPa, C44=4.38GPa, C66=6.57GPa;密度ρ=2.1g/cm3.但是在观测坐标系下,TI介质的对称轴连续变化,对应的倾向角$\theta \text{=}\frac{x}{11.98}\pi $,0≤x≤11.98km, 走向角φ=0.我们可以看到,在地下介质两侧是典型的VTI介质,随着介质对称轴的连续变化,到介质的中央部位地下介质逐步变成HTI介质.

图 7给出了模型4在T=1.62s时刻的波场快照.从快照中我们发现,对称轴的连续变化显著改变了波的传播特征.

图 7 对称轴连续变化模型4中u1u2u3弹性波场快照(T =1.62 s) (a, b, c)在:xoy截面上;(d, e, f)在:xoz截面上;(g, h,f)在yoz截面上. Fig. 7 The wave-field snapshots of displacement component u1,u2 and u31 in the model 4 (T=1.62 s) (a, b, c) In xoy-plane; (d, e, f) In xoz-plane; (g, h, i) In yoz-plane.

xoy截面,qP 波的最大速度方向发生在x方向,且qP波的u3 分量在xoy截面上能量较弱;在三个分量的快照中,均发现了三叉区和横波分离现象,但是图 8中的三叉区更加明显.

图 8 具有倾斜界面的三层介质模型(模型5) Fig. 8 Three-layer model with tilted discontinue surface (model 5)

xoz截面,qP 波在x方向和z方向均达到了最大速度,这使得qP 的波形呈现出近似菱形的形态,qP 波的u2 分量能量较弱;在三个截面上快qS波和慢qS波在各个方向上均发生了明显的横波分离,快qS波和慢qS波的波形构成了一个环形.

yoz截面,由于这个截面恰好是x轴最大坐标的中点,而在我们设计的模型中,此时介质为HTI介质,所以波在此截面上表现出明显的HTI特征:qP波的最大波速方向发生在z方向,最小速度方向发生在y方向,其qP 波的u1 分量能量较弱;快慢qS波在z方向发生横波分离;三个分量上都观测到了三叉区.

4.5 模型5 具有倾斜界面的三层介质模型

因为有限差分方法的网格不能准确刻画倾斜间断面,所以使用有限差分方法模拟地震波场时,在倾斜间断面处容易产生数值频散和人工绕射.模型5(图 8)用来说明并行WNAD 算法能较好地压制弹性波在具有较大速度反差的倾斜间断面附近的人工噪声.图 8所给出的模型分为三个不同的区域,每个区域内的介质速度如表 1所示.从表 1知,不同区域之间最大的层速度反差为2.17.震源位于介质中心,震源函数同模型1.

表 1 模型5参数 Table 1 The parameters of model 5

计算时选择空间步长Δxyz=20m, 时间步长Δt=1 ms, 利用并行WNAD 算法得到T=0.9s时的弹性波在xoz面的波场快照图 9.从图 9中可以看到弹性波在倾斜界面发生发射、折射和波形转换.整个波场中波形比较复杂,但是每个波形都十分清晰,在间断面没有产生明显的人工绕射和数值频散.说明并行WNAD 算法可以有效压制弹性波在倾斜间断面处的人工绕射和数值频散,得到清晰的波场快照.

图 9 模型5中u1(a)、u2(b)、u3(c)在xoz截面上T = 0.9 s时弹性波场快照 Fig. 9 The wave-field snapshots of displacement component u1(a),u2(b) and u3(c) in xoz-plane, respectively, in themodel 5 (T=0.9 s)
5 讨论与结论

本文给出了求解三维弹性波方程的并行加权近似解析离散化(WNAD)算法.这种方法利用位移及其各个梯度分量的加权组合共同近似位移的高阶偏导数.相对于一般高阶差分格式,NAD 类算法保留了更多的波场信息,使得这类算法可以在使用较大的空间网格步长前提下仍然可以得到清晰的、无可见数值频散的波场快照,极大地提高了弹性波场的模拟效率.数值算例表明,在同样的计算条件下,要达到消除数值频散的目的,WNAD 方法的计算速度比交错网格方法快6.38 倍,存储量只需交错网格方法的14.85%;利用并行计算机进行计算时,WNAD 方法的通信量只需交错网格方法的14.81%;同时,WNAD 算法的并行加速比分析表明,WNAD 算法的计算时间正比于进程数,比率约为0.5,说明WNAD 算法可以使用较少的计算资源得到更好的并行加速比,更适合在大规模并行机上进行计算.

使用并行的WNAD 方法,我们模拟了弹性波在TI介质中的传播.对于均匀VTI介质,我们观测到u1 分量在xoz截面,u2 分量在yoz截面的波场快照中可观测到横波分离,且横波分裂方向为裂隙发育方向.在对称轴发生变化的双层介质的波场快照中,在xoy截面,由于该截面是各向同性面,直达波、反射波和转换波表现出明显的各向同性特征,波前为一系列的圆形;在xoz面和yoz面,由于对称轴发生了变化,所以这个截面中的波形也随之发生了变化,将每种波的波前变的不规则,波场更加复杂.最后我们考察了对称轴连续变动的情形,在这种情况下,波场的各向异性特征更加明显,xoz截面上每个方向都观测到了横波分裂.xoy面和yoz面的快照中观测到了更加复杂的三叉区,并且两个截面中的qP 波和qS 波的最大速度方向刚好相反.需要特别说明的是,利用本文提出的并行WNAD 方法模拟大规模地震波场时,数值频散能够得到有效压制,波场非常清晰,数值结果表明并行WNAD 方法是一种十分有效的正演模拟技术.

参考文献
[1] Thomsen L. Weak elastic anisotropy. Geophysics , 1986, 51(10): 1954-1966. DOI:10.1190/1.1442051
[2] 牛滨华, 何樵登, 孙春岩. 六方各向异性介质方位矢量波动方程及其相速度. 石油物探 , 1994, 33(1): 19–29. Niu B H, He Q D, Sun C Y. Azimuth vector wave equation and its phase velocity in hexagonal anisotropy medium. Geophys. Prospect. Petrol. (in Chinese) , 1994, 33(1): 19-29.
[3] 吴国忱, 梁楷, 印兴耀. TTI介质弹性波相速度与偏振特征分析. 地球物理学报 , 2010, 53(8): 1914–1923. Wu G C, Liang K, Yin X Y. The analysis of phase velocity and polarization feature for elastic wave in TTI media. Chinese J. Geophys. (in Chinese) , 2010, 53(8): 1914-1923.
[4] 郝重涛, 姚陈. 任意空间取向TI介质中体波速度特征. 地球物理学报 , 2007, 50(2): 546–555. Hao C T, Yao C. A theoretical study on the effect of ozone below cloud on total ozone retrieved by TOMS and a new inversion algorithm. Chinese J. Geophys. (in Chinese) , 2007, 50(2): 546-555.
[5] 杜向东, 翁斌, 刘军荣, 等. TI介质偏移速度建模研究. 地球物理学报 , 2008, 51(2): 538–545. Du X D, Weng B, Liu J R, et al. Migration velocity modeling strategies of TI media. Chinese J. Geophys. (in Chinese) , 2008, 51(2): 538-545.
[6] Crampin S. An introduction to wave propagation in 2D anisotropic media. Geophys. J. Roy. Astron. Soc. , 1984, 76(1): 17-28. DOI:10.1111/j.1365-246X.1984.tb05018.x
[7] Carcione J M, Kosloff D, Behle A, et al. A spectral scheme for wave propagation simulation in 3-D elastic-anisotropic media. Geophysics , 1992, 57(12): 1593-1607. DOI:10.1190/1.1443227
[8] Igel H, Mora P, Riollet B. Anisotropic wave propagation through finite-difference grids. Geophysics , 1995, 60(4): 1203-1261. DOI:10.1190/1.1443849
[9] Juhlin C. Finite-difference elastic wave propagation in 2D heterogeneous transversely isotropic media. Geophys. Prospect. , 1995, 43(6): 843-858. DOI:10.1111/gpr.1995.43.issue-6
[10] 侯安宁, 何樵登. 各向异性介质中弹性波动高阶差分法及其稳定性的研究. 地球物理学报 , 1995, 38(2): 243–251. Hou A N, He Q D. Study of a elastic wave high-order difference method and its stability in anisotropic media. Chinese J. Geophys. (in Chinese) , 1995, 38(2): 243-251.
[11] 吴国忱. 各向异性介质地震波传播与成像. 东营: 中国石油大学出版社, 2006 . Wu G C. Seismic Wave Propagation and Imaging in Anisotropic Media (in Chinese). Dongying: China University of Petroleum Press, 2006 .
[12] 牟永光, 裴正林. 三维复杂介质地震数值模拟. 北京: 石油工业出版社, 2005 . Mou Y G, Pei Z L. Seismic Numerical Modeling for 3-D Complex Media (in Chinese). (版本). Beijing: Petrolum Industry Press, 2005 .
[13] 郝奇. 三维TTI介质地震波场正演模拟. 长春: 吉林大学, 2006 Hao Q. Seismic wavefield modeling in 3D TTI media (in Chinese). Changchun: Jilin University, 2006.
[14] Dablain M A. The application of high-order differencing to the scale wave equation Laser. Geophysics , 1986, 1(51): 54-66.
[15] Virieux J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics , 1986, 51(4): 889-901. DOI:10.1190/1.1442147
[16] Yang D H, Teng J W, Zhang Z J, et al. A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media. Bull. Seis. Soc. Am. , 2003, 93(2): 882-890. DOI:10.1785/0120020125
[17] 王磊, 杨顶辉, 邓小英. 非均匀介质中地震波应力场的WNAD方法及其数值模拟. 地球物理学报 , 2009, 52(6): 1526–1535. Wang L, Yang D H, Deng X Y. A WNAD method for seismic stress-field modeling in heterogeneous media. Chinese J. Geophys. (in Chinese) , 2009, 52(6): 1526-1535.
[18] Yang D H, Song G J, Hua B L, et al. Simulation of acoustic wavefields in heterogeneous media: A robust method for automatic suppression of numerical dispersion. Geophysics , 2010, 75(3): 99-110. DOI:10.1190/1.3428483
[19] 宋国杰, 杨顶辉, 陈亚丽, 等. 基于WNAD方法的非一致网格算法及其弹性波场模拟. 地球物理学报 , 2010, 53(8): 1985–1992. Song G J, Yang D H, Chen Y L, et al. Non-uniform grid algorithm based on the WNAD method and elastic wave-field simulations. Chinese J. Geophys. (in Chinese) , 2010, 53(8): 1985-1992.
[20] http://www.top500.org/lists/2010/11/press-release.
[21] Olsen K B, Archuleta R J, Matarese J R. 3-dimensional simulation of a magnitude 7.75 earthquake on the San Andreas fault. Sciences , 1995, 270(5242): 1628-1632. DOI:10.1126/science.270.5242.1628
[22] Bohlen T. Parallel 3-D viscoelastic finite difference seismic modelling. Comput. Geosci-UK , 2002, 28(8): 887-899. DOI:10.1016/S0098-3004(02)00006-7
[23] Komatitsch D, Ritsema J, Tromp J. The spectral-element method, Beowulf computing, and global seismology. Science , 2002, 298(5599): 1737-1742. DOI:10.1126/science.1076024
[24] Gao H W, Zhang J F. Parallel 3-D simulation of seismic wave propagation in heterogeneous anisotropic media: a grid method approach. Geophys. J. Int. , 2006, 165(3): 875-888. DOI:10.1111/gji.2006.165.issue-3
[25] Tape C, Liu QY, Maggi A, et al. Adjoint tomography of the Southern California crust. Science , 2009, 325(5943): 988-992. DOI:10.1126/science.1175298
[26] Ting T C T. Anisotropic Elasticity: Theory and Applications. New York: Oxford University Press, 1996 .
[27] Yang D H, Song G J, Lu M. Optimally accurate nearly analytic discrete scheme for wave-field simulation in 3D anisotropic media. Bull. Seis. Soc. Am. , 2007, 97(5): 1557-1569. DOI:10.1785/0120060209
[28] Grama A, Gupta A, Karypis G, et al. Introduction to Parallel Computing. Beijing: China Machine Press, 2003 .