CN111190224A - 一种基于三维地震波反向照明的动态采样全波形反演系统及方法 - Google Patents
一种基于三维地震波反向照明的动态采样全波形反演系统及方法 Download PDFInfo
- Publication number
- CN111190224A CN111190224A CN202010019988.1A CN202010019988A CN111190224A CN 111190224 A CN111190224 A CN 111190224A CN 202010019988 A CN202010019988 A CN 202010019988A CN 111190224 A CN111190224 A CN 111190224A
- Authority
- CN
- China
- Prior art keywords
- energy
- model
- formula
- residual
- velocity
- 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.)
- Granted
Links
- 238000005286 illumination Methods 0.000 title claims abstract description 110
- 238000000034 method Methods 0.000 title claims abstract description 35
- 238000005070 sampling Methods 0.000 title claims abstract description 26
- 238000004364 calculation method Methods 0.000 claims abstract description 16
- SYHGEUNFJIGTRX-UHFFFAOYSA-N methylenedioxypyrovalerone Chemical compound C=1C=C2OCOC2=CC=1C(=O)C(CCC)N1CCCC1 SYHGEUNFJIGTRX-UHFFFAOYSA-N 0.000 claims description 10
- 150000001875 compounds Chemical class 0.000 claims description 9
- 238000001914 filtration Methods 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000002945 steepest descent method Methods 0.000 claims description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 2
- 238000004519 manufacturing process Methods 0.000 abstract description 2
- 238000003384 imaging method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 238000009826 distribution Methods 0.000 description 4
- 238000013508 migration Methods 0.000 description 4
- 230000005012 migration Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000007789 gas Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000006467 substitution reaction Methods 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/303—Analysis for determining velocity profiles or travel times
-
- 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/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
- G01V2210/622—Velocity, density or impedance
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)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种基于三维地震波反向照明的动态采样全波形反演系统及方法。本发明通过输入观测系统文件、高密度炮记录、偏移速度模型、密度模型，建立规则采样的观测系统，求取三维全波形反演速度，计算确定浅层速度模型残差能量最强、最弱区域，进行反向照明得到地表反向照明能量最强、最弱区域，改变炮点数建立新的观测系统，再次计算三维全波形反演速度及浅层速度模型残差，判断是否满足误差条件，若不满足则继续更新速度模型及观测系统；若满足则采用相同方法依次对地表至中层、深层模型进行反向照明，直至整个模型满足误差条件，输出三维全波形反演得到的速度模型。本发明大幅度提高了计算效率，推动了三维全波形反演在生产上的应用。
Description
技术领域
本发明涉及石油勘探领域，具体涉及一种基于三维地震波反向照明的动态采样全波形反演系统及方法。
背景技术
当今世界对石油、天然气等能源的依赖不言而喻，地球物理勘探已成为能源发展的一个重要环节。为了更好地提高勘探的成功概率，针对地下构造目标层成像方法应用的研究，一直都是勘探工作者的重点研究目标。从最原始的射线理论的偏移技术，再到基于波动方程的偏移，可以说偏移成像方法变得越来越完善。但是随着偏移技术的改进以及勘探开发的深入，勘探地质目标日益复杂，对速度的要求也越来越严格。因此，高精度速度模型对高质量地震成像来说至关重要。
全波形反演利用全波形信息，对地震原始炮记录与模型正演炮记录的差值建立目标泛函，利用最优化思想计算出梯度方向，通过不断更新地下介质参数使地震记录差值越来越小，既而实现优化介质参数的目的，实现过程不仅考虑了地震波传播的旅行时等运动学信息，还加入了振幅、相位等动力学信息，能够适应强横向变速介质和各向异性介质的速度反演。作为目前精度最高的一种速度反演工具，全波形反演通过迭代反演得到高精度的地下构造，能够满足目前勘探开发日益复杂的需求，为叠前成像技术提供更准确的速度模型，最终提高油气勘探的成功概率。
发明内容
针对现有技术中高密度海量三维采集地震数据导致全波形反演计算量巨大的现状，本发明提出了一种基于三维地震波反向照明的动态采样全波形反演系统及方法，通过面向地质目标的地震反向照明，建立动态观测系统，减少炮点数量的同时保证了反演效果和精度，大幅度提高了计算效率。
为了实现上述目的，本发明采用如下技术方案：
一种基于三维地震波反向照明的动态采样全波形反演系统，该系统包括输入模块、观测系统建立模块、反向照明模块、全波形反演模块、判断模块、输出模块；
输入模块，被配置为用于输入观测系统文件、高密度炮记录、偏移速度模型、密度模型；
观测系统建立模块，被配置为用于建立规则采样的观测系统及更新观测系统；
反向照明模块，被配置为用于在地层的浅、中、深层速度模型残差能量最强和最弱区域放置炮点进行反向照明；
全波形反演模块，被配置为用于求取三维全波形反演速度及更新速度模型；
判断模块，被配置为用于判断地表反向照明能量最强区域A1和地表反向照明能量最弱区域B1及求取模型残差，判断是否满足误差条件，如果满足则输出反演速度，如果不满足继续更新观测系统，进行迭代，更新反演速度；
输出模块，被配置为用于输出三维全波形反演得到的速度模型。
本发明还提出了一种基于三维地震波反向照明的动态采样全波形反演方法，该方法采用如上所述的基于三维地震波反向照明的动态采样全波形反演系统，包括如下步骤：
步骤1：输入观测系统文件、高密度炮记录、偏移速度模型、密度模型；
步骤2：建立规则采样的观测系统：
xs(x,y,z＝0)＝seekzero(mod(i-1,4),mod(j-1,4),0),i＝1,nx,j＝1,ny (1)；
式中，xs表示震源坐标，x、y、z表示三维坐标，seekzero()表示求取函数值等于零的位置，i、j表示炮点纵横测线的计数单位，mod()为求余函数，表示取余数值；
步骤3：求取三维全波形反演速度：
采用最速下降法的梯度公式：
式中，g表示梯度，v表示背景速度，Tmax表示最大计算时间，t为计算时间，xs为震源坐标，p为正向延拓的三维背景波场，p*表示逆时延拓的残差波场值；
利用时域多尺度方法从低波数到高波数进行速度反演，使用维纳滤波进行多尺度分解，维纳滤波公式为：
式中f表示维纳滤波值，Wt表示目标波形，ω表示角频率，ε为一个小的常数，Wo *(ω)表示波形的共轭转置；
采用如下公式更新速度场：
v(k)＝v(k-1)-αg(k) (4)；
式中v(k+1)和v(k)分别表示第k+1和k次迭代的速度，g(k)表示第k次迭代的梯度，α表示根据抛物线拟合方法计算得到的最优步长；
步骤4：利用如下所示公式，求取地表至浅层速度模型残差能量：
式中，δ(v)1为速度模型残差能量，nz1为浅层深度的网格点数；
步骤5：在速度模型残差能量最强区域放置炮点进行反向照明，得到地表反向照明能量最强区域A1：
式中，
表示地表反向照明能量最强的区域A1，max()表示求取该函数的最大值，seeklarger()表示求取大于该函数值的位置，Ig1表示在速度模型残差能量最强区域放置炮点得到的反向照明能量值，Ig1(x,y,z＝0)表示速度模型残差能量最强区域地表位置上的反向照明能量值；
步骤6：在速度模型残差能量最弱区域放置炮点进行反向照明，得到地表反向照明能量最弱区域B1：
式中，
表示地表反向照明能量最弱区域B1，min()表示求取该函数最小值，seekless()表示求取小于该函数值的位置，Ig2表示在速度模型残差能量最弱区域放置炮点得到的反向照明能量值；Ig2(x,y,z＝0)表示速度模型残差能量最弱区域地表位置上的反向照明能量值；
步骤7：将地表反向照明能量最强区域A1的炮点增加一倍，地表反向照明能量最弱区域B1的炮点减少一半，建立新的观测系统：
式中，
和
分别表示地表反向照明能量最强区域A1和地表反向照明能量最弱区域B1的炮点震源坐标；
步骤8：再次计算三维全波形反演速度；
步骤9：求取地表至浅层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，重复步骤2至步骤8，继续更新观测系统，迭代更新反演速度直至满足误差条件；
如满足误差条件，利用公式(10)求取地表至中层速度模型残差能量，重复步骤5至步骤8，再次计算三维全波形反演速度后，进入步骤10；
式中，δ(v)1为速度模型残差能量，nz2为中层深度的网格点数；
步骤10：求取地表至中层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，重复步骤2至步骤8，继续更新观测系统，迭代更新反演速度直至满足误差条件；
如满足误差条件，利用公式(11)求取地表至深层速度模型残差能量，重复步骤5至步骤8，再次计算三维全波形反演速度后，进入步骤11；
式中，δ(v)1为速度模型残差能量，nz为整个模型垂直深度的网格点数；
步骤11：求取地表至深层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，继续重复步骤2至步骤8，更新观测系统，迭代更新反演速度，直至地表至深层速度残差满足误差条件，得到反演速度，进入步骤12；
如满足误差条件，直接得到反演速度，进入步骤12；
步骤12：输出三维全波形反演得到的速度模型。
优选地，在步骤5中，利用下式求取残差能量最小值：
xmin(x,y,z)＝Seek min(δ(v)1) (12)；
式中，Seek min()表示求取该函数最小值的位置，xmin为δ(v)1最小值所在的坐标点；
将震源点置于xmin处，并采用如下所示公式进行波场延拓：
式中，vx、vy、vz分别为x分量、y分量、z分量的质点速度，p为质点应力，ρ为密度，f为纵波震源；
利用下式求取反向照明能量：
式中，
表示在模型残差能量最强区域放置炮点得到的反向照明能量值，得到地表反向照明能量最强区域A1：
优选地，在步骤6中，利用下式求取残差能量最大值：
xmax(x,y,z)＝Seek max(δ(v)1) (16)；
式中，Seek max()表示求取该函数最大值的位置，xmax为δ(v)1最大值所在的坐标点；
将震源点置于xmax处，并采用如下所示公式进行波场延拓：
利用下式求取反向照明能量：
式中，Ig2表示在模型残差能量最弱区域放置炮点得到的反向照明能量值，得到地表反向照明能量最弱区域B1：
优选地，在步骤7中，对于地表反向照明能量最强区域A1、地表反向照明能量最弱区域B1的炮点震源数量：
如果下式成立，则地表反向照明能量最强区域A1不再增加炮点：
xsA(xs,ys,zs＝0)＝seekzero(mod(i-1,1),mod(j-1,1),0),i＝1,nx,j＝1,ny (18)；
式中，xsA表示需要增加炮点区域的震源坐标；
如果下式成立，则地表反向照明能量最弱区域B1不再减少炮点：
xsB(xs,ys,zs＝0)＝seekzero(mod(i-1,8),mod(j-1,8),0),i＝1,nx,j＝1,ny (19)；
式中，xsB表示需要减少炮点区域的震源坐标。
本发明所带来的有益技术效果：
本发明为了克服高密度海量三维采集地震数据导致全波形反演计算量巨大的难题，提出了一种基于三维地震波反向照明的动态采样全波形反演系统及方法，通过面向地质目标的地震反向照明，建立动态观测系统，在保证反演效果的同时减少炮点数量，并达到了与高密度三维数据全波形反演相同的精度，大幅度提高了计算效率，推动了三维全波形反演生产应用化，为建立更精准的三维速度场提供了地震成像基础。
附图说明
图1为本发明的一种基于三维地震波反向照明的动态采样全波形反演系统的结构示意图。
图2为本发明的一种基于三维地震波反向照明的动态采样全波形反演方法的流程图。
图3为三维断熔体实际速度模型。
图4为初始速度模型。
图5为中层速度模型残差能量最强区域的反向照明图。
图6为中层速度模型残差能量最强区域放置炮点进行反向照明得到的地表能量分布特征图。
图7为中层速度模型残差能量最弱区域的反向照明图。
图8为中层速度模型残差能量最弱区域放置炮点进行反向照明得到的地表能量分布特征图。
图9为需要增加炮点的区域A及需要减少炮点的区域B。
图10为本发明迭代15次得到的地表至浅层反演速度结果。
图11为本发明迭代15次得到的地表至中层反演速度结果。
图12为本发明迭代15次得到的最终反演速度结果。
图13为常规高密度全波形反演方法得到的反演速度结果。
图14为本发明方法与常规高密度全波形反演方法计算时间对比图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
实施例1
一种基于三维地震波反向照明的动态采样全波形反演系统，如图1所示，该系统包括输入模块、观测系统建立模块、反向照明模块、全波形反演模块、判断模块、输出模块；
输入模块，被配置为用于输入观测系统文件、高密度炮记录、偏移速度模型、密度模型；
观测系统建立模块，被配置为用于建立规则采样的观测系统及更新观测系统；
反向照明模块，被配置为用于在地层的浅、中、深层速度模型残差能量最强和最弱区域放置炮点进行反向照明；
全波形反演模块，被配置为用于求取三维全波形反演速度及更新速度模型；
判断模块，被配置为用于判断地表反向照明能量最强区域A1和地表反向照明能量最弱区域B1及求取模型残差，判断是否满足误差条件，如果满足则输出反演速度，如果不满足继续更新观测系统，进行迭代，更新反演速度；
输出模块，被配置为用于输出三维全波形反演得到的速度模型。
实施例2
在上述实施例的基础上，本发明还提出了一种基于三维地震波反向照明的动态采样全波形反演方法，如图2所示，具体包括如下步骤：
步骤1：输入观测系统文件、高密度炮记录、偏移速度模型、密度模型。
步骤2：建立规则采样的观测系统：
xs(xs,ys,zs＝0)＝seekzero(mod(i-1,4),mod(j-1,4),0),i＝1,nx,j＝1,ny (1)；
式中，xs表示震源坐标，x、y、z表示三维坐标，seekzero()表示求取函数值等于零的位置，i、j表示炮点纵横测线的计数单位，mod()为求余函数，表示取余数值。
步骤3：求取三维全波形反演速度：
采用最速下降法的梯度公式：
式中，g表示梯度，v表示背景速度，Tmax表示最大计算时间，t为计算时间，xs为震源坐标，p为正向延拓的三维背景波场，p*表示逆时延拓的残差波场值；
利用时域多尺度方法从低波数到高波数进行速度反演，使用维纳滤波进行多尺度分解，维纳滤波公式为：
其中f表示维纳滤波值，Wt表示目标波形，ω表示角频率，ε为一个小的常数，Wo *(ω)表示波形的共轭转置；
采用如下公式更新速度场：
v(k)＝v(k-1)-αg(k) (4)；
式中v(k+1)和v(k)分别表示第k+1和k次迭代的速度，g(k)表示第k次迭代的梯度，α表示根据抛物线拟合方法计算得到的最优步长。
步骤4：利用如下所示公式，求取地表至浅层速度模型残差能量：
式中，δ(v)1为速度模型残差能量，nz1为浅层深度的网格点数。
步骤5：在速度模型残差能量最强区域放置炮点进行反向照明，得到地表反向照明能量最强区域A1：
利用下式求取残差能量最小值：
xmin(x,y,z)＝Seek min(δ(v)1) (12)；
式中，Seek min()表示求取该函数最小值的位置，xmin为δ(v)1最小值所在的坐标点；
将震源点置于xmin处，并采用如下所示公式进行波场延拓：
式中，vx、vy、vz分别为x分量、y分量、z分量的质点速度，p为质点应力，ρ为密度，f为纵波震源；
利用下式求取反向照明能量：
式中，Ig1表示在模型残差能量最强区域放置炮点得到的反向照明能量值，得到地表反向照明能量最强区域A1：
式中，xA1(x,y,z＝0)表示地表反向照明能量最强的区域A1，max()表示求取该函数的最大值，seeklarger()表示求取大于该函数值的位置，Ig1表示在速度模型残差能量最强区域放置炮点得到的反向照明能量值，Ig1(x,y,z＝0)表示速度模型残差能量最强区域地表位置上的反向照明能量值。
步骤6：在模型残差能量最弱区域放置炮点进行反向照明，得到地表反向照明能量最弱区域B1：
利用下式求取残差能量最大值：
xmax(x,y,z)＝Seek max(δ(v)1) (15)；
式中，Seek max()表示求取该函数最大值的位置，xmax为δ(v)1最大值所在的坐标点；
将震源点置于xmax处，并采用如下所示公式进行波场延拓：
利用如下所示公式求取反向照明能量：
式中，Ig2表示在模型残差能量最弱区域放置炮点得到的反向照明能量值，得到地表反向照明能量最弱区域B1：
其中，
表示地表反向照明能量最强的区域B1，min()表示求取该函数最小值，seekless()表示求取小于该函数值的位置，Ig2表示反向照明能量值。
步骤7：将地表反向照明能量最强区域A1的炮点增加一倍，地表反向照明能量最弱区域B1的炮点减少一半，建立新的观测系统：
式中，
和
分别表示地表反向照明能量最强区域A1和地表反向照明能量最弱区域B1的炮点震源坐标；
如果下式成立，则地表反向照明能量最强区域A1不再增加炮点：
xsA(xs,ys,zs＝0)＝seekzero(mod(i-1,1),mod(j-1,1),0),i＝1,nx,j＝1,ny (18)；
式中，xsA表示需要增加炮点区域的震源坐标；
如果下式成立，则地表反向照明能量最弱区域B1不再减少炮点：
xsB(xs,ys,zs＝0)＝seekzero(mod(i-1,8),mod(j-1,8),0),i＝1,nx,j＝1,ny (19)；
式中，xsB表示需要减少炮点区域的震源坐标。
步骤8：再次计算三维全波形反演速度。
步骤9：求取地表至浅层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，重复步骤2至步骤8，继续更新观测系统，迭代更新反演速度直至满足误差条件；
如满足误差条件，利用如下公式求取地表至中层速度模型残差能量：
式中，δ(v)1为速度模型残差能量，nz2为中层深度的网格点数；
重复步骤5至步骤8，再次计算三维全波形反演速度后，进入步骤10。
步骤10：求取地表至中层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，重复步骤2至步骤8，继续更新观测系统，迭代更新反演速度直至满足误差条件；
如满足误差条件，利用如下公式求取地表至深层速度模型残差能量：
式中，δ(v)1为速度模型残差能量，nz为整个模型垂直深度的网格点数；
重复步骤5至步骤8，再次计算三维全波形反演速度后，进入步骤11。
步骤11：求取地表至深层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，继续重复步骤2至步骤8，更新观测系统，迭代更新反演速度，直至地表至深层速度残差满足误差条件，得到反演速度，进入步骤12；
如满足误差条件，直接得到反演速度，进入步骤12。
步骤12：输出三维全波形反演得到的速度模型。
应用实验
本发明一种基于三维地震波反向照明的动态采样全波形反演系统及方法，应用于三维断熔体实际速度模型数据，取得了理想的计算效果。图3所示为三维断熔体实际速度模型，输入观测系统文件、高密度炮记录、偏移速度模型、密度模型，建立规则采样的观测系统，求取得到如图4所示的初始速度模型，以求取地层中层反演速度为例，利用公式求取地表至中层速度模型残差能量，确定中层速度模型残差能量最强区域A和中层速度模型残差能量最弱区域B，通过对中层速度模型残差能量最强区域放置炮点进行反向照明，得到中层速度模型残差能量最强区域的反向照明图(如图5所示)和地表能量分布特征(如图6所示)；通过对中层速度模型残差能量最弱区域放置炮点进行反向照明，得到中层模型残差能量最弱区域的反向照明图(如图7所示)和地表能量分布特征(如图8所示)；将地表反向照明能量最强区域A炮点数量增加一倍，地表反向照明能量最弱区B炮点数量减少一半，如图9所示，建立新的观测系统；再次计算三维全波形反演速度；求取地表至中层速度模型残差并判断是否满足误差条件，经过15次迭代，得到如图11所示的地表至中层反演速度结果。运用相同方法对浅层速度模型进行反演，经过多次迭代计算，得到地表至浅层反演速度结果，如图10所示；运用相同方法对深层速度模型进行反演，经过多次迭代计算，得到最终反演速度结果，如图12所示。作为对比，计算得到如图13所示的常规高密度全波形反演方法得到的反演速度结果，从图中对比可得，本发明方法得到的反演速度结果与常规全波形反演方法结果类似，但是本实施例使用的炮点数仅为常规高密度方法炮点数的1/4，同时，图14所示为本发明方法与常规高密度全波形反演方法计算时间对比图，对比得到本发明方法的计算时间约为常规方法的1/3，计算速度大幅度提升。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。
Claims (5)
1.一种基于三维地震波反向照明的动态采样全波形反演系统，其特征在于，该系统包括输入模块、观测系统建立模块、反向照明模块、全波形反演模块、判断模块、输出模块；
输入模块，被配置为用于输入观测系统文件、高密度炮记录、偏移速度模型、密度模型；
观测系统建立模块，被配置为用于建立规则采样的观测系统及更新观测系统；
反向照明模块，被配置为用于在地层的浅、中、深层速度模型残差能量最强和最弱区域放置炮点进行反向照明；
全波形反演模块，被配置为用于求取三维全波形反演速度及更新速度模型；
判断模块，被配置为用于判断地表反向照明能量最强区域A1和地表反向照明能量最弱区域B1及求取模型残差，判断是否满足误差条件，如果满足则输出反演速度，如果不满足继续更新观测系统，进行迭代，更新反演速度；
输出模块，被配置为输出三维全波形反演得到的速度模型。
2.一种基于三维地震波反向照明的动态采样全波形反演方法，其特征在于，采用如权利要求1所述的基于三维地震波反向照明的动态采样全波形反演系统，包括如下步骤：
步骤1：输入观测系统文件、高密度炮记录、偏移速度模型、密度模型；
步骤2：建立规则采样的观测系统：
xs(x,y,z＝0)＝seekzero(mod(i-1,4),mod(j-1,4),0),i＝1,nx,j＝1,ny (1)；
式中，xs表示震源坐标，x、y、z表示三维坐标，seekzero()表示求取函数值等于零的位置，i、j表示炮点纵横测线的计数单位，mod()为求余函数，表示取余数值；
步骤3：求取三维全波形反演速度：
采用最速下降法的梯度公式：
式中，g表示梯度，v表示背景速度，Tmax表示最大计算时间，t为计算时间，xs为震源坐标，p为正向延拓的三维背景波场，p*表示逆时延拓的残差波场值；
利用时域多尺度方法从低波数到高波数进行速度反演，使用维纳滤波进行多尺度分解，维纳滤波公式为：
式中f表示维纳滤波值，Wt表示目标波形，ω表示角频率，ε为一个小的常数，
表示波形的共轭转置；
采用如下所示公式更新速度场：
v(k)＝v(k-1)-αg(k) (4)；
式中v(k+1)和v(k)分别表示第k+1和k次迭代的速度，g(k)表示第k次迭代的梯度，α表示根据抛物线拟合方法计算得到的最优步长；
步骤4：利用如下所示公式，求取地表至浅层速度模型残差能量：
式中，δ(v)1为速度模型残差能量，nz1为浅层深度的网格点数；
步骤5：在速度模型残差能量最强区域放置炮点进行反向照明，得到地表反向照明能量最强区域A1：
式中，
表示地表反向照明能量最强的区域A1，max()表示求取该函数的最大值，seeklarger()表示求取大于该函数值的位置，Ig1表示在速度模型残差能量最强区域放置炮点得到的反向照明能量值，Ig1(x,y,z＝0)表示速度模型残差能量最强区域地表位置上的反向照明能量值；
步骤6：在速度模型残差能量最弱区域放置炮点进行反向照明，得到地表反向照明能量最弱区域B1：
式中，
表示地表反向照明能量最弱区域B1，min()表示求取该函数最小值，seekless()表示求取小于该函数值的位置，Ig2表示在速度模型残差能量最弱区域放置炮点得到的反向照明能量值；Ig2(x,y,z＝0)表示速度模型残差能量最弱区域地表位置上的反向照明能量值；
步骤7：将地表反向照明能量最强区域A1的炮点增加一倍，地表反向照明能量最弱区域B1的炮点减少一半，建立新的观测系统：
式中，
和
分别表示地表反向照明能量最强区域A1和地表反向照明能量最弱区域B1的炮点震源坐标；
步骤8：再次计算三维全波形反演速度；
步骤9：求取地表至浅层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，重复步骤2至步骤8，继续更新观测系统，迭代更新反演速度直至满足误差条件；
如满足误差条件，利用公式(10)求取地表至中层速度模型残差能量，重复步骤5至步骤8，再次计算三维全波形反演速度后，进入步骤10；
式中，δ(v)1为速度模型残差能量，nz2为中层深度的网格点数；
步骤10：求取地表至中层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，重复步骤2至步骤8，继续更新观测系统，迭代更新反演速度直至满足误差条件；
如满足误差条件，利用公式(11)求取地表至深层速度模型残差能量，重复步骤5至步骤8，再次计算三维全波形反演速度后，进入步骤11；
式中，δ(v)1为速度模型残差能量，nz为整个模型垂直深度的网格点数；
步骤11：求取地表至深层速度模型残差，判断速度模型残差是否满足误差条件：
如不满足误差条件，继续重复步骤2至步骤8，更新观测系统，迭代更新反演速度，直至地表至深层速度残差满足误差条件，得到反演速度，进入步骤12；
如满足误差条件，直接得到反演速度，进入步骤12；
步骤12：输出三维全波形反演得到的速度模型。
3.根据权利要求2所述的一种基于三维地震波反向照明的动态采样全波形反演方法，其特征在于，在步骤5中，利用下式求取残差能量最小值：
xmin(x,y,z)＝Seekmin(δ(v)1) (12)；
式中，Seekmin()表示求取该函数最小值的位置，xmin为δ(v)1最小值所在的坐标点；
将震源点置于xmin处，并采用如下所示公式进行波场延拓：
式中，vx、vy、vz分别为x分量、y分量、z分量的质点速度，p为质点应力，ρ为密度，f为纵波震源；
利用下式求取反向照明能量：
式中，
表示在模型残差能量最强区域放置炮点得到的反向照明能量值，得到地表反向照明能量最强区域A1：
4.根据权利要求2所述的一种基于三维地震波反向照明的动态采样全波形反演方法，其特征在于，在步骤6中，利用下式求取残差能量最大值：
xmax(x,y,z)＝Seekmax(δ(v)1) (15)；
式中，Seekmax()表示求取该函数最大值的位置，xmax为δ(v)1最大值所在的坐标点；
将震源点置于xmax处，并采用如下所示公式进行波场延拓：
利用下式求取反向照明能量：
式中，Ig2表示在模型残差能量最弱区域放置炮点得到的反向照明能量值，得到地表反向照明能量最弱区域B1：
5.根据权利要求2所述的一种基于三维地震波反向照明的动态采样全波形反演方法，其特征在于，在步骤7中，对于地表反向照明能量最强区域A1、地表反向照明能量最弱区域B1的炮点数量：
如果下式成立，则地表反向照明能量最强区域A1不再增加炮点：
xsA(xs,ys,zs＝0)＝seekzero(mod(i-1,1),mod(j-1,1),0),i＝1,nx,j＝1,ny (18)；
式中，xsA表示需要增加炮点区域的震源坐标；
如果下式成立，则地表反向照明能量最弱区域B1不再减少炮点：
xsB(xs,ys,zs＝0)＝seekzero(mod(i-1,8),mod(j-1,8),0),i＝1,nx,j＝1,ny (19)；
式中，xsB表示需要减少炮点区域的震源坐标。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010019988.1A CN111190224B (zh) | 2020-01-09 | 2020-01-09 | 一种基于三维地震波反向照明的动态采样全波形反演系统及方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010019988.1A CN111190224B (zh) | 2020-01-09 | 2020-01-09 | 一种基于三维地震波反向照明的动态采样全波形反演系统及方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111190224A true CN111190224A (zh) | 2020-05-22 |
| CN111190224B CN111190224B (zh) | 2022-03-25 |
Family
ID=70705267
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010019988.1A Active CN111190224B (zh) | 2020-01-09 | 2020-01-09 | 一种基于三维地震波反向照明的动态采样全波形反演系统及方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111190224B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113341461A (zh) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | 地震速度预测方法、装置及服务器 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103499835A (zh) * | 2013-10-13 | 2014-01-08 | 中国石油集团西北地质研究所 | 初至波波形反演近地表速度模型方法 |
| CN105572728A (zh) * | 2014-10-16 | 2016-05-11 | 中国石油化工股份有限公司 | 基于最小二乘目标泛函的反向照明速度反演方法 |
| CN106526674A (zh) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | 一种三维全波形反演能量加权梯度预处理方法 |
| CN107783190A (zh) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | 一种最小二乘逆时偏移梯度更新方法 |
| US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
| US20210223424A1 (en) * | 2018-06-08 | 2021-07-22 | Total Se | Method for generating an image of a subsurface of an area of interest from seismic data |
-
2020
- 2020-01-09 CN CN202010019988.1A patent/CN111190224B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103499835A (zh) * | 2013-10-13 | 2014-01-08 | 中国石油集团西北地质研究所 | 初至波波形反演近地表速度模型方法 |
| CN105572728A (zh) * | 2014-10-16 | 2016-05-11 | 中国石油化工股份有限公司 | 基于最小二乘目标泛函的反向照明速度反演方法 |
| CN106526674A (zh) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | 一种三维全波形反演能量加权梯度预处理方法 |
| US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
| CN107783190A (zh) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | 一种最小二乘逆时偏移梯度更新方法 |
| US20210223424A1 (en) * | 2018-06-08 | 2021-07-22 | Total Se | Method for generating an image of a subsurface of an area of interest from seismic data |
Non-Patent Citations (2)
| Title |
|---|
| 庞新明等: "基于坡印廷矢量的地震照明分析", 《热带海洋学报》 * |
| 庞新明等: "照明分析方法的研究进展综述", 《地球物理学进展》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113341461A (zh) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | 地震速度预测方法、装置及服务器 |
| CN113341461B (zh) * | 2021-06-10 | 2023-09-01 | 中国石油大学(北京) | 地震速度预测方法、装置及服务器 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111190224B (zh) | 2022-03-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113740901B (zh) | 基于复杂起伏地表的陆上地震数据全波形反演方法及装置 | |
| CN107894618B (zh) | 一种基于模型平滑算法的全波形反演梯度预处理方法 | |
| CN102236108B (zh) | 一种磁性地表三维地形改正方法 | |
| CN102243678A (zh) | 一种基于沉积动力学反演的储集砂体分析方法 | |
| CN113031068B (zh) | 一种基于反射系数精确式的基追踪叠前地震反演方法 | |
| CN105549080A (zh) | 一种基于辅助坐标系的起伏地表波形反演方法 | |
| CN111290019B (zh) | 一种应用于最小二乘逆时偏移的l-bfgs初始矩阵求取方法 | |
| CN102901985A (zh) | 一种适用于起伏地表的深度域层速度修正方法 | |
| CN110531410A (zh) | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 | |
| CN106291682A (zh) | 一种基于基追踪方法的叠后声波阻抗反演方法 | |
| CN111580163B (zh) | 一种基于非单调搜索技术的全波形反演方法及系统 | |
| CN110737018B (zh) | Vsp地震资料各向异性建模方法 | |
| CN110361788B (zh) | 一种空-地联合三维重力数据特征分析及密度反演方法 | |
| CN111190224B (zh) | 一种基于三维地震波反向照明的动态采样全波形反演系统及方法 | |
| CN109239776B (zh) | 一种地震波传播正演模拟方法和装置 | |
| CN104597489A (zh) | 一种震源子波优化设置方法和装置 | |
| CN111273346B (zh) | 去除沉积背景的方法、装置、计算机设备及可读存储介质 | |
| CN111859251A (zh) | 一种基于pde的磁测数据等效源上延拓与下延拓方法 | |
| CN111175822B (zh) | 改进直接包络反演与扰动分解的强散射介质反演方法 | |
| Bychkov et al. | Near-surface correction on seismic and gravity data | |
| CN109901221B (zh) | 一种基于动校正速度参数的地震资料各向异性建模方法 | |
| CN115880455A (zh) | 基于深度学习的三维智能插值方法 | |
| Martyshko et al. | On solutions of forward and inverse problem for potential geophysical fields: Gravity inversion for Urals region | |
| CN113970789A (zh) | 全波形反演方法、装置、存储介质及电子设备 | |
| CN110568497A (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 |