CN109164488B - 一种梯形网格有限差分地震波场模拟方法 - Google Patents

一种梯形网格有限差分地震波场模拟方法 Download PDF

Info

Publication number
CN109164488B
CN109164488B CN201811178312.6A CN201811178312A CN109164488B CN 109164488 B CN109164488 B CN 109164488B CN 201811178312 A CN201811178312 A CN 201811178312A CN 109164488 B CN109164488 B CN 109164488B
Authority
CN
China
Prior art keywords
coordinate system
trapezoidal
under
time
difference
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811178312.6A
Other languages
English (en)
Other versions
CN109164488A (zh
Inventor
吴帮玉
孟德林
俞嘉旭
刘乃豪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201811178312.6A priority Critical patent/CN109164488B/zh
Publication of CN109164488A publication Critical patent/CN109164488A/zh
Application granted granted Critical
Publication of CN109164488B publication Critical patent/CN109164488B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
    • G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface

Abstract

一种梯形网格有限差分地震波场模拟方法，首先确定梯形坐标系坐标与直角坐标系坐标之间的转换关系，给出卷积完美匹配层吸收边界条件下的二维波动方程在该梯形坐标系表达形式。对于梯形坐标系波动方程中出现空间混合偏导，利用坐标变换将混合偏导化解为非混合偏导后求解梯形坐标系下的波动方程。本发明提出的方法可避免不同大小网格区域过渡所产生的虚假反射，能有效减少传统固定网格有限差分对深部高速区域的过采样。Marmousi模型测试表明，与传统的固定网格有限差分算法相比，本发明提高有限差分地震波正演效率的同时可显著降低内存使用量。

Description

一种梯形网格有限差分地震波场模拟方法
技术领域
本发明属于地震勘探技术领域，涉及一种地震波场模拟方法，尤其是一种梯形网格有限差分地震波场模拟方法。
背景技术
如何精确并有效地求解波动方程已成为地球物理学的一个重要课题，尤其在地震勘探的许多应用中都涉及这一问题，例如地震数据解释、观测系统设计、高精度成像和建立模型等(Etgen and O’Brien,2007)。用于模拟地震波场的常见方法主要有以下几种：有限差分法、有限元法、体积元法、伪谱法和粒子法等(黎昌成和张美根，2018；Hu and Jia,2016)。其中，有限差分法有着计算量小、内存需求少和易于实现等优点，因而在实际应用中被广泛采用(Etgen and O’Brien,2007；黄超和董良国，2009；Yao et al.,2016)。
在传统有限差分方法中，模型中的空间采样步长是固定的，即用固定网格对介质进行离散(Dablain，1986)。与固定网格有限差分相比，变网格有限差分可以在兼顾地震波模拟精度的同时提升模拟的效率(Moczo，1989；朱生旺等，2007；赵海波和王秀明，2007)。变网格有限差分在低速区域采用细网格剖分，而在高速区则采用粗网格剖分。该方法的关键在于如何处理粗细网格区域的过渡，根据处理方法的不同可以分为基于插值的算法(Jastram and Behle,1992；Hayashi and Burns,1999；Wang et al.,2001；Pasalic andMcGarry,2010)、基于邻域差分近似的算法(黄超和董良国，2009；孙小东等，2012；Takekawaet al.,2015；Li et al.,2017)和基于坐标变换的算法(Jastram and Tessmer,1994；Atleet al.,2011；Chen and Xu,2012)。其中，基于坐标变换的变网格有限差分可以有效避免不同大小网格区域过渡产生的虚假反射，以及较小的附加计算量和网格大小可随速度相应变化的优点。Chen and Xu(2012)提出梯形网格剖分，并用伪谱法对弹性波波场进行了数值模拟，结果表明梯形网格剖分在计算量和内存消耗上比固定网格剖分有显著降低。
发明内容
针对上述问题，本发明的目的在于针对地震波传播速度由浅入深整体逐渐增大的特点，提出一种基于卷积完美匹配层吸收边界的梯形网格有限差分地震波场模拟方法，该方法利用梯形坐标系设计耦合速度由浅入深逐渐增大的变化，采用卷积完美匹配层吸收边界条件来消除虚假的边界反射，并且在时间效率和电脑内存需求上均显著优于传统的固定网格有限差分算法。
本发明的目的是通过以下技术方案来解决的：
一种梯形网格有限差分地震波场模拟方法，包括以下步骤：
(1)对给定的直角坐标系下的速度模型v(x0,z0)，确定梯形坐标系形状参数γ为：
Figure BDA0001824347340000021
式中，
Figure BDA0001824347340000022
表示最深层的纵向深度坐标位置，
Figure BDA0001824347340000023
是深度z0处的横向采样间隔，
Figure BDA0001824347340000024
是深度z0处的最大横向采样间隔，
Figure BDA0001824347340000025
是最深层最大横向采样间隔；其中，vmin(z0)是深度z0的最小速度，f0为震源主频率，NG为一个波长内采样点数；
(2)对给定的直角坐标系下的速度模型v(x0,z0)，进行坐标变换将其转换到梯形坐标系的速度模型v(x,z)，坐标变换表达式为：
Figure BDA0001824347340000026
式中，x和z分别是梯形坐标系下的横向与纵向坐标，α是梯形坐标横向位置参数，在地震波模拟时与震源横向位置相同，γ是梯形坐标系形状参数；
(3)根据步骤(2)中的坐标变换，将直角坐标系(x0,z0)下的CPML条件下二维时域各向同性声波方程变换为梯形坐标系(x,z)下的声波方程：
Figure BDA0001824347340000031
式中，u＝u(x,z)是梯形坐标系下的压力场，t是记录时间，xs和zs分别是震源在梯形坐标系下的横纵坐标位置，f(t)是震源的时域波形函数，δ(·)表示狄拉克函数，辅助变量ψi和ζi(i∈{x,z})的定义如下：
Figure BDA0001824347340000032
Figure BDA0001824347340000033
Figure BDA0001824347340000034
Figure BDA0001824347340000035
式中j表示第j个时刻，系数ai和bi(i∈{x,z})表达式如下：
Figure BDA0001824347340000036
Figure BDA0001824347340000037
其中Δt是时间间隔，
Figure BDA0001824347340000038
R是理论边界反射系数，vmax是模型最大速度，Li是CPML随着i变化的吸收厚度，
Figure BDA0001824347340000039
是吸收区域到计算区域的距离，
Figure BDA00018243473400000310
αmax＝πf0，f0是震源主频率；
(4)利用坐标变换，将步骤(3)中的梯形坐标系(x,z)下的声波方程转化为
Figure BDA0001824347340000041
其中坐标变换定义为：
Figure BDA0001824347340000042
其中(x',z')是变换后的坐标，θ为沿x方向的网格线与对角线的夹角；
(5)利用高阶有限差分对步骤(4)中的声波方程进行离散，得到CPML边界条件下梯形坐标系中的时间-空间域二维声波方程的差分方程，并利用此差分方程迭代求解，得到各个时刻的波场：
Figure BDA0001824347340000043
式中
Figure BDA0001824347340000044
表示空间位置(xm,zn)＝(x+(m-1)Δx,z+(n-1)Δz)和第j个计算时刻点tj＝t0+(j-1)Δt处的波场，Δx和Δz是梯形坐标系横向和纵向空间采样间隔，t0为震源时延；cp和ηp(p＝1,2,…,N)表示空间一阶和二阶偏导差分系数，N是对应于空间差分精度的差分点数；其中Δt为时间采样间隔。
本发明进一步的改进在于，步骤(2)的具体过程如下：
采用如下的坐标变换将直角坐标系转换为梯形坐标系：
Figure BDA0001824347340000051
式中，(x0,z0)是直角坐标系坐标，(x,z)是梯形坐标系坐标；α是梯形坐标横向位置参数，在地震波模拟中与震源横向位置相同；γ是梯形坐标系形状参数；
若梯形坐标系下横向采样间隔为Δx，则其对应物理直角坐标系采样间隔Δx0和Δz0表示为
Δx0(z)＝(1+γg(z))Δx
Δz0(z)＝g(z+Δz)-g(z)
若直角坐标系下的速度模型v(x0,z0)中每个深度最小速度表示为vmin(z0)，且保持震源主频率下波数空间采样点数相同，则每个深度最大横向采样间隔表示为
Figure BDA0001824347340000052
其中f0为震源主频率，NG为一个波长内采样点数。
本发明进一步的改进在于，利用三次样条插值来建立函数g(z)。
本发明进一步的改进在于，步骤(3)的具体过程如下：
在给定的梯形坐标系与直角坐标系转换关系下，利用偏导数的求导规则，将波场关于直角坐标系(x0,z0)的偏导转化为梯形坐标系(x,z)下的偏导；将直角坐标系压力场u(x0,z0)和梯形坐标系压力场u(x,z)分别简化为u0和u，则二者的一阶和二阶偏导关系为：
Figure BDA0001824347340000053
Figure BDA0001824347340000054
Figure BDA0001824347340000061
Figure BDA0001824347340000062
对于CPML边界条件下直角坐标系中的二维各向同性声波方程：
Figure BDA0001824347340000063
其中t为记录时间，δ(x)表示狄拉克函数，(x0s,z0s)是震源在直角坐标系中的横纵坐标，f(t)为时域震源函数，辅助变量ψi0和ζi0(i∈{x,z})的结构如下：
Figure BDA0001824347340000064
Figure BDA0001824347340000065
在辅助变量ψi0和ζi0(i∈{x,z})的表达式中，j为第j个计算时刻，参数
Figure BDA0001824347340000066
和
Figure BDA0001824347340000067
表达式如下：
Figure BDA0001824347340000068
Figure BDA0001824347340000069
其中Δt是时间间隔，
Figure BDA00018243473400000610
R是理论边界反射系数，vmax是模型最大速度，
Figure BDA00018243473400000611
是CPML随着i0变化的吸收厚度，
Figure BDA00018243473400000612
是吸收区域到计算区域的距离，
Figure BDA00018243473400000613
αmax＝πf0，f0是震源主频率；
根据两个坐标系间的偏导关系，对于上述CPML边界条件下直角坐标系中的二维各向同性声波方程，其在梯形坐标系下的声波方程为：
Figure BDA0001824347340000071
式中，u＝u(x,z)是梯形坐标系下的压力场，t是记录时间，xs和zs分别是震源在梯形坐标系下的横纵坐标位置，f(t)是震源的时域波形函数，δ(·)表示狄拉克函数，辅助变量ψi和ζi(i∈{x,z})的结构如下：
Figure BDA0001824347340000072
Figure BDA0001824347340000073
Figure BDA0001824347340000074
Figure BDA0001824347340000075
系数ai和bi(i∈{x,z})表达式如下：
Figure BDA0001824347340000076
Figure BDA0001824347340000077
其中Δt是时间间隔，
Figure BDA0001824347340000078
R是理论边界反射系数，vmax是模型最大速度，Li是卷积完美匹配层随着i变化的吸收厚度，
Figure BDA0001824347340000079
是吸收区域到计算区域的距离，
Figure BDA00018243473400000710
αmax＝πf0，其中f0是震源主频率。
本发明进一步的改进在于，步骤(4)的具体过程如下：
通过坐标变换将梯形坐标系中CPML边界条件下二维时域各向同性声波方程中出现的空间混合偏导项
Figure BDA0001824347340000081
转换为二阶非混合偏导，定义如下的坐标变换：
Figure BDA0001824347340000082
式中，(x',z')是变换后的坐标，θ为沿x方向的网格线与对角线的夹角，且满足0<θ<π/2；在此坐标变换下，混合偏导项
Figure BDA0001824347340000083
表示为：
Figure BDA0001824347340000084
将此式代入梯形坐标系下的声波方程，得到用于波场模拟的声波方程为
Figure BDA0001824347340000085
本发明进一步的改进在于，步骤(5)的具体过程为：
对步骤(4)中得到的声波方程，空间偏导采用高阶有限差分格式：
Figure BDA0001824347340000086
Figure BDA0001824347340000087
Figure BDA0001824347340000088
Figure BDA0001824347340000091
其中
Figure BDA0001824347340000092
表示空间位置(xm,zn)＝(x+(m-1)Δx,z+(n-1)Δz)和时刻tj＝t0+(j-1)Δt处的波场，Δx和Δz是梯形坐标系横向和纵向空间采样间隔，Δt为时间采样间隔，t0为震源时延；cp和ηp(p＝1,2,…,N)表示空间一阶和二阶偏导差分系数，N是对应于空间差分精度的差分点数；若梯形坐标系中横向和纵向空间采样间隔相同，则坐标系(x',z')中采样有Δx'2＝Δz'2＝Δx2+Δz2，所以
Figure BDA0001824347340000093
Figure BDA0001824347340000094
时间偏导数采用2阶，代入并化简后，梯形坐标系中CPML边界条件下时间空间域二维声波方程差分方程为
Figure BDA0001824347340000095
式中，
Figure BDA0001824347340000096
表示当且仅当m＝m0时，
Figure BDA0001824347340000097
其他时候为0；
Figure BDA0001824347340000098
同理；辅助变量ψi和ζi(i∈{x,z})及他们的偏导也采用高阶有限差分形式。
本发明进一步的改进在于，当Δx＝Δz时θ＝π/4，则差分方程简化为
Figure BDA0001824347340000101
利用差分方程进行波场迭代时，满足有限差分的稳定条件；以角度φ传播的平面波定义如下：
Figure BDA0001824347340000102
其中ω是角频率，
Figure BDA0001824347340000103
是平面波的振幅；将上式带入差分方程，且只考虑最大波数，得到Δt满足下面的稳定性条件：
Figure BDA0001824347340000104
与现有技术相比，本发明具有以下有益的效果：
(1)由于重力引起的岩石压实效应，地震波的传播速度通常会随着深度的增加而加快。基于这一事实，梯形坐标系设计可以耦合地下介质速度整体变化，该坐标系中均匀网格采样所对应的物理直角坐标系网格由浅入深逐渐增大，即浅部低速区对应细网格，深部高速区对应粗网格。因此，本发明在梯形坐标系下用有限差分求解波动方程，对应笛卡尔直角坐标系中为变网格有限差分，从而模拟地震波。
(2)本发明应用的梯形坐标系较直角坐标系更为一般，对应物理模型水平和深度方向采样间隔均随深度增加而增大，是一种变网格有限差分方法，且不会产生由于网格大小变化引起的虚假反射，能够有效减少传统固定网格有限差分对深部高速区域的过采样。
(3)本发明采用了卷积完美匹配层吸收边界条件来消除虚假的边界反射；
(4)本发明在计算效率上显著优于传统的固定网格有限差分算法，同时对电脑内存的需求大为降低。
附图说明
图1是直角坐标系和梯形坐标系的对应关系示意图，其中，(a)为直角坐标系下的梯形网格，(b)为梯形坐标系下的固定网格；图中黑色区域表示相同的采样位置。
图2是坐标变换中角度定义示意图。
图3是两层模型CPML边界条件下和无CPML边界条件下的波场快照：其中，(a)为带CPML边界条件，(b)为无CPML边界条件。黑色虚线表示两层的交界面，黑色五角星表示震源位置。
图4是Marmousi速度模型：其中，(a)为地震波场模拟的实际区域，(b)为梯形坐标系Marmousi模型，(c)为相应的纵向伸缩函数，(d)为纵向采样间隔从浅部到深部的变化图。
图5是Marmousi模型模拟波场快照：其中，(a)为梯形坐标系下固定网格有限差分的波场快照，(b)为梯形坐标系波场((a)图)经坐标变换和插值到直角坐标系的波场快照，(c)为直角坐标系下传统固定网格有限差分的波场快照。
图6是Marmousi模型波场快照在不同横向位置深度波形对比图：其中，(a)为3.5km，(b)为4.6875km，(c)为5.75km。实线和虚线分别表示传统和梯形坐标系固定网格有限差分的模拟结果。
具体实施方式
本发明针对地震波传播速度随深度逐渐增加这一普遍特征，提出一种二维梯形坐标系波动方程有限差分波场模拟方法，并采用了卷积完美匹配层吸收边界条件来消除虚假的边界反射。本发明基于卷积完美匹配层吸收边界，在梯形坐标系下用有限差分求解波动方程，实现一种物理横向和纵向网格采样间隔随深度增加的有限差分地震波模拟方法。首先确定梯形坐标系与直角坐标系下坐标位置之间的转换关系，然后推导卷积完美匹配层吸收边界条件下的二维波动方程在该梯形坐标系表达形式。最后，由于梯形坐标系波动方程出现空间混合偏导，本发明提出利用坐标变换将混合偏导化解为非混合偏导后，可利用常规均匀网格有限差分算法求解梯形坐标系下的波动方程，从而得到模拟波场。
下面结合附图和实例对本发明进行详细的描述。
本发明提出的一种基于卷积完美匹配层吸收边界的梯形网格有限差分地震波场正演模拟方法，其具体步骤如下：
(1)对给定的直角坐标系下的速度模型v(x0,z0)，确定梯形坐标系形状参数γ为：
Figure BDA0001824347340000121
式中，
Figure BDA0001824347340000122
表示最深层的纵向深度坐标位置，
Figure BDA0001824347340000123
是深度z0处的横向采样间隔，
Figure BDA0001824347340000124
是深度z0处的最大横向采样间隔，
Figure BDA0001824347340000125
是最深层最大横向采样间隔。其中，vmin(z0)是深度z0的最小速度，f0为震源主频率，NG为一个波长内采样点数。
(2)对给定的直角坐标系下的速度模型v(x0,z0)，进行坐标变换将其转换到梯形坐标系的速度模型v(x,z)，坐标变换表达式为：
Figure BDA0001824347340000126
式中，x和z分别是梯形坐标系下的横向与纵向坐标，α是梯形坐标横向位置参数，在地震波模拟时与震源横向位置相同，γ是梯形坐标系形状参数。纵向比例函数g(z)与速度有关，为使之后的有限差分有效进行，该函数必须是一阶和二阶连续的。在实际操作中，本发明利用三次样条插值来建立函数g(z)。
本发明中步骤(2)的具体过程如下：
采用如下的坐标变换将直角坐标系转换为梯形坐标系：
Figure BDA0001824347340000131
式中，(x0,z0)是直角坐标系坐标，(x,z)是梯形坐标系坐标。α是梯形坐标横向位置参数，在地震波模拟中一般与震源横向位置相同。γ是梯形坐标系形状参数，与模型速度随深度变化相关。纵向比例函数g(z)与速度有关，为使之后的有限差分有效进行，该函数必须为一阶和二阶连续。在本发明中，利用三次样条插值来建立函数g(z)。
若梯形坐标系下横向采样间隔为Δx，则其对应物理直角坐标系采样间隔Δx0和Δz0可表示为
Δx0(z)＝(1+γg(z))Δx
Δz0(z)＝g(z+Δz)-g(z)
若直角坐标系下的速度模型v(x0,z0)中每个深度最小速度表示为vmin(z0)，且保持震源主频率下波数空间采样点数相同，则每个深度最大横向采样间隔可表示为
Figure BDA0001824347340000132
其中f0为震源主频率，NG为一个波长内采样点数，可根据对模拟精度需求设定。
其中，
Figure BDA0001824347340000133
表示最深层最大横向采样间隔；那么，梯形坐标系形状参数γ可由如下表达式确定
Figure BDA0001824347340000141
图1所示为直角坐标系和梯形坐标系中对同一区域离散采样示意图。图1中(a)为直角坐标系一梯形区域，则在梯形坐标系下为一矩形区域(图1中(b))。二图中总采样点数相同，但图1中(a)中横向和纵向采样间隔随深度逐渐增大，而图1中(b)中为均匀采样。两个网格覆盖区域完全相同，黑色网格表示相同采样点。
(3)根据步骤(2)中的坐标变换，将直角坐标系(x0,z0)下的CPML(卷积完美匹配层吸收边界)条件下二维时域各向同性声波方程变换为梯形坐标系(x,z)下的声波方程：
Figure BDA0001824347340000142
式中，u＝u(x,z)是梯形坐标系下的压力场，t是时间，xs和zs分别是震源在梯形坐标系下的横纵坐标位置，f(t)是震源的时域波形函数，δ(·)表示狄拉克函数，辅助变量ψi和ζi(i∈{x,z})的定义如下：
Figure BDA0001824347340000143
Figure BDA0001824347340000144
Figure BDA0001824347340000145
Figure BDA0001824347340000146
式中j表示第j个时刻，系数ai和bi(i∈{x,z})表达式如下：
Figure BDA0001824347340000151
Figure BDA0001824347340000152
其中Δt是时间间隔，
Figure BDA0001824347340000153
R是理论边界反射系数，vmax是模型最大速度，Li是CPML随着i变化的吸收厚度，
Figure BDA0001824347340000154
是吸收区域到计算区域的距离，
Figure BDA0001824347340000155
αmax＝πf0，其中f0是震源主频率。
步骤(3)的具体过程如下：
在给定的梯形坐标系与直角坐标系转换关系下，利用偏导数的求导规则，可将波场关于直角坐标系(x0,z0)的偏导转化为梯形坐标系(x,z)下的偏导。将直角坐标系压力场u(x0,z0)和梯形坐标系压力场u(x,z)分别简化为u0和u，则二者的一阶和二阶偏导关系为：
Figure BDA0001824347340000156
Figure BDA0001824347340000157
Figure BDA0001824347340000158
Figure BDA0001824347340000159
对于CPML边界条件下直角坐标系中的二维各向同性声波方程：
Figure BDA00018243473400001510
其中t为记录时间，δ(x)表示狄拉克函数，(x0s,z0s)是震源在直角坐标系中的横纵坐标，f(t)为时域震源函数，辅助变量
Figure BDA0001824347340000161
和
Figure BDA0001824347340000162
的结构如下：
Figure BDA0001824347340000163
Figure BDA0001824347340000164
在辅助变量的表达式中，j为第j个计算时刻，参数
Figure BDA0001824347340000165
和
Figure BDA0001824347340000166
表达式如下：
Figure BDA0001824347340000167
Figure BDA0001824347340000168
其中Δt是时间间隔，
Figure BDA0001824347340000169
R是理论边界反射系数，vmax是模型最大速度，
Figure BDA00018243473400001610
是CPML随着i0变化的吸收厚度，
Figure BDA00018243473400001611
是吸收区域到计算区域的距离，
Figure BDA00018243473400001612
αmax＝πf0，其中f0是震源主频率。
那么根据两个坐标系间的偏导关系，对于上述CPML边界条件下直角坐标系中的二维各向同性声波方程，其在梯形坐标系下的声波方程为：
Figure BDA00018243473400001613
式中，u＝u(x,z)是梯形坐标系下的压力场，t是时间，xs和zs分别是震源在梯形坐标系下的横纵坐标位置，f(t)是时域震源函数，δ(·)表示狄拉克函数，辅助变量ψi和ζi(i∈{x,z})的结构如下：
Figure BDA00018243473400001614
Figure BDA0001824347340000171
Figure BDA0001824347340000172
Figure BDA0001824347340000173
系数ai和bi(i∈{x,z})表达式如下：
Figure BDA0001824347340000174
Figure BDA0001824347340000175
其中Δt是时间间隔，
Figure BDA0001824347340000176
R是理论边界反射系数，vmax是模型最大速度，Li是卷积完美匹配层随着i变化的吸收厚度，
Figure BDA0001824347340000177
是吸收区域到计算区域的距离，
Figure BDA0001824347340000178
αmax＝πf0，其中f0是震源主频率。
(4)利用坐标变换，将步骤(3)中的梯形坐标系(x,z)下的声波方程转化为
Figure BDA0001824347340000179
其中坐标变换定义为：
Figure BDA00018243473400001710
其中(x',z')是变换后的坐标，θ为沿x方向的网格线与对角线的夹角，且满足0<θ<π/2。
步骤(4)的具体过程如下：
对于梯形坐标系中CPML边界条件下二维时域各向同性声波方程中出现的空间混合偏导项
Figure BDA0001824347340000181
其较非空间混合偏导计算量大，可以通过坐标变换将其转换为二阶非混合偏导，如图2所示，定义如下的坐标变换：
Figure BDA0001824347340000182
式中，(x',z')是变换后的坐标，θ为沿x方向的网格线与对角线的夹角，且满足0<θ<π/2。那么，在此坐标变换下，混合偏导项可以表示为：
Figure BDA0001824347340000183
将此式代入梯形坐标系下的声波方程，得到最终用于波场模拟的声波方程为
Figure BDA0001824347340000184
(5)利用高阶有限差分对步骤(4)中的声波方程进行离散，得到CPML边界条件下梯形坐标系中的时间-空间域二维声波方程的差分方程，并利用此差分方程迭代求解，得到各个时刻的波场：
Figure BDA0001824347340000185
式中
Figure BDA0001824347340000191
表示空间位置(xm,zn)＝(x+(m-1)Δx,z+(n-1)Δz)和第j个计算时刻点tj＝t0+(j-1)Δt处的波场，Δx和Δz是梯形坐标系横向和纵向空间采样间隔，t0为震源时延。cp和ηp(p＝1,2,…,N)表示空间一阶和二阶偏导差分系数，N是对应于空间差分精度的差分点数。其中Δt为时间采样间隔，需满足如下的稳定性条件：
Figure BDA0001824347340000192
步骤(5)的具体过程为：
对步骤(4)中得到的声波方程，为了提高模拟精度，空间偏导采用高阶有限差分格式：
Figure BDA0001824347340000193
Figure BDA0001824347340000194
Figure BDA0001824347340000195
Figure BDA0001824347340000196
其中
Figure BDA0001824347340000197
表示空间位置(xm,zn)＝(x+(m-1)Δx,z+(n-1)Δz)和时刻tj＝t0+(j-1)Δt处的波场，Δx和Δz是梯形坐标系横向和纵向空间采样间隔，Δt为时间采样间隔，t0为震源时延。cp和ηp(p＝1,2,…,N)表示空间一阶和二阶偏导差分系数，N是对应于空间差分精度的差分点数。若梯形坐标系中横向和纵向空间采样间隔相同(即Δx＝Δz)，则坐标系(x',z')中采样有关系Δx'2＝Δz'2＝Δx2+Δz2，所以有
Figure BDA0001824347340000201
Figure BDA0001824347340000202
时间偏导数采用2阶，代入并化简后，则梯形坐标系中CPML边界条件下时间空间域二维声波方程差分方程为
Figure BDA0001824347340000203
式中，
Figure BDA0001824347340000204
表示当且仅当m＝m0时，
Figure BDA0001824347340000205
其他时候为0；
Figure BDA0001824347340000206
同理。辅助变量ψi和ζi(i∈{x,z})及他们的偏导也采用高阶有限差分形式。
特别地，当Δx＝Δz时θ＝π/4，则差分方程可进一步简化为
Figure BDA0001824347340000211
利用上述差分方程进行波场迭代时，需要满足有限差分的稳定条件。以角度φ传播的平面波定义如下：
Figure BDA0001824347340000212
其中ω是角频率，
Figure BDA0001824347340000213
是平面波的振幅。将上式带入差分方程，且只考虑最大波数(带奎斯特波数)，得到Δt必须满足下面的稳定性条件：
Figure BDA0001824347340000214
数值仿真结果
下面给出不同速度模型的实验结果来验证本发明所提方法的有效性。首先在一个两层速度模型上展示了本算法的计算结果，随后在Marmousi模型下进行数值运算，并与传统固定网格有限差分的结果进行比较。
数值实验1：两层速度模型数值模拟
该两层模型第一层速度为2000m/s，第二层速度为4000m/s，用主频为20Hz的ricker子波为源放置在顶层正中间进行激发。横向采样间隔Δx0从第一层的4m增加至最底层的16m(对应γ为2.4834*10-4)。梯形坐标系下有限差分时间采样间隔为1.783ms。图3中(a)为采用卷积完美匹配层边界条件的波场快照，图3(b)为无CPML吸收边界条件下的波场快照。可以看出，采用CPML边界条件很好地避免了边界反射。
数值实验2：Marmousi模型数值模拟
该原始模型横向和纵向采样点分别为256和425，横向采样间隔Δx0为12.5m，纵向采样间隔Δz0为4m，用主频为12Hz的ricker子波为源放置在顶层正中间进行激发。梯形坐标系下有限差分时间采样间隔为1.404ms。图4中(a)展示了实际有限差分模拟区域，该区域在梯形坐标系下显示如图4(b)所示，相应的纵向比例函数如图4(c)所示。图4(d)表明实际的纵向采样间隔从第一层的8.333m增加至最底层的18.333m。图5中(a)展示了梯形坐标系下有限差分波场快照，图5(b)显示的是将图5(a)内插至与原始速度模型相同采样间隔波场。作为对比，图5(c)显示基于原始速度模型采样的固定网格传统有限差分相同时刻波场快照，与图5(b)相比二者整体上一致。图6展示了三个不同横向位置的深度域波形对比，虚线代表梯形网格有限差分的结果，实线代表常规固定网格有限差分的结果。图6中(a)、图6中(b)和图6中(c)分别距离震源3.5km、4.6875km和5.75km。可以看出这两种方法的结果是一致的。
本发明公开利用梯形坐标系可以耦合地下介质速度整体变化的优点，推导了梯形坐标系下波动方程的表达式并应用有限差分方法求解。数值算法上，采用常规均匀网格有限差分算法；而在物理上，则是对应变网格有限差分，即浅部低速区对应细网格，深部高速区对应粗网格，本发明可避免不同大小网格区域过渡所产生的虚假反射。此外，本发明采用卷积完美匹配层吸收边界条件来消除虚假的边界反射。与传统的固定网格有限差分算法相比，Marmousi模型测试结果表明，本发明提出的方法快2.9倍，同时节省了80％的电脑内存。
本发明通过首先确定梯形坐标系坐标与直角坐标系坐标之间的转换关系，给出卷积完美匹配层吸收边界条件下的二维波动方程在该梯形坐标系表达形式。对于梯形坐标系波动方程中出现空间混合偏导，利用坐标变换将混合偏导化解为非混合偏导后求解梯形坐标系下的波动方程。本发明提出的方法可避免不同大小网格区域过渡所产生的虚假反射，能有效减少传统固定网格有限差分对深部高速区域的过采样。Marmousi模型测试表明，与传统的固定网格有限差分算法相比，本发明提高有限差分地震波正演效率的同时可显著降低内存使用量。

Claims (7)

1.一种梯形网格有限差分地震波场模拟方法，其特征在于，包括以下步骤：
(1)对给定的直角坐标系下的速度模型v(x0,z0)，确定梯形坐标系形状参数γ为：
Figure FDA0002280669060000011
式中，
Figure FDA0002280669060000012
表示最深层的纵向深度坐标位置，
Figure FDA0002280669060000013
是深度z0处的横向采样间隔，
Figure FDA0002280669060000017
是深度z0处的最大横向采样间隔，
Figure FDA0002280669060000014
是最深层最大横向采样间隔；其中，vmin(z0)是深度z0的最小速度，f0为震源主频率，NG为一个波长内采样点数；
(2)对给定的直角坐标系下的速度模型v(x0,z0)，进行坐标变换将其转换到梯形坐标系的速度模型v(x,z)，坐标变换表达式为：
Figure FDA0002280669060000015
式中，x和z分别是梯形坐标系下的横向与纵向坐标，α是梯形坐标横向位置参数，在地震波模拟时与震源横向位置相同，γ是梯形坐标系形状参数，g(z)为纵向比例函数；
(3)根据步骤(2)中的坐标变换，将直角坐标系(x0,z0)下的CPML条件下二维时域各向同性声波方程变换为梯形坐标系(x,z)下的声波方程：
Figure FDA0002280669060000016
式中，u＝u(x,z)是梯形坐标系下的压力场，t是记录时间，xs和zs分别是震源在梯形坐标系下的横纵坐标位置，f(t)是震源的时域波形函数，δ(·)表示狄拉克函数，辅助变量ψi和ζi的定义如下：
Figure FDA0002280669060000021
Figure FDA0002280669060000022
Figure FDA0002280669060000023
Figure FDA0002280669060000024
式中j表示第j个时刻，i∈{x,z}，系数ai和bi表达式如下：
Figure FDA0002280669060000025
Figure FDA0002280669060000026
其中，Δt是时间间隔，
Figure FDA0002280669060000027
R是理论边界反射系数，vmax是模型最大速度，Li是CPML随着i变化的吸收厚度，
Figure FDA0002280669060000028
是吸收区域到计算区域的距离，
Figure FDA0002280669060000029
αmax＝πf0，f0是震源主频率；
(4)利用坐标变换，将步骤(3)中的梯形坐标系(x,z)下的声波方程转化为
Figure FDA00022806690600000210
其中坐标变换定义为：
Figure FDA0002280669060000031
其中(x',z')是变换后的坐标，θ为沿x方向的网格线与对角线的夹角，xs表示震源函数横坐标位置，zs表示震源函数纵坐标位置；
(5)利用高阶有限差分对步骤(4)中的声波方程进行离散，得到CPML边界条件下梯形坐标系中的时间-空间域二维声波方程的差分方程，并利用此差分方程迭代求解，得到各个时刻的波场：
Figure FDA0002280669060000032
式中
Figure FDA0002280669060000033
表示空间位置(xm,zn)＝(x+(m-1)Δx,z+(n-1)Δz)和第j个计算时刻点tj＝t0+(j-1)Δt处的波场，Δx和Δz是梯形坐标系横向和纵向空间采样间隔，t0为震源时延；cp和ηp表示空间一阶和二阶偏导差分系数，p＝1,2,…,N，N是对应于空间差分精度的差分点数；其中Δt为时间采样间隔，
Figure FDA0002280669060000034
表示当且仅当m＝m0时，
Figure FDA0002280669060000035
其他时候为0；
Figure FDA0002280669060000036
表示当且仅当n＝n0时，
Figure FDA0002280669060000037
其他时候为0。
2.根据权利要求1所述的一种梯形网格有限差分地震波场模拟方法，其特征在于，步骤(2)的具体过程如下：
采用如下的坐标变换将直角坐标系转换为梯形坐标系：
Figure FDA0002280669060000038
式中，(x0,z0)是直角坐标系坐标，(x,z)是梯形坐标系坐标；α是梯形坐标横向位置参数，在地震波模拟中与震源横向位置相同；γ是梯形坐标系形状参数；
若梯形坐标系下横向采样间隔为Δx，则其对应物理直角坐标系采样间隔Δx0和Δz0表示为
Δx0(z)＝(1+γg(z))Δx
Δz0(z)＝g(z+Δz)-g(z)
若直角坐标系下的速度模型v(x0,z0)中每个深度最小速度表示为vmin(z0)，且保持震源主频率下波数空间采样点数相同，则每个深度最大横向采样间隔表示为
Figure FDA0002280669060000041
其中f0为震源主频率，NG为一个波长内采样点数。
3.根据权利要求2所述的一种梯形网格有限差分地震波场模拟方法，其特征在于，利用三次样条插值来建立函数g(z)。
4.根据权利要求1所述的一种梯形网格有限差分地震波场模拟方法，其特征在于，步骤(3)的具体过程如下：
在给定的梯形坐标系与直角坐标系转换关系下，利用偏导数的求导规则，将波场关于直角坐标系(x0,z0)的偏导转化为梯形坐标系(x,z)下的偏导；将直角坐标系压力场u(x0,z0)和梯形坐标系压力场u(x,z)分别简化为u0和u，则二者的一阶和二阶偏导关系为：
Figure FDA0002280669060000042
Figure FDA0002280669060000043
Figure FDA0002280669060000044
Figure FDA0002280669060000051
对于CPML边界条件下直角坐标系中的二维各向同性声波方程：
Figure FDA0002280669060000052
其中t为记录时间，δ(x)表示狄拉克函数，
Figure FDA0002280669060000053
是震源在直角坐标系中的横纵坐标，f(t)为时域震源函数，辅助变量
Figure FDA0002280669060000054
和
Figure FDA0002280669060000055
的结构如下：
Figure FDA0002280669060000056
Figure FDA0002280669060000057
在辅助变量
Figure FDA0002280669060000058
和
Figure FDA0002280669060000059
的表达式中，j为第j个计算时刻，i∈{x,z}，参数
Figure FDA00022806690600000510
和
Figure FDA00022806690600000511
表达式如下：
Figure FDA00022806690600000512
Figure FDA00022806690600000513
其中Δt是时间间隔，
Figure FDA00022806690600000514
R是理论边界反射系数，vmax是模型最大速度，
Figure FDA00022806690600000515
是CPML随着i0变化的吸收厚度，
Figure FDA00022806690600000516
是吸收区域到计算区域的距离，
Figure FDA00022806690600000517
αmax＝πf0，f0是震源主频率；
根据两个坐标系间的偏导关系，对于上述CPML边界条件下直角坐标系中的二维各向同性声波方程，其在梯形坐标系下的声波方程为：
Figure FDA0002280669060000061
式中，u＝u(x,z)是梯形坐标系下的压力场，t是记录时间，xs和zs分别是震源在梯形坐标系下的横纵坐标位置，f(t)是震源的时域波形函数，δ(·)表示狄拉克函数，辅助变量ψi和ζi的结构如下：
Figure FDA0002280669060000062
Figure FDA0002280669060000063
Figure FDA0002280669060000064
Figure FDA0002280669060000065
系数ai和bi表达式如下：
Figure FDA0002280669060000066
Figure FDA0002280669060000067
其中，i∈{x,z}，Δt是时间间隔，
Figure FDA0002280669060000068
R是理论边界反射系数，vmax是模型最大速度，Li是卷积完美匹配层随着i变化的吸收厚度，
Figure FDA0002280669060000069
是吸收区域到计算区域的距离，
Figure FDA00022806690600000610
αmax＝πf0，其中f0是震源主频率。
5.根据权利要求4所述的一种梯形网格有限差分地震波场模拟方法，其特征在于，步骤(4) 的具体过程如下：
通过坐标变换将梯形坐标系中CPML边界条件下二维时域各向同性声波方程中出现的空间混合偏导项
Figure FDA0002280669060000071
转换为二阶非混合偏导，定义如下的坐标变换：
Figure FDA0002280669060000072
式中，(x',z')是变换后的坐标，θ为沿x方向的网格线与对角线的夹角，且满足0<θ<π/2；在此坐标变换下，混合偏导项
Figure FDA0002280669060000073
表示为：
Figure FDA0002280669060000074
将此式代入梯形坐标系下的声波方程，得到用于波场模拟的声波方程为
Figure FDA0002280669060000075
6.根据权利要求1所述的一种梯形网格有限差分地震波场模拟方法，其特征在于，步骤(5)的具体过程为：
对步骤(4)中得到的声波方程，空间偏导采用高阶有限差分格式：
Figure FDA0002280669060000076
Figure FDA0002280669060000077
Figure FDA0002280669060000081
Figure FDA0002280669060000082
其中
Figure FDA0002280669060000083
表示空间位置(xm,zn)＝(x+(m-1)Δx,z+(n-1)Δz)和时刻tj＝t0+(j-1)Δt处的波场，Δx和Δz是梯形坐标系横向和纵向空间采样间隔，Δt为时间采样间隔，t0为震源时延；cp和ηp表示空间一阶和二阶偏导差分系数，p＝1,2,…,N，N是对应于空间差分精度的差分点数；若梯形坐标系中横向和纵向空间采样间隔相同，则坐标系(x',z')中采样有Δx'2＝Δz'2＝Δx2+Δz2，所以
Figure FDA0002280669060000084
Figure FDA0002280669060000085
时间偏导数采用2阶，代入并化简后，梯形坐标系中CPML边界条件下时间空间域二维声波方程差分方程为
Figure FDA0002280669060000091
式中，
Figure FDA0002280669060000092
表示当且仅当m＝m0时，
Figure FDA0002280669060000093
其他时候为0；
Figure FDA0002280669060000094
同理；辅助变量ψi和ζi及他们的偏导也采用高阶有限差分形式，i∈{x,z}。
7.根据权利要求6所述的一种梯形网格有限差分地震波场模拟方法，其特征在于，当Δx＝Δz时θ＝π/4，则差分方程简化为
Figure FDA0002280669060000095
利用差分方程进行波场迭代时，满足有限差分的稳定条件；以角度φ传播的平面波定义如下：
Figure FDA0002280669060000096
其中ω是角频率，
Figure FDA0002280669060000097
Figure FDA0002280669060000098
是平面波的振幅；将上式带入差分方程，且只考虑最大波数，得到Δt满足下面的稳定性条件：
Figure FDA0002280669060000101
CN201811178312.6A 2018-10-10 2018-10-10 一种梯形网格有限差分地震波场模拟方法 Active CN109164488B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811178312.6A CN109164488B (zh) 2018-10-10 2018-10-10 一种梯形网格有限差分地震波场模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811178312.6A CN109164488B (zh) 2018-10-10 2018-10-10 一种梯形网格有限差分地震波场模拟方法

Publications (2)

Publication Number Publication Date
CN109164488A CN109164488A (zh) 2019-01-08
CN109164488B true CN109164488B (zh) 2020-03-17

Family

ID=64877914

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811178312.6A Active CN109164488B (zh) 2018-10-10 2018-10-10 一种梯形网格有限差分地震波场模拟方法

Country Status (1)

Country Link
CN (1) CN109164488B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112649859B (zh) * 2019-10-12 2024-03-22 中国石油化工股份有限公司 一种地震波速度自适应无网格场节点建立方法及系统
CN112379422A (zh) * 2020-10-30 2021-02-19 中国石油天然气集团有限公司 垂变网格地震波场外推方法及装置

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5999488A (en) * 1998-04-27 1999-12-07 Phillips Petroleum Company Method and apparatus for migration by finite differences
US8494778B2 (en) * 2009-10-05 2013-07-23 Chevron U.S.A. Inc. Variable grid for finite difference computation
CN103149586B (zh) * 2013-02-04 2016-01-20 西安交通大学 一种倾斜层状粘弹性介质中波场正演模拟方法
CN104990853B (zh) * 2015-06-30 2017-07-21 中国石油大学(北京) 多孔介质全阶渗透率张量的预测方法
CN107526105A (zh) * 2017-08-09 2017-12-29 西安交通大学 一种波场模拟交错网格有限差分方法
CN107657137B (zh) * 2017-11-09 2021-08-20 吉林大学 一种有理函数逼近的分数阶电磁反常扩散三维模拟方法
CN108051855B (zh) * 2017-12-13 2019-02-15 国家深海基地管理中心 一种基于拟空间域声波方程的有限差分计算方法

Also Published As

Publication number Publication date
CN109164488A (zh) 2019-01-08

Similar Documents

Publication Publication Date Title
Rawlinson et al. Seismic ray tracing and wavefront tracking in laterally heterogeneous media
CN104122585B (zh) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
CN109885945B (zh) 一种半空间环境下的边界元法近场声全息变换方法
JP7142968B2 (ja) フルウェーブフォームインバージョン方法、装置及び電子設備
CN110109177B (zh) 基于旋转时空双变网格有限差分法的地震波正演模拟方法
CN110456417B (zh) 一种地震数据多次波压制方法
CN110133644B (zh) 基于插值尺度函数法的探地雷达三维正演方法
CN109164488B (zh) 一种梯形网格有限差分地震波场模拟方法
CN111766628A (zh) 一种预条件的时间域弹性介质多参数全波形反演方法
Fang et al. A hybrid approach to solve the high-frequency Helmholtz equation with source singularity in smooth heterogeneous media
Muratov et al. Grid-characteristic method on unstructured tetrahedral meshes
US9928315B2 (en) Re-ordered interpolation and convolution for faster staggered-grid processing
Hu et al. Numerical modeling of seismic waves using frequency-adaptive meshes
CN116520418A (zh) 一种弹性波角度域共成像点道集高效提取方法
Wu et al. A dispersion minimizing subgridding finite difference scheme for the Helmholtz equation with PML
Wu et al. Trapezoid coordinate finite difference modeling of acoustic wave propagation using the CPML boundary condition
CN105319594A (zh) 一种基于最小二乘参数反演的傅里叶域地震数据重构方法
CN113221392A (zh) 一种非均匀粘声波在无限域内传播模型的构建方法
CN104750954A (zh) 一种在复杂各向异性介质中模拟地震波的方法及装置
Ramos-Martínez et al. Full-waveform inversion by pseudo-analytic extrapolation
CN110609325B (zh) 弹性波场数值模拟方法及系统
Li et al. The numerical modeling of frequency domain scalar wave equation based on the rhombus stencil
CN114460640A (zh) 有限差分模拟弹性波全波形反演方法和装置
CN112444849A (zh) 一种基于交错网格低秩有限差分的弹性逆时偏移成像方法
CN110764146A (zh) 一种基于声波算子的空间互相关弹性波反射波形反演方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant