CN111781647B - Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well - Google Patents

Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well Download PDF

Info

Publication number
CN111781647B
CN111781647B CN202010670672.9A CN202010670672A CN111781647B CN 111781647 B CN111781647 B CN 111781647B CN 202010670672 A CN202010670672 A CN 202010670672A CN 111781647 B CN111781647 B CN 111781647B
Authority
CN
China
Prior art keywords
dzmig
imaging
depth
shot
wave
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.)
Active
Application number
CN202010670672.9A
Other languages
Chinese (zh)
Other versions
CN111781647A (en
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 CN202010670672.9A priority Critical patent/CN111781647B/en
Publication of CN111781647A publication Critical patent/CN111781647A/en
Application granted granted Critical
Publication of CN111781647B publication Critical patent/CN111781647B/en
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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • 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

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

The invention discloses a method and a device for imaging a free surface multiple of a shot-inspection mobile VSP (vertical seismic profiling) of a highly deviated well, wherein the method comprises the following steps of: s1, inputting a shot-inspection mobile VSP full wavefield common detection point gather, a 2-dimensional geological model and offset parameters of a highly inclined shaft; s2, establishing coordinates of the shot point and the demodulator probe in the 2-dimensional geological model according to input data; s3, performing fast Fourier transform on the input common detection wave point gather to a frequency wave number domain, and setting a frequency domain seismic source; s4, free surface multiple single-pass wave pre-stack depth migration imaging of a single detection point; s5, completing shot detection moving VSP free surface multiple pre-stack depth migration imaging of all detection points; and S6, stacking the common imaging gathers to obtain shot-inspection mobile VSP free surface multiple pre-stack depth migration stacked imaging. The method effectively utilizes rich information carried by the free surface multiples, realizes shot-inspection moving VSP free surface multiples prestack depth migration imaging, and provides an accurate data basis for geological research and well drilling guidance.

Description

一种大斜井炮检移动VSP自由表面多次波成像方法和装置A method and device for multiple wave imaging of mobile VSP free surface for high-inclined well inspection

技术领域technical field

本发明涉及地球物理勘探中地震数据偏移成像，特别是涉及一种大斜井炮检移动VSP自由表面多次波成像方法和装置。The invention relates to seismic data migration imaging in geophysical exploration, in particular to a method and device for multiple wave imaging of a mobile VSP free surface for heavy-inclined well inspection.

背景技术Background technique

炮检移动VSP是针对大斜井设计的观测系统，即炮线布设于地表、位于井中检波器的正上方，炮线随着检波器移动。炮检移动VSP的优势是覆盖次更加均匀。这种观测系统主要用于斜井钻进过程中的导向入靶。通常，自由表面多次波作为干扰波被衰减了。与反射相比，自由表面多次波具有易识别、照明宽、反射角小、传播路径长的特点。Shot detection mobile VSP is an observation system designed for highly deviated wells, that is, the shot line is arranged on the surface, just above the geophone in the well, and the shot line moves with the geophone. The advantage of offsetting mobile VSPs is more uniform coverage. This kind of observation system is mainly used for guiding and targeting during the drilling of inclined wells. Usually, free surface multiples are attenuated as interference waves. Compared with reflection, free surface multiples have the characteristics of easy identification, wide illumination, small reflection angle and long propagation path.

吴世萍等研究了《Walkaway VSP多次波成像技术研究》，文献将地震相干成像用于Walkaway VSP多次波成像。李建国等研究了《Walkaway VSP自由表面多次波叠前深度偏移成像》，文献采用单程波延拓实现了变偏移距VSP的自由表面多次波成像；但是由于针对炮检移动VSP的自由表面多次波成像的研究还不成熟，在地址研究、钻进导向等领域还存在着诸多不便。Shiping Wu et al. studied "Research on Walkaway VSP Multiple Wave Imaging Technology", and the literature uses seismic coherence imaging for Walkaway VSP multiple wave imaging. Li Jianguo et al. studied "Walkaway VSP Free Surface Multiple Wave Prestack Depth Migration Imaging". The literature uses one-way wave continuation to realize free surface multiple imaging of variable offset VSP; The research on surface multiple wave imaging is still immature, and there are still many inconveniences in the fields of site research and drilling guidance.

发明内容SUMMARY OF THE INVENTION

本发明的目的在于克服现有技术的不足，提供一种大斜井炮检移动VSP自由表面多次波成像方法和装置，有效利用自由表面多次波所携带的丰富信息，实现了炮检移动VSP自由表面多次波叠前深度偏移成像，为地质研究、钻井导向提供了准确的数据基础。The purpose of the present invention is to overcome the deficiencies of the prior art, and to provide a method and device for imaging mobile VSP free surface multiples in a large-inclined well, which effectively utilizes the rich information carried by the free surface multiples to realize the mobile inspection. VSP free surface multiple wave prestack depth migration imaging provides an accurate data basis for geological research and drilling steering.

本发明的目的是通过以下技术方案来实现的：一种大斜井炮检移动VSP自由表面多次波成像方法，包括以下步骤：The object of the present invention is to be achieved through the following technical solutions: a method for imaging mobile VSP free surface multiples of large deviated wells, comprising the following steps:

S1.输入大斜井炮检移动VSP全波场共检波点道集、2维地质模型、偏移参数；S1. Input the high-inclined well shot detection mobile VSP full-wave field co-detection point gather, 2D geological model, and migration parameters;

S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；S2. According to the input data, establish the coordinates of the shot point and the detection point in the 2D geological model;

S3.对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源；S3. Perform fast Fourier transform on the input common detection point gathers to the frequency wavenumber domain, and set the frequency domain source;

S4.将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动VSP共检波点全波场、震源波场从自由表面向下延拓，实现单检波点的自由表面多次波单程波叠前深度偏移成像；S4. Extend the seismic source in the frequency domain from the depth of the receiver point upward to the free surface, the frequency wavenumber domain offsets the full wave field of the moving VSP common receiver point, and the source wave field extends downward from the free surface to realize the free surface of a single receiver point. Sub-wave one-way pre-stack depth migration imaging;

S5.循环执行步骤S4，完成所有检波点的炮检移动VSP自由表面多次波叠前深度偏移成像；S5. Execute step S4 in a loop to complete the detection of all the detection points and move the VSP free surface multiple wave pre-stack depth migration imaging;

S6.将步骤S5的成像重排成共成像道集，共成像道集叠加，得到炮检移动VSP自由表面多次波叠前深度偏移叠加成像。S6. Rearrange the imaging in step S5 into common imaging gathers, and superimpose the common imaging gathers to obtain multiple pre-stack depth migration stacking imaging of the offset moving VSP free surface.

进一步地，所述步骤S1包括：Further, the step S1 includes:

S101.输入大斜井炮检移动VSP全波场共检波点道集{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；S101. Input the mobile VSP full-wave field co-detection point gather {VSPData}, read the shot coordinates {SX, SY}, the shot depth SZ, the receiver coordinates {RX, RY}, the receiver depth RZ, geodetic coordinates of wellhead {WMX, WMY};

S102.输入纵波2维网格层速度地质模型{V}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；S102. Input the longitudinal wave 2-dimensional grid layer velocity geological model {V}, model abscissa MX, model abscissa MZ, wellhead model coordinate MWMX information;

S103.输入成像算子最大倾角dip、成像终止深度zmig2、成像深度步长dzmig、波场{VSPData}最小频率fl、波场{VSPData}最大频率参数fh。S103. Input the maximum dip angle dip of the imaging operator, the imaging termination depth zmig2, the imaging depth step size dzmig, the minimum frequency fl of the wave field {VSPData}, and the maximum frequency parameter fh of the wave field {VSPData}.

所述步骤S2包括：The step S2 includes:

S201.利用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：S201. Using the shot coordinates, shot depth, detection point coordinates, detection point depth, wellhead geodetic coordinates, and wellhead model coordinates in step S1 to establish the coordinates of the shot point and the detection point in the 2-dimensional geological model:

炮点在2维地质模型中的坐标为：The coordinates of the shot point in the 2D geological model are:

Figure GDA0002618028090000021
Figure GDA0002618028090000021

MSZi＝SZi MSZ i =SZ i

其中，MSXi是第i个炮点在2维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于2维地质模型的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；Among them, MSX i is the coordinate of the ith shot in the 2D geological model, SX i and SY i are the geodetic coordinates of the ith shot, and SZ i is the ith shot relative to the 2D geological model. Depth, WMX, WMY are the geodetic coordinates of the wellhead, MWMX is the wellhead model coordinates, and sign is the sign function;

检波点在2维地质模型中的坐标为：The coordinates of the detection point in the 2D geological model are:

Figure GDA0002618028090000022
Figure GDA0002618028090000022

MRZj＝RZj MRZj = RZj

其中，MRXj是第j个检波点在2维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。Among them, MRX j is the coordinate of the j-th detection point in the 2D geological model, RX j and RY j are the geodetic coordinates of the j-th detection point, RZ j is the depth of the j-th detection point, and WMX and WMY are the wellhead Geodetic coordinates, MWMX is the wellhead model coordinates, and sign is the sign function.

所述步骤S3包括：The step S3 includes:

S301.将步骤S1输入的第j个检波点的时间空间域全波场VSPDataj做快速傅里叶变换转换到频率波数域FVSPj；S301. the time-space domain full-wave field VSPData j of the j-th detection point input in step S1 is done fast Fourier transform and converted to frequency wavenumber domain FVSP j ;

S302.在步骤S2的坐标系下，将频率域震源FSouj设置在第j个检波点处。S302. In the coordinate system of step S2, set the frequency domain source FSou j at the jth detection point.

所述步骤S4包括：The step S4 includes:

S401.计算频率域震源向上延拓网格：S401. Calculate the upward extension grid of the frequency domain source:

kzmig＝(k-1)·dzmigkzmig=(k-1)·dzmig

Figure GDA0002618028090000031
Figure GDA0002618028090000031

其中，MRZj是第j个检波点的模型深度，dzmig是成像深度步长，kzmig是第k个向上延拓的深度；Among them, MRZ j is the model depth of the j-th detection point, dzmig is the imaging depth step, and kzmig is the depth of the k-th upward extension;

S402.频率域震源

Figure GDA00026180280900000312
用向上延拓一个深度步长：S402. Frequency Domain Source
Figure GDA00026180280900000312
Extend one depth step upwards:

Figure GDA0002618028090000032
Figure GDA0002618028090000032

f＝[fl:df:fh]f=[f l :df:f h ]

Figure GDA0002618028090000033
Figure GDA0002618028090000033

kk＝f/Vk=f/V

gs0＝1g s0 =1

gs1＝-i·π·dzmig/kkg s1 =-i·π·dzmig/kk

gs2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4 g s2 =(-i·π·dzmig·kk/4-π 2 ·dzmig 2 ·kk 2 /2)/kk 4

gs3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6 g s3 =(-i·π·dzmig·kk/8+π 2 ·dzmig 2 ·kk 2 /4-i·π 3 ·dzmig 3 ·kk 3 /6)/kk 6

gs4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8 g s4 =(-i·π · dzmig · kk · 5/64-π2·dzmig2·kk2·5/32+i·π3· dzmig3 · kk3 / 8 + π4 · dzmig4 ·kk 4/24 )/kk 8

Figure GDA0002618028090000034
Figure GDA0002618028090000034

Figure GDA0002618028090000035
Figure GDA0002618028090000035

Figure GDA0002618028090000036
Figure GDA0002618028090000036

Figure GDA0002618028090000037
Figure GDA0002618028090000037

Figure GDA0002618028090000038
Figure GDA0002618028090000038

Figure GDA0002618028090000039
Figure GDA0002618028090000039

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，表示频率从fl到fh间隔为df，V是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mutes是波数切除因子，kk、gs0-gs4、ts0-ts4是中间变量，

Figure GDA00026180280900000310
是
Figure GDA00026180280900000311
向上延拓一个深度步长的波场；where dip is the maximum inclination angle, fl is the minimum frequency, fh is the maximum frequency, f=[f l :df:f h ] is the frequency vector, indicating that the interval from f l to f h is df, V is the speed, k x is the wave number, dzmig is the imaging depth step, i is the imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function, km is the maximum wave number corresponding to the maximum dip angle, mutes is the wave number cutoff factor, kk , g s0 -g s4 , t s0 -t s4 are intermediate variables,
Figure GDA00026180280900000310
Yes
Figure GDA00026180280900000311
Extend the wavefield up by one depth step;

S403.循环执行步骤S401～S402实现频率域震源向上延拓至自由表面；S403. Steps S401-S402 are executed cyclically to realize the upward extension of the frequency domain source to the free surface;

S404.计算向下延拓成像网格：S404. Calculate the downward extension imaging grid:

kzmig＝(k-1)·dzmigkzmig=(k-1)·dzmig

Figure GDA0002618028090000041
Figure GDA0002618028090000041

其中，zmig2是成像终止深度，dzmig是成像深度步长，kzmig是第k个成像深度；Among them, zmig2 is the imaging termination depth, dzmig is the imaging depth step, and kzmig is the kth imaging depth;

S405.将频率波数域炮检移动VSP共检波点全波场

Figure GDA0002618028090000042
用向下延拓一个深度步长：S405. Detect the full wave field of the moving VSP co-detection point in the frequency wavenumber domain
Figure GDA0002618028090000042
Use downward continuation one depth step:

Figure GDA0002618028090000043
Figure GDA0002618028090000043

f＝[fl:df:fh]f=[f l :df:f h ]

Figure GDA0002618028090000044
Figure GDA0002618028090000044

kk＝f/Vk=f/V

g0＝1g 0 =1

g1＝-i·π·dzmig/kkg 1 =-i·π·dzmig/kk

g2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4 g 2 =(-i·π·dzmig·kk/4-π 2 ·dzmig 2 ·kk 2 /2)/kk 4

g3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6 g 3 =(-i·π·dzmig·kk/8+π 2 ·dzmig 2 ·kk 2 /4-i·π 3 ·dzmig 3 ·kk 3 /6)/kk 6

g4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8 g 4 =(-i·π·dzmig·kk·5/64-π 2 ·dzmig 2 ·kk 2 ·5/32+i·π 3 ·dzmig 3 ·kk 3 /8+π 4 ·dzmig 4 ·kk 4/24 )/kk 8

Figure GDA0002618028090000045
Figure GDA0002618028090000045

Figure GDA0002618028090000046
Figure GDA0002618028090000046

Figure GDA0002618028090000047
Figure GDA0002618028090000047

Figure GDA0002618028090000048
Figure GDA0002618028090000048

Figure GDA0002618028090000049
Figure GDA0002618028090000049

Figure GDA00026180280900000410
Figure GDA00026180280900000410

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，V是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mute是波数切除因子，kk、g0-g4、t0-t4是中间变量，

Figure GDA0002618028090000051
是
Figure GDA0002618028090000052
向下延拓一个深度步长的波场；where dip is the maximum dip, fl is the minimum frequency, fh is the maximum frequency, f=[ fl :df: fh ] is the frequency vector, V is the velocity, kx is the wavenumber, dzmig is the imaging depth step, and i is Imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function, km is the maximum wave number corresponding to the maximum dip angle dip, mute is the wave number cutoff factor, kk, g 0 -g 4 , t 0 -t 4 is the intermediate variable,
Figure GDA0002618028090000051
Yes
Figure GDA0002618028090000052
Extend the wavefield of a depth step down;

S406.将频率域震源

Figure GDA0002618028090000053
向下延拓一个深度步长：S406. The frequency domain source
Figure GDA0002618028090000053
Continue down one depth step:

Figure GDA0002618028090000054
Figure GDA0002618028090000054

f＝[fl:df:fh]f=[f l :df:f h ]

Figure GDA0002618028090000055
Figure GDA0002618028090000055

kk＝f/Vk=f/V

gs0＝1g s0 =1

gs1＝i·π·dzmig/kkg s1 =i·π·dzmig/kk

gs2＝(i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4 g s2 = (i·π·dzmig·kk/4-π 2 ·dzmig 2 ·kk 2 /2)/kk 4

gs3＝(i·π·dzmig·kk/8+π2·dzmig2·kk2/4+i·π3·dzmig3·kk3/6)/kk6 g s3 =(i·π·dzmig·kk/8+π 2 ·dzmig 2 ·kk 2 /4+i·π 3 ·dzmig 3 ·kk 3 /6)/kk 6

gs4＝(i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32-i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8 g s4 =(i·π · dzmig · kk · 5 /64-π2·dzmig2·kk2·5/32-i· π3 ·dzmig3·kk3/8+ π4 · dzmig4 · kk4 /24)/kk 8

Figure GDA0002618028090000056
Figure GDA0002618028090000056

Figure GDA0002618028090000057
Figure GDA0002618028090000057

Figure GDA0002618028090000058
Figure GDA0002618028090000058

Figure GDA0002618028090000059
Figure GDA0002618028090000059

Figure GDA00026180280900000510
Figure GDA00026180280900000510

Figure GDA00026180280900000511
Figure GDA00026180280900000511

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，V是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，

Figure GDA00026180280900000512
是
Figure GDA00026180280900000513
向下延拓一个深度步长的波场；where dip is the maximum dip, fl is the minimum frequency, fh is the maximum frequency, f=[ fl :df: fh ] is the frequency vector, V is the velocity, kx is the wavenumber, dzmig is the imaging depth step, and i is Imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function,
Figure GDA00026180280900000512
Yes
Figure GDA00026180280900000513
Extend the wavefield of a depth step down;

S407.相关成像条件提取成像值：S407. Extract imaging values from relevant imaging conditions:

Figure GDA00026180280900000514
Figure GDA00026180280900000514

其中，

Figure GDA00026180280900000515
是
Figure GDA00026180280900000516
向下延拓一个深度步长的波场，
Figure GDA00026180280900000517
是
Figure GDA00026180280900000518
向下延拓一个深度步长的波场，
Figure GDA00026180280900000519
是提取的kzmig+dzmig深度的成像值，conj是复共轭函数；in,
Figure GDA00026180280900000515
Yes
Figure GDA00026180280900000516
Extending the wavefield down a depth step,
Figure GDA00026180280900000517
Yes
Figure GDA00026180280900000518
Extending the wavefield down a depth step,
Figure GDA00026180280900000519
is the imaging value of the extracted kzmig+dzmig depth, and conj is the complex conjugate function;

S408.循环执行步骤S404～S407，至zmig2是成像终止深度，完成第j个检波点的炮检移动VSP自由表面多次波叠前深度偏移成像。S408. Steps S404 to S407 are performed cyclically until zmig2 is the imaging termination depth, and the offset imaging of the jth detection point is completed before moving the VSP free surface multiple wave stack depth migration.

一种大斜井炮检移动VSP自由表面多次波成像装置，包括：A mobile VSP free surface multiple wave imaging device for high-inclined well inspection, comprising:

数据输入模块，用于输入大斜井炮检移动VSP全波场共检波点道集、2维地质模型、偏移参数；The data input module is used to input the mobile VSP full-wave field co-detection point gather, 2D geological model, and migration parameters for the high-inclined well shot detection;

坐标建立模块，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；The coordinate establishment module is used to establish the coordinates of the shot point and the receiver point in the 2D geological model according to the input data;

多次波叠前成像模块，用于对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源；将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动VSP共检波点全波场、震源波场从自由表面向下延拓，实现单检波点的自由表面多次波单程波叠前深度偏移成像；并按照同样的方式完成所有检波点的炮检移动VSP自由表面多次波叠前深度偏移成像；The multiple wave prestack imaging module is used to perform fast Fourier transform on the input common receiver gathers to the frequency wavenumber domain, and set the frequency domain source; extend the frequency domain source from the depth of the receiver point upward to the free surface, and Wavenumber domain offsets the full wave field of the moving VSP common receiver point, and the source wave field extends downward from the free surface to realize the free surface multiple wave one-way wave prestack depth migration imaging of a single receiver point; and complete all the steps in the same way. Offsetting mobile VSP free surface multiple pre-stack depth migration imaging of the detection point;

叠加成像模块，用于将得到的多次波叠前深度偏移成像重排成共成像道集，共成像道集叠加，得到炮检移动VSP自由表面多次波叠前深度偏移叠加成像。The stacking imaging module is used for rearranging the obtained multiple pre-stack depth migration imaging into common imaging gathers, and stacking the co-imaging gathers to obtain the multiple pre-stack depth migration stacking imaging of the offset mobile VSP free surface.

本发明的有益效果是：本发明输入炮检移动VSP共检波点全波场、成像速度模型，给定频率域震源波场从检波点处向上延拓至自由表面，再从自由表面向下延拓共检波点全波场、上述的震源波场，采用相关成像条件提取成像值，实现炮检移动VSP自由表面多次波叠前深度偏移成像。The beneficial effects of the present invention are as follows: the present invention inputs the full wave field and imaging velocity model of the mobile VSP co-detection point, and the source wave field of a given frequency domain extends upward from the detection point to the free surface, and then extends downward from the free surface. Extend the full wave field of the common receiver point and the above-mentioned source wave field, and extract the imaging value by using the relevant imaging conditions to realize the multiple wave pre-stack depth migration imaging of the mobile VSP free surface.

附图说明Description of drawings

图1为本发明的方法流程图；Fig. 1 is the method flow chart of the present invention;

图2为实施例中大斜井炮检移动VSP全波场共检波点道集示意图；Fig. 2 is a schematic diagram of the collective detection point gathers of the mobile VSP full-wave field for artillery inspection of the highly deviated well in the embodiment;

图3为实施例中大斜井炮检移动VSP的检波点分布示意图；3 is a schematic diagram of the distribution of detection points of the mobile VSP in the high-inclined well shot detection in the embodiment;

图4为实施例中大斜井炮检移动VSP的炮点分布示意图；Fig. 4 is the shot distribution schematic diagram of the mobile VSP in the high-inclined well shot detection in the embodiment;

图5为实施例中成像速度模型示意图；5 is a schematic diagram of an imaging velocity model in an embodiment;

图6为实施例中炮检移动VSP自由表面多次波叠前深度偏移成像示意图。FIG. 6 is a schematic diagram of depth-migration imaging before multiple wave stacking of multiple waves of an offset mobile VSP free surface in an embodiment.

图7为实施例中炮检移动VSP自由表面多次波叠前深度偏移叠加成像示意图；7 is a schematic diagram of depth-migration stacking imaging before multiple wave stacking of mobile VSP free surface multiple waves in the embodiment;

图8为本发明的装置原理框图。FIG. 8 is a schematic block diagram of the apparatus of the present invention.

具体实施方式Detailed ways

下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。The technical solutions of the present invention are further described in detail below with reference to the accompanying drawings, but the protection scope of the present invention is not limited to the following.

如图1所示，一种大斜井炮检移动VSP自由表面多次波成像方法，包括以下步骤：As shown in Fig. 1, a method for multiple wave imaging of mobile VSP free surface for high-inclined well inspection includes the following steps:

S1.输入大斜井炮检移动VSP全波场共检波点道集、2维地质模型、偏移参数：S1. Input the high-inclined well shot detection mobile VSP full-wave field co-detection point gather, 2D geological model, and migration parameters:

S101.输入大斜井炮检移动VSP全波场共检波点道集{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；S101. Input the mobile VSP full-wave field co-detection point gather {VSPData}, read the shot coordinates {SX, SY}, the shot depth SZ, the receiver coordinates {RX, RY}, the receiver depth RZ, geodetic coordinates of wellhead {WMX, WMY};

S102.输入纵波2维网格层速度地质模型{V}，模型横坐标MX、模型横纵坐标MZ、井口模型坐标MWMX信息；S102. Input the longitudinal wave 2-dimensional grid layer velocity geological model {V}, model abscissa MX, model abscissa MZ, wellhead model coordinate MWMX information;

S103.输入成像算子最大倾角dip、成像终止深度zmig2、成像深度步长dzmig、波场{VSPData}最小频率fl、波场{VSPData}最大频率参数fh。S103. Input the maximum dip angle dip of the imaging operator, the imaging termination depth zmig2, the imaging depth step size dzmig, the minimum frequency fl of the wave field {VSPData}, and the maximum frequency parameter fh of the wave field {VSPData}.

在本申请的实施例中，输入的大斜井炮检移动VSP全波场共检波点道集如图2所示，图2中横坐标为道号，纵坐标为时间(单位：毫秒)；输入的大斜井炮检移动VSP的检波点分布如图3所示，图3中三角形为检波点所在位置，横坐标为长度(单位：米)，纵坐标为深度(单位：米)；输入的大斜井炮检移动VSP的炮点分布如图4所示，图4中，点线为炮线位置；横坐标为长度(单位：米)；纵坐标表示第几次激发；输入的成像速度模型如图5所示，图5中，横坐标为长度(单位：米)；纵坐标为深度(单位：米)。In the embodiment of the present application, the input high-inclined well shot detection mobile VSP full-wave field co-detection point gather is shown in Figure 2, where the abscissa in Figure 2 is the track number, and the ordinate is the time (unit: millisecond); Figure 3 shows the distribution of the detection points of the input high-inclined well shot detection mobile VSP. In Figure 3, the triangle is the location of the detection point, the abscissa is the length (unit: m), and the ordinate is the depth (unit: m); input The shot distribution of the mobile VSP in the high-inclined well inspection is shown in Fig. 4. In Fig. 4, the dotted line is the position of the shot line; the abscissa is the length (unit: m); the ordinate indicates the number of excitations; the input imaging The velocity model is shown in Figure 5. In Figure 5, the abscissa is the length (unit: m); the ordinate is the depth (unit: m).

S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；S2. According to the input data, establish the coordinates of the shot point and the detection point in the 2D geological model;

S201.利用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：S201. Using the shot coordinates, shot depth, detection point coordinates, detection point depth, wellhead geodetic coordinates, and wellhead model coordinates in step S1 to establish the coordinates of the shot point and the detection point in the 2-dimensional geological model:

炮点在2维地质模型中的坐标为：The coordinates of the shot point in the 2D geological model are:

Figure GDA0002618028090000071
Figure GDA0002618028090000071

MSZi＝SZi MSZ i =SZ i

其中，MSXi是第i个炮点在2维地质模型中的坐标，SXi、SYi是第i个炮点的大地坐标，SZi是第i个炮点的相对于2维地质模型的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；Among them, MSX i is the coordinate of the ith shot in the 2D geological model, SX i and SY i are the geodetic coordinates of the ith shot, and SZ i is the ith shot relative to the 2D geological model. Depth, WMX, WMY are the geodetic coordinates of the wellhead, MWMX is the wellhead model coordinates, and sign is the sign function;

检波点在2维地质模型中的坐标为：The coordinates of the detection point in the 2D geological model are:

Figure GDA0002618028090000072
Figure GDA0002618028090000072

MRZj＝RZj MRZj = RZj

其中，MRXj是第j个检波点在2维地质模型中坐标，RXj、RYj是第j个检波点的大地坐标，RZj是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。Among them, MRX j is the coordinate of the j-th detection point in the 2D geological model, RX j and RY j are the geodetic coordinates of the j-th detection point, RZ j is the depth of the j-th detection point, and WMX and WMY are the wellhead Geodetic coordinates, MWMX is the wellhead model coordinates, and sign is the sign function.

S3.对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源：S3. Perform fast Fourier transform on the input common-detection point gathers to the frequency wavenumber domain, and set the frequency domain source:

所述步骤S3包括：The step S3 includes:

S301.将步骤S1输入的第j个检波点的时间空间域全波场VSPDataj做快速傅里叶变换转换到频率波数域FVSPj；S301. the time-space domain full-wave field VSPData j of the j-th detection point input in step S1 is done fast Fourier transform and converted to frequency wavenumber domain FVSP j ;

S302.在步骤S2的坐标系下，将频率域震源FSouj设置在第j个检波点处。S302. In the coordinate system of step S2, set the frequency domain source FSou j at the jth detection point.

S4.将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动VSP共检波点全波场、震源波场从自由表面向下延拓，延拓算子为ω-x单程波算子，相关成像条件提取成像值，实现单检波点的自由表面多次波单程波叠前深度偏移成像：S4. Extend the seismic source in the frequency domain from the depth of the receiver point upward to the free surface, the frequency wavenumber domain offsets the full wave field of the moving VSP common receiver point, and the source wave field extends downward from the free surface, and the extension operator is ω-x The one-way wave operator, extracts the imaging value from the relevant imaging conditions, and realizes the free surface multiple wave one-way wave pre-stack depth migration imaging of a single detection point:

S401.计算频率域震源向上延拓网格：S401. Calculate the upward extension grid of the frequency domain source:

kzmig＝(k-1)·dzmigkzmig=(k-1)·dzmig

Figure GDA0002618028090000081
Figure GDA0002618028090000081

其中，MRZj是第j个检波点的模型深度，dzmig是成像深度步长，kzmig是第k个向上延拓的深度；Among them, MRZ j is the model depth of the j-th detection point, dzmig is the imaging depth step, and kzmig is the depth of the k-th upward extension;

S402.频率域震源

Figure GDA0002618028090000082
用向上延拓一个深度步长：S402. Frequency Domain Source
Figure GDA0002618028090000082
Extend one depth step upwards:

Figure GDA0002618028090000083
Figure GDA0002618028090000083

f＝[fl:df:fh]f=[f l :df:f h ]

Figure GDA0002618028090000084
Figure GDA0002618028090000084

kk＝f/Vk=f/V

gs0＝1g s0 =1

gs1＝-i·π·dzmig/kkg s1 =-i·π·dzmig/kk

gs2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4 g s2 =(-i·π·dzmig·kk/4-π 2 ·dzmig 2 ·kk 2 /2)/kk 4

gs3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6 g s3 =(-i·π·dzmig·kk/8+π 2 ·dzmig 2 ·kk 2 /4-i·π 3 ·dzmig 3 ·kk 3 /6)/kk 6

gs4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8 g s4 =(-i·π · dzmig · kk · 5/64-π2·dzmig2·kk2·5/32+i·π3· dzmig3 · kk3 / 8 + π4 · dzmig4 ·kk 4/24 )/kk 8

Figure GDA0002618028090000085
Figure GDA0002618028090000085

Figure GDA0002618028090000086
Figure GDA0002618028090000086

Figure GDA0002618028090000091
Figure GDA0002618028090000091

Figure GDA0002618028090000092
Figure GDA0002618028090000092

Figure GDA0002618028090000093
Figure GDA0002618028090000093

Figure GDA0002618028090000094
Figure GDA0002618028090000094

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，表示频率从fl到fh间隔为df，V是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mutes是波数切除因子，kk、gs0-gs4、ts0-ts4是中间变量，

Figure GDA0002618028090000095
是
Figure GDA0002618028090000096
向上延拓一个深度步长的波场；where dip is the maximum inclination angle, fl is the minimum frequency, fh is the maximum frequency, f=[f l :df:f h ] is the frequency vector, indicating that the interval from f l to f h is df, V is the speed, k x is the wave number, dzmig is the imaging depth step, i is the imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function, km is the maximum wave number corresponding to the maximum dip angle, mutes is the wave number cutoff factor, kk , g s0 -g s4 , t s0 -t s4 are intermediate variables,
Figure GDA0002618028090000095
Yes
Figure GDA0002618028090000096
Extend the wavefield up by one depth step;

S403.循环执行步骤S401～S402实现频率域震源向上延拓至自由表面；S403. Steps S401-S402 are executed cyclically to realize the upward extension of the frequency domain source to the free surface;

S404.计算向下延拓成像网格：S404. Calculate the downward extension imaging grid:

kzmig＝(k-1)·dzmigkzmig=(k-1)·dzmig

Figure GDA0002618028090000097
Figure GDA0002618028090000097

其中，zmig2是成像终止深度，dzmig是成像深度步长，kzmig是第k个成像深度；Among them, zmig2 is the imaging termination depth, dzmig is the imaging depth step, and kzmig is the kth imaging depth;

S405.将频率波数域炮检移动VSP共检波点全波场

Figure GDA0002618028090000098
用向下延拓一个深度步长：S405. Detect the full wave field of the moving VSP co-detection point in the frequency wavenumber domain
Figure GDA0002618028090000098
Use downward continuation one depth step:

Figure GDA0002618028090000099
Figure GDA0002618028090000099

f＝[fl:df:fh]f=[f l :df:f h ]

Figure GDA00026180280900000910
Figure GDA00026180280900000910

kk＝f/Vk=f/V

g0＝1g 0 =1

g1＝-i·π·dzmig/kkg 1 =-i·π·dzmig/kk

g2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4 g 2 =(-i·π·dzmig·kk/4-π 2 ·dzmig 2 ·kk 2 /2)/kk 4

g3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6 g 3 =(-i·π·dzmig·kk/8+π 2 ·dzmig 2 ·kk 2 /4-i·π 3 ·dzmig 3 ·kk 3 /6)/kk 6

g4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8 g 4 =(-i·π·dzmig·kk·5/64-π 2 ·dzmig 2 ·kk 2 ·5/32+i·π 3 ·dzmig 3 ·kk 3 /8+π 4 ·dzmig 4 ·kk 4/24 )/kk 8

Figure GDA0002618028090000101
Figure GDA0002618028090000101

Figure GDA0002618028090000102
Figure GDA0002618028090000102

Figure GDA0002618028090000103
Figure GDA0002618028090000103

Figure GDA0002618028090000104
Figure GDA0002618028090000104

Figure GDA0002618028090000105
Figure GDA0002618028090000105

Figure GDA0002618028090000106
Figure GDA0002618028090000106

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，V是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，km是最大倾角dip对应的最大波数，mute是波数切除因子，kk、g0-g4、t0-t4是中间变量，

Figure GDA0002618028090000107
是
Figure GDA0002618028090000108
向下延拓一个深度步长的波场；where dip is the maximum dip, fl is the minimum frequency, fh is the maximum frequency, f=[ fl :df: fh ] is the frequency vector, V is the velocity, kx is the wavenumber, dzmig is the imaging depth step, and i is Imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function, km is the maximum wave number corresponding to the maximum dip angle dip, mute is the wave number cutoff factor, kk, g 0 -g 4 , t 0 -t 4 is the intermediate variable,
Figure GDA0002618028090000107
Yes
Figure GDA0002618028090000108
Extend the wavefield of a depth step down;

S406.将频率域震源

Figure GDA0002618028090000109
向下延拓一个深度步长：S406. The frequency domain source
Figure GDA0002618028090000109
Continue down one depth step:

Figure GDA00026180280900001010
Figure GDA00026180280900001010

f＝[fl:df:fh]f=[f l :df:f h ]

Figure GDA00026180280900001011
Figure GDA00026180280900001011

kk＝f/Vk=f/V

gs0＝1g s0 =1

gs1＝i·π·dzmig/kkg s1 =i·π·dzmig/kk

gs2＝(i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4 g s2 = (i·π·dzmig·kk/4-π 2 ·dzmig 2 ·kk 2 /2)/kk 4

gs3＝(i·π·dzmig·kk/8+π2·dzmig2·kk2/4+i·π3·dzmig3·kk3/6)/kk6 g s3 =(i·π·dzmig·kk/8+π 2 ·dzmig 2 ·kk 2 /4+i·π 3 ·dzmig 3 ·kk 3 /6)/kk 6

gs4＝(i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32-i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8 g s4 =(i·π · dzmig · kk · 5 /64-π2·dzmig2·kk2·5/32-i· π3 ·dzmig3·kk3/8+ π4 · dzmig4 · kk4 /24)/kk 8

Figure GDA00026180280900001012
Figure GDA00026180280900001012

Figure GDA00026180280900001013
Figure GDA00026180280900001013

Figure GDA00026180280900001014
Figure GDA00026180280900001014

Figure GDA00026180280900001015
Figure GDA00026180280900001015

Figure GDA0002618028090000111
Figure GDA0002618028090000111

Figure GDA0002618028090000112
Figure GDA0002618028090000112

其中，dip是最大倾角，fl是最小频率，fh是最大频率，f＝[fl:df:fh]是频率向量，V是速度，kx是波数，dzmig是成像深度步长，i是虚数单位，fft是傅里叶变换函数，ifft是傅里叶逆变换函数，

Figure GDA0002618028090000113
是
Figure GDA0002618028090000114
向下延拓一个深度步长的波场；where dip is the maximum dip, fl is the minimum frequency, fh is the maximum frequency, f=[ fl :df: fh ] is the frequency vector, V is the velocity, kx is the wavenumber, dzmig is the imaging depth step, and i is Imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function,
Figure GDA0002618028090000113
Yes
Figure GDA0002618028090000114
Extend the wavefield of a depth step down;

S407.相关成像条件提取成像值：S407. Extract imaging values from relevant imaging conditions:

Figure GDA0002618028090000115
Figure GDA0002618028090000115

其中，

Figure GDA0002618028090000116
是
Figure GDA0002618028090000117
向下延拓一个深度步长的波场，
Figure GDA0002618028090000118
是
Figure GDA0002618028090000119
向下延拓一个深度步长的波场，
Figure GDA00026180280900001110
是提取的kzmig+dzmig深度的成像值，conj是复共轭函数；in,
Figure GDA0002618028090000116
Yes
Figure GDA0002618028090000117
Extending the wavefield down a depth step,
Figure GDA0002618028090000118
Yes
Figure GDA0002618028090000119
Extending the wavefield down a depth step,
Figure GDA00026180280900001110
is the imaging value of the extracted kzmig+dzmig depth, and conj is the complex conjugate function;

S408.循环执行步骤S404～S407，至zmig2是成像终止深度，完成第j个检波点的炮检移动VSP自由表面多次波叠前深度偏移成像。S408. Steps S404 to S407 are performed cyclically until zmig2 is the imaging termination depth, and the offset imaging of the jth detection point is completed before moving the VSP free surface multiple wave stack depth migration.

S5.循环执行步骤S4，完成所有检波点的炮检移动VSP自由表面多次波叠前深度偏移成像；S5. Execute step S4 in a loop to complete the detection of all the detection points and move the VSP free surface multiple wave pre-stack depth migration imaging;

在本申请的实施例中，炮检移动VSP自由表面多次波叠前深度偏移成像如图6所示，与图2对应，横坐标为道号；纵坐标为深度(单位：米)；In the embodiment of the present application, the pre-stack depth migration imaging of multiple waves of the offset mobile VSP free surface is shown in Figure 6, corresponding to Figure 2, the abscissa is the track number; the ordinate is the depth (unit: meter);

S6.将步骤S5的成像重排成共成像道集，共成像道集叠加，得到炮检移动VSP自由表面多次波叠前深度偏移叠加成像。S6. Rearrange the imaging in step S5 into common imaging gathers, and superimpose the common imaging gathers to obtain multiple pre-stack depth migration stacking imaging of the offset moving VSP free surface.

在本申请的实施例中，炮检移动VSP自由表面多次波叠前深度偏移叠加成像，如图7所示，横坐标为长度(单位：米)；纵坐标为深度(单位：米)。In the embodiment of the present application, the multiple wave pre-stack depth migration stack imaging of the mobile VSP free surface is detected, as shown in Figure 7, the abscissa is the length (unit: m); the ordinate is the depth (unit: m) .

如图8所示，一种大斜井炮检移动VSP自由表面多次波成像装置，包括：As shown in Figure 8, a mobile VSP free surface multiple wave imaging device for high-inclined well inspection includes:

数据输入模块，用于输入大斜井炮检移动VSP全波场共检波点道集、2维地质模型、偏移参数；The data input module is used to input the mobile VSP full-wave field co-detection point gather, 2D geological model, and migration parameters for the high-inclined well shot detection;

坐标建立模块，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；The coordinate establishment module is used to establish the coordinates of the shot point and the receiver point in the 2D geological model according to the input data;

多次波叠前成像模块，用于对输入的共检波点道集做快速傅里叶变换到频率波数域，设置频率域震源；将频率域震源从检波点深度向上延拓至自由表面，频率波数域炮检移动VSP共检波点全波场、震源波场从自由表面向下延拓，实现单检波点的自由表面多次波单程波叠前深度偏移成像；并按照同样的方式完成所有检波点的炮检移动VSP自由表面多次波叠前深度偏移成像；The multiple wave prestack imaging module is used to perform fast Fourier transform on the input common receiver gathers to the frequency wavenumber domain, and set the frequency domain source; extend the frequency domain source from the depth of the receiver point upward to the free surface, and Wavenumber domain offsets the full wave field of the moving VSP common receiver point, and the source wave field extends downward from the free surface to realize the free surface multiple wave one-way wave prestack depth migration imaging of a single receiver point; and complete all the steps in the same way. Offsetting mobile VSP free surface multiple pre-stack depth migration imaging of the detection point;

叠加成像模块，用于将得到的多次波叠前深度偏移成像重排成共成像道集，共成像道集叠加，得到炮检移动VSP自由表面多次波叠前深度偏移叠加成像。The stacking imaging module is used for rearranging the obtained multiple pre-stack depth migration imaging into common imaging gathers, and stacking the co-imaging gathers to obtain the multiple pre-stack depth migration stacking imaging of the offset mobile VSP free surface.

综上，本发明输入炮检移动VSP共检波点全波场、成像速度模型，给定频率域震源波场从检波点处向上延拓至自由表面，再从自由表面向下延拓共检波点全波场、上述的震源波场，采用相关成像条件提取成像值，实现炮检移动VSP自由表面多次波叠前深度偏移成像。To sum up, the present invention inputs the full wave field and imaging velocity model of the offset moving VSP common detection point, and the source wave field of the given frequency domain extends upward from the detection point to the free surface, and then extends the common detection point downward from the free surface. The full wave field, the above-mentioned source wave field, and the relevant imaging conditions are used to extract the imaging values to realize the multiple wave pre-stack depth migration imaging of the offset mobile VSP free surface.

以上所述是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应该看作是对其他实施例的排除，而可用于其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。The above are preferred embodiments of the present invention, it should be understood that the present invention is not limited to the form disclosed herein, should not be regarded as an exclusion of other embodiments, but can be used in other combinations, modifications and environments, and can be used herein Within the scope of the stated concept, modifications can be made through the above teachings or skill or knowledge in the relevant field. However, modifications and changes made by those skilled in the art do not depart from the spirit and scope of the present invention, and should all fall within the protection scope of the appended claims of the present invention.

Claims (5)

1. A highly deviated well shot-inspection mobile VSP free surface multiple imaging method is characterized by comprising the following steps: the method comprises the following steps:
s1, inputting a shot-inspection mobile VSP full wavefield common-inspection point gather, a 2-dimensional geological model and offset parameters of a highly deviated well;
s2, establishing coordinates of the shot point and the wave detection point in the 2-dimensional geological model according to input data;
s3, performing fast Fourier transform on the input common detection wave point gather to a frequency wave number domain, and setting a frequency domain seismic source;
s4, extending a frequency domain seismic source upwards from the depth of a wave detection point to the free surface, extending a frequency wave number domain shot detection moving VSP common detection wave point full wave field and a seismic source wave field downwards from the free surface, and realizing free surface multiple-order single wave prestack depth migration imaging of a single detection point;
the step S4 includes:
s401, calculating an upward continuation grid of a seismic source in a frequency domain:
kzmig＝(k-1)·dzmig
Figure FDA0003583716160000011
wherein MRZjIs the model depth of the jth demodulator probe, dzmig is the imaging depth step, kzmig is the depth of the kth upward continuation;
s402. frequency domain seismic source
Figure FDA0003583716160000012
One depth step is extended upwards:
Figure FDA0003583716160000013
f＝[fl:df:fh]
Figure FDA0003583716160000014
kk＝f/V
gs0＝1
gs1＝-i·π·dzmig/kk
gs2＝(-i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4
gs3＝(-i·π·dzmig·kk/8+π2·dzmig2·kk2/4-i·π3·dzmig3·kk3/6)/kk6
gs4＝(-i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32+i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8
Figure FDA0003583716160000015
Figure FDA0003583716160000016
Figure FDA0003583716160000021
Figure FDA0003583716160000022
Figure FDA0003583716160000023
Figure FDA0003583716160000024
where dip is the maximum tilt angle, fl is the minimum frequency, fh is the maximum frequency, and f ═ fl:df:fh]Is a frequency vector representing the frequency from flTo fhInterval df, V velocity, kxIs wavenumber, dzmig is imaging depth step size, i is imaginaryNumber units, fft is the Fourier transform function, ifft is the inverse Fourier transform function, kmIs the maximum wave number corresponding to the maximum dip, mutes is the wave number cut-off factor, kk, gs0-gs4、ts0-ts4Is the intermediate variable that is the variable between,
Figure FDA0003583716160000025
is that
Figure FDA0003583716160000026
Extending a wave field of a depth step upwards;
s403, circularly executing the steps S401-S402 to extend the frequency domain seismic source upwards to the free surface;
s404, calculating a downward continuation imaging grid:
kzmig＝(k-1)·dzmig
Figure FDA0003583716160000027
where zmig2 is the imaging stop depth, dzmig is the imaging depth step, kzmig is the kth imaging depth;
s405, carrying out shot detection on the frequency wave number domain to move a VSP (vertical seismic profiling) common detection wave point full wave field
Figure FDA0003583716160000028
Extend one depth step down:
Figure FDA0003583716160000029
f＝[fl:df:fh]
Figure FDA00035837161600000210
kk＝f/V
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 FDA0003583716160000031
Figure FDA0003583716160000032
Figure FDA0003583716160000033
Figure FDA0003583716160000034
Figure FDA0003583716160000035
Figure FDA0003583716160000036
where dip is the maximum tilt angle, fl is the minimum frequency, fh is the maximum frequency, and f ═ fl:df:fh]Is a frequency vector, V is velocity, kxIs the wave number of the wave, and,dzmig is the imaging depth step, i is the imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function, kmIs the maximum wave number corresponding to the maximum dip, mute is the wave number cut-off factor, kk, g0-g4、t0-t4Is the intermediate variable that is the variable between,
Figure FDA0003583716160000037
is that
Figure FDA0003583716160000038
Extending a depth step wave field downwards;
s406, seismic source of frequency domain
Figure FDA0003583716160000039
Continue one depth step down:
Figure FDA00035837161600000310
f＝[fl:df:fh]
Figure FDA00035837161600000311
kk＝f/V
gs0＝1
gs1＝i·π·dzmig/kk
gs2＝(i·π·dzmig·kk/4-π2·dzmig2·kk2/2)/kk4
gs3＝(i·π·dzmig·kk/8+π2·dzmig2·kk2/4+i·π3·dzmig3·kk3/6)/kk6
gs4＝(i·π·dzmig·kk·5/64-π2·dzmig2·kk2·5/32-i·π3·dzmig3·kk3/8+π4·dzmig4·kk4/24)/kk8
Figure FDA00035837161600000312
Figure FDA00035837161600000313
Figure FDA00035837161600000314
Figure FDA00035837161600000315
Figure FDA0003583716160000041
Figure FDA0003583716160000042
where dip is the maximum tilt angle, fl is the minimum frequency, fh is the maximum frequency, and f ═ fl:df:fh]Is a frequency vector, V is velocity, kxIs the wavenumber, dzmig is the imaging depth step, i is the imaginary unit, fft is the Fourier transform function, ifft is the inverse Fourier transform function,
Figure FDA0003583716160000043
is that
Figure FDA0003583716160000044
Downwardly extending a depth step wave field;
s407, extracting an imaging value according to the relevant imaging condition:
Figure FDA0003583716160000045
wherein,
Figure FDA0003583716160000046
is that
Figure FDA0003583716160000047
The wavefield is extended down by one depth step,
Figure FDA0003583716160000048
is that
Figure FDA0003583716160000049
The wavefield is extended down by one depth step,
Figure FDA00035837161600000410
is the imaging value of the extracted kzmig + dzmig depth, conj is the complex conjugate function;
s408, circularly executing the steps S404 to S407 until the imaging termination depth zmig2 is extended, and finishing shot detection moving VSP free surface multiple pre-stack depth migration imaging of the jth wave detection point;
s5, circularly executing the step S4 to finish shot-inspection moving VSP free surface multiple prestack depth migration imaging of all inspection points;
s6, rearranging the imaging in the step S5 into a common imaging gather, and superposing the common imaging gather to obtain shot-geophone test moving VSP free surface multiple pre-stack depth migration superposition imaging.
2. The method of claim 1, wherein the method comprises the following steps: the step S1 includes:
s101, inputting a shot moving VSP full wave field common geophone gather { VSPDATA }, and reading shot coordinates { SX, SY }, shot depth SZ, geophone coordinates { RX, RY }, geophone depth RZ and well-head geodetic coordinates { WMX, WMY };
s102, inputting longitudinal wave 2-dimensional grid layer speed geological model { V }, model horizontal coordinate MX, model horizontal and vertical coordinate MZ and wellhead model coordinate MWMX information;
s103, inputting a maximum inclination dip of an imaging operator, an imaging termination depth zmig2, an imaging depth step dzmig, a wave field { VSPDATA } minimum frequency fl and a wave field { VSPDATA } maximum frequency parameter fh.
3. The method of claim 1, wherein the method comprises the following steps: the step S2 includes:
s201, establishing coordinates of the shot point and the wave detection point in the 2-dimensional geological model by using the shot point coordinates, the shot point depth, the wave detection point coordinates, the wave detection point depth, the wellhead geodetic coordinates and the wellhead model coordinates in the step S1:
the coordinates of the shot point in the 2-dimensional geological model are:
Figure FDA0003583716160000051
MSZi＝SZi
wherein, MSXiIs the coordinate of the ith shot point in the 2-dimensional geological model, SXi、SYiIs the geodetic coordinate of the ith shot point, SZiThe depth of the ith shot point relative to the 2-dimensional geological model, WMX and WMY are the geodetic coordinates of a wellhead, MWMX is the coordinates of the wellhead model, and sign is a symbolic function;
the coordinates of the demodulator probe in the 2-dimensional geological model are:
Figure FDA0003583716160000052
MRZj＝RZj
wherein MRXjIs the coordinates of the jth receiver point in the 2-dimensional geological model, RXj、RYjIs the geodetic coordinate of the jth detection point, RZjIs the jthAnd the depths of the wave detection points, WMX and WMY are geodetic coordinates of a wellhead, MWMX is a wellhead model coordinate, and sign is a symbolic function.
4. The method of claim 1, wherein the method comprises the following steps: the step S3 includes:
s301, inputting the time-space domain full wavefield VSPDATA of the j-th wave detection point input in the step S1jFast Fourier transform is carried out to convert the fast Fourier transform into a frequency wave number domain FVSPj；
S302, under the coordinate system of the step S2, the frequency domain seismic source FSou is processedjIs set at the j-th detection point.
5. A highly deviated well shot-inspection moving VSP free surface multiple imaging device adopting the method of any one of claims 1-4, characterized in that: the method comprises the following steps:
the data input module is used for inputting a shot-inspection mobile VSP full wavefield common detection point gather, a 2-dimensional geological model and an offset parameter of the highly inclined shaft;
the coordinate establishing module is used for establishing the coordinates of the shot point and the wave detection point in the 2-dimensional geological model according to the input data;
the multiple pre-stack imaging module is used for performing fast Fourier transform on the input common detection point gather to a frequency wave number domain and setting a frequency domain seismic source; extending a frequency domain seismic source from the depth of a demodulator probe upwards to the free surface, extending a frequency wave number domain shot-detection moving VSP common-detector full wave field and a seismic source wave field from the free surface downwards, and realizing free surface multiple-order single-pass wave prestack depth migration imaging of a single-detector point; completing shot-inspection moving VSP free surface multiple prestack depth migration imaging of all inspection points in the same way;
and the superposition imaging module is used for rearranging the obtained multiple prestack depth migration imaging into a common imaging gather, and superposing the common imaging gather to obtain shot-inspection mobile VSP free surface multiple prestack depth migration superposition imaging.
CN202010670672.9A 2020-07-13 2020-07-13 Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well Active CN111781647B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010670672.9A CN111781647B (en) 2020-07-13 2020-07-13 Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010670672.9A CN111781647B (en) 2020-07-13 2020-07-13 Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well

Publications (2)

Publication Number Publication Date
CN111781647A CN111781647A (en) 2020-10-16
CN111781647B true CN111781647B (en) 2022-05-20

Family

ID=72767595

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010670672.9A Active CN111781647B (en) 2020-07-13 2020-07-13 Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well

Country Status (1)

Country Link
CN (1) CN111781647B (en)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
NO883738D0 (en) * 1987-11-16 1988-08-19 Western Atlas Int Inc PROCEDURES FOR VERTICAL SEISMIC PROFILING (VSP).
US4870580A (en) * 1983-12-30 1989-09-26 Schlumberger Technology Corporation Compressional/shear wave separation in vertical seismic profiling
US4894809A (en) * 1985-05-23 1990-01-16 Mobil Oil Corporation Method for bin, moveout correction and stack of offset vertical seismic profile data in media with dip
US5696735A (en) * 1994-10-19 1997-12-09 Exxon Production Research Company Seismic migration using offset checkshot data
CN103454680A (en) * 2013-08-27 2013-12-18 中国石油集团川庆钻探工程有限公司地球物理勘探公司 Method for calculating vertical coverage times of Walk-away VSP observing system
CN104853822A (en) * 2014-09-19 2015-08-19 杨顺伟 Method for evaluating shale gas reservoir and searching sweet spot region
CN105911587A (en) * 2016-04-22 2016-08-31 中国地质大学（北京） Two-way wave pre-stack depth migration method through one-way wave operator
CN109991662A (en) * 2019-05-15 2019-07-09 中油奥博(成都)科技有限公司 Apparatus and method for measuring and calculating two-dimensional or three-dimensional elastic parameters of shallow formation
CN209946406U (en) * 2019-05-15 2020-01-14 中油奥博(成都)科技有限公司 Device for measuring and calculating two-dimensional or three-dimensional elastic parameters of shallow stratum

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10267937B2 (en) * 2014-04-17 2019-04-23 Saudi Arabian Oil Company Generating subterranean imaging data based on vertical seismic profile data and ocean bottom sensor data

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4870580A (en) * 1983-12-30 1989-09-26 Schlumberger Technology Corporation Compressional/shear wave separation in vertical seismic profiling
US4894809A (en) * 1985-05-23 1990-01-16 Mobil Oil Corporation Method for bin, moveout correction and stack of offset vertical seismic profile data in media with dip
NO883738D0 (en) * 1987-11-16 1988-08-19 Western Atlas Int Inc PROCEDURES FOR VERTICAL SEISMIC PROFILING (VSP).
US5696735A (en) * 1994-10-19 1997-12-09 Exxon Production Research Company Seismic migration using offset checkshot data
CN103454680A (en) * 2013-08-27 2013-12-18 中国石油集团川庆钻探工程有限公司地球物理勘探公司 Method for calculating vertical coverage times of Walk-away VSP observing system
CN104853822A (en) * 2014-09-19 2015-08-19 杨顺伟 Method for evaluating shale gas reservoir and searching sweet spot region
CN105911587A (en) * 2016-04-22 2016-08-31 中国地质大学（北京） Two-way wave pre-stack depth migration method through one-way wave operator
CN109991662A (en) * 2019-05-15 2019-07-09 中油奥博(成都)科技有限公司 Apparatus and method for measuring and calculating two-dimensional or three-dimensional elastic parameters of shallow formation
CN209946406U (en) * 2019-05-15 2020-01-14 中油奥博(成都)科技有限公司 Device for measuring and calculating two-dimensional or three-dimensional elastic parameters of shallow stratum

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Walkaway VSP处理技术及应用;臧胜涛等;《石油地球物理勘探》;20171230;全文 *
Walkaway VSP多次波成像技术研究;吴世萍等;《石油物探》;20110325(第02期);全文 *
Walkaway VSP深度域叠前偏移成像;王欣等;《石油地球物理勘探》;20160815(第04期);第745-750页 *
Walkaway VSP自由表面多次波叠前深度偏移成像;李建国等;《石油地球物理勘探》;20181230;第77-82页 *

Also Published As

Publication number Publication date
CN111781647A (en) 2020-10-16

Similar Documents

Publication Publication Date Title
CN101907727B (en) Multi-component converted wave static correction method by using surface waves
CN102053261B (en) Method for processing seismic data
CN103713323B (en) Omnibearing aeolotropy amplitude-preservation imaging and gather extracting method
CN107526101A (en) A kind of collection for obtaining earthquake reflected wave and processing method
CN102540250A (en) Azimuth fidelity angle domain imaging-based fractured oil and gas reservoir seismic exploration method
CN106526678B (en) A kind of wave field separation method and device of reflected acoustic wave well logging
CN103454681B (en) Evaluate the method and apparatus of 3 D seismic observation system imaging effect
CN113640881B (en) Multi-offset two-dimensional lateral high-resolution transient surface wave detection method
Brodic et al. Three-component seismic land streamer study of an esker architecture through S-and surface-wave imaging
CN117826248B (en) Surface wave dispersion extraction method based on multi-scale observation background noise bunching
CN102053260B (en) Method for acquiring azimuth velocity of primary wave and method for processing earthquake data
CN111751876B (en) A method and device for variable offset VSP converted shear wave one-way wave prestack depth migration
CN105093318B (en) A kind of adaptive wave equation wave field extrapolation static correcting method
CN111781647B (en) Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well
CN107340537A (en) A kind of method of P-SV converted waves prestack reverse-time depth migration
Schapper et al. Anisotropic velocities and offset vector tile prestack-migration processing of the Durham Ranch 3D, Northwest Colorado
Lei et al. The application of ambient noise and reflection seismic exploration in an urban active fault survey
Panea et al. Analysis of crooked-line 2D seismic reflection data recorded in areas with complex surface and subsurface conditions
Almholt et al. High resolution 2D reflection seismic land streamer survey for groundwater mapping: Case study from south east Denmark
WO2020180857A1 (en) Seismic surveys using two-way virtual source redatuming
Wang et al. Research on the anisotropy of gas hydrate reservoirs in South China Sea
CN103576188B (en) A kind of seismic source location method of release rate error impact
Li et al. High-density analysis of surface wave profile imaging based on a multiple coverage common-midpoint signal-couple array
CN105242308A (en) Superposition method and apparatus of wide-azimuth earthquake data
O'Connell et al. Bullwiukle: A unique 3-D experiment

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