CN111751876A - 一种变偏移距vsp转换横波单程波叠前深度偏移方法和装置 - Google Patents

一种变偏移距vsp转换横波单程波叠前深度偏移方法和装置 Download PDF

Info

Publication number
CN111751876A
CN111751876A CN202010643476.2A CN202010643476A CN111751876A CN 111751876 A CN111751876 A CN 111751876A CN 202010643476 A CN202010643476 A CN 202010643476A CN 111751876 A CN111751876 A CN 111751876A
Authority
CN
China
Prior art keywords
wave
depth
dzmig
imaging
frequency
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
CN202010643476.2A
Other languages
English (en)
Other versions
CN111751876B (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 CN202010643476.2A priority Critical patent/CN111751876B/zh
Publication of CN111751876A publication Critical patent/CN111751876A/zh
Application granted granted Critical
Publication of CN111751876B publication Critical patent/CN111751876B/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/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
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • 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

  • 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

本发明公开了一种变偏移距VSP转换横波单程波叠前深度偏移方法和装置，所述方法包括以下步骤：S1.输入变偏移距VSP共检波点反射转换横波波场、纵波2维地质模型、横波2维地质模型、偏移参数；S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；S3.对于反射转换横波时间空间域波场，进行快速傅里叶变换到频率波数域，并设置频率域震源；S4.实现单检波点的单程波叠前深度偏移成像；S5.完成所有检波点的变偏移距VSP转换横波单程波叠前深度偏移；S6.共成像道集叠加，得到变偏移VSP转换横波单程波叠前深度偏移叠加成像。本发明实现了变偏移距VSP共检波点单程波叠前深度偏移成像，并得到变偏移距VSP叠加成像，为井旁构造研究、储层预测提供了数据支撑。

Description

一种变偏移距VSP转换横波单程波叠前深度偏移方法和装置
技术领域
本发明涉及地震数据偏移成像，特别是涉及一种变偏移距VSP转换横波单程波叠前深度偏移方法和装置。
背景技术
转换横波勘探通常由纵波震源激发，多分量检波器接收，同时获得纵波资料与转换横波资料。转换横波资料主要用于构造成像、岩性辨别、弹性参数反演及裂缝发育预测等。
变偏移距VSP是指多个检波器组成的阵列陈放在观测井中，地表布设一条人工激发的多个炮点组成的炮线，是一种VSP的二维观测方式，在炮点人工激发地震数据的条件下，实现变偏移距VSP井中数据采集，采集系统立体图如图1所示；变偏移距VSP采集，其实就是在每一个炮点激发时，由各个检波点的检波器进行数据采集；变偏移距VSP可以观测到反射纵波、反射转换横波，利用偏移成像研究井旁构造、储层预测。
宋勇研究了《Kirchhoff积分转换波叠前深度偏移方法研究》，文献采用最短路径法射线追踪计算走时，实现了地面地震转换波叠前深度偏移成像。岳玉波等研究了《转换波Kirchhoff叠前时间偏移的成像优化方案》，文献研究了偏移振幅加权函数、反假频算子，实现了VTI介质地面地震转换波Kirchhoff叠前时间偏移。但是变偏移距VSP偏移成像研究主要集中于时间域纵波、深度域纵波，转换横波很少。
发明内容
本发明的目的在于克服现有技术的不足，提供一种变偏移距VSP转换横波单程波叠前深度偏移方法和装置，实现了变偏移距VSP共检波点单程波叠前深度偏移成像，并将多个检波点分别成像后叠加，得到变偏移距VSP叠加成像，为井旁构造研究、储层预测提供了数据支撑。
本发明的目的是通过以下技术方案来实现的：一种变偏移距VSP转换横波单程波叠前深度偏移方法，包括以下步骤：
S1.输入变偏移距VSP共检波点反射转换横波波场、纵波2维地质模型、横波2维地质模型、偏移参数；
S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；
S3.对于第j个检波点的反射转换横波时间空间域波场，进行快速傅里叶变换到频率波数域，并设置频率域震源；
S4.对于第j个检波点，频率域震源从检波点深度用横波速度向下延拓，频率波数域变偏移距VSP共检波点反射转换横波波场用纵波速度从地表向下延拓，延拓算子为ω-x单程波算子，并实现单检波点的单程波叠前深度偏移成像；
S5.对每一个检波点，重复执行步骤S3～S4，完成所有检波点的变偏移距VSP转换横波单程波叠前深度偏移；
S6.将步骤S5中得到的叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP转换横波单程波叠前深度偏移叠加成像。
进一步地，所述步骤S1包括：
S101.输入变偏移距VSP共检波点反射转换横波波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；
S202.输入纵波2维网格层速度地质模型{Vp}，输入横波2维网格层速度地质模型{Vs}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；
S203.输入成像算子的最大倾角dip、成像终止深度zmig2、成像深度步长dzmig、变偏移距VSP转换横波波场的最小频率fl、变偏移距VSP转换横波波场的最大频率参数fh。
进一步地，所述步骤S2包括：
利用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：
炮点在2维地质模型中的坐标为：
Figure BDA0002572148540000021
MSZi＝SZi
其中，MSXi是第i个炮点在2维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于2维地质模型的深度、表示起伏地表，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；
检波点在2维地质模型中的坐标为：
Figure BDA0002572148540000022
MRZj＝RZj
其中，MRXj是第j个检波点在2维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。
进一步地，所述步骤S3包括：
S301.将步骤S1中第j个检波点的反射转换横波时间空间域波场VSPDataj做快速傅里叶变换转换到频率波数域FVSPj；
S302.在步骤S2的坐标系下，将频率域震源FSouj设置在第j个检波点处。
进一步地，所述步骤S4中，对于第j个检波点，完成变偏移距VSP转换横波单程波叠前深度偏移的过程包括：
S401.计算深度成像网格：
kzmig＝(k-1)·dzmig
Figure BDA0002572148540000031
其中，zmig2是成像终止深度，dzmig是成像深度步长，kzmig是第k个成像深度；
S402.频率波数域变偏移距VSP共检波点反射转换横波波场
Figure BDA0002572148540000032
用纵波速度向下延拓一个深度步长：
Figure BDA0002572148540000033
f＝[fl:df:fh]
Figure BDA0002572148540000034
kk＝f/Vp
g0＝1
g1＝-i·π·dzmig/kk
g2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4
g3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6
g4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8
Figure BDA0002572148540000035
Figure BDA0002572148540000036
Figure BDA0002572148540000037
Figure BDA0002572148540000041
Figure BDA0002572148540000042
Figure BDA0002572148540000043
其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，表示频率从fl到fh间隔为df，Vp纵波速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mute是波数切除因子，kk、g0-g4、t0-t4是中间变量，
Figure BDA0002572148540000044
是
Figure BDA0002572148540000045
向下延拓一个深度步长的波场；
S403.若kzmig＜＝MRZj，循环执行步骤S401～S402；
若kzmig＞MRZj，则跳转至步骤S404；
S404.频率域震源
Figure BDA0002572148540000046
用横波速度向下延拓一个深度步长：
Figure BDA0002572148540000047
f＝[fl:df:fh]
Figure BDA0002572148540000048
kks＝f/Vs
gs0＝1
gs1＝i·π·dzmig/kks
gs2＝(i·π·dzmig·kks/4-π2·dzmig2·kks2/2)/kks4
gs3＝(i·π·dzmig·kks/8+π2·dzmig2·kks2/4+i·π3·dzmig3·kks3/6)/kks6
gs4＝(i·π·dzmig·kks·5/64-π2·dzmig2·kks2·5/32-i·π3·dzmig3·kks3/8+π4·dzmig4·kks4/24)/kks8
Figure BDA0002572148540000049
Figure BDA00025721485400000410
Figure BDA00025721485400000411
Figure BDA00025721485400000412
Figure BDA0002572148540000051
Figure BDA0002572148540000052
其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，Vs横波速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，
Figure BDA0002572148540000053
是
Figure BDA0002572148540000054
向下延拓一个深度步长的波场；
S405.相关成像条件提取成像值：
Figure BDA0002572148540000055
其中，
Figure BDA0002572148540000056
是
Figure BDA0002572148540000057
向下延拓一个深度步长的波场，
Figure BDA0002572148540000058
是
Figure BDA0002572148540000059
向下延拓一个深度步长的波场，
Figure BDA00025721485400000510
是提取的kzmig+dzmig深度的成像值，conj是复共轭函数；
S406.循环执行步骤S401～S405，直至zmig2是成像终止深度，完成第j个检波点的变偏移距VSP转换横波单程波叠前深度偏移。
一种变偏移距VSP转换横波单程波叠前深度偏移装置，包括：
数据输入模块，用于输入变偏移距VSP共检波点反射转换横波波场、纵波2维地质模型、横波2维地质模型和偏移参数；
坐标建立模块，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；
单程波叠前深度偏移成像模块，用于对每一个检波点的反射转换横波时间空间域波场，进行快速傅里叶变换到频率波数域，并设置频率域震源，将频率域震源从检波点深度用横波速度向下延拓，频率波数域变偏移距VSP共检波点反射转换横波波场用纵波速度从地表向下延拓，并实现每一个检波点的单程波叠前深度偏移；
叠加成像模块，用于将得到的叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP转换横波单程波叠前深度偏移叠加成像。
本发明的有益效果是：本发明输入变偏移距VSP共检波点反射转换横波波场，将震源置于井中检波点处用横波速度向下延拓，将变偏移距VSP共检波点反射转换横波波场用纵波速度从地表向下延拓，延拓算子为频率波数域单程波算子，相关成像条件提取成像值，实现变偏移距VSP共检波点单程波叠前深度偏移成像；多个检波点分别成像后叠加，得到变偏移距VSP叠加成像，为井旁构造研究、储层预测提供数据支撑。
附图说明
图1为变偏移距VSP采集系统的立体图。
图2为本发明的方法流程图；
图3为实施例中输入的纵波2维地质模型示意图；
图4为实施例中输入的横波2维地质模型示意图；
图5为实施例中输入的共检波点变偏移距VSP转换横波波场示意图；
图6为实施例中单程波叠前深度偏移成像结果示意图；
图7为实施例中多个共检波点转换横波单程波叠前深度偏移成像的叠加成像结果示意图；
图8为本发明装置的原理框图。
具体实施方式
下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。
如图2所示，一种变偏移距VSP转换横波单程波叠前深度偏移方法，包括以下步骤：
S1.输入变偏移距VSP共检波点反射转换横波波场、纵波2维地质模型、横波2维地质模型、偏移参数：
S101.输入变偏移距VSP共检波点反射转换横波波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；
S202.输入纵波2维网格层速度地质模型{Vp}，输入横波2维网格层速度地质模型{Vs}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；
S203.输入最大倾角dip、成像终止深度zmig2、成像深度步长dzmig、最小频率fl、最大频率参数fh。
在本申请的实施例中，输入的纵波2维地质模型如图3所示，横坐标为模型长度(单位：米)；纵坐标为深度(单位：米)；输入的横波2维地质模型，如图4所示，横坐标为模型长度(单位：米)；纵坐标为深度(单位：米)；输入的共检波点变偏移距VSP转换横波波场如图5所示，横坐标为道号；纵坐标为时间(单位：毫秒)。
S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标：
利用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：
炮点在2维地质模型中的坐标为：
Figure BDA0002572148540000061
MSZi＝SZi
其中，MSXi是第i个炮点在2维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于2维地质模型的深度、表示起伏地表，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；
检波点在2维地质模型中的坐标为：
Figure BDA0002572148540000071
MRZj＝RZj
其中，MRXj是第j个检波点在2维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。
S3.对于第j个检波点的反射转换横波时间空间域波场，进行快速傅里叶变换到频率波数域，并设置频率域震源；
S301.将步骤S1中第j个检波点的反射转换横波时间空间域波场VSPDataj做快速傅里叶变换转换到频率波数域FVSPj；
S302.在步骤S2的坐标系下，将频率域震源FSouj设置在第j个检波点处。
S4.对于第j个检波点，频率域震源从检波点深度用横波速度向下延拓，频率波数域变偏移距VSP共检波点反射转换横波波场用纵波速度从地表向下延拓，延拓算子为ω-x单程波算子，并实现单检波点的单程波叠前深度偏移成像；
所述步骤S4中，对于第j个检波点，完成变偏移距VSP转换横波单程波叠前深度偏移的过程包括：
S401.计算深度成像网格：
kzmig＝(k-1)·dzmig
Figure BDA0002572148540000072
其中，zmig2是成像终止深度，dzmig是成像深度步长，kzmig是第k个成像深度；
S402.频率波数域变偏移距VSP共检波点反射转换横波波场
Figure BDA0002572148540000073
用纵波速度向下延拓一个深度步长：
Figure BDA0002572148540000074
f＝[fl:df:fh]
Figure BDA0002572148540000081
kk＝f/Vp
g0＝1
g1＝-i·π·dzmig/kk
g2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4
g3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6
g4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8
Figure BDA0002572148540000082
Figure BDA0002572148540000083
Figure BDA0002572148540000084
Figure BDA0002572148540000085
Figure BDA0002572148540000086
Figure BDA0002572148540000087
其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，表示频率从fl到fh间隔为df，，Vp纵波速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mute是波数切除因子，kk、g0-g4、t0-t4是中间变量，
Figure BDA0002572148540000088
是
Figure BDA0002572148540000089
向下延拓一个深度步长的波场；
S403.若kzmig＜＝MRZj，循环执行步骤S401～S402；
若kzmig＞MRZj，则跳转至步骤S404；
S404.频率域震源
Figure BDA00025721485400000810
用横波速度向下延拓一个深度步长：
Figure BDA00025721485400000811
f＝[fl:df:fh]
Figure BDA0002572148540000091
kks＝f/Vs
gs0＝1
gs1＝i·π·dzmig/kks
gs2＝(i·π·dzmig·kks/4-π2·dzmig2·kks2/2)/kks4
gs3＝(i·π·dzmig·kks/8+π2·dzmig2·kks2/4+i·π3·dzmig3·kks3/6)/kks6
gs4＝(i·π·dzmig·kks·5/64-π2·dzmig2·kks2·5/32-i·π3·dzmig3·kks3/8+π4·dzmig4·kks4/24)/kks8
Figure BDA0002572148540000092
Figure BDA0002572148540000093
Figure BDA0002572148540000094
Figure BDA0002572148540000095
Figure BDA0002572148540000096
Figure BDA0002572148540000097
其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，Vs横波速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，
Figure BDA0002572148540000098
是
Figure BDA0002572148540000099
向下延拓一个深度步长的波场；
S405.相关成像条件提取成像值：
Figure BDA00025721485400000910
其中，
Figure BDA00025721485400000911
是
Figure BDA00025721485400000912
向下延拓一个深度步长的波场，
Figure BDA00025721485400000913
是
Figure BDA00025721485400000914
向下延拓一个深度步长的波场，
Figure BDA00025721485400000915
是提取的kzmig+dzmig深度的成像值，conj是复共轭函数；
S406.循环执行步骤S401～S405，直至zmig2是成像终止深度，完成第j个检波点的变偏移距VSP转换横波单程波叠前深度偏移。
S5.对每一个检波点，重复执行步骤S3～S4，完成所有检波点的变偏移距VSP转换横波单程波叠前深度偏移；
在本申请的实施例中，对于图5对应的共检波点变偏移距VSP转换横波波场，其单程波叠前深度偏移成像结果如图6所示，横坐标为道号；纵坐标为深度(单位：米)；
S6.将步骤S5中得到的叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP转换横波单程波叠前深度偏移叠加成像。在本申请的实施例中，多个共检波点转换横波单程波叠前深度偏移成像的叠加成像结果如图7所示，横坐标为道号；纵坐标为深度(单位：米)。
如图8所示，一种变偏移距VSP转换横波单程波叠前深度偏移装置，包括：
数据输入模块，用于输入变偏移距VSP共检波点反射转换横波波场、纵波2维地质模型、横波2维地质模型和偏移参数；
坐标建立模块，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；
单程波叠前深度偏移成像模块，用于对每一个检波点的反射转换横波时间空间域波场，进行快速傅里叶变换到频率波数域，并设置频率域震源，将频率域震源从检波点深度用横波速度向下延拓，频率波数域变偏移距VSP共检波点反射转换横波波场用纵波速度从地表向下延拓，并实现每一个检波点的单程波叠前深度偏移；
叠加成像模块，用于将得到的叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP转换横波单程波叠前深度偏移叠加成像。
以上所述是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应该看作是对其他实施例的排除，而可用于其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (6)

1.一种变偏移距VSP转换横波单程波叠前深度偏移方法，其特征在于：包括以下步骤：
S1.输入变偏移距VSP共检波点反射转换横波波场、纵波2维地质模型、横波2维地质模型、偏移参数；
S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；
S3.对于第j个检波点的反射转换横波时间空间域波场，进行快速傅里叶变换到频率波数域，并设置频率域震源；
S4.对于第j个检波点，频率域震源从检波点深度用横波速度向下延拓，频率波数域变偏移距VSP共检波点反射转换横波波场用纵波速度从地表向下延拓，延拓算子为ω-x单程波算子，并实现单检波点的单程波叠前深度偏移成像；
S5.对每一个检波点，重复执行步骤S3～S4，完成所有检波点的变偏移距VSP转换横波单程波叠前深度偏移；
S6.将步骤S5中得到的叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP转换横波单程波叠前深度偏移叠加成像。
2.根据权利要求1所述的一种变偏移距VSP转换横波单程波叠前深度偏移方法，其特征在于：所述步骤S1包括：
S101.输入变偏移距VSP共检波点反射转换横波波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；
S202.输入纵波2维网格层速度地质模型{Vp}，输入横波2维网格层速度地质模型{Vs}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；
S203.输入成像算子的最大倾角dip、成像终止深度zmig2、成像深度步长dzmig、变偏移距VSP转换横波波场的最小频率fl、变偏移距VSP转换横波波场的最大频率参数fh。
3.根据权利要求2所述的一种变偏移距VSP转换横波单程波叠前深度偏移方法，其特征在于：所述步骤S2包括：
利用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：
炮点在2维地质模型中的坐标为：
Figure FDA0002572148530000011
MSZi＝SZi
其中，MSXi是第i个炮点在2维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于2维地质模型的深度、表示起伏地表，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；
检波点在2维地质模型中的坐标为：
Figure FDA0002572148530000021
MRZj＝RZj
其中，MRXj是第j个检波点在2维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。
4.根据权利要求3所述的一种变偏移距VSP转换横波单程波叠前深度偏移方法，其特征在于：所述步骤S3包括：
S301.将步骤S1中第j个检波点的反射转换横波时间空间域波场VSPDataj做快速傅里叶变换转换到频率波数域FVSPj；
S302.在步骤S2的坐标系下，将频率域震源FSouj设置在第j个检波点处。
5.根据权利要求1所述的一种变偏移距VSP转换横波单程波叠前深度偏移方法，其特征在于：所述步骤S4中，对于第j个检波点，完成变偏移距VSP转换横波单程波叠前深度偏移的过程包括：
S401.计算深度成像网格：
kzmig＝(k-1)·dzmig
Figure FDA0002572148530000022
其中，zmig2是成像终止深度，dzmig是成像深度步长，kzmig是第k个成像深度；
S402.频率波数域变偏移距VSP共检波点反射转换横波波场
Figure FDA0002572148530000023
用纵波速度向下延拓一个深度步长：
Figure FDA0002572148530000024
f＝[fl:df:fh]
Figure FDA0002572148530000031
kk＝f/Vp
g0＝1
g1＝-i·π·dzmig/kk
g2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4
g3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6
g4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8
Figure FDA0002572148530000032
Figure FDA0002572148530000033
Figure FDA0002572148530000034
Figure FDA0002572148530000035
Figure FDA0002572148530000036
Figure FDA0002572148530000037
其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，表示频率从fl到fh间隔为df，Vp纵波速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mute是波数切除因子，kk、g0-g4、t0-t4是中间变量，
Figure FDA0002572148530000038
是
Figure FDA0002572148530000039
向下延拓一个深度步长的波场；
S403.若kzmig＜＝MRZj，循环执行步骤S401～S402；
若kzmig＞MRZj，则跳转至步骤S404；
S404.频率域震源
Figure FDA00025721485300000310
用横波速度向下延拓一个深度步长：
Figure FDA00025721485300000311
f＝[fl:df:fh]
Figure FDA0002572148530000041
kks＝f/Vs
gs0＝1
gs1＝i·π·dzmig/kks
gs2＝(i·π·dzmig·kks/4-π2·dzmig2·kks2/2)/kks4
gs3＝(i·π·dzmig·kks/8+π2·dzmig2·kks2/4+i·π3·dzmig3·kks3/6)/kks6
gs4＝(i·π·dzmig·kks·5/64-π2·dzmig2·kks2·5/32-i·π3·dzmig3·kks3/8+π4·dzmig4·kks4/24)/kks8
Figure FDA0002572148530000042
Figure FDA0002572148530000043
Figure FDA0002572148530000044
Figure FDA0002572148530000045
Figure FDA0002572148530000046
Figure FDA0002572148530000047
其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，Vs横波速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，
Figure FDA0002572148530000048
是
Figure FDA0002572148530000049
向下延拓一个深度步长的波场；
S405.相关成像条件提取成像值：
Figure FDA00025721485300000410
其中，
Figure FDA00025721485300000411
是
Figure FDA00025721485300000412
向下延拓一个深度步长的波场，
Figure FDA00025721485300000413
是
Figure FDA00025721485300000414
向下延拓一个深度步长的波场，
Figure FDA00025721485300000415
是提取的kzmig+dzmig深度的成像值，conj是复共轭函数；
S406.循环执行步骤S401～S405，直至zmig2是成像终止深度，完成第j个检波点的变偏移距VSP转换横波单程波叠前深度偏移。
6.一种变偏移距VSP转换横波单程波叠前深度偏移装置，采用如权利要求1～5中任意一项所述的方法，其特征在于：包括：
数据输入模块，用于输入变偏移距VSP共检波点反射转换横波波场、纵波2维地质模型、横波2维地质模型和偏移参数；
坐标建立模块，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；
单程波叠前深度偏移成像模块，用于对每一个检波点的反射转换横波时间空间域波场，进行快速傅里叶变换到频率波数域，并设置频率域震源，将频率域震源从检波点深度用横波速度向下延拓，频率波数域变偏移距VSP共检波点反射转换横波波场用纵波速度从地表向下延拓，并实现每一个检波点的单程波叠前深度偏移；
叠加成像模块，用于将得到的叠前深度偏移重排成共成像道集，共成像道集叠加，得到变偏移VSP转换横波单程波叠前深度偏移叠加成像。
CN202010643476.2A 2020-07-07 2020-07-07 一种变偏移距vsp转换横波单程波叠前深度偏移方法和装置 Active CN111751876B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010643476.2A CN111751876B (zh) 2020-07-07 2020-07-07 一种变偏移距vsp转换横波单程波叠前深度偏移方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010643476.2A CN111751876B (zh) 2020-07-07 2020-07-07 一种变偏移距vsp转换横波单程波叠前深度偏移方法和装置

Publications (2)

Publication Number Publication Date
CN111751876A true CN111751876A (zh) 2020-10-09
CN111751876B CN111751876B (zh) 2022-05-20

Family

ID=72679635

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010643476.2A Active CN111751876B (zh) 2020-07-07 2020-07-07 一种变偏移距vsp转换横波单程波叠前深度偏移方法和装置

Country Status (1)

Country Link
CN (1) CN111751876B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112782767A (zh) * 2020-12-26 2021-05-11 中油奥博(成都)科技有限公司 一种das采集vsp变偏移距波场径向补偿方法与装置

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5696735A (en) * 1994-10-19 1997-12-09 Exxon Production Research Company Seismic migration using offset checkshot data
US20040220743A1 (en) * 2003-04-30 2004-11-04 Conocophillips Company Method for determining shear-wave velocity model for depth migration of mode-converted data
CN101419292A (zh) * 2007-10-25 2009-04-29 中国石油天然气集团公司 采用纵波源多分量地震数据生成横波地震剖面的方法
CN101984367A (zh) * 2010-10-18 2011-03-09 中联煤层气国家工程研究中心有限责任公司 一种转换波叠前成像方法和装置
CN103487834A (zh) * 2013-09-09 2014-01-01 中国石油集团川庆钻探工程有限公司地球物理勘探公司 转换波共检波点叠加静校正方法
CN103630934A (zh) * 2012-08-23 2014-03-12 中国石油天然气集团公司 一种确定转换波检波点大的横波静校正量的方法
CN104155694A (zh) * 2014-08-28 2014-11-19 中国石油集团川庆钻探工程有限公司地球物理勘探公司 一种反射转换横波共检波点叠加剖面的剩余静校正方法
CN105005075A (zh) * 2015-06-25 2015-10-28 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于地震频率信息的多波匹配方法
CN105445785A (zh) * 2014-06-12 2016-03-30 中国石油化工股份有限公司 一种横波偏移速度建模方法
US20170371054A1 (en) * 2016-06-24 2017-12-28 Pgs Geophysical As Migrating a Horizontal Component of a Wavefield
CN108845351A (zh) * 2018-06-26 2018-11-20 中国石油大学(华东) 一种vsp地震资料转换波全波形反演方法
CN209821405U (zh) * 2019-06-14 2019-12-20 中油奥博(成都)科技有限公司 地面分布式光纤三分量地面地震数据采集系统
CN110749927A (zh) * 2019-11-19 2020-02-04 中油奥博(成都)科技有限公司 光纤声波传感正交偶极声波测井系统及其测量方法

Patent Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5696735A (en) * 1994-10-19 1997-12-09 Exxon Production Research Company Seismic migration using offset checkshot data
US20040220743A1 (en) * 2003-04-30 2004-11-04 Conocophillips Company Method for determining shear-wave velocity model for depth migration of mode-converted data
CN101419292A (zh) * 2007-10-25 2009-04-29 中国石油天然气集团公司 采用纵波源多分量地震数据生成横波地震剖面的方法
CN101984367A (zh) * 2010-10-18 2011-03-09 中联煤层气国家工程研究中心有限责任公司 一种转换波叠前成像方法和装置
CN103630934A (zh) * 2012-08-23 2014-03-12 中国石油天然气集团公司 一种确定转换波检波点大的横波静校正量的方法
CN103487834A (zh) * 2013-09-09 2014-01-01 中国石油集团川庆钻探工程有限公司地球物理勘探公司 转换波共检波点叠加静校正方法
CN105445785A (zh) * 2014-06-12 2016-03-30 中国石油化工股份有限公司 一种横波偏移速度建模方法
CN104155694A (zh) * 2014-08-28 2014-11-19 中国石油集团川庆钻探工程有限公司地球物理勘探公司 一种反射转换横波共检波点叠加剖面的剩余静校正方法
CN105005075A (zh) * 2015-06-25 2015-10-28 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于地震频率信息的多波匹配方法
US20170371054A1 (en) * 2016-06-24 2017-12-28 Pgs Geophysical As Migrating a Horizontal Component of a Wavefield
CN108845351A (zh) * 2018-06-26 2018-11-20 中国石油大学(华东) 一种vsp地震资料转换波全波形反演方法
CN209821405U (zh) * 2019-06-14 2019-12-20 中油奥博(成都)科技有限公司 地面分布式光纤三分量地面地震数据采集系统
CN110749927A (zh) * 2019-11-19 2020-02-04 中油奥博(成都)科技有限公司 光纤声波传感正交偶极声波测井系统及其测量方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
孙涛: ""分步傅里叶转换波叠前深度偏移方法研究"", 《中国优秀博硕士学位论文全文数据库(硕士)基础科学辑》 *
岳玉波等: "转换波Kirchhoff叠前时间偏移的成像优化方案", 《地球物理学报》 *
康亮: "起伏地表多波共炮记录叠前深度偏移", 《内蒙古石油化工》 *
李建国等: "Walkaway VSP自由表面多次波叠前深度偏移成像", 《石油地球物理勘探》 *
李建国等: "复杂地表变偏移距VSP角度域旅行时层析及叠前深度偏移（英文）", 《APPLIED GEOPHYSICS》 *
董海杰等: "巷道掘进中地震超前探测波场特征研究", 《煤矿安全》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112782767A (zh) * 2020-12-26 2021-05-11 中油奥博(成都)科技有限公司 一种das采集vsp变偏移距波场径向补偿方法与装置

Also Published As

Publication number Publication date
CN111751876B (zh) 2022-05-20

Similar Documents

Publication Publication Date Title
CN102053261B (zh) 一种地震数据处理方法
CN103995288B (zh) 一种高斯束叠前深度偏移方法及装置
US20020103606A1 (en) Method and system for deghosting
CN105182408A (zh) 一种合成地震记录的制作方法和装置
AU2012254102A1 (en) Two component source seismic acquisition and source de-ghosting
WO2010082126A2 (en) Processing seismic data
CN104533396A (zh) 一种远探测声波的处理方法
Colombo et al. Near-surface full-waveform inversion in a transmission surface-consistent scheme
CN111751876B (zh) 一种变偏移距vsp转换横波单程波叠前深度偏移方法和装置
CN105044779B (zh) 基于相控接收指向性的反射界面方位定量判定方法及装置
EP3004939B1 (en) Device and method for velocity function extraction from the phase of ambient noise
CN106257309A (zh) 叠后地震数据体处理方法及装置
Oren et al. An overview of reproducible 3D seismic data processing and imaging using Madagascar
CN115469362B (zh) 一种地震勘探中的能流密度矢量计算方法
CN110579799A (zh) 一种等旅行时间间隔的地震采集观测方法及系统
Lord et al. A source-synchronous filter for uncorrelated receiver traces from a swept-frequency seismic source
RU2747628C1 (ru) Способ определения углов наклона отражающих границ по данным МОГТ 2D
CN114721044A (zh) 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统
CN111781647B (zh) 一种大斜井炮检移动vsp自由表面多次波成像方法和装置
Panea et al. Analysis of crooked-line 2D seismic reflection data recorded in areas with complex surface and subsurface conditions
Wege et al. Field and Synthetic Waveform Tests on Using Large‐Offset Seismic Streamer Data to Derive Shallow Seabed Shear‐Wave Velocity and Geotechnical Properties
Ding et al. Reverse-time ray-tracing method for microseismic source localization
CN117826248B (zh) 基于多尺度观测背景噪声聚束的面波频散提取方法
Giustiniani et al. Imaging subsurface structures using wave equation datuming advanced seismic techniques
Bao et al. A new technique to support future energy exploration of continental sedimentary basin: BWH full-frequency fidelity and amplitude preserving processing technology

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