CN113933894A - Virtual source imaging method and device - Google Patents
Virtual source imaging method and device Download PDFInfo
- Publication number
- CN113933894A CN113933894A CN202010670466.8A CN202010670466A CN113933894A CN 113933894 A CN113933894 A CN 113933894A CN 202010670466 A CN202010670466 A CN 202010670466A CN 113933894 A CN113933894 A CN 113933894A
- Authority
- CN
- China
- Prior art keywords
- actual
- wave field
- virtual
- source
- time
- 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.)
- Pending
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 33
- 238000000034 method Methods 0.000 claims abstract description 48
- 230000009466 transformation Effects 0.000 claims description 27
- 230000006978 adaptation Effects 0.000 claims description 10
- 230000001131 transforming effect Effects 0.000 claims description 5
- 238000010586 diagram Methods 0.000 description 18
- 238000010606 normalization Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 4
- 230000002902 bimodal effect Effects 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000005284 excitation Effects 0.000 description 2
- 239000004576 sand Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005422 blasting Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000007935 neutral effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
- G01V2210/322—Trace stacking
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The application provides a method and a device for virtual source imaging, and belongs to the technical field of exploration of earth. The method comprises the following steps: determining an actual down-going wavefield received by the first virtual source and a first actual up-going wavefield received by the second detector; performing time domain convolution on the actual downlink wave field and the cross-correlation channel to obtain a first predicted uplink wave field; carrying out normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the earth surface seismic source; determining a first virtual channel corresponding to the first virtual seismic source according to the weights of the cross-correlation channel and the surface seismic source; a superimposed cross-sectional image of the survey area is determined based on the first virtual source and the at least one second virtual trace. Determining the weights of the surface seismic sources through the actual down-going wave field and the first actual up-going wave field, wherein the actual down-going wave field and the first actual up-going wave field are not influenced by the near-surface environment; therefore, the accuracy of the weight of the earth surface seismic source is improved, and the imaging quality of the superposed section image is improved.
Description
Technical Field
The application relates to the technical field of exploration earth, in particular to a method and a device for virtual source imaging.
Background
In the field of exploration for the earth, near-surface imaging is an important method of exploration for the earth. However, in near-surface imaging, the complex near-surface environment introduces multiples, surface waves, and scattering noise, which interfere with the imaging. Thus, the earth can be explored by means of virtual source imaging, thereby eliminating interference from complex near-surface environments.
In the related art, a plurality of earth surface seismic sources and a plurality of detectors are arranged in an exploration area; the earth surface seismic source is used for exciting seismic waves, and the geophone is used for detecting the seismic waves. For each earth surface seismic source, performing cross correlation on seismic waves excited by the earth surface seismic source and received by any two detectors to obtain a wave field of a virtual earth surface seismic source corresponding to the earth surface seismic source; determining the offset distance between the virtual earth surface seismic source and the earth surface seismic source, wherein the earth surface seismic source with small offset distance has large influence on the virtual earth surface seismic source, and the weight of the earth surface seismic source is large; conversely, the weight of the surface seismic source is small. And superposing the wave fields of the virtual earth surface seismic sources corresponding to the earth surface seismic sources according to the weight of each earth surface seismic source to obtain a superposed section image, thereby finishing the near-earth surface imaging.
However, the complex near-surface environment can seriously distort the path of the seismic wave, and the path of the seismic wave excited by the near-offset ground-surface seismic source to the virtual ground-surface seismic source may be lengthened, so that the influence on the virtual ground-surface seismic source is reduced; the path from the earth surface seismic source with the long offset distance to the virtual earth surface seismic source is possibly shortened, and the influence on the virtual earth surface seismic source is increased; as a result, the weights of the earth's surface seismic source determined from the offset are inaccurate, resulting in poor quality earth's surface imaging from the earth's surface seismic source.
Disclosure of Invention
The embodiment of the application provides a virtual source imaging method and device, which can improve the imaging quality of near-surface imaging. The technical scheme is as follows:
in one aspect, the present application provides a method of virtual source imaging, the method comprising:
determining a first detector and a second detector from a plurality of detectors in a survey area, the first detector being a first virtual source;
for each of a plurality of surface seismic sources in a survey area, determining an actual down-going wavefield of the surface seismic source excited seismic waves received by the first virtual source and a first actual up-going wavefield of the surface seismic source excited seismic waves received by the second detector;
performing time domain cross-correlation on the actual downlink wave field and the first actual uplink wave field to obtain a cross-correlation channel of the first virtual source and the second detector;
performing time domain convolution on the actual downlink wave field and the cross-correlation channel to obtain a first predicted uplink wave field;
performing normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the earth surface seismic source;
determining a first virtual gather corresponding to the first virtual seismic source according to the cross-correlation trace and the weight of the surface seismic source;
determining at least one second virtual gather corresponding to a second virtual source by taking at least one third detector except the first detector in the plurality of detectors in the exploration area as the second virtual source;
determining an overlay profile image of the survey area from the first set of virtual traces and at least one of the second set of virtual traces.
In one possible implementation, the determining, according to the cross-correlation trace and the weight of the surface seismic source, a first virtual gather corresponding to the first virtual seismic source includes:
determining a first virtual gather corresponding to the first virtual seismic source according to the weight of the cross-correlation trace and the weight of the surface seismic source by using a formula I;
the formula I is as follows:
wherein,a weight representing the earth's surface seismic source,represents a time domain cross-correlation; r isA,rBAnd rSRespectively representing the spatial location of the first detector, the spatial location of the second detector and the spatial location of the surface seismic source; t represents a time shift amount; vw(rB|rA(ii) a t) represents the first virtual gather; d (r)A|rS(ii) a t) represents the actual down-going wavefield; u (r)B|rS(ii) a t) represents the first actual upgoing wavefield.
In another possible implementation manner, the performing time-domain convolution on the actual down-going wavefield and the cross-correlation trace to obtain a first predicted up-going wavefield includes:
performing time domain convolution on the actual downlink wave field and the cross-correlation channel, and obtaining a first predicted uplink wave field through a second formula;
the formula II is as follows:
wherein,representing the first predicted upgoing wavefield; denotes time domain convolution; t represents a time shift amount; d (r)A|rS(ii) a t) represents the actual down-going wavefield;representing the cross-correlation trace.
In another possible implementation manner, the performing a normalized cross-correlation process on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the surface seismic source includes:
respectively carrying out Fourier transform on the first predicted upgoing wave field and the first actual upgoing wave field to obtain a second predicted upgoing wave field and a second actual upgoing wave field;
according to the second predicted upgoing wave field and the second actual upgoing wave field, carrying out normalized cross correlation processing on the second predicted upgoing wave field and the second actual upgoing wave field through a third formula to obtain the weight of the earth surface seismic source;
the formula III is as follows:
wherein,representing a weight of the surface seismic source;representing a second predicted upgoing wavefield; u (r)B|rS(ii) a f) Representing the second actual upgoing wavefield and f representing the frequency of the fourier transform.
In another possible implementation manner, the performing a normalized cross-correlation process on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the surface seismic source includes:
acquiring time, frequency and wave number of a same-phase axis of the earth surface seismic source before target stacking;
according to the time, performing wavelet transformation on the first predicted uplink wave field to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field, and performing wavelet transformation on the first actual uplink wave field to obtain a second wavelet coefficient corresponding to the first actual uplink wave field;
converting the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data body and a second time-frequency space data body;
and carrying out normalized cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body to obtain the weight of the earth surface seismic source.
In another possible implementation manner, the performing wavelet transform on the first predicted upgoing wave field according to the time to obtain a first wavelet coefficient corresponding to the first predicted upgoing wave field includes:
according to the time, performing wavelet transformation on the first predicted uplink wave field through a fourth formula to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field;
the formula four is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing the first wavelet coefficients; is thatA seismic representation of the first predicted upgoing wavefieldRoad phi*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
In another possible implementation manner, the performing wavelet transform on the first actual upgoing wave field according to the time to obtain a second wavelet coefficient corresponding to the first predicted upgoing wave field includes:
according to the time, performing wavelet transformation on the first actual uplink wave field through a formula five to obtain a second wavelet coefficient corresponding to the first actual uplink wave field;
the formula five is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing the second wavelet coefficients; is U (t) a seismic trace, φ, representing the first actual upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
In another possible implementation manner, the performing normalized cross-correlation processing on the first time-frequency space data volume and the second time-frequency space data volume to obtain the weight of the surface seismic source includes:
performing normalized cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body through a sixth formula to obtain the weight of the earth surface seismic source;
formula six:
wherein,representing a weight of the surface seismic source;representing a time-frequency wavenumber domain wavefield;representing a first time-frequency spatial data volume;representing a second time-frequency space data volume; f represents the frequency; τ represents a local time; k is a radical ofrRepresenting the wavenumber.
In another possible implementation manner, the method according to the first virtual channel set and at least one of the second virtual channel sets includes:
determining a common center point between the first virtual gather and at least one second virtual gather to obtain at least one common center point;
and overlapping the at least one common central point to obtain an overlapped section image of the exploration area.
In another aspect, the present application provides an apparatus for virtual source imaging, the apparatus comprising:
a first determining module, configured to determine a first detector and a second detector from a plurality of detectors in a survey area, with the first detector as a first virtual source;
a second determination module for determining, for each of a plurality of earth seismic sources in a survey area, an actual down-going wavefield of the earth source-excited seismic waves received by the first virtual source and a first actual up-going wavefield of the earth source-excited seismic waves received by the second detector;
the time domain cross-correlation module is used for performing time domain cross-correlation on the actual downlink wave field and the first actual uplink wave field to obtain a cross-correlation channel of the first virtual source and the second detector;
the time domain convolution module is used for performing time domain convolution on the actual downlink wave field and the cross-correlation channel to obtain a first predicted uplink wave field;
the processing module is used for carrying out normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the earth surface seismic source;
the third determining module is used for determining a first virtual gather corresponding to the first virtual seismic source according to the cross-correlation trace and the weight of the surface seismic source;
a fourth determining module, configured to determine, by using at least one third detector, except for the first detector, of the multiple detectors in the exploration area as a second virtual source, at least one second virtual gather corresponding to the second virtual source;
a fifth determination module to determine an overlay profile image of the survey area based on the first set of virtual traces and at least one of the second set of virtual traces.
In a possible implementation manner, the third determining module is configured to determine, according to the cross-correlation trace and the weight of the surface seismic source, a first virtual gather corresponding to the first virtual seismic source by using the following formula one;
the formula I is as follows:
wherein,a weight representing the earth's surface seismic source,represents a time domain cross-correlation; r isA,rBAnd rSRespectively representing the spatial location of the first detector, the spatial location of the second detector and the spatial location of the surface seismic source; t represents a time shift amount; vw(rB|rA(ii) a t) represents the first virtual gather; d (r)A|rS(ii) a t) represents the actual down-going wavefield; u (r)B|rS(ii) a t) represents the secondAn actual up-going wavefield.
In another possible implementation manner, the time domain convolution module is configured to perform time domain convolution on the actual down-going wavefield and the cross-correlation trace, and obtain a first predicted up-going wavefield according to a second formula;
the formula II is as follows:
wherein,representing the first predicted upgoing wavefield; denotes time domain convolution; t represents a time shift amount; d (r)A|rS(ii) a t) represents the actual down-going wavefield;representing the cross-correlation trace.
In another possible implementation manner, the processing module is configured to perform fourier transform on the first predicted upgoing wave field and the first actual upgoing wave field respectively to obtain a second predicted upgoing wave field and a second actual upgoing wave field;
according to the second predicted upgoing wave field and the second actual upgoing wave field, carrying out normalized cross correlation processing on the second predicted upgoing wave field and the second actual upgoing wave field through a third formula to obtain the weight of the earth surface seismic source;
the formula III is as follows:
wherein,representing a weight of the surface seismic source;representing a second predicted upgoing wavefield; u (r)B|rS(ii) a f) Representing the second actual upgoing wavefield and f representing the frequency of the fourier transform.
In another possible implementation manner, the processing module includes:
the acquisition unit is used for acquiring the time, the frequency and the wave number of the in-phase axis of the earth surface seismic source before target stacking;
the transformation unit is used for performing wavelet transformation on the first predicted uplink wave field according to the time to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field, and performing wavelet transformation on the first actual uplink wave field to obtain a second wavelet coefficient corresponding to the first actual uplink wave field;
a conversion unit, configured to convert the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data volume and a second time-frequency space data volume;
and the normalization unit is used for performing normalized cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body to obtain the weight of the earth surface seismic source.
In another possible implementation manner, the transformation unit is configured to perform wavelet transformation on the first predicted upgoing wave field according to the time by using a following formula four to obtain a first wavelet coefficient corresponding to the first predicted upgoing wave field;
the formula four is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing the first wavelet coefficients; is thatA seismic trace, φ, representing said first predicted upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
In another possible implementation manner, the transforming unit is configured to perform wavelet transformation on the first actual uplink wave field according to the time by using a formula five below to obtain a second wavelet coefficient corresponding to the first actual uplink wave field;
the formula five is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing the second wavelet coefficients; is U (t) a seismic trace, φ, representing the first actual upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
In another possible implementation manner, the normalization unit is configured to perform normalized cross-correlation processing on the first time-frequency space data volume and the second time-frequency space data volume through a sixth following formula to obtain a weight of the earth surface seismic source;
formula six:
wherein,representing a weight of the surface seismic source;representing a time-frequency wavenumber domain wavefield;representing a first time-frequency spatial data volume;representing a second time-frequency space data volume; f represents the frequency; τ represents a local time; k is a radical ofrRepresenting the wavenumber.
In another possible implementation manner, the fifth determining module is configured to determine a common center point between the first virtual gather and at least one of the second virtual gathers, so as to obtain at least one common center point; and overlapping the at least one common central point to obtain an overlapped section image of the exploration area.
The technical scheme provided by the embodiment of the application has the following beneficial effects:
in the embodiment of the application, the weights of the surface seismic source are determined due to the actual downlink wave field of the seismic waves excited by the surface seismic source received by the first virtual source and the first actual uplink wave field of the seismic waves excited by the surface seismic source received by the second detector, and the actual downlink wave field and the first actual uplink wave field are not influenced by the near-surface environment; therefore, the accuracy of determining the weight of the earth surface seismic source through the actual down-going wave field and the first actual up-going wave field is high, and the imaging quality of the superposition section image of the exploration area determined according to the weight of the earth surface seismic source is improved.
Drawings
In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed to be used in the description of the embodiments are briefly introduced below, and it is obvious that the drawings in the following description are only some embodiments of the present application, and it is obvious for those skilled in the art to obtain other drawings based on these drawings without creative efforts.
Fig. 1 is a flow chart of a method of virtual source imaging provided in accordance with an embodiment of the present application;
FIG. 2 is a schematic diagram of a method of virtual source imaging provided in accordance with an embodiment of the present application;
FIG. 3 is a schematic flow chart diagram illustrating a method for virtual source imaging according to an embodiment of the present application;
FIG. 4 is a schematic diagram of a first actual upgoing wavefield and a first predicted upgoing wavefield received by a second detector according to an embodiment of the present application;
FIG. 5 is a schematic diagram of a first time-frequency spatial data volume corresponding to an original first actual upgoing wavefield according to an embodiment of the present application;
FIG. 6 is a diagram illustrating a second time-frequency space data volume corresponding to a first predicted upgoing wavefield, according to an embodiment of the present application;
fig. 7 is a schematic diagram of coherence coefficients of a first time-frequency space data volume and a second time-frequency space data volume according to an embodiment of the present application;
FIG. 8 is a graph of the weighting of seismic sources at a first detector A determined from seismic waves of different batches of surveys according to an embodiment of the present application;
FIG. 9 is a graph of a weighting profile of a seismic source of the earth's surface determined from seismic waves of different batches of surveys at another first geophone B provided in accordance with embodiments of the present application;
FIG. 10 is a superimposed sectional image of a survey area from different surveys provided by embodiments of the present application;
FIG. 11 is a schematic diagram of a comparison of 13 surveys of superimposed cross-sectional images at a common center point within the survey area, according to an embodiment of the present application;
FIG. 12 is a schematic diagram illustrating comparison of NRMS values of a target layer window of a superimposed sectional image according to an embodiment of the present application;
fig. 13 is a block diagram of an apparatus for virtual source imaging according to an embodiment of the present application.
Detailed Description
To make the objects, technical solutions and advantages of the present application more clear, embodiments of the present application will be described in further detail below with reference to the accompanying drawings.
Fig. 1 is a flowchart of a method for virtual source imaging according to an embodiment of the present application. Referring to fig. 1, the method includes:
101. the computer device determines a first detector and a second detector from a plurality of detectors in the survey area, the first detector being a first virtual source.
The exploration area comprises a plurality of detectors, and the detectors can directly receive seismic waves excited by an earth seismic source and also can receive seismic waves reflected by a target layer in the exploration area. Optionally, a plurality of detectors in the exploration area are uniformly distributed; wherein, the distance between any two detectors can be any value between 20m and 50 m; for example, the distance between any two detectors is 30 m. The number of first detectors may be one or more. In the embodiment of the present application, referring to fig. 2, an example of including a plurality of first detectors in a survey area is described.
In one possible implementation, the first detector and the second detector are any two different detectors of the plurality of detectors. Accordingly, the step of the computer device determining a first detector and a second detector from the plurality of detectors in the survey area is: the computer device selects a first detector from the plurality of detectors, and sequentially uses the other detectors except the first detector as second detectors.
In another possible implementation, the computer device selects a first detector from the plurality of detectors and determines a second detector based on the position of the first detector and the other detectors. Accordingly, the step of the computer device determining a first detector and a second detector from the plurality of detectors in the survey area is: the computer device selects a first detector from the plurality of detectors, determines first coordinate information of the first detector and coordinate information of other detectors other than the first detector, and sequentially selects a second detector closest to the first detector from the other detectors other than the first detector.
For example, the number of the plurality of detectors in the survey area is 5, 1 first detector is determined from the 5 detectors, and a second detector closest to the first detector is sequentially selected from the other 4 detectors except the first detector by the first detector.
It should be noted that the number of second detectors may be one or more less than the number of multiple detectors in the survey area. Wherein the higher the number of second receivers, the higher the accuracy of the survey.
In one possible implementation, the computer device randomly selects a first detector from the plurality of detectors.
In another possible implementation, the computer device selects a first detector from the plurality of detectors based on a distance of the detector from the target layer. Accordingly, the computer device determines the first detector by: the computer equipment acquires third coordinate information of the plurality of detectors and fourth coordinate information of the target layer; a first detector closest to the target layer is selected from the plurality of detectors.
In this step, the first detector is used as a first virtual source, the positions of the first virtual source and the first detector are the same, and the received seismic waves of the earth surface seismic source are the same.
102. The computer device determines, for each of a plurality of earth surface seismic sources in the survey area, an actual down-going wavefield of earth surface source-excited seismic waves received by the first virtual source and a first actual up-going wavefield of earth surface source-excited seismic waves received by the second detector.
The exploration area comprises a plurality of earth surface seismic sources; the surface seismic source may be a shot point located at the surface. Optionally, a plurality of earth surface seismic sources in the exploration area are uniformly distributed; wherein, the distance between any two earth surface seismic sources can be any value between 5m and 50 m; for example, the distance between any two surface seismic sources is 7.5 m.
In the embodiment of the application, for each geophone, the computer stores seismic waves of a plurality of earth surface seismic sources received by the geophone within a preset time period. Wherein the preset time length can be any value between 50ms and 5000 ms; such as 60ms, 600ms, etc. For each earth source, the geophone receives seismic waves including an actual down-going wavefield of seismic waves and a first actual up-going wavefield of seismic waves.
In one possible implementation, the computer device stores therein a correspondence between geophone coordinates, surface seismic source identification, and seismic waves. The coordinates of the first virtual source are the same as the coordinates of the first detector.
Correspondingly, for each of the plurality of earth surface seismic sources in the survey area, the step of determining the actual down-going wavefield of seismic waves excited by the earth surface seismic source received by the first virtual source and the first actual up-going wavefield of seismic waves excited by the earth surface seismic source received by the second detector is as follows: the computer equipment acquires an earth surface seismic source identifier of an earth surface seismic source and first coordinate information of a first detector, and determines an actual downlink wave field of seismic waves excited by the earth surface seismic source received by a first virtual source from a corresponding relation among detector coordinates, the earth surface seismic source identifier and the seismic waves stored in the computer equipment according to the first coordinate information and the earth surface seismic source identifier; and the computer equipment acquires second coordinate information of the second geophone, and determines a first actual uplink wave field of seismic waves excited by the earth surface seismic source received by the second geophone from the corresponding relation among the geophone coordinates, the earth surface seismic source identification and the seismic waves stored in the computer equipment according to the second coordinate information and the earth surface seismic source identification.
The surface seismic source identification is used for distinguishing different seismic sources; for example, referring to fig. 3, N surface seismic sources are included within the survey area, identified as shot 1, shot 2, … …, shot k, … …, N, respectively.
In another possible implementation manner, the computer device stores the geophone identifier, the surface seismic source identifier and the corresponding relation with the seismic wave. The first virtual source identification of the first virtual source is the same as the first detector identification of the first detector.
Correspondingly, for each of the plurality of earth surface seismic sources in the survey area, the step of determining the actual down-going wavefield of seismic waves excited by the earth surface seismic source received by the first virtual source and the first actual up-going wavefield of seismic waves excited by the earth surface seismic source received by the second detector is as follows: the method comprises the steps that computer equipment obtains an earth surface seismic source identifier of an earth surface seismic source and a first detector identifier of a first detector, and according to the first detector identifier and the earth surface seismic source identifier, an actual downlink wave field of seismic waves excited by the earth surface seismic source and received by a first virtual source is determined from a corresponding relation among the detector identifiers, the earth surface seismic source identifier and the seismic waves stored in the computer equipment; and the computer equipment acquires a second detector identifier of the second detector, and determines a first actual uplink wave field of seismic waves excited by the earth surface seismic source received by the second detector from the corresponding relation among the detector coordinates, the earth surface seismic source identifier and the seismic waves stored in the computer equipment according to the second detector identifier and the earth surface seismic source identifier. Wherein the detector identification is used for distinguishing different detectors; for example, the detector identification may be A, B, C, etc.
103. And the computer equipment performs time domain cross-correlation on the actual downlink wave field and the first actual uplink wave field to obtain a cross-correlation channel of the first virtual source and the second detector.
In a possible implementation manner, the computer device performs time domain cross-correlation on the actual down-going wave field and the first actual up-going wave field through a first prediction formula to obtain a cross-correlation trace of the first virtual source and the second detector. A first predictive formula is stored within the computer device.
Optionally, the computer device stores a first prediction formula:
wherein t represents the amount of time shift; g (r) represents a weight operator;represents a time domain cross-correlation; r isA,rBAnd rSRespectively representing the spatial position of the first detector, the spatial position of the second detector and the spatial position of the earth surface seismic source; t represents a time shift amount;a cross-correlation trace representing the first virtual source and the second detector; d (r)A|rS(ii) a t) represents the actual down-going wavefield; u (r)B|rS(ii) a t) represents the first actual upgoing wavefield.
It should be noted that the weight operator in the first prediction formula is related to the offset between the earth's surface source and the first detector. Optionally, the weight operator is a gaussian function centered at the first detector.
Optionally, the gaussian function corresponding to the weight operator is stored in the computer device and is:
wherein g (r) represents a weight operator; r represents the radius of the area where the plurality of surface seismic sources are located, and R represents the offset distance between the surface seismic sources and the first detector. Where g (r) typically ranges between 0 and 1, with a maximum of 1 directly above the virtual source and 0 at the maximum offset from the surface seismic source. For example, referring to fig. 2, g (r) is a bell curve.
104. And the computer equipment performs time domain convolution on the actual downlink wave field and the cross-correlation trace to obtain a first predicted uplink wave field.
In one possible implementation, a first prediction function is stored within the computer device. The first prediction function is the following equation two. Correspondingly, the steps can be as follows: and the computer equipment performs time domain convolution on the first actual up-going wave field and the cross-correlation trace to obtain a first predicted up-going wave field through the following formula II.
The formula II is as follows:
wherein,representing a first predicted upgoing wavefield; denotes time domain convolution; t represents a time shift amount; d (r)A|rS(ii) a t) represents the actual down-going wavefield;representing the cross-correlation trace.
105. And the computer equipment performs normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the earth surface seismic source.
In one possible implementation, the method includes the following steps: the computer equipment respectively performs Fourier transform on the first predicted upgoing wave field and the first actual upgoing wave field to obtain a second predicted upgoing wave field and a second actual upgoing wave field; according to the second predicted upgoing wave field and the second actual upgoing wave field, performing normalized cross-correlation processing on the second predicted upgoing wave field and the second actual upgoing wave field through a third formula to obtain the weight of the earth surface seismic source;
the formula III is as follows:
wherein,representing the weight of the earth surface seismic source;representing a second predicted upgoing wavefield; u (r)B|rS(ii) a f) Representing the second actual upgoing wavefield and f representing the frequency of the fourier transform.
The fourier transform equation is:
wherein f represents the frequency of the fourier transform; t represents the amount of time shift.
It should be noted that the computer device converts the first predicted upgoing wavefield from a time domain function to a second predicted upgoing wavefield from a frequency domain function by a fourier transform formula; the first actual upgoing wavefield is converted from a time-domain function to a second actual upgoing wavefield of a frequency-domain function.
It should be noted that the weight of the earth seismic source determined by the above formula three is a two-dimensional function, including a time domain dimension and a frequency domain dimension. The more the dimensionality of the weight function is, the higher the accuracy of obtaining the weight of the earth surface seismic source through the weight function is.
In another possible implementation, to increase the dimensionality of the initial weight function, the computer device increases the wavenumber domain dimension in the weights of the surface seismic sources. Correspondingly, the step of carrying out normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field by the computer equipment to obtain the weight of the earth surface seismic source comprises the following steps: the computer equipment acquires the time, frequency and wave number of a ground surface seismic source in a same phase axis before target stacking; according to time, performing wavelet transformation on the first predicted upgoing wave field to obtain a first wavelet coefficient corresponding to the first predicted upgoing wave field, and performing wavelet transformation on the first actual upgoing wave field to obtain a second wavelet coefficient corresponding to the first actual upgoing wave field; converting the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data body and a second time-frequency space data body; and carrying out normalized cross correlation on the first time-frequency space data body and the second time-frequency space data body to obtain the weight of the earth surface seismic source.
In a possible implementation manner, the wavelet transforming, by the computer device, the first predicted upgoing wave field according to time to obtain a first wavelet coefficient corresponding to the first predicted upgoing wave field includes: according to time, the computer equipment performs wavelet transformation on the first predicted uplink wave field through the following formula IV to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field;
the formula four is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing a first wavelet coefficient; is thatSeismic traces, phi, representing a first predicted upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
In a possible implementation manner, the wavelet transforming, by the computer device, the first actual upgoing wave field according to time to obtain a second wavelet coefficient corresponding to the first actual upgoing wave field includes: according to time, the computer equipment performs wavelet transformation on the first actual uplink wave field through the following formula V to obtain a second wavelet coefficient corresponding to the first actual uplink wave field;
the formula five is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing second wavelet coefficients; is U (t) a seismic trace, φ, representing a first actual upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
The computer equipment converts the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data body and a second time-frequency space data body, and the steps are as follows: and the computer equipment performs Fourier transform on the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data body and a second time-frequency space data body.
In a possible implementation manner, the step of performing, by the computer device, normalized cross-correlation processing on the first time-frequency space data volume and the second time-frequency space data volume to obtain the weight of the earth surface seismic source is as follows: the computer equipment performs normalized cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body through a sixth formula to obtain the weight of the earth surface seismic source;
formula six:
wherein,representing the weight of the earth surface seismic source;representing a time-frequency wavenumber domain wavefield;representing a first time-frequency spatial data volume;representing a second time-frequency space data volume; f represents a frequency; τ represents a local time; k is a radical ofrRepresenting the wave number.
It should be noted that, with reference to fig. 3, a surface seismic source is taken as a shot point for explanation. After determining the weight of K shots, the computer device determines whether K shots are the last shots. If K shots are the last shots, execution continues with step 106. If K shot is not the last shot, a new shot K +1 is determined, and steps 102-105 continue for the new shot K +1 until the weights of all shots are determined, and step 106 continues.
106. And the computer equipment determines a first virtual gather corresponding to the first virtual seismic source according to the cross-correlation trace and the weight of the earth surface seismic source.
And the computer equipment determines a first virtual gather corresponding to the first virtual seismic source according to the weights of the cross-correlation trace and the earth surface seismic source by the following formula I.
The formula I is as follows:
wherein t represents the amount of time shift; vw(rB|rA(ii) a t) represents a first virtual gather;the weight of the earth's seismic source is represented,represents a time domain cross-correlation; r isA,rBAnd rSRespectively representing the spatial position of the first detector, the spatial position of the second detector and the spatial position of the earth surface seismic source; d (r)A|rS(ii) a t) represents the actual down-going wavefield; u (r)B|rS(ii) a t) represents the first actual upgoing wavefield.
It should be noted that, with continued reference to fig. 3, after determining the first virtual gather corresponding to the first virtual source, the computer device determines whether the first virtual source is the last virtual source. If the first virtual source shot is the last shot, step 108 is performed directly. If the first virtual source is not the last shot, step 107 is performed to determine a new second virtual source, and for the second virtual source, steps 101 through 106 are continued until all second virtual sources determine corresponding second virtual gathers, and step 108 is not continued.
107. The computer device takes at least one third detector except the first detector in the plurality of detectors in the exploration area as at least one second virtual source, and determines at least one second virtual gather corresponding to the at least one second virtual source.
The number of the third detectors is one or more less than the number of the plurality of detectors in the survey area. In this step, the second virtual source and the third detector are located at the same position, and the received seismic waves of the earth surface seismic source are the same.
In one possible implementation, the number of third detectors is one less than the number of the plurality of detectors in the survey area.
In another possible implementation, the computer device may also select the third detector from among the plurality of detectors at intervals in order to reduce the number of third detectors and improve imaging efficiency. Optionally, the computer device selects a third detector from the plurality of detectors at intervals of a first predetermined number of detectors. Wherein the first preset number may be any number between 1 and 10; e.g., 1, 3, 5, etc.
Wherein the computer device determines at least one second virtual gather corresponding to the at least one second virtual source. For each second virtual source, the step of determining, by the computer device, the second virtual gather corresponding to the second virtual source is the same as the step of determining, by the computer device, the first virtual gather corresponding to the first virtual source in step 106, and details are not repeated here.
It should be noted that the number of third detectors may be one or more less than the number of the plurality of detectors in the survey area. The more the third detectors are, the more the second virtual sources corresponding to the third detectors are, and the higher the exploration precision is.
108. The computer device determines a superimposed profile image of the survey area from the first set of virtual traces and the at least one second set of virtual traces.
In one possible implementation, the superimposed cross-sectional image is constructed from a plurality of common neutral points within the survey area. Correspondingly, the method comprises the following steps: the computer device determines at least one common center point of the first virtual gather and the at least one second virtual gather, and superimposes the at least one common center point to obtain a superimposed profile image of the survey area.
Next, the imaging quality of the virtual source imaging method in the present application is actually tested by time-lapse seismic data in a desert environment. Wherein the step of testing the imaging quality comprises the following steps (1) to (5):
(1) the computer equipment acquires seismic waves of a plurality of detectors in the desert area.
The desert area was subjected to 13 repeated two-dimensional surveys within 19 months, the first 3 (S1-S6) completed within 3 months, and the last 7 surveys within a week after 17 months of interruption (S7-S13). The 13 surveys all used Mertz 26 model controlled surface seismic sources manufactured by MERTZ, USA. Wherein, the precision of the Mertz 26 type controllable earth surface seismic source repeated excitation reaches 1 m.
The desert area contains 80 detectors, the average interval of the 80 detectors is 30m, and the depth of the detectors is 50 m; 80 detectors form a two-dimensional side line. The desert area comprises 9 gun lines. Wherein, the shot points of 9 gun lines are excited by adopting a dense three-dimensional surface, and the gun distance between the shot points is 7.5m multiplied by 7.5 m. It should be noted that the shot point uses dense three-dimensional surface excitation to effectively remove linear and scattered noise. Wherein, the offset distance between the detector and the earth surface seismic source is 0 to 2400 m.
The target layer is located about 2000 meters underground, and a plurality of sets of strata with large speed differences are covered on the target layer. The near-surface is a non-uniformly distributed sand layer, and the thickness of the sand layer is any value between several meters and dozens of meters. Below the near surface is a simple layered structure with the angles of inclination of the reflecting interfaces of the layered structure all less than 5 deg.
(2) And the computer equipment performs noise reduction processing on the seismic waves.
The computer equipment directly superposes 9 gun lines and suppresses blasting noise, random noise and scattered energy with strong amplitude. Wherein, a fifth gun line positioned in the center is obtained and is overlapped with a two-dimensional side line consisting of 80 detectors. And (4) integrating shot points into a common channel set by computer equipment, and carrying out f-k filtering to remove linear noise.
(3) The computer equipment performs imaging through the virtual source imaging method in the application.
301. The first detector is the same as the plurality of detectors in the desert area. The computer equipment sequentially determines first geophones from a plurality of geophones in the desert area, determines the angle of seismic waves received by the first geophones for each first geophone, and determines a second geophone through which the path of the seismic waves passes according to the angle of the seismic waves.
302. The computer device determines an actual down-going wavefield of the seismic waves received by the first geophone and a first actual up-going wavefield of the seismic waves received by the second geophone. Wherein the actual downlink wave field of the seismic waves comprises the first arrival waves and the direct arrival waves received by the first detector within 60 ms. The first actual upgoing wavefield of seismic waves includes first arrival waves and direct arrival waves received by the second geophone within 60 ms.
303. The computer equipment performs time domain cross correlation on the actual downlink wave field and the first actual uplink wave field to obtain a cross correlation channel of the earth surface seismic source at the first detector; and performing time domain convolution on the actual downlink wave field and the cross-correlation channel to obtain a first predicted uplink wave field. A schematic of the first actual upgoing wavefield and the first predicted upgoing wavefield received by the second detector is shown in figure 4. Where the graph a) in figure 4 is a schematic representation of the first actual upgoing wavefield received by the second detector and the graph b) in figure 4 is a schematic representation of the first predicted upgoing wavefield received by the second detector.
304. The computer equipment performs wavelet transformation on the first predicted uplink wave field to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field and performs wavelet transformation on the first actual uplink wave field to obtain a second wavelet coefficient corresponding to the first actual uplink wave field according to time; and converting the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data body and a second time-frequency space data body.
Wherein the original first actual upgoing wave field corresponds to the first time-frequency space data volumeAs shown in fig. 5. Wherein the second time-frequency space data body corresponding to the first prediction upgoing wave fieldAs shown in fig. 6. The coherence coefficients of the first time-frequency space data volume and the second time-frequency space data volume are shown in fig. 7.
305. And the computer equipment performs normalized cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body to obtain the weight of the earth surface seismic source.
It should be noted that the weight of the earth surface seismic source is determined according to the seismic wave of each exploration, and in the 13 exploration processes, the seismic waves received by the same detector in different batches of exploration are different. Correspondingly, the weights of the same earth surface seismic source at the same detector are different in 13 exploration processes.
FIG. 8 is a graph of the weight distribution of the earth's surface seismic source determined by the computer device at the first receiver A from seismic waves of different batches of surveys. Wherein, a) in FIG. 8 is a weight distribution diagram of the earth surface seismic source determined by the first geophone A according to the seismic waves of the 1 st exploration; FIG. 8 b) is a graph of the weight distribution of the earth's source determined by the first geophone A from the seismic waves of the 11 th survey; fig. 8, c) is a graph of the weight distribution of the earth's source determined by the first detector a according to prior art methods. Wherein the first detector A is a virtual seismic source A.
FIG. 9 is a graph of the weight distribution of the earth's surface seismic source determined by the computer device at another first detector B from seismic waves of a different survey batch. Wherein, a) in FIG. 9 is a weight distribution diagram of the earth surface seismic source determined by the first geophone B according to the seismic waves of the 1 st exploration; FIG. 9B) is a graph of the weight distribution of the earth's source determined by the first geophone B from the seismic waves of the 11 th survey; fig. 9, c) is a graph of the weight distribution of the earth's source determined by the first detector B according to prior art methods. Wherein the first detector B is a virtual seismic source B.
306. And the computer equipment determines the superposed section image of the exploration area according to the common center point of the virtual seismic source and the second detector.
(4) And comparing the signal-to-noise ratio of the superposed sectional image obtained by the application with the signal-to-noise ratio of the superposed sectional image obtained by the related technology by the computer equipment.
FIG. 10 is a cross-sectional view of a computer device surveying the survey area in different batches. Wherein, a) in FIG. 10 is a superimposed sectional image of the survey area obtained by the method of the related art for the 1 st survey; FIG. 10 b) is a superimposed sectional image of survey area obtained by the method of the present application for survey 1; FIG. 10 c) is a superimposed sectional image of the survey area obtained by the method of the related art for the 11 th survey; FIG. 10 d) is a superimposed sectional image of the survey area obtained by the method of the present application for the 11 th survey. It should be noted that the signal-to-noise ratio at the arrow of the superimposed sectional image of the exploration area obtained by the method in the related art is low; the signal-to-noise ratio at the arrow of the superimposed sectional image of the exploration area obtained by the method in the application is high.
(5) The computer device compares the repeatability of the overlay cross-sectional image obtained by the present application with that obtained in the related art.
In one possible implementation, the computer device compares the repeatability by overlapping profile images at some concentric point within the survey area.
See, for example, fig. 11; wherein, the a) diagram in fig. 11 is a comparison schematic diagram of the superposed section images obtained at a certain concentric point in the exploration area through a non-virtual source method in 13 exploration. Fig. 11 b) is a diagram showing a comparison of 13 surveys of the superimposed sectional images obtained at a certain concentric point in the survey area by the method of the related art. Fig. 11 c) is a schematic diagram of a comparison of 13 surveys with superimposed sectional images obtained at a concentric point within the survey area by the method of the present application.
It should be noted that the repeatability of the overlapping profile images obtained at a certain concentric point in the exploration area through the non-virtual source method and the method in the related art has the problem of obvious wavelet amplitude and waveform inconsistency in the comparison schematic diagrams of 13 exploration. The repeatability of the superimposed sectional image obtained by the method of the application at a certain concentric point in the exploration area is enhanced in the comparison schematic diagram of 13 exploration.
In another possible implementation, the computer device determines an NRMS (Normalized Root Mean Square) value of a target layer window of the overlay cross-sectional image, and determines a repeatability of the overlay cross-sectional image based on the NRMS value. Wherein, the smaller the NRMS value is, the stronger the reproducibility of the superimposed sectional image is.
Referring to fig. 12, the computer sets NRMS values for the target layer windows of the superimposed sectional images for 13 determinations, respectively. The non-virtual source method determines that the NRMS value of a target layer window of the overlay section image is in bimodal distribution, the peak values respectively appear at 20% and 65%, and the NRMS value is 47%. The traditional virtual source method determines that the NRMS value of the target layer window of the overlay section image is in a bimodal distribution, the peak values occur at about 33% and 52%, and the NRMS value is 42%. The method determines that the NRMS value of a target layer window of the superposed section image is slightly bimodal, and the NRMS value is 37%. The NRMS value of the overlapped sectional image determined by the method is the minimum, so that the repeatability of the overlapped sectional image determined by the method is the strongest.
In the embodiment of the application, the weights of the surface seismic source are determined due to the actual downlink wave field of the seismic waves excited by the surface seismic source received by the first virtual source and the first actual uplink wave field of the seismic waves excited by the surface seismic source received by the second detector, and the actual downlink wave field and the first actual uplink wave field are not influenced by the near-surface environment; therefore, the accuracy of determining the weight of the earth surface seismic source through the actual down-going wave field and the first actual up-going wave field is high, and the imaging quality of the superposition section image of the exploration area determined according to the weight of the earth surface seismic source is improved.
Fig. 13 is a block diagram of an apparatus for virtual source imaging according to an embodiment of the present application. Referring to fig. 13, the apparatus includes:
a first determining module 1301, configured to determine a first detector and a second detector from a plurality of detectors in a survey area, with the first detector as a first virtual source;
a second determining module 1302 for determining, for each of a plurality of earth seismic sources in the survey area, an actual down-going wavefield of seismic waves excited by the earth seismic source received by the first virtual source and a first actual up-going wavefield of seismic waves excited by the earth seismic source received by the second detector;
the time domain cross-correlation module 1303 is used for performing time domain cross-correlation on the actual downlink wave field and the first actual uplink wave field to obtain a cross-correlation channel between the first virtual source and the second detector;
a time domain convolution module 1304, configured to perform time domain convolution on the actual downlink wave field and the cross-correlation trace to obtain a first predicted uplink wave field;
the processing module 1305 is configured to perform normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field to obtain a weight of the earth surface seismic source;
a third determining module 1306, configured to determine, according to the cross-correlation trace and the weight of the earth surface seismic source, a first virtual gather corresponding to the first virtual seismic source;
a fourth determining module 1307, configured to determine, by using at least one third detector, except for the first detector, in the plurality of detectors in the exploration area as a second virtual source, at least one second virtual gather corresponding to the second virtual source;
a fifth determination module 1308 for determining an overlay profile image of the survey area based on the first set of virtual traces and the at least one second set of virtual traces.
In a possible implementation manner, the third determining module 1306 is configured to determine, according to the cross-correlation trace and the weight of the earth surface seismic source, a first virtual gather corresponding to the first virtual seismic source by the following formula one;
the formula I is as follows:
wherein,the weight of the earth's seismic source is represented,represents a time domain cross-correlation; r isA,rBAnd rSRespectively representing the spatial position of the first detector, the spatial position of the second detector and the spatial position of the earth surface seismic source; t represents a time shift amount; vw(rB|rA(ii) a t) represents a first virtual gather; d(rA|rS(ii) a t) represents the actual down-going wavefield; u (r)B|rS(ii) a t) represents the first actual upgoing wavefield.
In another possible implementation manner, the time domain convolution module is configured to perform time domain convolution on the actual downlink wavefield and the cross-correlation trace, and obtain a first predicted uplink wavefield according to a second formula;
the formula II is as follows:
wherein,representing a first predicted upgoing wavefield; denotes time domain convolution; t represents a time shift amount; d (r)A|rS(ii) a t) represents the actual down-going wavefield;representing the cross-correlation trace.
In another possible implementation manner, the processing module 1305 is configured to perform fourier transform on the first predicted upgoing wave field and the first actual upgoing wave field respectively to obtain a second predicted upgoing wave field and a second actual upgoing wave field;
according to the second predicted upgoing wave field and the second actual upgoing wave field, performing normalized cross-correlation processing on the second predicted upgoing wave field and the second actual upgoing wave field through a third formula to obtain the weight of the earth surface seismic source;
the formula III is as follows:
wherein,representing the weight of the earth surface seismic source;representing a second predicted upgoing wavefield; u (r)B|rS(ii) a f) Representing the second actual upgoing wavefield and f representing the frequency of the fourier transform.
In another possible implementation, the processing module 1305 includes: the acquisition unit is used for acquiring the time, the frequency and the wave number of a ground surface seismic source in a same phase axis before target stacking; the transformation unit is used for performing wavelet transformation on the first predicted uplink wave field to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field and performing wavelet transformation on the first actual uplink wave field to obtain a second wavelet coefficient corresponding to the first actual uplink wave field according to time; the conversion unit is used for converting the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data body and a second time-frequency space data body; and the normalization unit is used for performing normalization cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body to obtain the weight of the earth surface seismic source.
In another possible implementation manner, the transformation unit is configured to perform wavelet transformation on the first predicted upgoing wave field according to time by using the following formula four to obtain a first wavelet coefficient corresponding to the first predicted upgoing wave field;
the formula four is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing a first wavelet coefficient; is thatSeismic traces, phi, representing a first predicted upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
In another possible implementation manner, the transformation unit is configured to perform wavelet transformation on the first actual upgoing wave field according to time by using the following formula five to obtain a second wavelet coefficient corresponding to the first actual upgoing wave field;
the formula five is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing second wavelet coefficients; is U (t) a seismic trace, φ, representing a first actual upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
In another possible implementation manner, the normalization unit is configured to perform normalized cross-correlation processing on the first time-frequency space data volume and the second time-frequency space data volume by using a formula six below, so as to obtain a weight of the earth surface seismic source;
formula six:
wherein,representing the weight of the earth surface seismic source;representing a time-frequency wavenumber domain wavefield;representing a first time-frequency spatial data volume;representing a second time-frequency space data volume; f represents a frequency; τ represents a local time; k is a radical ofrRepresenting the wave number.
In another possible implementation manner, the fifth determining module 1208 is configured to determine a common center point between the first virtual channel set and the at least one second virtual channel set, so as to obtain at least one common center point; and overlapping at least one common central point to obtain an overlapped section image of the exploration area.
In the embodiment of the application, the weights of the surface seismic source are determined due to the actual downlink wave field of the seismic waves excited by the surface seismic source received by the first virtual source and the first actual uplink wave field of the seismic waves excited by the surface seismic source received by the second detector, and the actual downlink wave field and the first actual uplink wave field are not influenced by the near-surface environment; therefore, the accuracy of determining the weight of the earth surface seismic source through the actual down-going wave field and the first actual up-going wave field is high, and the imaging quality of the superposition section image of the exploration area determined according to the weight of the earth surface seismic source is improved.
The above description is only exemplary of the present application and should not be taken as limiting, as any modification, equivalent replacement, or improvement made within the spirit and principle of the present application should be included in the protection scope of the present application.
Claims (10)
1. A method of virtual source imaging, the method comprising:
determining a first detector and a second detector from a plurality of detectors in a survey area, the first detector being a first virtual source;
for each of a plurality of surface seismic sources in a survey area, determining an actual down-going wavefield of the surface seismic source excited seismic waves received by the first virtual source and a first actual up-going wavefield of the surface seismic source excited seismic waves received by the second detector;
performing time domain cross-correlation on the actual downlink wave field and the first actual uplink wave field to obtain a cross-correlation channel of the first virtual source and the second detector;
performing time domain convolution on the actual downlink wave field and the cross-correlation channel to obtain a first predicted uplink wave field;
performing normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the earth surface seismic source;
determining a first virtual gather corresponding to the first virtual seismic source according to the cross-correlation trace and the weight of the surface seismic source;
determining at least one second virtual gather corresponding to a second virtual source by taking at least one third detector except the first detector in the plurality of detectors in the exploration area as the second virtual source;
determining an overlay profile image of the survey area from the first set of virtual traces and at least one of the second set of virtual traces.
2. The method of claim 1, wherein determining the first virtual gather corresponding to the first virtual source according to the cross-correlation trace and the weight of the surface source comprises:
determining a first virtual gather corresponding to the first virtual seismic source according to the weight of the cross-correlation trace and the weight of the surface seismic source by using a formula I;
the formula I is as follows:
wherein,a weight representing the earth's surface seismic source,represents a time domain cross-correlation; r isA,rBAnd rSRespectively representing the spatial location of the first detector, the spatial location of the second detector and the spatial location of the surface seismic source; t represents a time shift amount; vw(rB|rA(ii) a t) represents the first virtual gather; d (r)A|rS(ii) a t) represents the actual down-going wavefield; u (r)B|rS(ii) a t) represents the first actual upgoing wavefield.
3. The method of claim 1, wherein time-domain convolving the actual down-going wavefield with the cross-correlation traces to obtain a first predicted up-going wavefield comprises:
performing time domain convolution on the actual downlink wave field and the cross-correlation channel, and obtaining a first predicted uplink wave field through a second formula;
the formula II is as follows:
4. The method of claim 1, wherein the normalized cross-correlation of the first predicted upgoing wavefield and the first actual upgoing wavefield to obtain the weights for the seismic surface source comprises:
respectively carrying out Fourier transform on the first predicted upgoing wave field and the first actual upgoing wave field to obtain a second predicted upgoing wave field and a second actual upgoing wave field;
according to the second predicted upgoing wave field and the second actual upgoing wave field, carrying out normalized cross correlation processing on the second predicted upgoing wave field and the second actual upgoing wave field through a third formula to obtain the weight of the earth surface seismic source;
the formula III is as follows:
5. The method of claim 1, wherein the normalized cross-correlation of the first predicted upgoing wavefield and the first actual upgoing wavefield to obtain the weights for the seismic surface source comprises:
acquiring time, frequency and wave number of a same-phase axis of the earth surface seismic source before target stacking;
according to the time, performing wavelet transformation on the first predicted uplink wave field to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field, and performing wavelet transformation on the first actual uplink wave field to obtain a second wavelet coefficient corresponding to the first actual uplink wave field;
converting the first wavelet coefficient and the second wavelet coefficient according to the frequency and the wave number to obtain a first time-frequency space data body and a second time-frequency space data body;
and carrying out normalized cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body to obtain the weight of the earth surface seismic source.
6. The method of claim 5, wherein said wavelet transforming the first predicted upgoing wavefield based on the time to obtain first wavelet coefficients corresponding to the first predicted upgoing wavefield comprises:
according to the time, performing wavelet transformation on the first predicted uplink wave field through a fourth formula to obtain a first wavelet coefficient corresponding to the first predicted uplink wave field;
the formula four is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing the first wavelet coefficients; is thatA seismic trace, φ, representing said first predicted upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
7. The method of claim 5, wherein said wavelet transforming said first actual upgoing wavefield based on said time to obtain second wavelet coefficients corresponding to said first predicted upgoing wavefield comprises:
according to the time, performing wavelet transformation on the first actual uplink wave field through a formula five to obtain a second wavelet coefficient corresponding to the first actual uplink wave field;
the formula five is as follows:
wherein t represents the amount of time shift; alpha represents a wavelet scale; τ represents a local time;representing the second wavelet coefficients; is U (t) a seismic trace, φ, representing the first actual upgoing wavefield*Representing the parent wavelet that satisfies the adaptation condition and has a mean value of zero.
8. The method of claim 5, wherein the performing normalized cross-correlation on the first time-frequency space data volume and the second time-frequency space data volume to obtain the weights of the earth seismic source comprises:
performing normalized cross-correlation processing on the first time-frequency space data body and the second time-frequency space data body through a sixth formula to obtain the weight of the earth surface seismic source;
formula six:
wherein,representing a weight of the surface seismic source;representing a time-frequency wavenumber domain wavefield;representing a first time-frequency spatial data volume;representing a second time-frequency space data volume; f represents the frequency; τ represents a local time; k is a radical ofrRepresenting the wavenumber.
9. The method of claim 1, wherein said determining from said first virtual gather and at least one of said second virtual gathers comprises:
determining a common center point between the first virtual gather and at least one second virtual gather to obtain at least one common center point;
and overlapping the at least one common central point to obtain an overlapped section image of the exploration area.
10. An apparatus for virtual source imaging, the apparatus comprising:
a first determining module, configured to determine a first detector and a second detector from a plurality of detectors in a survey area, with the first detector as a first virtual source;
a second determination module for determining, for each of a plurality of earth seismic sources in a survey area, an actual down-going wavefield of the earth source-excited seismic waves received by the first virtual source and a first actual up-going wavefield of the earth source-excited seismic waves received by the second detector;
the time domain cross-correlation module is used for performing time domain cross-correlation on the actual downlink wave field and the first actual uplink wave field to obtain a cross-correlation channel of the first virtual source and the second detector;
the time domain convolution module is used for performing time domain convolution on the actual downlink wave field and the cross-correlation channel to obtain a first predicted uplink wave field;
the processing module is used for carrying out normalized cross-correlation processing on the first predicted upgoing wave field and the first actual upgoing wave field to obtain the weight of the earth surface seismic source;
the third determining module is used for determining a first virtual gather corresponding to the first virtual seismic source according to the cross-correlation trace and the weight of the surface seismic source;
a fourth determining module, configured to determine, by using at least one third detector, except for the first detector, of the multiple detectors in the exploration area as a second virtual source, at least one second virtual gather corresponding to the second virtual source;
a fifth determination module to determine an overlay profile image of the survey area based on the first set of virtual traces and at least one of the second set of virtual traces.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010670466.8A CN113933894A (en) | 2020-07-13 | 2020-07-13 | Virtual source imaging method and device |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010670466.8A CN113933894A (en) | 2020-07-13 | 2020-07-13 | Virtual source imaging method and device |
Publications (1)
Publication Number | Publication Date |
---|---|
CN113933894A true CN113933894A (en) | 2022-01-14 |
Family
ID=79273667
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010670466.8A Pending CN113933894A (en) | 2020-07-13 | 2020-07-13 | Virtual source imaging method and device |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113933894A (en) |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110130967A1 (en) * | 2008-01-11 | 2011-06-02 | Andrey Victorovich Bakulin | Method of correcting amplitudes in virtual source imaging of seismic data |
US20190187107A1 (en) * | 2017-12-15 | 2019-06-20 | Regents Of The University Of Minnesota | Methods for ultrasonic non-destructive testing using analytical reverse time migration |
CN111257939A (en) * | 2020-03-26 | 2020-06-09 | 中国石油大学(北京) | Time-lapse seismic virtual source bidirectional wave field reconstruction method and system |
CN111257938A (en) * | 2020-03-25 | 2020-06-09 | 中国石油大学(北京) | Time-lapse seismic virtual source wave field reconstruction method and system based on wavelet cross-correlation |
-
2020
- 2020-07-13 CN CN202010670466.8A patent/CN113933894A/en active Pending
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110130967A1 (en) * | 2008-01-11 | 2011-06-02 | Andrey Victorovich Bakulin | Method of correcting amplitudes in virtual source imaging of seismic data |
US20190187107A1 (en) * | 2017-12-15 | 2019-06-20 | Regents Of The University Of Minnesota | Methods for ultrasonic non-destructive testing using analytical reverse time migration |
CN111257938A (en) * | 2020-03-25 | 2020-06-09 | 中国石油大学(北京) | Time-lapse seismic virtual source wave field reconstruction method and system based on wavelet cross-correlation |
CN111257939A (en) * | 2020-03-26 | 2020-06-09 | 中国石油大学(北京) | Time-lapse seismic virtual source bidirectional wave field reconstruction method and system |
Non-Patent Citations (1)
Title |
---|
YANG ZHAO ET AL.: ""Target-oriented diversity stacking for virtual-source imaging and monitoring"", 《THE LEADING EDGE》, pages 766 - 773 * |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US5995905A (en) | Source signature determination and multiple reflection reduction | |
US6058074A (en) | Method and system for detecting hydrocarbon reservoirs using amplitude-versus-offset analysis with improved measurement of background statistics | |
EP2786176B1 (en) | Separation of simultaneous source data | |
CN101556339B (en) | Method for deghosting marine seismic streamer data with irregular receiver positions | |
US5621700A (en) | Method for attenuation of reverberations using a pressure-velocity bottom cable | |
WO2006054181A1 (en) | Method for processing at least two sets of seismic data | |
GB2448415A (en) | Method for prediction of surface related multiples from marine towed dual sensor seismic streamer data | |
AU2014280832B2 (en) | Seismic data spectrum restoring and broadening | |
CA2741888A1 (en) | 4d seismic signal analysis | |
CN112327362B (en) | Submarine multiple prediction and tracking attenuation method in velocity domain | |
US4363113A (en) | Seismic exploration with simulated plane waves | |
US11442186B2 (en) | Source-receiver position estimation using direct arrival modeling and inversion | |
Lu et al. | 3D supervirtual refraction interferometry | |
CN110579799B (en) | Seismic acquisition observation method and system with equal travel time intervals | |
CN113933894A (en) | Virtual source imaging method and device | |
CN115061197A (en) | Two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal and flow measurement equipment | |
CN112505782B (en) | Interference reference plane reconstruction method and system for radiation mode correction in four-dimensional earthquake | |
Rabinovich et al. | Location technology for construction of seismic images | |
GB2368911A (en) | Computing a stacked seismic line by interpolation between known stacks | |
US4682307A (en) | Underwater seismic testing | |
CN105572742A (en) | Method and device for determining depth of seawater | |
CN111665549A (en) | Inversion method of stratum acoustic wave attenuation factor | |
CN111665546A (en) | Acoustic parameter acquisition method for combustible ice detection | |
CN116609834B (en) | Data processing method based on ocean vertical cable seismic exploration | |
AU2018325477B2 (en) | Source-receiver position estimation using direct arrival modeling and inversion |
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 |