CN108303736B - 各向异性ti介质最短路径射线追踪正演方法 - Google Patents
各向异性ti介质最短路径射线追踪正演方法 Download PDFInfo
- Publication number
- CN108303736B CN108303736B CN201711287053.6A CN201711287053A CN108303736B CN 108303736 B CN108303736 B CN 108303736B CN 201711287053 A CN201711287053 A CN 201711287053A CN 108303736 B CN108303736 B CN 108303736B
- Authority
- CN
- China
- Prior art keywords
- nodes
- group velocity
- node
- medium
- travel time
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 35
- 230000014509 gene expression Effects 0.000 claims description 29
- 239000013598 vector Substances 0.000 claims description 22
- 238000001514 detection method Methods 0.000 claims description 17
- 239000011435 rock Substances 0.000 claims description 15
- 239000000523 sample Substances 0.000 claims description 14
- 238000007689 inspection Methods 0.000 claims description 4
- 238000012986 modification Methods 0.000 claims description 4
- 230000004048 modification Effects 0.000 claims description 4
- 239000013078 crystal Substances 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 2
- 238000006467 substitution reaction Methods 0.000 claims description 2
- 238000009826 distribution Methods 0.000 description 16
- 238000010586 diagram Methods 0.000 description 6
- 238000005452 bending Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 3
- 241000287196 Asthenes Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 230000002250 progressing effect Effects 0.000 description 1
- 230000001681 protective effect Effects 0.000 description 1
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. for interpretation or for event detection
- 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
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Image Generation (AREA)
Abstract
本发明提供一种各向异性TI介质最短路径射线追踪正演方法,其包括如下步骤:S1、建立各向异性介质的群速度模型;S2、通过各向异性介质最短路径射线追踪算法进行各向异性介质最短路径射线追踪。
Description
技术领域
本发明涉及各向异性介质射线追踪正演技术领域,特别涉及一种各向异性TI介质最短路径射线追踪正演方法。
背景技术
现有技术对各向异性介质射线追踪正演进行研究时,大多数采用各向异性介质的简化形式,即对称轴倾角垂直的介质,也称为VTI介质。选择VTI介质进行研究有两大理由:首先,自然界很多岩石矿物的结构特点类似于VTI介质的组成结构,例如沉积岩中的泥岩、成层性较好的岩石和具有水平裂隙的岩石。其次,VTI介质包含五个弹性模量参数使各向异性介质在一定的程度得到了简化,各向异性介质射线追踪算法里可以避免复杂的几何关系式。Qin与Schuster(1993)根据惠更斯原理提出了二维VTI介质的初至走时计算方法;Ruger与Alkhalifah(1996)针对二维VTI介质提出了一种有效的射线追踪方法;Qian与Symes(2002)针对二维、三维VTI介质提出了利用有限差分法求解体波走时的计算方法。自然界的岩石沉积之后往往会伴随一系列地壳运动,这会导致地层的对称轴倾角存在不同情况,而不仅仅只有垂直对称轴倾角(VTI介质情形)。这里简介的射线追踪方法就是针对具有任意对称轴倾向的地层,这种方法是各向同性介质最短路径射线追踪方法向各向异性介质情形的一种拓展应用。
实际上,将地层假设为弱各向异性介质进行射线追踪具有不合理性,因为自然界的岩石和地层具有任意各向异性强度。Thomsen(1986)指出Thomsen参数既可以描述弱各向异性介质,也可以描述强各向异性介质,差别在于定义的五个参数不一样,表征任意强度各向异性介质的参数为(α0,β0,ε,δ*,γ),表征弱各向异性介质的参数为(α0,β0,ε,δ,γ)。虽然这两组参数大部分相同,但是利用它们表达的群速度和相速度表达式差别比较大。弱各向异性介质地震射线追踪的参数模型为vp0,ε,δ。然而,任意强度各向异性介质地震地震射线追踪的参数模型为α0,β0,ε,δ*,γ,或者对应的弹性模量参数。
发明内容
有鉴于此,本发明提出一种各向异性TI介质最短路径射线追踪正演方法。
一种各向异性TI介质最短路径射线追踪正演方法,其包括如下步骤:
S1、建立各向异性介质的群速度模型;
S2、通过各向异性介质最短路径射线追踪算法进行各向异性介质最短路径射线追踪。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中,
二维情形的群速度表达式{Um=(x,θ0,n),m=1,2,3}与三维情形的群速度表达式其中m=1,2,3分别对应三种波的模式qP,qSV和qSH,三种波模式的群速度都是空间坐标x、岩石对称轴倾角θ0与方位角以及相慢度方向n的函数。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中,
所述步骤S2包括:
S21、模型参数步骤;
S22、对走时场进行初始化;
S23、计算一个晶格内两个节点之间的最小走时;
S24、求取所有炮点的走时和射线;
S25、求取反射波与衍射波的射线路径。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中,
所述步骤S21包括:
以一定的网格单元大小,将二维各向异性介质模型m(x)划分成由主节点组成的网格模型
{m(xk),k=1,2,...,N1,N1=(Nx+1)(Nz+1)},其中Nx、Nz分别代表x和z方向的晶格数目;二维各向异性介质模型由多个矩形晶格组成,每个矩形晶格的四个角被定义为主节点;采用较大的晶格单元,同时在矩形晶格的每条边增加次一级节点来提高计算的精度与效率;将网格化的各向异性介质模型转化成为群速度模型Um(xk,θ0,n),m=1,2,3分别对应qP、qSV和qSH三种波模式,而k∈N1;群速度写成射线方向r0的函数Um(xk,θ0,r0),选择群速度的最小值从而确保群速度曲线的一致连续性;群速度模型晶格内部的群速度根据晶格四个顶点所在节点的群速度值,利用双线性插值函数或者拉格朗日插值函数进行求取
在模型参数化过程中增加次节点密度;如果次一级的节点数为N2,那么整个模型的节点数为N1+N2;此外,如果炮点、检波点位置和模型网格节点的位置不重合,在炮点和检波点位置增加额外的节点,因此整个节点网络包含三类节点:主节点、次节点、炮点和炮检点位置的节点,整个节点网络用节点序号l和空间坐标{xl=(xl,zl),l=1,2,...,N}来加以识别。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中,
所述步骤S22包括:
对炮点所在节点的走时赋初值"0",其它网格节点的走时值赋为无穷大,从炮点所在节点开始,计算炮点周围节点的走时,再从相邻节点计算其它节点的走时,并记录相邻节点的走时最小值。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中,
所述步骤S23包括:
将一个晶格内两个节点之间的走时通过如下公式近似求取,即利用节点i和节点j之间的距离除以两个节点之间的平均群速度
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中,
所述步骤S24包括:
当所有检波点的最小走时计算完成之后,拾取检点所在节点的走时,并根据入射节点的序号以反向追踪的形式,求取连接检波点与炮点之间的射线路径,得到共炮点道集的走时和射线路径;继续求取下一个炮点对应检波点的走时和射线路径,反复利用这一技术路线,以完成所有炮点的走时和射线求取。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中,
所述步骤S25包括:
同构一系列密集的节点来代表界面,把界面节点与所述的主节点、次节点、炮点和炮检点位置的节点组合起来,得到一个新的网络,包括主节点、次节点、界面节点和炮、检点节点;每个界面对应的反射波或者衍射波走时τ的计算方法如下:
τ=min{τS(xint)+τR(xint),xint∈Ωint},其中,τS(xint)和τR(xint)分别是炮点S与检波点R到界面节点位置xint的初至走时值,Ωint代表界面节点集合;分别假设炮点S与检波点R为震源点位置,进行两次正演走时计算;τS(xint)在第一次最短路径射线追踪计算时已经有值了,剩下步骤将检波点R作为震源点位置再做一次最短路径射线追踪计算求取走时τR(xint),然后根据步骤S25中公式寻找相应的界面反射节点,以及反射波的走时与射线路径。
实施本发明提供的各向异性TI介质最短路径射线追踪正演方法与现有技术相比具有以下有益效果:(1)针对的研究对象各向异性介质可以是任意强度各向异性介质情形,即包含强各向异性性质的介质和弱各向异性性质的介质。(2)各向异性介质的倾角可以包含VTI、HTI和TTI三种情形。(3)地震射线追踪算法是基于非均匀节点网格速度模型,这种网格模型的方案是采用较大的晶格单元,同时在矩形晶格的每条边增加次一级节点,这样可以提高射线追踪运算的精度与效率。
附图说明
图1是本发明实施例的各向异性TI介质最短路径射线追踪正演方法流程图;
图2是qP-qP反射波在倾斜界面模型中的射线路径(a)与走时(b)分布图;
图3是qP-qSV转换波在倾斜界面模型中的射线路径(a)与走时(b)分布图;
图4是qP-qSH转换波在倾斜界面模型中的射线路径(a)与走时(b)分布图;
图5是qP-qP反射波在弯曲界面模型中的射线路径(a)与走时(b)分布图;
图6是qP-qSV转换波在弯曲界面模型中的射线路径(a)与走时(b)分布图;
图7是qP-qSH转换波在弯曲界面模型中的射线路径(a)与走时(b)分布图。
具体实施方式
如图1所示,本发明提供一种各向异性TI介质最短路径射线追踪正演方法,其包括如下步骤:
S1、建立各向异性介质的群速度模型;
S2、通过各向异性介质最短路径射线追踪算法进行各向异性介质最短路径射线追踪。
如图2-7所示,以下对本发明实施例的原理进行进一步说明:
1.各向异性介质的群速度模型
如果将最短路径射线追踪方法从各向同性介质拓展至各向异性介质,首先需要组成各向异性介质的五个弹性模量参数{a11,a13,a33,a44,a66}(量纲为N/m2,牛顿每平方米)与对称轴倾角θ0(x),并且进一步研究它们和三个群速度公式{Um(x,θ0,r0),m=1,2,3}之间的关系,三个群速度分别对应qP波、qSV波和qSH波的群速度,可以看出群速度表达式不仅是空间位置x的函数,而且还是射线方向r0与介质对称轴倾角θ0(x)的函数。根据群速度模型,可以利用最短路径射线追踪方法计算各向异性介质qP波、qSV波和qSH波的初至走时和射线路径,通过进一步的修改还可以计算反射波、衍射波和转换波的走时和射线路径。VTI介质的相速度表达式,即克里斯托弗矩阵的特征值如下:
式中P和Q的表达式为
式中Q1、Q2和Q3的表达式为
上述关系式中,角度为相慢度向量和各向异性介质对称轴倾角之间的夹角,{a11,a13,a33,a44,a66}是各向异性介质的五个弹性模量参数,方程(4.1)的下标1、2和3分别代表qP波、qSV波和qSH波。计算群速度向量大小的具体表达式如下:
为了避免m与水平、垂直分量标记h,z相混淆,式中波的模式m写成上标形式。
为了使公式(4.1)-(4.7)能够应用于一般TI介质,可以利用标准三角关系(坐标旋转)得到入射角的余弦表达式。入射角是指相慢度方向与TI介质对称轴的夹角,而相慢度方向通常是由计算坐标系给出,不是由对称轴坐标系给出。如果知道相慢度方向和各向异性介质的对称轴倾角那么入射角可以表示为然后,才能应用公式(4.1)计算三个波的相速度,所以要由和来求出
代入ex,ez,可以得到U(m)的表达式:
可以采用如下表达式,对二维与三维各向异性介质模型进行描述
m(x)={a11(x),a13(x),a33(x),a44(x),a66(x),θ0(x)},2D case (4.14)
模型向量的所有分量都随着空间坐标(x)的变化而变化,将公式(4.14)和(4.15)代入公式(4.1)和(4.8)可以得到二维情形的群速度表达式{Um=(x,θ0,n),m=1,2,3}与三维情形的群速度表达式(m=1,2,3分别对应三种波的模式qP,qSV和qSH),三种波模式的群速度都是空间坐标x、岩石对称轴倾角θ0与方位角以及相慢度方向n的函数。将公式(4.14)和(4.15)定义的弹性参数模型代入二维与三维群速度表达式就可以得到弹性参数与群速度之间的关系式,这种关系式是后续进行射线追踪算法的基础。群速度可以用相慢度方向n或者射线方向r0表示,利用公式(4.13)可以计算出二维情形的射线角
因此射线方向可以写成r0=(sinθ′,cosθ′),群速度{Um=(x,θ0,r0),m=1,2,3}随射线方向r0的变化而变化。
2各向异性介质最短路径射线追踪算法
2.1各向异性介质最短路径射线追踪的原理
如果已知射线路径,那么可以沿着射线路径利用线性积分计算走时,计算表达式如下
式中ds=|dx|是射线段的长度,x是局部射线向量;由公式(4.17)可知走时计算需要确定弹性参数模型m(x)的群速度U(x,θ0,r0)与射线路径R(x)。上一节介绍了怎样建立弹性参数模型m(x)与群速度U(x,θ0,r0)之间的关系,这一节将描述如何确定射线路径段R(x),并进一步利用群速度U(x,θ0,r0)计算走时τ。
根据费马原理,射线路径R(x)是使积分公式(4.17)不变(δτ=0)的一条轨迹,这意味着在空间位置xB处的走时计算表达式如下
式中ΩB是模型内xB点的邻域,r0是射线方向从xA点到xB点的单位向量,公式(4.18)为各向同性介质的射线追踪拓展至各向异性介质提供了可能性。可以明显地看出这种扩展过程需要考虑到基于射线方向r0的群速度与搜索从xA点到xB点最小走时的有效方法。这使得我们可以将Nakanishi和Yamaguchi提出的最短路径方法从各向同性介质扩展至各向异性介质。
2.2各向异性介质最短路径射线追踪的步骤
最短路径射线追踪方法遵循惠更斯原理,根据这种原理各向异性介质的最短路径方法需要包括以下几个步骤。
(1)模型的参数化:以一定的网格单元大小,将二维各向异性介质模型m(x)划分成由主节点组成的网格模型
{m(xk),k=1,2,...,N1,N1=(Nx+1)(Nz+1)} (4.19)
式中Nx、Nz分别代表x和z方向的晶格数目。二维各向异性介质模型由许多矩形晶格组成,每个矩形晶格的四个角被定义为主节点,根据这种规则网格方法,在实际走时计算时须将晶格划分的非常小,但是我们可以采用非规则网格方法:采用较大的晶格单元,同时在矩形晶格的每条边增加次一级节点来提高计算的精度与效率(Bai[82],2007)。将网格化的各向异性介质模型转化成为群速度模型Um(xk,θ0,n),m=1,2,3分别对应qP、qSV和qSH三种波模式,而k∈N1;群速度可以写成射线方向r0的函数Um(xk,θ0,r0),但是qSV波的群速度存在多值现象,在处理qSV波的多值问题时,我们选择群速度的最小值从而确保群速度曲线的一致连续性。群速度模型晶格内部的群速度可以根据晶格四个顶点所在节点的群速度值,利用双线性插值函数或者拉格朗日插值函数进行求取
式中矩形晶格四个节点的位置与网格模型的节点位置{xk,k∈N1}相重合,这些节点被称为初始节点或者主节点。为了提高计算的效率和精度,模型参数化过程中需要增加次节点密度。如果次一级的节点数为N2,那么整个模型的节点数为N1+N2。另外,如果炮点、检波点位置和模型网格节点的位置不重合,那么必须在炮点和检波点位置增加额外的节点,因此整个节点网络包含三类节点:主节点、次节点、炮点和炮检点位置的节点,整个节点网络需要用节点序号l和空间坐标{xl=(xl,zl),l=1,2,...,N}来加以识别。
(2)初始化走时场:对炮点所在节点的走时赋初值"0",其它网格节点的走时值赋为无穷大,例如赋初值"107",从炮点所在节点开始,计算炮点周围节点的走时,再从相邻节点计算其它节点的走时,并记录相邻节点的走时最小值。
(3)最小走时计算方法:一个晶格内两个节点之间的走时可以利用如下公式近似求取,即利用节点i和节点j之间的距离除以两个节点之间的平均群速度
式中Nj是未计算节点j周围的已知走时节点数目,群速度和均由公式(4.20)进行计算求得;由公式(4.20)可以看出如果模型是各向同性介质,则群速度与射线方向无关,公式(4.21)可以简化为标准的最短射线路径走时计算方程。
(4)射线路径的求取:当所有检波点的最小走时计算完成之后,可以拾取检点所在节点的走时,并根据入射节点的序号以反向追踪的形式,求取连接检波点与炮点之间的射线路径,这样就可以得到共炮点道集的走时和射线路径;根据这种走时和路径的求取方法,可以继续求取下一个炮点对应检波点的走时和射线路径,反复利用这一技术路线,便可以完成所有炮点的走时和射线求取。
(5)反射波与衍射波的射线路径求取:如果复杂的各向异性介质模型存在界面,那么地震波在界面就会产生反射波与衍射波。我们可以用一系列密集的节点来代表界面,把界面节点与前面所述的三类节点组合起来,可以得到一个新的网络,它包括主节点、次节点、界面节点和炮、检点节点。网络中大部分晶格的节点分布是一样的,除了含有界面节点的那些晶格,但是它们只占整个模型晶格数目的一小部分。每个界面对应的反射波或者衍射波走时τ的计算方法如下
τ=min{τS(xint)+τR(xint),xint∈Ωint} (4.22)
式中,τS(xint)和τR(xint)分别是炮点S与检波点R到界面节点位置xint的初至走时值,Ωint代表界面节点集合。分别假设炮点S与检波点R为震源点位置,进行两次正演走时计算。τS(xint)在第一次最短路径射线追踪计算时就已经有值了,剩下步骤是将检波点R作为震源点位置再做一次最短路径射线追踪计算求取走时τR(xint),然后根据公式(4.22)寻找相应的界面反射节点,以及反射波的走时与射线路径。值得注意的问题是:在计算反射波与衍射波走时的时候,群速度模型只与界面以上的部分有关,也就是说在计算τS(xint)和τR(xint)时,只截取了模型界面以上部分的节点,这样可以大大提高计算效率。如果一个晶格包含界面节点,那么需要将界面以下节点的群速度值用界面以上节点的群速度值来替代,再用群速度计算公式(4.20)来求取τS(xint)和τR(xint)。同理,这种计算方法也可以运用于转换波走时与射线路径的求取。
3各向异性介质最短路径射线追踪数值模拟
上述内容介绍了各向异性介质最短路径射线追踪算法的基本原理与实现步骤,这里建立倾斜界面模型和弯曲界面模型,分别研究qP-qP、qP-qSV与qP-qSH波在两种模型中的反射波射线路径和走时分布。
3.1倾斜界面模型的射线追踪数值模拟
倾斜界面模型水平方向的长度为800米,垂直方向的深度为500米,第一层介质的弹性模量参数大小为a11=9.08,a13=2.98,a33=7.53,a44=2.27,a66=3.84,相应的对称轴倾角和方位角分别为θ0=90°,第二层介质的弹性模量参数大小为a11=20.3,a13=9.58,a33=22.3,a44=8.35,a66=11.35,相应的对称轴倾角和方位角分别为θ0=90°,第三层介质的弹性模量参数大小为a11=13.86,a13=4.31,a33=10.93,a44=3.31,a66=4.34,相应的对称轴倾角和方位角分别为θ0=0°,炮点位于地表400米处,共有33个检波点,检波点间距为25米。水平与垂直方向的网格间距为5米。下面依次研究qP-qP、qP-qSV与qP-qSH波在倾斜界面模型的反射波射线路径和走时分布。
(1)qP-qP反射波在倾斜界面模型中的射线路径与走时分布,如图2所示。
(2)qP-qSV转换波在倾斜界面模型中的射线路径与走时分布,如图3所示。
(3)qP-qSH转换波在倾斜界面模型中的射线路径与走时分布,如图4所示。
3.2弯曲界面模型的射线追踪数值模拟
弯曲界面模型水平方向的长度为800米,垂直方向的深度为500米,第一层介质的弹性模量参数大小为a11=9.08,a13=2.98,a33=7.53,a44=2.27,a66=3.84,相应的对称轴倾角和方位角分别为θ0=90°,第二层介质的弹性模量参数大小为a11=20.3,a13=9.58,a33=22.3,a44=8.35,a66=11.35,相应的对称轴倾角和方位角分别为θ0=90°,第三层介质的弹性模量参数大小为a11=13.86,a13=4.31,a33=10.93,a44=3.31,a66=4.34,相应的对称轴倾角和方位角分别为θ0=0°,炮点位于地表400米处,共有33个检波点,检波点间距为25米。水平与垂直方向的网格间距为5米。下面依次研究qP-qP、qP-qSV与qP-qSH波在弯曲界面模型的反射波射线路径和走时分布。
(1)qP-qP反射波在弯曲界面模型中的射线路径与走时分布,如图5所示。
(2)qP-qSV转换波在弯曲界面模型中的射线路径与走时分布,如图6所示。
(3)qP-qSH转换波在弯曲界面模型中的射线路径与走时分布,如图7所示。
可以理解的是,对于本领域的普通技术人员来说,可以根据本发明的技术构思做出其它各种相应的改变与变形,而所有这些改变与变形都应属于本发明权利要求的保护范围。
Claims (5)
1.一种各向异性TI介质最短路径射线追踪正演方法,其特征在于,其包括如下步骤:
S1、建立各向异性介质的群速度模型;
S2、通过各向异性介质最短路径射线追踪算法进行各向异性介质最短路径射线追踪;其中,
所述步骤S1包括:
根据群速度模型,可以利用最短路径射线追踪方法计算各向异性介质qP波、qSV波和qSH波的初至走时和射线路径,通过进一步的修改还可以计算反射波、衍射波和转换波的走时和射线路径,TI介质的相速度表达式,即克里斯托弗矩阵的特征值如下:
式中P和Q的表达式为
P=0.5(Q1+Q2)
Q=Q1Q2-Q3 (2)
式中Q1、Q2和Q3的表达式为
上述关系式中,入射角为第一相慢度向量和各向异性介质对称轴倾角之间的夹角,a11,a13,a33,a44,a66是各向异性介质的五个弹性模量参数,方程(1)的c1、c2以及c3分别代表qP波、qSV波和qSH波,计算群速度向量大小的具体表达式如下:
为了避免m与水平、垂直分量标记h,z相混淆,式中波的模式m写成上标形式;
代入ex,ez,可以得到U(m)的表达式:
其中,b、c、d为三维坐标系的三维单位向量,
可以采用如下表达式,对二维与三维各向异性介质模型进行描述
m(x)={a11(x),a13(x),a33(x),a44(x),a66(x),θ0(x)}, (14)
模型向量的所有分量都随着空间坐标(x)的变化而变化,将公式(14)和(15)代入公式(1)和(8)得到二维情形的群速度表达式Um=(x,θ0,n),m=1,2,3与三维情形的群速度表达式m=1,2,3分别对应三种波的模式qP,qSV和qSH,三种波模式的群速度都是空间坐标x、岩石对称轴倾角θ0与方位角以及相慢度方向n的函数;将公式(14)和(15)定义的弹性参数模型代入二维与三维群速度表达式就可以得到弹性参数与群速度之间的关系式;群速度用相慢度方向n或者射线方向r0表示,利用公式(13)可以计算出二维情形的射线角
射线方向写成r0=(sinθ′,cosθ′),群速度{Um=(x,θ0,r0),m=1,2,3}随射线方向r0的变化而变化;
所述步骤S2包括:
S21、模型参数步骤;
S22、对走时场进行初始化;
S23、计算一个晶格内两个节点之间的最小走时;
S24、求取所有炮点的走时和射线;
S25、求取反射波与衍射波的射线路径;
所述步骤S21包括:
以一定的网格单元大小,将二维各向异性介质模型划m(x)分成由主节点组成的网格模型{m(xk),k=1,2,...,N1,N1=(Nx+1)(Nz+1)},其中Nx、Nz分别代表x和z方向的晶格数目;二维各向异性介质模型由多个矩形晶格组成,每个矩形晶格的四个角被定义为主节点;采用较大的晶格单元,同时在矩形晶格的每条边增加次一级节点来提高计算的精度与效率;将网格化的各向异性介质模型转化成为群速度模型Um(xk,θ0,n),m=1,2,3分别对应qP、qSV和qSH三种波模式,而;群速度写成射线方向r0的函数Um(xk,θ0,r0),选择群速度的最小值从而确保群速度曲线的一致连续性;群速度模型晶格内部的群速度根据晶格四个顶点所在节点的群速度值,利用双线性插值函数或者拉格朗日插值函数进行求取其中矩形晶格四个节点的位置与网格模型的节点位置xk相重合,这些节点被称为初始节点或者主节点,分别为矩形晶格四个节点的位置的两个插值,i,j为插值序号,为晶格四个顶点所在节点的群速度值;
在模型参数化过程中增加次节点密度;如果次一级的节点数为N2,那么整个模型的节点数为N1+N2;
此外,如果炮点、检波点位置和模型网格节点的位置不重合,在炮点和检波点位置增加额外的节点,因此整个节点网络包含三类节点:主节点、次节点、炮点和炮检点位置的节点,整个节点网络用节点序号l和空间坐标xl=(xl,zl),l=1,2,...,N来加以识别。
2.如权利要求1所述的各向异性TI介质最短路径射线追踪正演方法,其特征在于,
所述步骤S22包括:
对炮点所在节点的走时赋初值"0",其它网格节点的走时值赋为无穷大,从炮点所在节点开始,计算炮点周围节点的走时,再从相邻节点计算其它节点的走时,并记录相邻节点的走时最小值。
4.如权利要求3所述的各向异性TI介质最短路径射线追踪正演方法,其特征在于,
所述步骤S24包括:
当所有检波点的最小走时计算完成之后,拾取检波点所在节点的走时,并根据入射节点的序号以反向追踪的形式,求取连接检波点与炮点之间的射线路径,得到共炮点道集的走时和射线路径;继续求取下一个炮点对应检波点的走时和射线路径,反复利用这一技术路线,以完成所有炮点的走时和射线求取。
5.如权利要求3所述的各向异性TI介质最短路径射线追踪正演方法,其特征在于,
所述步骤S25包括:
同构一系列密集的节点来代表界面,把界面节点与所述的主节点、次节点、炮点和炮检点位置的节点组合起来,得到一个新的网络,包括主节点、次节点、界面节点和炮、检点节点;每个界面对应的反射波或者衍射波走时τ的计算方法如下:
τ=min{τS(xint)+τR(xint),xint∈Ωint},其中,τS(xint)和τR(xint)分别是炮点S与检波点R到界面节点位置xint的初至走时值,Ωint代表界面节点集合;分别假设炮点S与检波点R为震源点位置,进行两次正演走时计算;τS(xint)在第一次最短路径射线追踪计算时已经有值了,剩下步骤将检波点R作为震源点位置再做一次最短路径射线追踪计算求取走时τR(xint),然后根据步骤S25中公式寻找相应的界面反射节点,以及反射波的走时与射线路径。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711287053.6A CN108303736B (zh) | 2017-12-07 | 2017-12-07 | 各向异性ti介质最短路径射线追踪正演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711287053.6A CN108303736B (zh) | 2017-12-07 | 2017-12-07 | 各向异性ti介质最短路径射线追踪正演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108303736A CN108303736A (zh) | 2018-07-20 |
CN108303736B true CN108303736B (zh) | 2020-11-17 |
Family
ID=62870251
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201711287053.6A Expired - Fee Related CN108303736B (zh) | 2017-12-07 | 2017-12-07 | 各向异性ti介质最短路径射线追踪正演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108303736B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108983290B (zh) * | 2018-09-12 | 2019-05-21 | 中国石油大学(华东) | 一种三维横向各向同性介质中旅行时确定方法及系统 |
CN110261896B (zh) * | 2019-04-26 | 2021-07-20 | 中国石油化工股份有限公司 | 稳定的各向异性ti介质正演模拟方法 |
CN114428292B (zh) * | 2020-09-22 | 2023-06-02 | 中国石油化工股份有限公司 | 近地表速度模型的构建方法和存储介质 |
CN112926235B (zh) * | 2021-01-27 | 2022-04-26 | 浙江大学 | 一种可指定晶格各向异性性能的晶格结构设计方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104360396A (zh) * | 2014-12-04 | 2015-02-18 | 中国海洋石油总公司 | 一种海上井间tti介质三种初至波走时层析成像方法 |
CN104391327A (zh) * | 2014-12-04 | 2015-03-04 | 中国海洋石油总公司 | 一种海上斜井井间地震叠前逆时深度偏移成像方法 |
-
2017
- 2017-12-07 CN CN201711287053.6A patent/CN108303736B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104360396A (zh) * | 2014-12-04 | 2015-02-18 | 中国海洋石油总公司 | 一种海上井间tti介质三种初至波走时层析成像方法 |
CN104391327A (zh) * | 2014-12-04 | 2015-03-04 | 中国海洋石油总公司 | 一种海上斜井井间地震叠前逆时深度偏移成像方法 |
Non-Patent Citations (3)
Title |
---|
Seismic ray tracing in anisotropic media: A modified Newton algorithm for solving highly nonlinear systems;Yanghua Wang;《GEOPHYSICS》;20140228;第79卷(第1期);第T1-T7页 * |
各向异性介质中的最短路径射线追踪;刘媛等;《中国地球物理第二十一届年会论文集》;20050831;第552页 * |
非均匀节点网格TI介质反射波射线追踪研究;黄光南等;《石油物探》;20160131;第55卷(第1期);第26-28页 * |
Also Published As
Publication number | Publication date |
---|---|
CN108303736A (zh) | 2018-07-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108303736B (zh) | 各向异性ti介质最短路径射线追踪正演方法 | |
EP2869096B1 (en) | Systems and methods of multi-scale meshing for geologic time modeling | |
US9846260B2 (en) | Systems and methods to build sedimentary attributes | |
US9759826B2 (en) | System and method for generating an implicit model of geological horizons | |
EP3293552B1 (en) | System and method for editing geological models by switching between volume-based models and surface-based structural models augmented with stratigraphic fiber bundles | |
US7379854B2 (en) | Method of conditioning a random field to have directionally varying anisotropic continuity | |
Cameron et al. | Seismic velocity estimation from time migration | |
CN107843922A (zh) | 一种基于地震初至波和反射波走时联合的层析成像方法 | |
WO2019242045A9 (zh) | 一种虚拟震源二维波前构建地震波走时计算方法 | |
CN109444955B (zh) | 三维地震射线追踪的双线性走时扰动插值方法 | |
CN111948708B (zh) | 一种浸入边界起伏地表地震波场正演模拟方法 | |
CN105549077B (zh) | 基于多级多尺度网格相似性系数计算的微震震源定位方法 | |
CN104360396B (zh) | 一种海上井间tti介质三种初至波走时层析成像方法 | |
CN109459787B (zh) | 基于地震槽波全波形反演的煤矿井下构造成像方法及系统 | |
CN105353406B (zh) | 一种生成角道集的方法和装置 | |
Waheed et al. | An iterative fast sweeping based eikonal solver for tilted orthorhombic media | |
CN109709602A (zh) | 一种远探测声波偏移成像方法、装置及系统 | |
CN106842314A (zh) | 地层厚度的确定方法 | |
Yu et al. | A minimum traveltime ray tracing global algorithm on a triangular net for propagating plane waves | |
CN108983290B (zh) | 一种三维横向各向同性介质中旅行时确定方法及系统 | |
Sun et al. | Joint 3D traveltime calculation based on fast marching method and wavefront construction | |
CN106353801A (zh) | 三维Laplace域声波方程数值模拟方法及装置 | |
CN102841374A (zh) | 基于扫描面正演的伪三维快速微地震正演方法 | |
CN107918144B (zh) | 各向异性介质初至波射线追踪方法及系统 | |
CN115358087A (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 | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20201117 Termination date: 20211207 |