CN113933894A

CN113933894A - Virtual source imaging method and device

Info

Publication number: CN113933894A
Application number: CN202010670466.8A
Authority: CN
Inventors: 刘文平; 赵杨; 周阳田丽; 雍锐; 左佳卉; 吴建发; 杨学锋; 石学文; 文山师; 王星皓
Original assignee: Petrochina Co Ltd
Current assignee: Petrochina Co Ltd
Priority date: 2020-07-13
Filing date: 2020-07-13
Publication date: 2022-01-14

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

Virtual source imaging method and device

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 is_A，r_BAnd r_SRespectively 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; v_w(r_B|r_A(ii) a t) represents the first virtual gather; d (r)_A|r_S(ii) a t) represents the actual down-going wavefield; u (r)_B|r_S(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|r_S(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|r_S(ii) a f) Representing the second actual upgoing wavefield and f representing the frequency of the fourier transform.

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 that

A 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:

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 of_rRepresenting 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 is_A，r_BAnd r_SRespectively 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; v_w(r_B|r_A(ii) a t) represents the first virtual gather; d (r)_A|r_S(ii) a t) represents the actual down-going wavefield; u (r)_B|r_S(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 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;

the formula III is as follows:

wherein,

representing a weight of the surface seismic source;

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:

representing the first wavelet coefficients; is that

A 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:

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;

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 is_A，r_BAnd r_SRespectively 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|r_S(ii) a t) represents the actual down-going wavefield; u (r)_B|r_S(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|r_S(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;

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:

representing a first wavelet coefficient; is that

Seismic 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:

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 of_rRepresenting 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; v_w(r_B|r_A(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 is_A，r_BAnd r_SRespectively 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|r_S(ii) a t) represents the actual down-going wavefield; u (r)_B|r_S(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 volume

As shown in fig. 5. Wherein the second time-frequency space data body corresponding to the first prediction upgoing wave field

As 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.

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 is_A，r_BAnd r_SRespectively 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; v_w(r_B|r_A(ii) a t) represents a first virtual gather; d(r_A|r_S(ii) a t) represents the actual down-going wavefield; u (r)_B|r_S(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 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;

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:

representing a first wavelet coefficient; is that

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:

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;

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.

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

1. A method of virtual source imaging, the method comprising:

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:

the formula I is as follows:

wherein,

a weight representing the earth's surface seismic source,

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:

the formula II is as follows:

wherein,

representing the cross-correlation trace.

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:

the formula III is as follows:

wherein,

representing a weight of the surface seismic source;

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:

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:

the formula four is as follows:

representing the first wavelet coefficients; is that

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:

the formula five is as follows:

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:

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;

9. The method of claim 1, wherein said determining from said first virtual gather and at least one of said second virtual gathers comprises:

10. An apparatus for virtual source imaging, the apparatus comprising: