高阶弹性波数值模拟及逆时偏移成像研究

薛 凡，张 智，李学来，李 奇，李喆祥

(桂林理工大学 地球科学学院，广西 桂林 541006)

0 引言

地震数据正演、逆时偏移及全波形反演等方法均需要以高精度的地震波场数值模拟算法为基础。有限差分法因其计算效率高、实现简单而被广泛应用。Alterman和Karal[1]在对均匀介质进行正演模拟时首先使用了有限差分法；Boore[2]用二阶弹性波有限差分法实现了地震波在非均匀介质条件下的正演模拟；Kelly[3]等使用有限差分法制作合成人工地震记录；随着研究的不断发展，Madariaga[4]提出了一种较为先进的交错网格有限差分法；Virieux[5]首先使用交错网格差分法求解一阶速度-应力方程，相对于常规网格解法，精度提高了4倍，成功模拟了转换横波在非均匀介质中的传播；Dablain[6]为解决计算精度与运算效率这一矛盾，提出了声波方程高阶差分解法，其优点是可以采用较大的空间步长来提高计算效率，同时计算精度也得到保障；董良国[7]将时间三阶导数通过一阶速度-应力方程转化为空间导数，解决了在不增加储存空间的前提下提高计算精度的问题；印兴耀[8]等提出了旋转交错网格与紧致有限差分相结合的方法，并基于模拟退火算法进行全局优化，有效地压制数值频散，提高了模拟精度。

逆时偏移的思想是1983年Whitmore[9]在SEG年会上提出的，利用逆时偏移结果拾取层位更新速度模型，通过迭代求取速度；同年McMechan[10]采用时间和空间域二阶精度有限差分法求解全声波方程，实现了逆时偏移；Chang[11-12]等先后采用有限差分实现了声波、弹性波的叠前逆时偏移；张会星[13]与杜启振[14]利用高阶有限差分法实现标量和矢量波场的叠前逆时偏移；刘伊克[15]提出利用地震多次波对盐丘下部成像的新方法，为地震成像研究开拓了新方向；杨江蜂[16]等对碳酸盐岩缝-洞储层进行叠前逆时偏移成像研究。逆时偏移是基于全波动方程，可以处理多次波现象，且不存在倾角限制，是对复杂介质成像非常理想的偏移方法。

本文在前人研究成果的基础上，选取不同地质构造模型，对正演模拟中的各模型进行逆时偏移成像，模型试算结果表明，交错网格有限差分法可以准确的模拟地震波在介质中的传播，分析偏移成像结果，验证了逆时偏移成像的准确性和有效性。

1 正演模拟

1.1 基本方程

本文研究的介质为各向同性介质，选用弹性波本构方程：

(1)

将柯西方程带入公式(1)中，与运动平衡微分方程联立，再将Vx=∂ux/∂x、Vz=∂uz/∂z代入，可以得到二维情况下的各向同性介质弹性波一阶速度-应力方程[公式(2)]：

式中：Vx、Vz表示质点速度，σxx、σzz、σxz表示应力向量，ρ为密度，λ、μ为拉梅常数。

(2)

同纵、横波速度之间的关系：

1.2 离散公式

Vx、Vz、σxx、σzz、σxz的离散值分别以U、V、R、T、H代表，i、j分别表示x、z方向的网格节点，k表示时间节点。通过对时间导数采用二阶近似，得到弹性波一阶速度-应力的交错网格差分格式(即离散形式)[公式(3)]：

(3)

在设置交错网格时，将介质的波场条件和相应的波场参数定义在图1的主体网格(图1a)和交错网格(图1b)上。表1网格点的数量表示波场的弹性分量和弹性参数在网格中的空间位置。

1.3 PML边界条件

大量研究者将Berenger[17]提出的完全匹配层(PML)边界条件应用到地震学研究中[18]。在二维情况下，PML吸收边界条件的主要思想是将边界中的未知量分解成垂直方向和平行界面方向两个部分，在不同的边界匹配层按照不同的方向进行衰减。用方程表示：

图1 主体网格(a)与交错网格设计图(b)Fig.1 Main grid design map(a)andstaggered grid design map(b)

表1 弹性波场分量及弹性参数空间位置Table 1 Spatial location of elastic wave fieldcomponent and elastic parameter

(4)

其中：衰减因子为

(5)

式中：Vp为最大纵波速度，δ为匹配层宽度，R为理想的边界反射系数(一般取10-6)，当dx、dz不等于0时表示衰减，等于0时表示不衰减。创建完全匹配层吸收边界条件基本的做法是在计算区域四周加入匹配层，图2b中左右两侧匹配层4、6区域中在x方向衰减，在z方向不衰减，表示为dx≠0，dz=0；在上下匹配层的2、8区域中在z方向衰减，x方向不衰减，表示为dx=0，dz≠0；在四周1、3、7、9区域内在x方向和z方向都衰减，表示为dx≠0，dz≠0。

图2 衰减边界(a)和PML边界设计方法(b)Fig.2 Attenuation boundary design method(a)andPML boundary design method(b)

为了验证PML边界条件的吸收效果和适用性，建立大小为300 m×300 m的均匀介质模型，模型参数设置为Vp=3000 m/s，Vs=2500 m/s，ρ=2000 kg/m2。网格间距Δx=Δz=1 m，时间采样间隔为Δt=0.1 ms，震源采用雷克子波，震源主频为35 Hz，震源位置在(150 m，150 m)处。设定PML边界吸收宽度为20 m。

由图3可见，采用PML边界条件对均匀介质的数值模拟效果良好，且边界处无明显反射，适用于本文数值模拟方法。

2 叠前逆时偏移成像

2.1 基本原理

假设地层内存在一个反射界面，在界面位置下行波到达时间和上行波到达时间是一致的，在偏移时，将检波器接收到的叠前或叠后地震资料，以地下真实速度沿着逆时间轴方向，从记录道的末尾时刻延拓到零时刻，在空间波场的每一时刻的波场快照中，选择出适用于成像条件的信息，再将这些波场信息叠加，得到反射界面的偏移成像结果。叠前偏移处理是将震源正演模拟和检波器逆时外推同时进行，当两者处在同一时刻(TS=TR)得到成像结果(图4)。

图3 均匀介质模型在t=0.08 s时刻的波场快照

图4 时间一致性原理Fig.4 Time consistency principle

2.2 接收波场的逆时延拓

根据Whitmore提出的逆时偏移基本物理思想：若函数g(x，z，t)是波动方程的一个解，那么函数g(x，z，t-Δt)也是波动方程的一个解。意味着波动现象可以按照时间前进的方向来观察事件发生的过程和最终结果。同理反之，沿着时间倒退的方向来观察这一事件发生的过程。基于这种思想，弹性波逆时偏移可以视为弹性波正演的逆过程。接收波场的逆时延拓是地震波逆时传播的过程，本质上是正演模拟的逆时计算，求解波动方程边值问题。基于二维各向同性介质弹性波方程进行全波动方程逆时偏移，表示为

(6)

T代表的是接收点记录的最大时间，Rx、Rz分别代表接收点记录的x分量和z分量，xR、zR分别代表接收点的位置。实现过程可以解释如下：逆时延拓是从检波器所接收到的最终时间t=T开始直至t=0，逆时延拓就是随着设置的时间步长进行运算，每一次计算就代表着将该时刻接收点的波场值插入到接收点的位置上，将该值作为接收点位置的波场值，即完成这一时刻的计算，当计算完整个时间时，也就完成了逆时延拓过程。

2.3 互相关成像条件

成像结果的质量对后续解释有很大的影响，而成像结果的质量取决于成像条件的好坏，因此对成像条件的研究是逆时偏移方法中最重要的内容。自Claerbout[19]提出时间一致性原理后，经过不断地发展，目前叠前偏移成像所应用的成像条件主要分为三种类型：互相关成像条件、激发时间成像条件和上下行波振幅比成像条件。本文主要使用互相关成像条件。

震源正演模拟波场和接收点逆时延拓波场零延迟互相关是实现时间一致性成像原理基本的方法，用方程可以表示为

(7)

式中：S为震源波场，R为接收波场，T为地震记录的最大时间。对弹性波场而言，就是震源波场和接收波场的x分量或z分量。互相关成像振幅的单位是震源波场或接收波场振幅的平方，其强度取决于震源强度，没有物理上的意义。

3 复杂模型

本文基于Matlab语言编程建立几种典型的介质模型，对地震波场进行数值模拟计算和逆时偏移成像处理。模型大小取300 m×300 m，网格间距为Δx=Δz=1 m，采样间隔为Δt=0.1 ms，采样时间为3.0 s。采用雷克子波作为震源子波，震源主频为35 Hz，设定的PML边界吸收宽度为20 m，检波器均匀分布在z=125 m处。所有成像只取x方向(水平分量)结果。

3.1 单一裂缝模型

建立单一裂缝模型(图5a)，模型参数见表2。介质分为三层，第一层深度为160 m，第二层介质厚40 m，在第二层介质中间存在一个裂缝，裂缝大小为10 m×40 m。参数与第三层介质相同。图5b为单一裂缝模型的地震记录图，震源位置在(150 m，20 m)处，图中可见直达波，分界面反射P、PS波，透射P、PS波，由裂缝边缘的地形尖灭点产生的绕射P、PS波及多次反射波等；图5c为含单一裂缝的地质模型在地震记录中去除直达波和分界面反射波后的地震记录图；图5d为多炮逆时偏移成像结果；图5e为散射波成像结果。

分析可知，偏移结果对地层和裂缝有较好的还原，但是可以看到下层部位有低频干扰。分析散射波成像结果，在去除反射波之后，裂缝形态更加明显，与模型基本吻合，并且逆时偏移成像的位置清晰准确，偏移取得了很好的结果。

图5 单一裂缝模型正演记录及其偏移结果

表2 单一裂缝模型参数设置Table 2 Parameter setting of single fracture model

3.2 凹陷地质模型

建立凹陷模型(图6a)，参数见表3。第一层深度为100 m，第二层厚度为75 m，为了研究地震波在凹陷地质模型中的传播，第三层介质中嵌入一个凹陷模型，凹陷区长度从x=100 m到x=200 m，深度从z=175 m到z=225 m。介质参数与第二层介质相同。图6b为凹陷地质模型震源位置在(150 m，20 m)处时的地震记录图，根据第二分界面的上行反射波和凹陷区底层界面的反射波，可以判断出凹陷区的位置和大体范围；图6c为凹陷地质模型去除直达波和分界面反射波后的地震记录图；图6d为多炮逆时偏移成像结果；图6e为散射波成像结果。

偏移结果能够准确的反映出凹陷地形的分界面，与给定的凹陷地质模型完全一致，多炮逆时偏移成像对地层分界面有较好的成像，凹陷区分界线明显，但也存在低频噪声的干扰；散射波成像效果较好，不受反射波的干扰，对凹陷区域能够准确还原，位置与模型吻合。

3.3 断层地质模型

建立断层地质模型(图7a)，模型参数见表4。断层起始位置为(100 m，100 m)，终止位置为(200 m，200 m)，断层面倾角为45°。断层上、下盘中填充物质与最下方介质参数相同，宽度为5 m。图7b为断层地质模型震源位置在(150 m，20 m)处时的地震记录，图中包含直达波、断层上、下界面反射P、PS，透射P、PS波和尖灭点产生的绕射P、PS波等。分析各波形分布，可以看出上盘上界面反射波和绕射波混合且能量较小，不易分辨。图7c为去除直达波和压制反射波后的地震记录图；图7d为多炮逆时偏移成像结果；图7e为散射波成像结果。

各偏移结果对断层上、下盘都有很好的还原，对断层模型绕射点(地形尖灭点)也可以还原，但对下盘下界面的偏移效果不明显，存在低频噪声的干扰；散射波成像结果对绕射点还原得更好，对断层界面的还原也有一定的效果。与给定的断层模型位置基本吻合。

4 结语

高阶弹性波动方程正演数值模拟和偏移成像一直是地球物理学界的一个重要研究领域。既可以作为一种认知手段来提高研究者对各种复杂介质中地震波传播规律的认识，也可以作为一项检测工具，检验各种正演和偏移算法的应用效果。本文从几种速度模型入手，通过对各模型成像结果分析并结合前人的研究成果，从中得到以下认识：

图6 凹陷模型正演记录及偏移结果

图7 断层地质模型正演记录及其偏移结果Fig.7 Forward recording and migration results of geological fault model

表3 凹陷地质模型参数设置Table 3 Parameter setting of geological depression model

表4 断层地质模型参数设置Table 4 Parameter setting of geological fault model

1)交错网格高阶有限差分法是种计算精度较高的数值模拟方法，采用这种方法得到的地震波波场信息丰富，对研究地震波在地下介质中的传播规律和波场特征十分有利。其缺点是难以克服频散效应，要解决这一问题，必须加密计算网格，但这会导致计算量的增加，使效率下降；另外，当地表起伏较大或地质构造复杂时，有限差分法难以确保数值模拟精度，适应性和灵活性较差。

2)通过对几种绕射点模型进行叠前逆时偏移成像处理，验证了本文的叠前逆时偏移方法是切实有效的，该方法可以对理论模型的反射点、绕射点进行准确成像。对不同的模型进行正演模拟，再进行叠前逆时偏移成像，可以对整个正反演的过程有一个系统的了解。利用叠前逆时偏移方法结合速度分析以及其他的反演方法，可以解决对更复杂的介质模型成像，进而对实际工作中遇到的地质问题提供参考依据，为解决地球物理反演多解性的问题提供有效的手段。