第29卷第3期2017年6月
岩性油气藏
LITHOLOGIC RESERVOIRS
V ol.29No.3Jun.2017
收稿日期：2016-09-12；修回日期：2017-01-03基金项目：国家重点基础研究发展计划（973）项目“纵横波品质因子A VO 反演方法研究”（编号：41604108）资助作者简介：邓帅（1986-），男，中国地质大学（北京）在读博士研究生，研究方向为地震资料高保真处理技术、GPU 高性能计算等。地址：
（100083）北京市海淀区学院路29号中国地质大学（北京）。Email ：*******************。
文章编号：1673-8926（2017）03-0118-08
DOI ：10.3969/j.issn.1673-8926.2017.03.014
上覆水平界面对目的层地震波振幅的影响
邓帅，刘学伟，王祥春
（中国地质大学（北京）地球物理与信息技术学院，北京100083）
摘要：地震波振幅属性是识别目的层岩性的重要参数，但该属性受上覆界面的影响较大。为了弄清上覆水平界面对目的层地震波振幅的影响，从正演模型出发，利用双程波动方程进行模拟，研究地震波在不同速度、不同尺度的上覆地层中传播时振幅的变化，建立上覆界面与目的层地震波振幅之间的联系，通过构建包含完美匹配层（PML ）边界衰减系数的高阶交错网格有限差分公式，使用GPU 进行并行计算加速，最终获得三维层状水平速度模型目的层地震波振幅的变化规律。结果表明，当上覆水平界面的反射系数发生变化时，目的层地震波成像振幅的变化范围及变化趋势都会随之发生相应的改变，尤其当上覆水平界面的反射系数趋于一致时，目的层地震波成像振幅便会发生异常性的变化。该研究成果可为地震资料处理解释过程中正确识别目的层岩性提供一定的理论帮助。关键词：双程波动方程；正演模拟；反射系数；振幅属性；GPU 中图分类号：P315.0
文献标志码：A
Influence of overlying horizontal stratum on ismic amplitude of target zone
DENG Shuai ，LIU Xuewei ，WANG Xiangchun
（School of Geophysics and Information Technology ，China University of Geosciences ，Beijing 100083，China ）
Abstract ：Seismic amplitude is an important parameter to identify the lithology of target zone ，however ，it is
greatly affected by overlying stratum.Bad on the forward modeling ，the numerical simulation of two-way wave equation was applied to study the amplitude change of ismic wave when it travels through the overlying stra-tum with different velocity and scale ，and the relationship between the overlying stratum structure and ismic amplitude of the target zone was discusd.By building high ordered staggered-grid finite difference algorithm containing PML （perfect matching layer ）attenuation coefficient and using GPU parallel acceleration ，the ampli-tude change regularity of the target zone of 3D layered velocity model was defined.The experiment results show that the amplitude of the target zone will change accordingly with the change of reflection coefficient of the over-lying stratum ，and it will change abnormally especially when the reflection coefficient of overlying stratum tends to be consistent.This study results are favorable for correct lithology identification by ismic data.Key words ：two-way wave equation ；forward modeling ；reflection coefficient ；amplitude ；graphic process unit
引用：邓帅，刘学伟，王祥春.上覆水平界面对目的层地震波振幅的影响，2017，29（3）：118-125.Cite ：DENG S ，LIU X W ，WANG X C.Influence of overlying horizontal stratum on ismic a
mplitude of target zone.Lithologic
Rervoirs ，2017，29（3）：118-125.
2017年邓帅等：上覆水平界面对目的层地震波振幅的影响119
0引言
振幅属性是地震数据处理解释过程中的一个重要参数［1］，但是，受上覆地层的影响，从偏移成像剖面中提取的目的层地震波振幅信息，不一定能表征目的层的岩性变化，而可能是上覆地层岩性、构造等要素的反映。
一方面，通过改进地震勘探仪器的性能，地震反射信号的记录精度得到不断提高，为叠前反演和储层预测提供了更高质量的地震资料［2-3］。另一方面，随着勘探领域从构造圈闭到岩性地层圈闭的逐渐延伸，国内外众多学者［4-7］提出了不同的保幅方法来尽可能地保持地震波振幅与目的层界面反射系数之间合理的比例关系，也有学者［8-10］从不同角度评价了这些方法的优劣性。整体而言，这些方法的研究重点大多是如何消除上覆界面的影响以得到真实的反射振幅，但是，受众多因素的制约，现有的保幅方法尚不能做到完全消除上覆界面的影响，业内对上覆界面和目的地层振幅之间的关系仍缺乏研究。
笔者通过实验模拟了地震波在三维倾斜界面速度模型中的传播过程，分析了倾斜界面及其上覆圈闭对目的层地震波振幅的影响，发现受上覆圈闭的影响无法得到目的层岩性的真实反射信息，而目的层的振幅变化同上覆界面之间具有一定的规律
性［11］。针对这一现象，本次研究从正演模型出发，利用波动方程进行数值模拟，获取目的层的反射振幅，并通过研究上覆水平界面反射系数与目的层地震波振幅之间的关系，确定目的层地震波振幅随上覆界面反射系数变化而产生的变化规律，以期为利用地震资料正确识别目的层岩性提供理论依据。
1研究方法
通过模型正演，利用双程波动方程模拟目的层地震波振幅的分布，由于双程波动差分方程数值模拟的计算速度较慢，因此采用GPU 并行加速来提高差分方程的计算效率。1.1高阶交错网格有限差分
为了简化模型、方便计算，本次研究使用声波波动方程来模拟地震波传播。相比于单程波，双程波包含的波场信息更丰富，其模拟的地震波场更接
近于野外实际观测到的地震波场，因此选取双程声波方程来进行地震波场的数值模拟计算。前人［12-15］针对声波方程交错网格有限差分算法的研究已经相当成熟，在实际地震资料成像处理中已经获得
了很好的应用效果。公式（1）为波动方程时间二阶、空间十六阶的三维交错网格有限差分公式，其中，空间维度的阶数越高，计算的精度越高，但相应的计算效率越低。
ìíî
ïïï
ïïïïï
ïïï
ïïïïïv t +1/2x ()i +1/2,j ,k =v t -1/2x ()i +1/2,j ,k -Δt
Δx ∑n =18
C ()8n []p t i +n ,j ,k -p t i -n +1,j ,k v t +1/2y ()i ,j +1/2,k =v t -1/2y ()i ,j +1/2,k -Δt Δy ∑n =18
C ()8n []p t i ,j +n ,k -p t i ,j -n +1,k v t +1/2
z ()i ,j ,k +1/2=v t -1/2z ()
i ,j ,k +1/2-Δt Δz ∑n =18
C ()8n []p t i ,j ,k +n -p t i ,j ,k -n +1p t +1x ()i ,j ,k =p t
x ()i ,j ,k -Δt Δx v 2p ∑n =1
8
C ()8n [
]v t +1/2x []i +()2n -1/2,j ,k -v t +1/2x []i -()2n -1/2,j ,k p t +1y ()i ,j ,k =p t
y ()
i ,j ,k -Δt Δy v 2p ∑n =18
C ()8n []v t +1/2y []i ,j +()2n -1/2,k -v t +1/2y []i ,j -()2n -1/2,k p t +1z ()i ,j ,k =p t z ()i ,j ,k -Δt Δz v 2p ∑n =18
C ()
8n []
v t +1/2z []i ,j ,k +()2n -1/2-v t +1/2z []i ,j ,k -()2n -1/2p t +1i ,j ,k =p t +1x ()i ,j ,k +p t +1y ()i ,j ,k +p t +1z ()i ,j ,k
（1）
式中：v ()x,y ,z 为速度在3个维度上的分量，m /s ；v t
()
i ,j ,k 为t 时刻网格点（i ，j ，k ）处的速度值，m /s ；p ()x,y ,z 为应
力在3个维度上的分量，MPa ；p
t
()
i ,j ,k 为t 时刻网格
点（i ，j ，k ）处的应力值，MPa ；v p 为声波速度，m /s ；Δt 为时间步长，s ；
Δx ，Δy ，Δz 分别为3个维度上网格步长，m ；C 为对应的网格参数。式（1）的稳定性条件为
Δt·V max
≤1∑n =1
N ||C n （2）
式中：V max 为速度最大值，m /s ；N 为空间阶数的一
半，此处N =8（空间十六阶）。
在使用波动方程模拟地震记录时，模型的四周会产生边界反射，因此需要人为添加一定厚度的吸收边界来消除边界反射。完美匹配层（PML ）吸收
120岩性油气藏第29卷第3期
边界［16-19］是目前业内应用最广泛且吸收效果比较好的边界条件，具体做法是将衰减系数d 引入到式（1）中，得到对应的含PML 衰减系数的三维交错网
格有限差分公式
ìíî
ïïïï
ïïïïïï
ïïïïïïïïïï
ïïïïïï
ïïïïïïv t +1/2x ()i +1/2,j ,k =2-d x Δt 2+d x Δt v t -1/2x ()i +1/2,j ,k -2Δt 2+d x Δt 1Δx ∑n =1
8
C ()8n []p t i +n ,j ,k -p t i -n +1,j ,k v t +1/2
y ()i ,j +1/2,k =2-d y Δt 2+d y Δt v t -1/2y ()i ,j +1/2,k -2Δt 2+d y Δt 1Δy ∑n =18
C ()8n []p t i ,j +n ,k -p t i ,j -n +1,k v t +1/2z ()i ,j ,k +1/2=2-d z Δt 2+d z Δt v t -1/2z ()i ,j ,k +1/2-2Δt 2+d z Δt 1Δz ∑n =18
C ()8n []p t i ,j ,k +n -p t i ,j ,k -n +1p t +1x ()i ,j ,k =2-d x Δt 2+d x Δt p t x ()i ,j ,k -2Δt 2+d x Δt v 2p Δx ∑n =18C ()8n [
]v t +1/2x []i +()2n -1/2,j ,k -v t +1/2x []i -()2n -1/2,j ,k p t +1y ()i ,j ,k =2-d y Δt 2+d y Δt p t y ()i ,j ,k -2Δt 2+d y Δt v 2p Δy ∑n =18
C ()8n []v t +1/2y []i ,j +()2n -1/2,k -v t +1/2y []i ,j -()2n -1/2,k p t +1z ()i ,j ,k =2-d z Δt 2+d z Δt p t z ()i ,j ,k -2Δt 2+d z Δt v 2p Δz ∑n =18
C ()8n [
]
v t +1/2z []i ,j ,k +()2n -1/2-v t +1/2z []i ,j ,k -()2n -1/2p t +1i ,j ,k =p t +1x ()i ,j ,k +p t +1y ()i ,j ,k +p t +1z ()i ,j ,k
（3）衰减系数d 在非边界吸收区域时其值等于0，此时式（3）退化为式（1）；在边界吸收区域内，衰减系数
d 的取值由式（4）计算得到
d =2πA F æèöøi δ2
（4）
式中：δ为吸收边界的厚度，m ；A 为一个常系数；F 为子波主频，Hz 。1.2GPU 并行计算GPU 的硬件架构中包括较多的计算核心，且其逻辑判断单元较少，这样的架构决定了GPU 适合于进行高密集、低复杂度的计算。交错网格有限差分算法中包含大量的重复迭代计算，而且采用式（3）将边界和核心计算区域统一起来可以避免在PML 边界位置进行的大量逻辑判断，因此选用GPU 并行计算来实现交错网格有限差分算法。式（3）中每个维度的分量在时间域内计算时都需要前一个时刻的值参与，在空间域内计算时则需要周围点在该时刻的值参与。在进行任务划分时，只能使用串行代码执行时间域内的计算，而将空间
airasia com域的运算放到GPU 上去执行。但是，空间高阶交错网格有限差分方程计算时需要大量的显存读写，对于式（3）中三维十六阶差分方程来说，每计算一个网格点的值都需要读取网格点周围48个网格点
的数据，因此需要借助GPU 片内高速的共享存储器来降低显存读取冗余度［20-22］。图1是对一个三维均匀介质模型分别沿x ，y ，z
等3个方向使用GPU 模拟得到的波场快照。该模型为一个正立方体，每个方向的宽度均为1000m ，激
发点被放置在模型正中心点位置。图1（a ）为x =500m 时yz 平面的波场快照，图1（b ）为y =500m 时xz 平面的波场快照，图1（c ）为z =500m 时xy 平面的波场快照。3个方向的波场快照结果几乎一样，说明GPU 计算结果正确可靠，同时，波场快照
也揭示出PML 边界的吸收效果显著。此外，运行速度上，使用GPU 并行计算比传统CPU 数值模拟提速大约30倍。
图1GPU 模拟三维均匀介质模型x ，y ，z 等3个方向的波场快照
Fig.1Wave field snapshot in three directions of homogeneous medium model calculated by GPU
太原网络营销
（c ）
（a ）（b ）
2017年
邓帅等：上覆水平界面对目的层地震波振幅的影响121
1.3实验模型参数
三维实验模型基本参数如表1所列，各参数均满足差分方程稳定性条件［式（4）］。震源使用雷克子波，激发时，炮点放置在地表位置，接收点放置在地下深度为2580m 处的目的层上。
表1三维实验模型基本参数
Table 1Parameters of 3D experimental model
三维实验模型中设计有一速度异常体，该速度异常体为一厚度（540m ）恒定的圆柱体，如图2所示。实验时，只改变异常体所在蓝色区域的速度V 和半径R 的大小，不改变异常体厚度及围岩参数，围岩速度设计及俯视效果如图2（a ）、图2（b ）所示。
图2三维实验模型
Fig.2
规格英文
Diagrams of 3D experimental model
1.4目的层地震波振幅提取
在本模拟实验中，将激发点放置在地表处，接收点放置在地下目的层之上，模拟二分之一自激自
收时间的地震波场分布。激发时不选择单炮依次激发，而选择所有激发点同时激发，将全部点震源合并为一个面震源，GPU 数值模拟程序会自动保存检波点位置的波场值，随后再提取每个检波点记录长度内的最大振幅值并将其输出，便可得到目的层地震波振幅分布。对于图2所示的实验模型，地震波从地表传播到目的层的最长时间约为2550ms ，所以在数值模拟时，应将记录长度设置的比2550ms 还要长，才能使地震波传播到目的层并被检波点记录。
2上覆水平界面反射系数变化对目
的层地震波振幅的影响
由于模型中速度异常体位于中心位置，提取的目的层地震波振幅图及振幅曲线都呈中心对称，为了便于分析，只截取中心位置（y width =2550m ）的振
幅曲线来进行分析。图3展示的是当异常体速度V =2600m/s 、半径R =750m 时，目的层地震波的
振幅［图3（a ）］及振幅曲线［图3（b ）］。
图3目的层地震波振幅记录结果
Fig.3Seismic amplitude of the target zone
模型大小（长×宽×深）/m 35100×5100×2580
网格大小/m 15
炮间距/m 15
子波主频/Hz 30
采样率/ms 1
记录长度/ms 3000
异常体
厚度/m 540
（b ）三维实验模型俯视图
（a ）三维实验模型剖面图
x wid th /m
（b ）y width =2550m 处目的层地震波振幅曲线
（a ）目的层地震波振幅图
750
1290
2580
深度/m
3000m /s
zhilian3000m /s
5100
x wid th /m
5100
y w i d t h /m
5100
x wid th /m
5100
y w i d t h /m
5100
01.81.61.41.2
1.00.80.60.40.2E
750
1500
225030003750
turnon45005100
x wid th /m
V =2600m/s
122岩性油气藏
第29卷第3期
2.1
上覆异常体速度变化对目的层地震波振幅的
影响
将异常体速度V 的取值范围设置为2000～4000m /s ，每隔100m/s 记录一次目的层的振幅结果，振幅曲线汇总情况如图4所示。所有振幅曲线
关于x width =2550m 处的中心点对称，曲线发生剧烈
软件学校抖动区域的x width 为1800～3300m ，即速度异常体（蓝色区域）所在的区域，各曲线变化最明显之处是其中心峰值的大小。此外，振幅曲线的最小值出现在速度异常区域边界附近位置。
图4速度V 变化时，目的层地震波振幅曲线汇总
Fig.4
Seismic amplitude curves of the target zone when velocity changes
ouz>general图5是异常体速度V 分别为2500m /s 、3000m /s 和3500m /s 时记录的目的层地震波振幅曲线的对比图，结合图4可以看出，V =3000m /s 是一个分界点，此时模型退化为层状介质模型，目的层的振幅曲线为一条水平线。V =2500m /s 时的目的层地震波振幅曲线的峰值明显大于V =3500m /s 时的目的层地震波振幅曲线峰值。
图5目的层地震波振幅曲线对比
Fig.5
Seismic amplitude curves of the target zone 图6所示的变化曲线是由图4中所有振幅曲线的峰值（振幅曲线的中心点值）绘制而成的目的层地震波振幅峰值变化曲线，目的层地震波振幅曲线的峰值先随
着异常体的速度增大而增大，并在异常体速度V =2700m /s 时达到最大值，而后陡然变小，且在异常体速度V =3000m /s 时，峰值最小，之
图6目的层地震波振幅曲线峰值变化统计
Fig.6Change of amplitude peak values of the target zone
深度/m
3000m /s
3000m /s
六年级数学辅导750
1290
22502580075015002250
3000375045005100
750
1500
2250
3000
3750
4500
5100
2.2
2.01.81.61.41.21.00.80.60.40.2E
x wid th /mtoefl考试培训
750
1500
2250
3000
3750
45005100
x wid th /m
1.61.41.2
1.00.80.60.40.2E
V =2500m/s
V =3000m/s V =3500m/s
2.22.01.8
1.61.41.21.00.8
0.6E
V /（m ·s -1）
1900
23002700310035003900
