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
shot
depth
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

Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well
Technical Field
The invention relates to seismic data migration imaging in geophysical exploration, in particular to a free surface multiple imaging method and device for shot-picking mobile VSP (vertical seismic profiling) of a deviated well.
Background
The shot detection moving VSP is an observation system designed for a large inclined well, namely a shot line is arranged on the ground surface and positioned right above a detector in the well, and the shot line moves along with the detector. The advantage of shot-picking to move the VSP is that the coverage is more uniform. The observation system is mainly used for guiding target entering in the inclined shaft drilling process. In general, the free surface multiples are attenuated as interference waves. Compared with reflection, the free surface multiple has the characteristics of easy identification, wide illumination, small reflection angle and long propagation path.
Wu Shineh et al studied Walkaway VSP multiple imaging technology research, which uses seismic coherence imaging for Walkaway VSP multiple imaging. Li Jian et al studied Walkaway VSP free surface multiple prestack depth migration imaging, and the literature adopted single-pass wave continuation to realize free surface multiple imaging of variable offset VSP; however, since research on free surface multiple imaging of shot-scan VSP is not yet mature, there are many inconveniences in the fields of address research, drill guidance, and the like.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provides a free surface multiple imaging method and device for shot-inspection mobile VSP of a highly deviated well, which effectively utilize rich information carried by the free surface multiple, realize pre-stack depth migration imaging of the free surface multiple of the shot-inspection mobile VSP and provide an accurate data base for geological research and well drilling guidance.
The purpose of the invention is realized by the following technical scheme: a method for imaging a free surface multiple of a shot-picking mobile VSP (vertical seismic profiling) of a deviated well 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;
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.
Further, 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.
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 GDA0002618028090000021
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 GDA0002618028090000022
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 depth of the jth demodulator probe, WMX, WMY are the geodetic coordinates of the wellhead, MWMX are the coordinates of the wellhead model, sign is the sign function.
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.
The step S4 includes:
s401, calculating an upward continuation grid of a seismic source in a frequency domain:
kzmig＝(k-1)·dzmig
Figure GDA0002618028090000031
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 GDA00026180280900000312
Extending a depth step upwards:
Figure GDA0002618028090000032
f＝[fl:df:fh]
Figure GDA0002618028090000033
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 GDA0002618028090000034
Figure GDA0002618028090000035
Figure GDA0002618028090000036
Figure GDA0002618028090000037
Figure GDA0002618028090000038
Figure GDA0002618028090000039
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 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, 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 GDA00026180280900000310
is that
Figure GDA00026180280900000311
Extending a depth step wave field 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 GDA0002618028090000041
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 GDA0002618028090000042
With downward continuation by one depth step:
Figure GDA0002618028090000043
f＝[fl:df:fh]
Figure GDA0002618028090000044
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 GDA0002618028090000045
Figure GDA0002618028090000046
Figure GDA0002618028090000047
Figure GDA0002618028090000048
Figure GDA0002618028090000049
Figure GDA00026180280900000410
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, 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 GDA0002618028090000051
is that
Figure GDA0002618028090000052
Extending a depth step wave field downwards;
s406, seismic source of frequency domain
Figure GDA0002618028090000053
Extend one depth step down:
Figure GDA0002618028090000054
f＝[fl:df:fh]
Figure GDA0002618028090000055
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 GDA0002618028090000056
Figure GDA0002618028090000057
Figure GDA0002618028090000058
Figure GDA0002618028090000059
Figure GDA00026180280900000510
Figure GDA00026180280900000511
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 GDA00026180280900000512
is that
Figure GDA00026180280900000513
Extending a depth step wave field downwards;
s407, extracting an imaging value according to the relevant imaging condition:
Figure GDA00026180280900000514
wherein,
Figure GDA00026180280900000515
is that
Figure GDA00026180280900000516
The wavefield is extended down by one depth step,
Figure GDA00026180280900000517
is that
Figure GDA00026180280900000518
The wavefield is extended down by one depth step,
Figure GDA00026180280900000519
is the imaging value of the extracted kzmig + dzmig depth, conj is the complex conjugate function;
s408, circularly executing the steps S404 to S407, and completing shot moving VSP free surface multiple prestack depth migration imaging of the jth wave detection point until zmig2 is the imaging termination depth.
A highly deviated well shot-inspection moving VSP free surface multiple imaging device comprises:
the data input module is used for inputting a shot-inspection moving VSP full wavefield common-inspection point gather, a 2-dimensional geological model and offset parameters of the deviated well;
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.
The invention has the beneficial effects that: the method inputs a shot-inspection mobile VSP common-inspection wave point full wave field and an imaging speed model, a seismic source wave field of a given frequency domain extends upwards from a detection point to a free surface, then the common-inspection wave point full wave field and the seismic source wave field extend downwards from the free surface, and an imaging value is extracted by adopting related imaging conditions, so that the pre-stack depth migration imaging of the shot-inspection mobile VSP free surface multiple waves is realized.
Drawings
FIG. 1 is a flow chart of a method of the present invention;
FIG. 2 is a schematic diagram of a gather of common detection points of a full wavefield of a shot-detection moving VSP of a deviated well in an embodiment;
FIG. 3 is a schematic view of the distribution of the wave detection points of the shot detection moving VSP of the deviated well in the embodiment;
FIG. 4 is a schematic diagram of shot point distribution of a shot inspection mobile VSP of a deviated well in the embodiment;
FIG. 5 is a schematic diagram of an imaging velocity model in an embodiment;
FIG. 6 is a diagram of shot-migrated VSP free surface multiple prestack depth migration imaging in an example.
FIG. 7 is a schematic diagram of prestack depth migration superposition imaging of shot-geophone-receiver (VSP) moving VSP free surface multiples in an embodiment;
fig. 8 is a schematic block diagram of the apparatus of the present invention.
Detailed Description
The technical solutions of the present invention are further described in detail below with reference to the accompanying drawings, but the scope of the present invention is not limited to the following.
As shown in fig. 1, a method for shot-geophone moving VSP free surface multiple imaging in a deviated well comprises the following steps:
s1, inputting a shot-inspection mobile VSP full wavefield common-inspection point gather of a deviated well, a 2-dimensional geological model and offset parameters:
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.
In the embodiment of the application, the input steep-hole shot-detection moving VSP full wavefield common-detection-point gather is shown in FIG. 2, wherein the abscissa in FIG. 2 is the track number and the ordinate is the time (unit: millisecond); the distribution of the wave detection points of the input shot-inspection moving VSP of the deviated well is shown in figure 3, wherein a triangle in figure 3 is the position of the wave detection point, the abscissa is the length (unit: meter), and the ordinate is the depth (unit: meter); shot point distribution of the input shot-picking mobile VSP of the deviated well is shown in FIG. 4, wherein in FIG. 4, a dotted line is a shot line position; the abscissa is the length (unit: meter); the ordinate represents the number of shots; the input imaging velocity model is shown in fig. 5, the abscissa is the length (unit: meter); the ordinate is the depth (unit: meter).
S2, establishing coordinates of the shot point and the wave detection point in the 2-dimensional geological model according to input data;
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 GDA0002618028090000071
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 GDA0002618028090000072
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 depth of the jth demodulator probe, WMX, WMY are the geodetic coordinates of the wellhead, MWMX are the coordinates of the wellhead model, sign is the sign function.
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:
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.
S4, extending a frequency domain seismic source upwards from the depth of a demodulator probe to a free surface, extending a frequency wave number domain shot detection mobile VSP common demodulator probe full wave field and a seismic source wave field downwards from the free surface, wherein the extension operator is an omega-x single-pass operator, extracting an imaging value under related imaging conditions, and realizing free surface multiple single-pass wave prestack depth migration imaging of a single demodulator probe:
s401, calculating an upward continuation grid of a seismic source in a frequency domain:
kzmig＝(k-1)·dzmig
Figure GDA0002618028090000081
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 GDA0002618028090000082
Extending a depth step upwards:
Figure GDA0002618028090000083
f＝[fl:df:fh]
Figure GDA0002618028090000084
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 GDA0002618028090000085
Figure GDA0002618028090000086
Figure GDA0002618028090000091
Figure GDA0002618028090000092
Figure GDA0002618028090000093
Figure GDA0002618028090000094
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 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, 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 GDA0002618028090000095
is that
Figure GDA0002618028090000096
Extending a depth step wave field 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 GDA0002618028090000097
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 GDA0002618028090000098
With downward continuation by one depth step:
Figure GDA0002618028090000099
f＝[fl:df:fh]
Figure GDA00026180280900000910
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 GDA0002618028090000101
Figure GDA0002618028090000102
Figure GDA0002618028090000103
Figure GDA0002618028090000104
Figure GDA0002618028090000105
Figure GDA0002618028090000106
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, 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 GDA0002618028090000107
is that
Figure GDA0002618028090000108
Extending a depth step wave field downwards;
s406, seismic source of frequency domain
Figure GDA0002618028090000109
Extend one depth step down:
Figure GDA00026180280900001010
f＝[fl:df:fh]
Figure GDA00026180280900001011
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 GDA00026180280900001012
Figure GDA00026180280900001013
Figure GDA00026180280900001014
Figure GDA00026180280900001015
Figure GDA0002618028090000111
Figure GDA0002618028090000112
where dip is the maximum tilt angle, fl is the minimum frequency, fh is the maximum frequency, and f ═ fl:df:fh]Is the frequency vector, V is the 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 GDA0002618028090000113
is that
Figure GDA0002618028090000114
Extending a depth step wave field downwards;
s407, extracting an imaging value according to the relevant imaging condition:
Figure GDA0002618028090000115
wherein,
Figure GDA0002618028090000116
is that
Figure GDA0002618028090000117
The wavefield is extended down by one depth step,
Figure GDA0002618028090000118
is that
Figure GDA0002618028090000119
The wavefield is extended down by one depth step,
Figure GDA00026180280900001110
is the imaging value of the extracted kzmig + dzmig depth, conj is the complex conjugate function;
s408, circularly executing the steps S404 to S407, and completing shot moving VSP free surface multiple prestack depth migration imaging of the jth wave detection point until zmig2 is the imaging termination depth.
S5, circularly executing the step S4 to finish shot-inspection moving VSP free surface multiple prestack depth migration imaging of all inspection points;
in the embodiment of the application, shot detection moving VSP free surface multiple prestack depth migration imaging is shown in FIG. 6, which corresponds to FIG. 2, and the abscissa is the track number; the ordinate is the depth (unit: meter);
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.
In the embodiment of the application, shot detection moving VSP free surface multiple pre-stack depth migration superposition imaging is shown in FIG. 7, and the abscissa is the length (unit: meter); the ordinate is the depth (unit: m).
As shown in fig. 8, an oblique borehole shot-inspection moving VSP free-surface multiple imaging apparatus includes:
the data input module is used for inputting a shot-inspection moving VSP full wavefield common-inspection point gather, a 2-dimensional geological model and offset parameters of the deviated well;
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.
In conclusion, the invention inputs a shot detection moving VSP common detection wave point full wave field and imaging speed model, a given frequency domain seismic source wave field extends upwards from a detection point to a free surface, then the common detection wave point full wave field and the seismic source wave field extend downwards from the free surface, and imaging values are extracted by adopting related imaging conditions, so that the shot detection moving VSP free surface multiple wave prestack depth migration imaging is realized.
The foregoing is a preferred embodiment of the present invention, it is to be understood that the invention is not limited to the form disclosed herein, but is not to be construed as excluding other embodiments, and is capable of other combinations, modifications, and environments and is capable of changes within the scope of the inventive concept as expressed herein, commensurate with the above teachings, or the skill or knowledge of the relevant art. And that modifications and variations may be effected by those skilled in the art without departing from the spirit and scope of the invention as defined by the appended claims.

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 中油奥博(成都)科技有限公司 Shallow stratum two dimension or the device and method of three dimensional elasticity parameter measurement and calculating
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 中油奥博(成都)科技有限公司 Shallow stratum two dimension or the device and method of three dimensional elasticity parameter measurement and calculating
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
Satoh et al. S-wave velocity structure of the Taichung basin, Taiwan, estimated from array and single-station records of microtremors
He et al. Petroleum electromagnetic prospecting advances and case studies in China
CN102053261A (en) Method for processing seismic data
CN103713323A (en) Omnibearing aeolotropy amplitude-preservation imaging and gather extracting method
CN109765615A (en) Stratum quality factor inversion method and device
EP2593815B1 (en) Method for accentuating specular and non-specular seismic events from within shallow subsurface rock formations
Wang et al. Physical modelling studies of 3-DP-wave seismic for fracture detection
CN105093318B (en) A kind of adaptive wave equation wave field extrapolation static correcting method
CN104793237A (en) Method and device for acquiring broadband controllable seismic source scanning signal
Zhenwu et al. Practices and expectation of high-density seismic exploration technology in CNPC
CN111781647B (en) Method and device for imaging free surface multiple of VSP (vertical seismic profiling) in shot-inspection mobile process in steep well
CN111751876B (en) Method and device for shifting prestack depth of converted shear wave single-pass wave of variable offset VSP (vertical seismic profiling)
Schapper et al. Anisotropic velocities and offset vector tile prestack-migration processing of the Durham Ranch 3D, Northwest Colorado
Chen et al. Evidence of Complex Faulting near the Huangcheng‐Shuangta Fault, Gansu, China, from the 11 May 2012 M w 4.8 Sunan Earthquake
Li et al. Fractured reservoir characterization using VSP and surface data: a case study from the Paris basin
Panea et al. Analysis of crooked-line 2D seismic reflection data recorded in areas with complex surface and subsurface conditions
Chironi et al. Imaging of intrabasalt and subbasalt structure with full wavefield seismic tomography
Villamizar et al. Seismic imaging of crystalline structures: improving energy focusing and signal alignment with azimuthal binning and 2.5-D full-waveform inversion
WO2020180857A1 (en) Seismic surveys using two-way virtual source redatuming
CN111691876A (en) Method and device for imaging adjacent well by using acoustic logging and storage medium
Almholt et al. High resolution 2D reflection seismic land streamer survey for groundwater mapping: Case study from south east Denmark
Zeng* et al. Wide-azimuth seismic data processing technology for small fault and crack identification
Bao et al. A new technique to support future energy exploration of continental sedimentary basin: BWH full-frequency fidelity and amplitude preserving processing technology
CN103576188B (en) A kind of seismic source location method of release rate error impact
CN117826248B (en) Surface wave dispersion extraction method based on multiscale observation background noise bunching

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