CN108303736B - 各向异性ti介质最短路径射线追踪正演方法 - Google Patents

各向异性ti介质最短路径射线追踪正演方法 Download PDF

Info

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
Application number
CN201711287053.6A
Other languages
English (en)
Other versions
CN108303736A (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.)
East China Institute of Technology
Original Assignee
East China Institute of Technology
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 East China Institute of Technology filed Critical East China Institute of Technology
Priority to CN201711287053.6A priority Critical patent/CN108303736B/zh
Publication of CN108303736A publication Critical patent/CN108303736A/zh
Application granted granted Critical
Publication of CN108303736B publication Critical patent/CN108303736B/zh
Expired - Fee Related 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. 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介质最短路径射线追踪正演方法
技术领域
本发明涉及各向异性介质射线追踪正演技术领域，特别涉及一种各向异性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}与三维情形的群速度表达式
Figure BDA0001498654640000021
其中m＝1,2,3分别对应三种波的模式qP,qSV和qSH，三种波模式的群速度都是空间坐标x、岩石对称轴倾角θ0与方位角
Figure BDA0001498654640000022
以及相慢度方向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)，选择群速度的最小值从而确保群速度曲线的一致连续性；群速度模型晶格内部的群速度根据晶格四个顶点所在节点的群速度值，利用双线性插值函数或者拉格朗日插值函数进行求取
Figure BDA0001498654640000023
其中矩形晶格四个节点的位置
Figure BDA0001498654640000031
与网格模型的节点位置{xk,k∈N1}相重合，这些节点被称为初始节点或者主节点；
在模型参数化过程中增加次节点密度；如果次一级的节点数为N2，那么整个模型的节点数为N1+N2；此外，如果炮点、检波点位置和模型网格节点的位置不重合，在炮点和检波点位置增加额外的节点，因此整个节点网络包含三类节点：主节点、次节点、炮点和炮检点位置的节点，整个节点网络用节点序号l和空间坐标{xl＝(xl,zl),l＝1,2,...,N}来加以识别。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中，
所述步骤S22包括：
对炮点所在节点的走时赋初值"0"，其它网格节点的走时值赋为无穷大，从炮点所在节点开始，计算炮点周围节点的走时，再从相邻节点计算其它节点的走时，并记录相邻节点的走时最小值。
在本发明所述的各向异性TI介质最短路径射线追踪正演方法中，
所述步骤S23包括：
将一个晶格内两个节点之间的走时通过如下公式近似求取，即利用节点i和节点j之间的距离除以两个节点之间的平均群速度
Figure BDA0001498654640000032
其中Nj是未计算节点j周围的已知走时节点数目，群速度
Figure BDA0001498654640000033
和
Figure BDA0001498654640000034
均由步骤S21中公式求得。
在本发明所述的各向异性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介质的相速度表达式，即克里斯托弗矩阵的特征值如下：
Figure BDA0001498654640000051
Figure BDA0001498654640000052
式中P和Q的表达式为
Figure BDA0001498654640000053
式中Q1、Q2和Q3的表达式为
Figure BDA00014986546400000512
Figure BDA00014986546400000513
Figure BDA00014986546400000514
上述关系式中，角度
Figure BDA00014986546400000515
为相慢度向量
Figure BDA0001498654640000054
和各向异性介质对称轴倾角
Figure BDA0001498654640000055
之间的夹角，{a11,a13,a33,a44,a66}是各向异性介质的五个弹性模量参数，方程(4.1)的下标1、2和3分别代表qP波、qSV波和qSH波。计算群速度向量大小的具体表达式如下：
Figure BDA0001498654640000056
式中m＝1,2,3分别对应qP波、qSV波和qSH波。群速度向量的水平和垂直分量
Figure BDA0001498654640000057
和
Figure BDA0001498654640000058
可以写成如下形式
Figure BDA0001498654640000059
为了避免m与水平、垂直分量标记h,z相混淆，式中波的模式m写成上标形式。
公式(4.4)和(4.5)的偏导数项
Figure BDA00014986546400000510
可以利用公式(4.1)-(4.3)对角度
Figure BDA00014986546400000516
进行求偏导数得到，它们的具体形式如下
Figure BDA00014986546400000511
Figure BDA0001498654640000061
式中
Figure BDA0001498654640000062
和
Figure BDA0001498654640000063
的表达式为
Figure BDA0001498654640000064
Figure BDA0001498654640000065
由公式(4.1)-(4.7)可以看出VTI介质的群速度与五个弹性模量参数{a11,a13,a33,a44,a66}以及相慢度向量
Figure BDA00014986546400000622
有关。
为了使公式(4.1)-(4.7)能够应用于一般TI介质，可以利用标准三角关系(坐标旋转)得到入射角的余弦表达式。入射角
Figure BDA00014986546400000623
是指相慢度方向与TI介质对称轴的夹角，而相慢度方向通常是由计算坐标系给出，不是由对称轴坐标系给出。如果知道相慢度方向
Figure BDA0001498654640000066
和各向异性介质的对称轴倾角
Figure BDA0001498654640000067
那么入射角
Figure BDA00014986546400000624
可以表示为
Figure BDA00014986546400000625
然后，才能应用公式(4.1)计算三个波的相速度，所以要由
Figure BDA0001498654640000068
和
Figure BDA0001498654640000069
来求出
Figure BDA00014986546400000626
Figure BDA00014986546400000627
Figure BDA00014986546400000628
Figure BDA00014986546400000610
Figure BDA00014986546400000611
将岩石的对称轴从垂直状态旋转至
Figure BDA00014986546400000612
可以得到一般TI介质群速度向量的三个分量形式，它的具体计算方法为：已知岩石对称轴的方向向量为
Figure BDA00014986546400000613
由它可以得到与ez垂直的x轴方向向量ex的表达式
Figure BDA00014986546400000614
然后，由ey＝(ez×ex)/|(ez×ex)|得到
Figure BDA00014986546400000615
程序里就是根据这个原理来得到岩石坐标系的。因此我们可以得到群速度的矢量形式为
Figure BDA00014986546400000616
代入ex,ez，可以得到U(m)的表达式：
Figure BDA00014986546400000617
这样，就得到了
Figure BDA00014986546400000618
在计算机坐标系的表达式
Figure BDA00014986546400000619
Figure BDA00014986546400000620
Figure BDA00014986546400000621
二维各向异性介质模型观测方位角
Figure BDA0001498654640000071
与介质的对称轴方位角
Figure BDA0001498654640000072
均为零，而且介质的对称轴倾角θ0、相慢度方向倾角θ与射线入射角
Figure BDA00014986546400000710
之间的关系可以进一步简化为
Figure BDA00014986546400000711
因此二维各向异性介质群速度的表达式可以写成
Figure BDA0001498654640000073
可以采用如下表达式，对二维与三维各向异性介质模型进行描述
m(x)＝{a11(x),a13(x),a33(x),a44(x),a66(x),θ0(x)}，2D case (4.14)
Figure BDA0001498654640000074
模型向量的所有分量都随着空间坐标(x)的变化而变化，将公式(4.14)和(4.15)代入公式(4.1)和(4.8)可以得到二维情形的群速度表达式{Um＝(x,θ0,n),m＝1,2,3}与三维情形的群速度表达式
Figure BDA0001498654640000075
(m＝1,2,3分别对应三种波的模式qP,qSV和qSH)，三种波模式的群速度都是空间坐标x、岩石对称轴倾角θ0与方位角
Figure BDA0001498654640000076
以及相慢度方向n的函数。将公式(4.14)和(4.15)定义的弹性参数模型代入二维与三维群速度表达式就可以得到弹性参数与群速度之间的关系式，这种关系式是后续进行射线追踪算法的基础。群速度可以用相慢度方向n或者射线方向r0表示，利用公式(4.13)可以计算出二维情形的射线角
Figure BDA0001498654640000077
因此射线方向可以写成r0＝(sinθ′,cosθ′)，群速度{Um＝(x,θ0,r0),m＝1,2,3}随射线方向r0的变化而变化。
2各向异性介质最短路径射线追踪算法
2.1各向异性介质最短路径射线追踪的原理
如果已知射线路径，那么可以沿着射线路径利用线性积分计算走时，计算表达式如下
Figure BDA0001498654640000078
式中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处的走时计算表达式如下
Figure BDA0001498654640000079
式中Ω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波的多值问题时，我们选择群速度的最小值从而确保群速度曲线的一致连续性。群速度模型晶格内部的群速度可以根据晶格四个顶点所在节点的群速度值，利用双线性插值函数或者拉格朗日插值函数进行求取
Figure BDA0001498654640000081
式中矩形晶格四个节点的位置
Figure BDA0001498654640000082
与网格模型的节点位置{xk,k∈N1}相重合，这些节点被称为初始节点或者主节点。为了提高计算的效率和精度，模型参数化过程中需要增加次节点密度。如果次一级的节点数为N2，那么整个模型的节点数为N1+N2。另外，如果炮点、检波点位置和模型网格节点的位置不重合，那么必须在炮点和检波点位置增加额外的节点，因此整个节点网络包含三类节点：主节点、次节点、炮点和炮检点位置的节点，整个节点网络需要用节点序号l和空间坐标{xl＝(xl,zl),l＝1,2,...,N}来加以识别。
(2)初始化走时场：对炮点所在节点的走时赋初值"0"，其它网格节点的走时值赋为无穷大，例如赋初值"107"，从炮点所在节点开始，计算炮点周围节点的走时，再从相邻节点计算其它节点的走时，并记录相邻节点的走时最小值。
(3)最小走时计算方法：一个晶格内两个节点之间的走时可以利用如下公式近似求取，即利用节点i和节点j之间的距离除以两个节点之间的平均群速度
Figure BDA0001498654640000091
式中Nj是未计算节点j周围的已知走时节点数目，群速度
Figure BDA0001498654640000092
和
Figure BDA0001498654640000093
均由公式(4.20)进行计算求得；由公式(4.20)可以看出如果模型是各向同性介质，则群速度与射线方向
Figure BDA0001498654640000094
无关，公式(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°,
Figure BDA0001498654640000101
第二层介质的弹性模量参数大小为a11＝20.3,a13＝9.58,a33＝22.3,a44＝8.35,a66＝11.35，相应的对称轴倾角和方位角分别为θ0＝90°,
Figure BDA0001498654640000102
第三层介质的弹性模量参数大小为a11＝13.86,a13＝4.31,a33＝10.93,a44＝3.31,a66＝4.34，相应的对称轴倾角和方位角分别为θ0＝0°,
Figure BDA0001498654640000103
炮点位于地表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°,
Figure BDA0001498654640000104
第二层介质的弹性模量参数大小为a11＝20.3,a13＝9.58,a33＝22.3,a44＝8.35,a66＝11.35，相应的对称轴倾角和方位角分别为θ0＝90°,
Figure BDA0001498654640000105
第三层介质的弹性模量参数大小为a11＝13.86,a13＝4.31,a33＝10.93,a44＝3.31,a66＝4.34，相应的对称轴倾角和方位角分别为θ0＝0°,
Figure BDA0001498654640000106
炮点位于地表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包括:
所述群速度模型包括二维情形和三维情形的群速度,二维情形的群速度表达式Um＝(x，θ0，n)，m＝1，2，3与三维情形的群速度表达式
Figure FDA0002510397910000011
其中，m＝1，2，3分别对应三种波的模式qP，qSV和qSH，三种波模式的群速度都是空间坐标x、岩石对称轴倾角θ0与方位角
Figure FDA0002510397910000012
以及相慢度方向n的函数；
根据群速度模型，可以利用最短路径射线追踪方法计算各向异性介质qP波、qSV波和qSH波的初至走时和射线路径，通过进一步的修改还可以计算反射波、衍射波和转换波的走时和射线路径,TI介质的相速度表达式，即克里斯托弗矩阵的特征值如下：
Figure FDA0002510397910000013
Figure FDA0002510397910000014
式中P和Q的表达式为
P＝0.5(Q1+Q2)
Q＝Q1Q2-Q3 (2)
式中Q1、Q2和Q3的表达式为
Figure FDA0002510397910000015
Figure FDA0002510397910000016
Figure FDA0002510397910000017
上述关系式中，入射角
Figure FDA00025103979100000110
为第一相慢度向量
Figure FDA0002510397910000018
和各向异性介质对称轴倾角
Figure FDA0002510397910000019
之间的夹角，a11，a13，a33，a44，a66是各向异性介质的五个弹性模量参数，方程(1)的c1、c2以及c3分别代表qP波、qSV波和qSH波，计算群速度向量大小的具体表达式如下：
Figure FDA0002510397910000021
式中m＝1,2,3分别对应qP波、qSV波和qSH波，群速度向量的水平和垂直分量
Figure FDA0002510397910000022
和
Figure FDA0002510397910000023
如下
Figure FDA0002510397910000024
为了避免m与水平、垂直分量标记h,z相混淆，式中波的模式m写成上标形式；
公式(4)和(5)的偏导数项
Figure FDA0002510397910000025
利用公式(1)-(3)对角度
Figure FDA0002510397910000026
进行求偏导数得到，具体形式如下
Figure FDA0002510397910000027
Figure FDA0002510397910000028
式中
Figure FDA0002510397910000029
和
Figure FDA00025103979100000210
的表达式为
Figure FDA00025103979100000211
Figure FDA00025103979100000212
由公式(1)-(7)得出TI介质的群速度与五个弹性模量参数a11，a13，a33，a44，a66以及第二相慢度向量
Figure FDA00025103979100000213
有关；
入射角
Figure FDA00025103979100000214
是指相慢度方向与TI介质对称轴的夹角，第一相慢度向量
Figure FDA00025103979100000215
和各向异性介质的对称轴倾角
Figure FDA00025103979100000216
那么入射角
Figure FDA00025103979100000217
可以表示为
Figure FDA00025103979100000218
并应用公式(1)计算三个波的相速度，所以要由
Figure FDA00025103979100000219
和
Figure FDA00025103979100000220
来求出
Figure FDA00025103979100000221
Figure FDA00025103979100000222
Figure FDA00025103979100000223
Figure FDA00025103979100000224
Figure FDA00025103979100000225
将岩石的对称轴从垂直状态旋转至
Figure FDA00025103979100000226
得到一般TI介质群速度向量的三个分量形式，其具体计算方法为：已知岩石对称轴的方向向量为
Figure FDA00025103979100000227
可以得到与ez垂直的x轴方向向量ex的表达式
Figure FDA00025103979100000228
Figure FDA0002510397910000031
然后，由ey＝(ez×ex)/|(ez×ex)|得到
Figure FDA0002510397910000032
得到的群速度的矢量形式为
Figure FDA0002510397910000033
代入ex,ez，可以得到U(m)的表达式：
Figure FDA0002510397910000034
其中，b、c、d为三维坐标系的三维单位向量，
得到
Figure FDA0002510397910000035
的坐标系表达式
Figure FDA0002510397910000036
Figure FDA0002510397910000037
Figure FDA0002510397910000038
二维各向异性介质模型观测方位角
Figure FDA0002510397910000039
与介质的对称轴方位角
Figure FDA00025103979100000310
均为零，而且介质的对称轴倾角θ0、相慢度方向倾角θ与射线入射角
Figure FDA00025103979100000311
之间的关系可以进一步简化为
Figure FDA00025103979100000312
二维各向异性介质群速度的表达式可以写成
Figure FDA00025103979100000313
可以采用如下表达式，对二维与三维各向异性介质模型进行描述
m(x)＝{a11(x),a13(x),a33(x),a44(x),a66(x),θ0(x)}， (14)
Figure FDA00025103979100000314
模型向量的所有分量都随着空间坐标(x)的变化而变化，将公式(14)和(15)代入公式(1)和(8)得到二维情形的群速度表达式Um＝(x，θ0，n)，m＝1，2，3与三维情形的群速度表达式
Figure FDA00025103979100000315
m＝1,2,3分别对应三种波的模式qP,qSV和qSH，三种波模式的群速度都是空间坐标x、岩石对称轴倾角θ0与方位角
Figure FDA00025103979100000316
以及相慢度方向n的函数；将公式(14)和(15)定义的弹性参数模型代入二维与三维群速度表达式就可以得到弹性参数与群速度之间的关系式；群速度用相慢度方向n或者射线方向r0表示，利用公式(13)可以计算出二维情形的射线角
Figure FDA00025103979100000317
射线方向写成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)，选择群速度的最小值从而确保群速度曲线的一致连续性；群速度模型晶格内部的群速度根据晶格四个顶点所在节点的群速度值，利用双线性插值函数或者拉格朗日插值函数进行求取
Figure FDA0002510397910000041
其中矩形晶格四个节点的位置
Figure FDA0002510397910000042
与网格模型的节点位置xk相重合，这些节点被称为初始节点或者主节点,
Figure FDA0002510397910000043
分别为矩形晶格四个节点的位置的两个插值，i,j为插值序号，
Figure FDA0002510397910000044
为晶格四个顶点所在节点的群速度值；
在模型参数化过程中增加次节点密度；如果次一级的节点数为N2，那么整个模型的节点数为N1+N2；
此外，如果炮点、检波点位置和模型网格节点的位置不重合，在炮点和检波点位置增加额外的节点，因此整个节点网络包含三类节点：主节点、次节点、炮点和炮检点位置的节点，整个节点网络用节点序号l和空间坐标xl＝(xl，zl)，l＝1，2，...，N来加以识别。
2.如权利要求1所述的各向异性TI介质最短路径射线追踪正演方法，其特征在于，
所述步骤S22包括：
对炮点所在节点的走时赋初值"0"，其它网格节点的走时值赋为无穷大，从炮点所在节点开始，计算炮点周围节点的走时，再从相邻节点计算其它节点的走时，并记录相邻节点的走时最小值。
3.如权利要求2所述的各向异性TI介质最短路径射线追踪正演方法，其特征在于，
所述步骤S23包括：
将一个晶格内两个节点之间的走时通过如下公式近似求取，即利用节点g和节点h之间的距离除以两个节点之间的平均群速度
Figure FDA0002510397910000051
其中Nh是未计算节点h周围的已知走时节点数目，群速度
Figure FDA0002510397910000052
和
Figure FDA0002510397910000053
均由步骤S21中公式求得，其中，xg，xh分别是节点g,h的节点坐标值，
Figure FDA0002510397910000054
分别节点g,h的对称轴倾角，
Figure FDA0002510397910000055
为节点g,h的射线方向值，
Figure FDA0002510397910000056
分别是利用节点g和节点h到参考节点之间的走时。
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中公式寻找相应的界面反射节点，以及反射波的走时与射线路径。
CN201711287053.6A 2017-12-07 2017-12-07 各向异性ti介质最短路径射线追踪正演方法 Expired - Fee Related CN108303736B (zh)

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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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 中国海洋石油总公司 一种海上斜井井间地震叠前逆时深度偏移成像方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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