CN111624648A - 一种起伏地表变偏移距vsp积分叠前深度偏移方法 - Google Patents

一种起伏地表变偏移距vsp积分叠前深度偏移方法 Download PDF

Info

Publication number
CN111624648A
CN111624648A CN202010502932.1A CN202010502932A CN111624648A CN 111624648 A CN111624648 A CN 111624648A CN 202010502932 A CN202010502932 A CN 202010502932A CN 111624648 A CN111624648 A CN 111624648A
Authority
CN
China
Prior art keywords
imaging
point
coordinate
depth
ktr
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
Application number
CN202010502932.1A
Other languages
English (en)
Other versions
CN111624648B (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.)
Optical Science and Technology Chengdu Ltd of CNPC
Original Assignee
Optical Science and Technology Chengdu Ltd of CNPC
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 Optical Science and Technology Chengdu Ltd of CNPC filed Critical Optical Science and Technology Chengdu Ltd of CNPC
Priority to CN202010502932.1A priority Critical patent/CN111624648B/zh
Publication of CN111624648A publication Critical patent/CN111624648A/zh
Application granted granted Critical
Publication of CN111624648B publication Critical patent/CN111624648B/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. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • 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/308—Time lapse or 4D effects, e.g. production related effects to the formation
    • 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/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G01V2210/512—Pre-stack
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/70—Other details related to processing
    • G01V2210/74—Visualisation of seismic data

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (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)
  • Fluid Mechanics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供一种起伏地表变偏移距VSP积分叠前深度偏移方法，包括输入变偏移距VSP共检波点数据、二维地质模型和偏移参数；根据输入数据，建立炮点、检波点在二维地质模型中的坐标；采用两点射线追踪计算散射点到炮点的走时，保存各个散射点的走时表；采用两点射线追踪计算散射点到检波点的走时，保存各个检波点的走时表；实现起伏地表变偏移距VSP共检波点积分叠前深度偏移成像；积分叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP积分叠前深度偏移叠加成像。本发明避免了走时表重复计算，算效率高；直接从起伏地表偏移，避免了静校正对波场的影响；深度域成像，有利于井旁构造研究、储层描述，有利于指导油气开发。

Description

一种起伏地表变偏移距VSP积分叠前深度偏移方法
技术领域
本发明涉及地球物理勘探中地震数据成像方法，属于垂直地震处理技术范畴，是一种起伏地表变偏移距VSP积分叠前深度偏移方法和装置。
背景技术
近年来，变偏移距VSP(即Walkaway VSP)已实现工业化应用。主要是通过变偏移距VSP成像，研究井旁构造、储层描述。
刘少勇等研究了《起伏地表Kirchhoff积分法叠前深度偏移方法研究与应用》，文献研究了三维地面地震的最短路径走时计算方法，实现了三维地面地震起伏地表Kirchhoff积分法叠前深度偏移。
现有的变偏移距VSP积分成像，多为时间域、平地表的，起伏地表的变偏移距VSP波场需要静校正处理。本发明提出了一种起伏地表变偏移距VSP积分叠前深度偏移方法。无需静校正处理；采用两点射线追踪，计算散射点到起伏地表炮点的走时存为走时表，计算检波点到散射点走时存为走时表；将孔径内的散射点到炮点的走时加上检波点到散射点的走时得到炮检对的共散射走时，将共散射走时对应的波场映射为散射点深度对应的波场，孔径内积分求和得到积分叠前深度偏移成像。
发明内容
本发明目的在于实现起伏地表变偏移距VSP积分叠前深度偏移成像，为深度域井旁构造研究、储层描述提供数据。
具体技术方案为：
一种起伏地表变偏移距VSP积分叠前深度偏移方法，包括以下步骤：
(1)输入变偏移距VSP共检波点数据、二维地质模型和偏移参数；
(1.1)输入变偏移距VSP共检波点波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；
(1.2)输入二维网格层速度地质模型{ModVel}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；
(1.3)输入偏移孔径aper、成像起始坐标xmig1、成像结束坐标xmig2、成像网格dxmig、成像终止深度zmig2、成像深度步长dzmig参数。
(2)根据输入数据，建立炮点、检波点在二维地质模型中的坐标；
用步骤(1)中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在二维地质模型中的坐标：
炮点在二维地质模型中的坐标为：
Figure BDA0002525422470000021
MSZi＝SZi
其中，MSXi是第i个炮点在二维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于二维地质模型的深度、表示起伏地表，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；
检波点在二维地质模型中的坐标为：
Figure BDA0002525422470000022
MRZj＝RZj
其中，MRXj是第j个检波点在二维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。
(3)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到炮点的走时，保存各个散射点的走时表；
(3.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：
Figure BDA0002525422470000023
Figure BDA0002525422470000024
其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长。
(3.2)两点射线追踪第kxtr列第kztr行散射点到第i个炮点的走时：
Figure BDA0002525422470000025
Figure BDA0002525422470000026
其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MSXi、MSZi是第i个炮点在二维地质模型中的坐标，ModVelxtr,ztr是xtr、ztr处的速度，max(ModVelxtr,k·dzmig)是k∈[(ztr-MSZi)/dzmig,ktr]中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到。
Figure BDA0002525422470000031
其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标。
Figure BDA0002525422470000032
其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVelxtr,ztr是xtr、ztr处的速度，MSXi、MSZi是第i个炮点在二维地质模型中的坐标，
Figure BDA0002525422470000033
是第i个炮点的xtr、ztr处的走时。
(3.3)循环步骤(3.2)计算第kxtr列的散射点到第i个炮点的走时
Figure BDA0002525422470000034
(3.4)循环步骤(3.2)～(3.3)，计算第kxtr列的散射点到全部炮点的走时tskxtr，保存为第kxtr列散射点的走时表。
(3.5)循环步骤(3.1)～(3.4)，计算全部列散射点到全部炮点的走时表。
(4)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到检波点的走时，保存各个检波点的走时表；
(4.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：
Figure BDA0002525422470000035
Figure BDA0002525422470000036
其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长。
(4.2)两点射线追踪第kxtr列第kztr行散射点到第j个检波点的走时：
Figure BDA0002525422470000037
Figure BDA0002525422470000041
其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MRXj、MRZj是第j个炮点在二维地质模型中的坐标，ModVelxtr,ztr是xtr、ztr处的速度，max(ModVelxtr,k·dzmig)是
Figure BDA0002525422470000042
中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到。
Figure BDA0002525422470000043
其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标。
Figure BDA0002525422470000044
其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVelxtr,ztr是xtr、ztr处的速度，MRXj、MRZj是第j个检波点在二维地质模型中的坐标，
Figure BDA0002525422470000045
是j个检波点的xtr、ztr处的走时。
(4.3)循环步骤(4.2)计算第kxtr列的散射点到第j个检波点的走时
Figure BDA0002525422470000046
(4.4)循环步骤(4.1)～(4.3)，计算全部列的散射点到第j个检波点的走时trj，保存为第j个检波点到散射点的走时表。
(5)在步骤(2)的模型坐标系中，用步骤(3)、步骤(4)计算的共散射走时，实现起伏地表变偏移距VSP共检波点积分叠前深度偏移成像。
(5.1)用步骤(1)的偏移参数定义成像坐标，选择步骤(1)输入数据中第j个检波点的波场数据VSPDataj
xmigktr＝xmig1+(ktr-1)·dxmig ktr＝[1,(xmig2-xmig1)/dxmig]
其中，xmigktr是第ktr道的成像坐标，xmig1是成像起始坐标，xmig2是成像结束坐标。
MigDataktr,j＝{0}
其中，MigDataktr,j第j个检波点第ktr道的成像数据，先置为零。
(5.2)用步骤(1)的偏移孔径、步骤(2)的炮点模型坐标、步骤(5.1)的成像坐标确定成像输入地震道范围：
xmigktr-aper≤MSXkaper,j≤xmigktr+aper
其中，xmigktr是第ktr道的成像坐标，aper是偏移孔径，MSXkaper,j是第j个检波点第ktr个成像道孔径内第kaper个炮点的模型坐标。
(5.3)选择步骤(3)第ktr个散射点的第kaper个炮点的走时加上步骤(4)第个j检波点第ktr个散射点的走时，即为xmigktr处kaper、j炮检对的共散射走时。
TSR＝tskaper,ktr+trj,ktr
其中，TSR是xmigktr处kaper、j炮检对的共散射走时，tskaper,ktr是第ktr个散射点的第kaper个炮点的走时，trj,ktr第个j检波点第ktr个散射点的走时。
(5.4)选择步骤(5.2)中第j个检波点第kaper个炮点的波场数据VSPDatakaper,j，把TSR对应的振幅映射为散射点深度，与xmigktr处的成像道MigDataktr,j积分求和。
(5.5)循环步骤(5.2)～(5.4)，实现第j个检波点第ktr个成像道孔径内全部地震道积分叠前深度偏移成像，得到MigDataktr,j。
(5.6)循环步骤(5.1)～(5.5)，实现第j个检波点全部成像道的积分叠前深度偏移成像得到MigDataj。
(5.7)并行循环步骤(5.1)～(5.6)，实现全部检波点全部成像道的积分叠前深度偏移成像得到MigData。
(6)步骤(5)的积分叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP积分叠前深度偏移叠加成像。
本发明提出了一种起伏地表变偏移距VSP积分叠前深度偏移方法，避免了走时表重复计算，算效率高；直接从起伏地表偏移，避免了静校正对波场的影响；深度域成像，有利于井旁构造研究、储层描述，有利于指导油气开发。
附图说明
图1为输入的二维地质模型；横坐标为模型长度(单位：米)；纵坐标为模型深度(单位：米)。
图2为输入的变偏移距VSP共检波点波场数据；横坐标为道数(单位：无)；纵坐标为时间(单位：毫秒)。
图3为某列的散射点到全部炮点的走时；横坐标为炮数(单位：无)；纵坐标为深度(单位：米)。
图4为某个检波点到所有列散射点的走时；横坐标散射点列数(单位：无)；纵坐标为深度(单位：米)。
图5为图4中检波点到图3中列散射点再到全部炮点的共散射走时；横坐标为炮数(单位：无)；纵坐标为深度(单位：米)。
图6为图5中偏移孔径内的走时；横坐标为炮数(单位：无)；纵坐标为深度(单位：米)。
图7为图2中波场的共检波点叠前深度偏移成像；横坐标为道号(单位：无)；纵坐标为深度(单位：米)。
图8为多个共检波点叠前深度偏移成像的叠加成像；横坐标为道号(单位：无)；纵坐标为深度(单位：米)。
具体实施方式
结合实施例说明本发明的具体技术方案。
(1)输入变偏移距VSP共检波点数据、二维地质模型和偏移参数；
(1.1)输入变偏移距VSP共检波点波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；
(1.2)输入二维网格层速度地质模型{ModVel}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；
(1.3)输入偏移孔径aper、成像起始坐标xmig1、成像结束坐标xmig2、成像网格dxmig、成像终止深度zmig2、成像深度步长dzmig参数。
输入的二维地质模型如图1所示，输入的变偏移距VSP共检波点波场数据如图2所示。
(2)根据输入数据，建立炮点、检波点在二维地质模型中的坐标；
用步骤(1)中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在二维地质模型中的坐标：
炮点在二维地质模型中的坐标为：
Figure BDA0002525422470000061
MSZi＝SZi
其中，MSXi是第i个炮点在二维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于二维地质模型的深度、表示起伏地表，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；
检波点在二维地质模型中的坐标为：
Figure BDA0002525422470000062
MRZj＝RZj
其中，MRXj是第j个检波点在二维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。
(3)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到炮点的走时，保存各个散射点的走时表；
(3.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：
Figure BDA0002525422470000071
Figure BDA0002525422470000072
其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长。
(3.2)两点射线追踪第kxtr列第kztr行散射点到第i个炮点的走时：
Figure BDA0002525422470000073
Figure BDA0002525422470000074
其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MSXi、MSZi是第i个炮点在二维地质模型中的坐标，ModVelxtr,ztr是xtr、ztr处的速度，max(ModVelxtr,k·dzmig)是k∈[(ztr-MSZi)/dzmig,ktr]中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到。
Figure BDA0002525422470000075
其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标。
Figure BDA0002525422470000076
其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVelxtr,ztr是xtr、ztr处的速度，MSXi、MSZi是第i个炮点在二维地质模型中的坐标，
Figure BDA0002525422470000077
是第i个炮点的xtr、ztr处的走时。
(3.3)循环步骤(3.2)计算第kxtr列的散射点到第i个炮点的走时
Figure BDA0002525422470000078
(3.4)循环步骤(3.2)～(3.3)，计算第kxtr列的散射点到全部炮点的走时tskxtr，保存为第kxtr列散射点的走时表。
(3.5)循环步骤(3.1)-(3.4)，计算全部列散射点到全部炮点的走时表。
某列的散射点到全部炮点的走时，如图3所示。
(4)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到检波点的走时，保存各个检波点的走时表；
(4.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：
Figure BDA0002525422470000081
Figure BDA0002525422470000082
其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长。
(4.2)两点射线追踪第kxtr列第kztr行散射点到第j个检波点的走时：
Figure BDA0002525422470000083
Figure BDA0002525422470000084
其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MRXj、MRZj是第j个炮点在二维地质模型中的坐标，ModVelxtr,ztr是xtr、ztr处的速度，max(ModVelxtr,k·dzmig)是
Figure BDA0002525422470000085
中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到。
Figure BDA0002525422470000086
其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标。
Figure BDA0002525422470000087
其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVelxtr,ztr是xtr、ztr处的速度，MRXj、MRZj是第j个检波点在二维地质模型中的坐标，
Figure BDA0002525422470000091
是j个检波点的xtr、ztr处的走时。
(4.3)循环(4.2)计算第kxtr列的散射点到第j个检波点的走时
Figure BDA0002525422470000092
(4.4)循环(4.1)～(4.3)，计算全部列的散射点到第j个检波点的走时trj，保存为第j个检波点到散射点的走时表。
某个检波点到所有列散射点的走时，如图4所示。
(5)在步骤(2)的模型坐标系中，用步骤(3)、步骤(4)计算的共散射走时，实现起伏地表变偏移距VSP共检波点积分叠前深度偏移成像。
(5.1)用步骤(1)的偏移参数定义成像坐标，选择步骤(1)输入数据中第j个检波点的波场数据VSPDataj
xmigktr＝xmig1+(ktr-1)·dxmig ktr＝[1,(xmig2-xmig1)/dxmig]
其中，xmigktr是第ktr道的成像坐标，xmig1是成像起始坐标，xmig2是成像结束坐标。
MigDataktr,j＝{0}
其中，MigDataktr,j第j个检波点第ktr道的成像数据，先置为零。
(5.2)用步骤(1)的偏移孔径、步骤(2)的炮点模型坐标、步骤(5.1)的成像坐标确定成像输入地震道范围
xmigktr-aper≤MSXkaper,j≤xmigktr+aper
其中，xmigktr是第ktr道的成像坐标，aper是偏移孔径，MSXkaper,j是第j个检波点第ktr个成像道孔径内第kaper个炮点的模型坐标。
(5.3)选择步骤(3)第ktr个散射点的第kaper个炮点的走时加上步骤(4)第个j检波点第ktr个散射点的走时，即为xmigktr处kaper、j炮检对的共散射走时。
TSR＝tskaper,ktr+trj,ktr
其中，TSR是xmigktr处kaper、j炮检对的共散射走时，tskaper,ktr是第第ktr个散射点的第kaper个炮点的走时，trj,ktr第个j检波点第ktr个散射点的走时。
(5.4)选择步骤(5.2)中第j个检波点第kaper个炮点的波场数据VSPDatakaper,j，把TSR对应的振幅映射为散射点深度，与xmigktr处的成像道MigDataktr,j积分求和。
(5.5)循环步骤(5.2)～(5.4)，实现第j个检波点第ktr个成像道孔径内全部地震道积分叠前深度偏移成像，得到MigDataktr,j。
(5.6)循环步骤(5.1)～(5.5)，实现第j个检波点全部成像道的积分叠前深度偏移成像得到MigDataj。
(5.7)并行循环步骤(5.1)～(5.6)，实现全部检波点全部成像道的积分叠前深度偏移成像得到MigData。
图4中检波点到图3中列散射点再到全部炮点的共散射走时，如图5所示；图5中偏移孔径内的共散射走时，如图6所示。图2中波场的叠前深度偏移成像，如图7所示。
(6)步骤(5)的积分叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP积分叠前深度偏移叠加成像。
为多个共检波点叠前深度偏移成像的叠加成像，如图8所示。

Claims (6)

1.一种起伏地表变偏移距VSP积分叠前深度偏移方法，其特征在于，包括以下步骤：
(1)输入变偏移距VSP共检波点数据、二维地质模型和偏移参数；
(2)根据输入数据，建立炮点、检波点在二维地质模型中的坐标；
(3)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到炮点的走时，保存各个散射点的走时表；
(4)在步骤(2)的模型坐标系中，采用两点射线追踪计算散射点到检波点的走时，保存各个检波点的走时表；
(5)在步骤(2)的模型坐标系中，用步骤(3)、步骤(4)计算的共散射走时，实现起伏地表变偏移距VSP共检波点积分叠前深度偏移成像；
(6)步骤(5)的积分叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP积分叠前深度偏移叠加成像。
2.根据权利要求1所述的一种起伏地表变偏移距VSP积分叠前深度偏移方法，其特征在于，步骤(1)具体包括以下子步骤：
(1.1)输入变偏移距VSP共检波点波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；
(1.2)输入二维网格层速度地质模型{ModVel}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；
(1.3)输入偏移孔径aper、成像起始坐标xmig1、成像结束坐标xmig2、成像网格dxmig、成像终止深度zmig2、成像深度步长dzmig参数。
3.根据权利要求1所述的一种起伏地表变偏移距VSP积分叠前深度偏移方法，其特征在于，步骤(2)具体包括以下过程：
用步骤(1)中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在二维地质模型中的坐标：
炮点在二维地质模型中的坐标为：
Figure FDA0002525422460000011
MSZi＝SZi
其中，MSXi是第i个炮点在二维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于二维地质模型的深度、表示起伏地表，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；
检波点在二维地质模型中的坐标为：
Figure FDA0002525422460000021
MRZj＝RZj
其中，MRXj是第j个检波点在二维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。
4.根据权利要求1所述的一种起伏地表变偏移距VSP积分叠前深度偏移方法，其特征在于，步骤(3)具体包括以下子步骤：
(3.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：
Figure FDA0002525422460000022
Figure FDA0002525422460000023
其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长；
(3.2)两点射线追踪第kxtr列第kztr行散射点到第i个炮点的走时：
Figure FDA0002525422460000024
Figure FDA0002525422460000025
其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MSXi、MSZi是第i个炮点在二维地质模型中的坐标，ModVelxtr,ztr是xtr、ztr处的速度，max(ModVelxtr,k·dzmig)是k∈[(ztr-MSZi)/dzmig,ktr]中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到；
Figure FDA0002525422460000026
其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标；
Figure FDA0002525422460000027
其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVelxtr,ztr是xtr、ztr处的速度，MSXi、MSZi是第i个炮点在二维地质模型中的坐标，
Figure FDA0002525422460000031
是第i个炮点的xtr、ztr处的走时；
(3.3)循环步骤(3.2)计算第kxtr列的散射点到第i个炮点的走时
Figure FDA0002525422460000032
(3.4)循环步骤(3.2)～(3.3)，计算第kxtr列的散射点到全部炮点的走时tskxtr，保存为第kxtr列散射点的走时表；
(3.5)循环步骤(3.1)～(3.4)，计算全部列散射点到全部炮点的走时表。
5.根据权利要求1所述的一种起伏地表变偏移距VSP积分叠前深度偏移方法，其特征在于，步骤(4)具体包括以下子步骤：
(4.1)用步骤(1)输入的成像起始坐标、成像结束坐标、成像网格、成像终止深度、成像深度步长，计算散射点坐标：
Figure FDA0002525422460000033
Figure FDA0002525422460000034
其中，xmig1是成像起始坐标、xmig2是成像结束坐标、dxmig是成像网格，kxtr是散射点序号，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，zmig2是成像终止深度，dzmig成像深度步长；
(4.2)两点射线追踪第kxtr列第kztr行散射点到第j个检波点的走时：
Figure FDA0002525422460000035
Figure FDA0002525422460000036
其中，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，MRXj、MRZj是第j个炮点在二维地质模型中的坐标，ModVelxtr,ztr是xtr、ztr处的速度，max(ModVelxtr,k·dzmig)是
Figure FDA0002525422460000038
中的速度的最大值，dzmig成像深度步长；L是未知数，通过求解非线性方程得到；
Figure FDA0002525422460000037
其中，p是射线参数，L是上式计算得到，dzmig成像深度步长，xtr是第kxtr列散射点坐标；
Figure FDA0002525422460000041
其中，p是上式计算的射线参数，dzmig成像深度步长，xtr是第kxtr列散射点坐标，ztr是第kztr行散射点深度，ModVelxtr,ztr是xtr、ztr处的速度，MRXj、MRZj是第j个检波点在二维地质模型中的坐标，
Figure FDA0002525422460000042
是j个检波点的xtr、ztr处的走时；
(4.3)循环步骤(4.2)计算第kxtr列的散射点到第j个检波点的走时
Figure FDA0002525422460000043
(4.4)循环步骤(4.1)～(4.3)，计算全部列的散射点到第j个检波点的走时trj，保存为第j个检波点到散射点的走时表。
6.根据权利要求1所述的一种起伏地表变偏移距VSP积分叠前深度偏移方法，其特征在于，步骤(5)具体包括以下子步骤：
(5.1)用步骤(1)的偏移参数定义成像坐标，选择步骤(1)输入数据中第j个检波点的波场数据VSPDataj：
xmigktr＝xmig1+(ktr-1)·dxmig ktr＝[1,(xmig2-xmig1)/dxmig]
其中，xmigktr是第ktr道的成像坐标，xmig1是成像起始坐标，xmig2是成像结束坐标；
MigDataktr,j＝{0}
其中，MigDataktr,j第j个检波点第ktr道的成像数据，先置为零；
(5.2)用步骤(1)的偏移孔径、步骤(2)的炮点模型坐标、步骤(5.1)的成像坐标确定成像输入地震道范围：
xmigktr-aper≤MSXkaper,j≤xmigktr+aper
其中，xmigktr是第ktr道的成像坐标，aper是偏移孔径，MSXkaper,j是第j个检波点第ktr个成像道孔径内第kaper个炮点的模型坐标；
(5.3)选择步骤(3)第ktr个散射点的第kaper个炮点的走时加上步骤(4)第个j检波点第ktr个散射点的走时，即为xmigktr处kaper、j炮检对的共散射走时；
TSR＝tskaper,ktr+trj,ktr
其中，TSR是xmigktr处kaper、j炮检对的共散射走时，tskaper,ktr是第ktr个散射点的第kaper个炮点的走时，trj,ktr第个j检波点第ktr个散射点的走时；
(5.4)选择步骤(5.2)中第j个检波点第kaper个炮点的波场数据VSPDatakaper,j，把TSR对应的振幅映射为散射点深度，与xmigktr处的成像道MigDataktr,j积分求和；
(5.5)循环步骤(5.2)～(5.4)，实现第j个检波点第ktr个成像道孔径内全部地震道积分叠前深度偏移成像，得到MigDataktr,j；
(5.6)循环步骤(5.1)～(5.5)，实现第j个检波点全部成像道的积分叠前深度偏移成像得到MigDataj；
(5.7)并行循环步骤(5.1)～(5.6)，实现全部检波点全部成像道的积分叠前深度偏移成像得到MigData。
CN202010502932.1A 2020-06-05 2020-06-05 一种起伏地表变偏移距vsp积分叠前深度偏移方法 Active CN111624648B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010502932.1A CN111624648B (zh) 2020-06-05 2020-06-05 一种起伏地表变偏移距vsp积分叠前深度偏移方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010502932.1A CN111624648B (zh) 2020-06-05 2020-06-05 一种起伏地表变偏移距vsp积分叠前深度偏移方法

Publications (2)

Publication Number Publication Date
CN111624648A true CN111624648A (zh) 2020-09-04
CN111624648B CN111624648B (zh) 2022-04-01

Family

ID=72260234

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010502932.1A Active CN111624648B (zh) 2020-06-05 2020-06-05 一种起伏地表变偏移距vsp积分叠前深度偏移方法

Country Status (1)

Country Link
CN (1) CN111624648B (zh)

Citations (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040041815A1 (en) * 2002-08-27 2004-03-04 Conoco Inc. Method of building and updating an anisotropic velocity model for depth imaging of seismic data
WO2005010797A2 (en) * 2003-07-23 2005-02-03 Lee Wook B Improved 3d veloctiy modeling, with calibration and trend fitting using geostatistical techniques, particularly advantageous for curved-ray prestack time migration and for such migration followed by prestack depth migration
CN102193109A (zh) * 2011-03-10 2011-09-21 中国科学院地质与地球物理研究所 起伏地表采集的三维地震资料的直接叠前时间偏移方法
EP2624015A2 (en) * 2012-02-02 2013-08-07 CGGVeritas Services SA Sweep design for seismic sources
CN104216009A (zh) * 2013-06-05 2014-12-17 中国石油天然气集团公司 一种斜井三维垂直地震剖面时间偏移的方法
CN104216007A (zh) * 2013-06-05 2014-12-17 中国石油天然气集团公司 二维垂直地震剖面和三维地面数据的空间同步成像处理方法
CA2926179A1 (en) * 2013-10-02 2015-04-09 Bp Corporation North America Inc. System and method for seismic adaptive optics
CN105629299A (zh) * 2015-12-19 2016-06-01 中国石油集团川庆钻探工程有限公司地球物理勘探公司 角度域叠前深度偏移的走时、角度表获取方法及成像方法
US20170038490A1 (en) * 2014-04-17 2017-02-09 Saudi Arabian Oil Company Generating subterranean imaging data based on vertical seismic profile data and ocean bottom sensor data
CN107229069A (zh) * 2016-03-24 2017-10-03 中国石油化工股份有限公司 一种对共检波点数据道集进行速度分析的方法
CN107561583A (zh) * 2017-08-01 2018-01-09 中国石油集团川庆钻探工程有限公司地球物理勘探公司 用于高斯束叠前深度偏移的局部角度计算方法及成像方法
WO2018075738A1 (en) * 2016-10-19 2018-04-26 Saudi Arabian Oil Company Generating subterranean imaging data based on vertical seismic profile data and ocean bottom sensor data
CN108363101A (zh) * 2018-02-02 2018-08-03 西安石油大学 一种斜井井间地震高斯束叠前深度偏移成像方法
CN108693560A (zh) * 2017-04-12 2018-10-23 中国石油化工股份有限公司 一种基于互相关道的散射波成像方法及系统
CN108919354A (zh) * 2018-09-27 2018-11-30 中国科学院地质与地球物理研究所 近地表q偏移方法及装置
EP3488075A1 (en) * 2016-09-27 2019-05-29 Halliburton Energy Services, Inc. Iterative migration velocity optimization for a vsp survey using semblance
CN110907995A (zh) * 2018-09-14 2020-03-24 中国石油天然气股份有限公司 井中vsp地震数据的逆时偏移方法及装置

Patent Citations (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040041815A1 (en) * 2002-08-27 2004-03-04 Conoco Inc. Method of building and updating an anisotropic velocity model for depth imaging of seismic data
WO2005010797A2 (en) * 2003-07-23 2005-02-03 Lee Wook B Improved 3d veloctiy modeling, with calibration and trend fitting using geostatistical techniques, particularly advantageous for curved-ray prestack time migration and for such migration followed by prestack depth migration
CN102193109A (zh) * 2011-03-10 2011-09-21 中国科学院地质与地球物理研究所 起伏地表采集的三维地震资料的直接叠前时间偏移方法
EP2624015A2 (en) * 2012-02-02 2013-08-07 CGGVeritas Services SA Sweep design for seismic sources
CN104216009A (zh) * 2013-06-05 2014-12-17 中国石油天然气集团公司 一种斜井三维垂直地震剖面时间偏移的方法
CN104216007A (zh) * 2013-06-05 2014-12-17 中国石油天然气集团公司 二维垂直地震剖面和三维地面数据的空间同步成像处理方法
US20160377755A1 (en) * 2013-10-02 2016-12-29 Bp Corporation North America Inc. System and method for seismic adaptive optics
CA2926179A1 (en) * 2013-10-02 2015-04-09 Bp Corporation North America Inc. System and method for seismic adaptive optics
US20170038490A1 (en) * 2014-04-17 2017-02-09 Saudi Arabian Oil Company Generating subterranean imaging data based on vertical seismic profile data and ocean bottom sensor data
CN105629299A (zh) * 2015-12-19 2016-06-01 中国石油集团川庆钻探工程有限公司地球物理勘探公司 角度域叠前深度偏移的走时、角度表获取方法及成像方法
CN107229069A (zh) * 2016-03-24 2017-10-03 中国石油化工股份有限公司 一种对共检波点数据道集进行速度分析的方法
EP3488075A1 (en) * 2016-09-27 2019-05-29 Halliburton Energy Services, Inc. Iterative migration velocity optimization for a vsp survey using semblance
WO2018075738A1 (en) * 2016-10-19 2018-04-26 Saudi Arabian Oil Company Generating subterranean imaging data based on vertical seismic profile data and ocean bottom sensor data
CN108693560A (zh) * 2017-04-12 2018-10-23 中国石油化工股份有限公司 一种基于互相关道的散射波成像方法及系统
CN107561583A (zh) * 2017-08-01 2018-01-09 中国石油集团川庆钻探工程有限公司地球物理勘探公司 用于高斯束叠前深度偏移的局部角度计算方法及成像方法
CN108363101A (zh) * 2018-02-02 2018-08-03 西安石油大学 一种斜井井间地震高斯束叠前深度偏移成像方法
CN110907995A (zh) * 2018-09-14 2020-03-24 中国石油天然气股份有限公司 井中vsp地震数据的逆时偏移方法及装置
CN108919354A (zh) * 2018-09-27 2018-11-30 中国科学院地质与地球物理研究所 近地表q偏移方法及装置

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李建国: "Walkaway VSP自由表面多次波叠前深度偏移成像", 《石油地球物理勘探》 *

Also Published As

Publication number Publication date
CN111624648B (zh) 2022-04-01

Similar Documents

Publication Publication Date Title
CN110133715B (zh) 一种基于初至时差和波形叠加的微地震震源定位方法
CN104570125B (zh) 一种利用井数据提高成像速度模型精度的方法
CN102282481B (zh) 基于地震能见度分析的数据采集和叠前偏移
CN102841379B (zh) 一种基于共散射点道集的叠前时间偏移与速度分析方法
CN102841375A (zh) 一种复杂条件下基于角度域共成像点道集的层析速度反演方法
CN101957455A (zh) 三维保幅叠前时间偏移方法
US8078406B2 (en) Processing seismic data in common group-center gathers
Wuestefeld et al. Benchmarking earthquake location algorithms: A synthetic comparison
Tan et al. Microseismic velocity model inversion and source location: The use of neighborhood algorithm and master station method
Zhang et al. Deep learning for efficient microseismic location using source migration‐based imaging
Dando et al. Relocating microseismicity from downhole monitoring of the Decatur CCS site using a modified double-difference algorithm
Talukdar et al. Sub-basalt imaging of hydrocarbon-bearing Mesozoic sediments using ray-trace inversion of first-arrival seismic data and elastic finite-difference full-wave modeling along Sinor–Valod profile of Deccan Syneclise, India
CN106054252B (zh) 一种叠前时间偏移的方法及装置
CN102053269A (zh) 一种对地震资料中速度分析方法
CN105137479A (zh) 一种面元覆盖次数的计算方法及装置
CN104597496B (zh) 一种二维地震资料中速度的三维空间归位方法
CN111624648B (zh) 一种起伏地表变偏移距vsp积分叠前深度偏移方法
CN104849751A (zh) 叠前地震资料成像的方法
CN102466818B (zh) 一种利用井间地震数据对各向异性介质成像的方法
CN104614762A (zh) 疏松砂岩气藏边界确定方法及装置
CN102778691B (zh) 一种计算检波器组内静校正时差的方法
Sylvester Automated multi‐well stratigraphic correlation and model building using relative geologic time
CN105259577B (zh) 一种确定地层界面的角度信息的方法及装置
Wang et al. Cross-related microseismic location based on improved particle swarm optimization and the double-difference location method of jointed coal rock mass
Saffarzadeh et al. Improving fault image by determination of optimum seismic survey parameters using ray-based modeling

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