CN111257939B

CN111257939B - Time-lapse seismic virtual source bidirectional wave field reconstruction method and system

Info

Publication number: CN111257939B
Application number: CN202010222440.7A
Authority: CN
Inventors: 赵杨; 周阳田丽; 左佳卉
Original assignee: China University of Petroleum Beijing
Current assignee: China University of Petroleum Beijing
Priority date: 2020-03-26
Filing date: 2020-03-26
Publication date: 2021-06-01
Anticipated expiration: 2040-03-26
Also published as: CN111257939A

Abstract

The invention relates to a time-lapse seismic virtual source bidirectional wave field reconstruction method and a time-lapse seismic virtual source bidirectional wave field reconstruction system, which comprise the following steps: s1, converting a pre-acquired seismic data signal into a frequency domain, reconstructing wave fields of a seismic source end and a receiving end, and reconstructing a seismic source and a receiver below the earth surface; s2, transforming the reconstructed seismic data signal back to a time domain; and S3, carrying out common center point superposition on the seismic data signals of the time domain to obtain a virtual seismic source excitation wave field. The seismic source and the detector are reconstructed to an observation system below the earth surface, reconstruction data are interpolated, the influence of a stratum structure on a wave field is reduced, and high-quality seismic wave field imaging is realized at lower cost.

Description

Time-lapse seismic virtual source bidirectional wave field reconstruction method and system

Technical Field

The invention relates to a time-lapse seismic virtual source bidirectional wave field reconstruction method and a time-lapse seismic virtual source bidirectional wave field reconstruction system, and belongs to the technical field of seismic data wave processing.

Background

The time-lapse seismic technology is a technology for monitoring a reservoir, perfecting an oil reservoir management scheme and improving oil and gas recovery ratio by using the difference of effective seismic information observed at different times. The repeatability of time lapse seismic is the most important influencing factor for long-term reservoir monitoring. The repeatability of the time-lapse earthquake is closely related to acquisition factors such as instrument factors, environment, background noise and the like in acquisition. The instrument factors are easy to control, and the environmental and background noises are often difficult to control, which is a main reason for unrepeatable monitoring results. Environmental, background noise is mainly due to time delays and phase and amplitude distortions caused by changes in the source, detector and near-surface conditions. In the prior art, the repeatability of time-lapse seismic is generally improved by burying the source and the detector in the ground. However, in desert environments, the imaging results of areas with a karst structure developed near the surface are very complex, and multiple waves, surface waves, scattering noise and interference of different wave types are introduced, so that the up-going wave field and the down-going wave field are distorted.

In order to solve the problems, a virtual seismic source method (VS) is introduced to detect the time-shifting earthquake, the virtual seismic source method does not depend on a near-surface velocity model, cross-correlation reconstruction is carried out on an up-going wave field and a down-going wave field, a ground seismic source is reconstructed to the position of an underground observation system, noise formed by near-surface instrument coupling and observation system change is suppressed, and repeatability is improved in reservoir monitoring and time-shifting earthquake processing. However, all interference-based wavefield reconstruction methods involve only one-way reconstruction, i.e., reconstructing one end of the source or detector to the position of the subsurface observation system, which can avoid the influence of the surface factors on the wavefield, but since the source or detector is located in the subsurface observation system and the opposite detector or source is located at the surface, the reconstructed wavefield still needs to pass through the stratigraphic structure. Moreover, in order to obtain a high-quality imaging effect, a high-density observation system is required, which is costly.

Disclosure of Invention

In view of the above-mentioned deficiencies of the prior art, the present invention aims to provide a time-lapse seismic virtual source bidirectional wave field reconstruction method and system, which reconstructs a seismic source and a detector to an observation system below the earth surface, interpolates reconstruction data, reduces the influence of a stratum structure on a wave field, and realizes high-quality seismic wave field imaging with low cost.

In order to achieve the aim, the invention provides a time-lapse seismic virtual source bidirectional wave field reconstruction method, which comprises the following steps: s1, converting a pre-acquired seismic data signal into a frequency domain, reconstructing wave fields of a seismic source end and a receiving end, and reconstructing a seismic source and a receiver below the earth surface; s2, transforming the reconstructed seismic data signal back to a time domain; s3, common center point superposition is carried out on the seismic data signals of the time domain to obtain a virtual seismic source excitation wave field

Further, the seismic data signals are divided into underground observation data and surface observation data, and the underground observation data comprises direct waves at a seismic source end and direct waves at a receiving end; the earth surface observation data comprises reflected waves for removing surface waves, and direct waves at the seismic source end, direct waves at the receiving end and reflected waves are converted into a frequency domain.

Further, in step S2, the wavefield reconstruction formula at the source end is:

x, Y and B are the space coordinates of the seismic source, the earth surface receiver and the underground receiver respectively; v (BY; omega) is a wave field received by the earth surface receiver Y when the underground receiver B is used as a vertical virtual seismic source; v^*(BY; ω) is the complex conjugate of V (BY; ω); d^*(X | B; omega) is the direct wavefield excited by the seismic source received by the underground receiver B; u (X | Y; omega) is the earth detector Y receptionA reflected up-going wave arriving from the seismic source X; n is_xIs the normal vector of the integration surface.

Further, in step S2, the wavefield reconstruction formula at the receiving end is:

where V (A | B; ω) is the reconstructed wavefield at another subsurface receiver B when the subsurface receiver A is considered a virtual source; k is the wave number emitted by the surface seismic source; d^*(X | A) is the wavefield excited by the surface source X received by the geophone A; d^*(BxX; ω) is the reconstructed wavefield received at source X, with subsurface receiver B considered as the virtual source.

Further, to D^*(X | A) and D^*And (B | X; omega) carrying out zero-phase whitening filtering, wherein the reconstructed wave field after filtering is as follows:

where H is the zero-phase whitening filter function and x is the complex conjugate.

Further, the formula for transforming the reconstructed seismic data signal back to the time domain in step S3 is:

V_H(A|B；t)＝F_ω→t{V_H(A|B；ω)}

wherein, F_ω→tRepresenting a frequency domain to time domain fourier transform.

Further, the receivers of the subsurface observation data are less than those of the surface observation data, and in order to compensate for the insufficient coverage times of the subsurface receivers, the seismic data signals need to be interpolated after being converted into the time domain in step S3.

Further, processing the seismic data signals by adopting a plane wave decomposition interpolation method, wherein the interpolation is written as a least square equation:

wherein, V_nIs the reconstructed seismic data signal of the nth trace, P_n,n-1Is the dip angle domain plane wave decomposition operator between the nth track and the (n-2) th track, and is determined by the formula:

and solving a least square equation.

Further, in the plane wave decomposition interpolation method, the reconstructed wave field obtained after interpolation is as follows:

wherein,

is smoothed to form an expanded tilt angle field

The length of (a) of (b),

is that

And reconstructing data after the strip difference value.

The invention also discloses a time-lapse seismic virtual source bidirectional wave field reconstruction system, which comprises: the bidirectional wave field reconstruction module is used for converting the pre-acquired seismic data signals into a frequency domain, reconstructing wave fields of a seismic source end and a receiving end and reconstructing a seismic source and a receiver below the earth surface; a time domain conversion module for converting the reconstructed seismic data signal back to the time domain; and the common-center-point superposition module is used for carrying out common-center-point superposition on the seismic data signals of the time domain to obtain a virtual seismic source excitation wave field.

Due to the adoption of the technical scheme, the invention has the following advantages:

1. the invention uses the sparse observation system to reconstruct the seismic source and the detector to the observation system below the earth surface, interpolates the reconstructed data, reduces the influence of the stratum structure on the wave field, and realizes high-quality seismic wave field imaging with lower cost.

2. The method expands the reconstruction method of the virtual source wave field from one direction to two directions, further suppresses the complex near-surface interference and improves the imaging quality and repeatability;

3. the dip-guided seismic channel interpolation method is applied to wavelength reconstruction, so that the defect of the coverage times of a sparse observation system is compensated, and the cost is reduced.

Drawings

FIG. 1 is a flow chart of a method for reconstructing a bidirectional wavefield of a time-lapse seismic virtual source according to an embodiment of the invention;

FIG. 2 is a schematic diagram of the principle of time-lapse seismic virtual wavefield detection in a complex near-surface context;

FIG. 3 is a seismic signal received by a surface receiver Y;

FIG. 4 is a seismic signal received by subsurface receivers A and B;

FIG. 5 is a schematic diagram of source-side wavefield reconstruction in a bi-directional wavefield reconstruction;

FIG. 6 is a schematic diagram of the wavefield reconstruction at the receiving end in a bi-directional wavefield reconstruction;

FIG. 7 is a schematic of a down-going direct wave D (Y | A) in the reconstructed received wavefield;

FIG. 8 is a single trace signal trace with the direct arrival D (X | A) from the seismic source X to the subsurface geophone point A arranged in a line offset;

FIG. 9(a) of FIG. 9 is a single shot signal record obtained after bi-directional wavefield reconstruction; FIG. 9(b) is a single shot signal record obtained after initial tilt estimation using a plane wave decomposition filter; FIG. 9(c) is a single shot signal record after interpolation;

FIG. 10 is a schematic diagram of the wavefield detection after bi-directional wavefield reconstruction and interpolation;

FIG. 11 is a layout view of a seismic survey data observation system;

FIG. 12 is a schematic view of a three-dimensional data volume acquired in accordance with an embodiment of the present invention;

FIG. 13 is a cross-sectional view of the common midpoint stacking of the earth surface observation system, FIG. 13 is a cross-sectional view of the common midpoint stacking of the subsurface observation system, and FIG. 13 is a cross-sectional view of the common midpoint stacking of the two-way reconstructed wave fields;

FIG. 14 is an overlay profile comparison of a first survey set (sets 1-6) and a second survey set (sets 7-13), with the first survey set on the left and the second survey set on the right, FIG. 14(a) being an overlay profile of the common center of the earth's surface survey system, FIG. 14(b) being an overlay profile of the common center of the subsurface survey system, and FIG. 14(c) being an overlay profile of the common center of the wavefield through two-way reconstruction;

FIG. 15 is a schematic of root mean square mean values for a surface observation system, a subsurface observation system, and a bidirectionally reconstructed wavefield system.

Detailed Description

The present invention is described in detail by way of specific embodiments in order to better understand the technical direction of the present invention for those skilled in the art. It should be understood, however, that the detailed description is provided for a better understanding of the invention only and that they should not be taken as limiting the invention. In describing the present invention, it is to be understood that the terminology used is for the purpose of description only and is not intended to be indicative or implied of relative importance.

Example one

The embodiment discloses a time-lapse seismic virtual source bidirectional wave field reconstruction method, as shown in fig. 1, comprising the following steps:

s1, acquiring seismic data signals;

s2, converting the seismic data signals into a frequency domain, reconstructing wave fields of a seismic source end and a receiving end, and reconstructing a seismic source and a receiver below the earth surface;

s3, transforming the reconstructed seismic data signal back to a time domain;

and S4, carrying out common center point superposition on the seismic data signals to obtain a virtual seismic source excitation wave field.

In the method in the embodiment, the seismic source and the detector are reconstructed to the observation system below the earth surface, and interpolation is performed on the reconstructed data, so that the influence of a stratum structure on a wave field is reduced, and high-quality seismic wave field imaging is realized at lower cost.

The seismic data signals are divided into underground observation data and surface observation data, and the underground observation data comprise direct waves at a seismic source end and direct waves at a receiving end; the earth surface observation data comprises reflected waves for removing surface waves, and direct waves at the seismic source end, direct waves at the receiving end and reflected waves are converted into a frequency domain.

The reflected wave seismic interference method is that different receivers receive seismic signals sent by the same seismic source, the seismic signals are subjected to cross correlation, and then each relevant seismic source is subjected to integration to obtain a seismic wave field diagram. Fig. 2 is a schematic diagram of the principle of time-lapse seismic virtual wavefield detection in a complex near-surface background, in which a star pattern represents a seismic source X and an inverted triangular pattern represents a surface receiver Y, as shown in fig. 2. Both a and B represent receivers disposed in the ground. The receiver here is typically a geophone for receiving seismic signals. The bottom circular pattern in the figure represents the target reflector. The method comprises the steps that a down direct wave D (X | B) generated by a seismic source is received and recorded by an underground receiver B, is converted into an up-going wave U (X | A) after being reflected by a target reflector and is received and recorded by the underground receiver A, and the up-going wave U (X | A) continues to propagate to the earth surface until being received and recorded by an earth surface receiver Y to generate an original earth surface record U (X | Y). FIG. 3 is a seismic signal received at a surface receiver Y. The horizontal axis in fig. 3 is the shot number, each trace in a single shot record is arranged by offset, and the vertical axis is time. FIG. 4 is a seismic signal received by subsurface receiver A and subsurface receiver B. Fig. 4 is similar to fig. 3, with the horizontal axis in fig. 4 representing shot number, each trace in a single shot record arranged by offset, and the vertical axis representing time. Comparing the seismic signal plot of FIG. 4 with the seismic signal plot of FIG. 3, it can be seen that the signals recorded by the subsurface receivers A and B have a higher signal-to-noise ratio and better reflection. Therein, the original near offset wavefield is considered to be the direct down-going wavefield (offset <30m, time <100ms) because the early wavefield consists mainly of direct P-waves. Similarly, wavefields arriving later in time (time >300ms) are considered upgoing reflected wavefields. Based on the analysis of fig. 3 and 4, the wavefields of the sources and receivers located at the surface need to be reconstructed below the surface to avoid the effects of surface conditions and stratigraphic structures on the seismic wavefield.

In the present embodiment, the principle of wave field reconstruction at the seismic source end in step S2 is shown in fig. 5, where seismic signal data U (X | Y) received by an earth surface receiver Y and a down-going direct wave D (X | B) received by an underground receiver B, which are emitted by an earth surface seismic source, are cross-correlated, and a plurality of single shot records are further superimposed to obtain a reconstructed wave field V (B | Y).

V (B | Y) satisfies the following formula:

where ω is the angular frequency and X, Y and B are the spatial coordinates of the seismic source, surface receivers and subsurface receivers, respectively; v (BY; omega) is a wave field received by the earth surface receiver Y when the underground receiver B is used as a vertical virtual seismic source; v^*(BY; ω) is the complex conjugate of V (BY; ω); d^*(X | B; omega) is the direct wavefield excited by the seismic source received by the underground receiver B; u (X | Y; ω) is the reflected up-going wave from source X received by surface detector Y; n is_xIs the normal vector of the integration surface.

Equation (1) can be simplified approximately to the following form in the case of far-field seismic:

V(B|Y；ω)-V^*(B|Y；ω)≈2ik∫_surfaceD^*(X|B；ω)U(X|Y；ω)d²X (2)

k is the wave number emitted by the surface seismic source. Far-field earthquakes are earthquakes with a seismic center distance of more than 1000 kilometers. Although source wavefield reconstruction reconstructs the subsurface receivers as virtual sources V (B | Y; ω), it is still not completely free of complex near-surface interference, as the seismic signals inevitably pass through the near-surface layer as they reach the surface receivers Y.

Therefore, the wavefield at the receiving end needs to be reconstructed into the subsurface receiver wavefield, as shown in fig. 6, fig. 6 is a schematic diagram of the wavefield reconstruction at the receiving end. And converting the signal V (B | Y) with the underground receiver B as a virtual seismic source into a signal V (Y | B) with the earth surface receiver Y as a virtual seismic source, then performing cross correlation with a direct wave D (Y | A) received by the underground detector A from the virtual seismic source Y, and superposing a plurality of single-shot signals with Y as the seismic source to generate a reconstructed wave field V (A | B).

The wave field reconstruction V (A | B; omega) at the receiving end meets the following formula:

V(A|B；ω)-V^*(A|B；ω)≈2ik∫_buriedD^*(Y|A；ω)U(Y|B；ω)d²Y (3)

where V (A | B; ω) is the reconstructed wavefield at another subsurface receiver B when the subsurface receiver A is considered a virtual source; v^*(A | B; ω) is the complex conjugate of V (A | B; ω); k is the wave number emitted by the surface seismic source; d^*(Y | A) is the wavefield excited by the virtual source Y received by the geophone A; u (Y | B; ω) is the reflected up-going wave from virtual source Y received by subsurface receiver B. (ii) a D^*(Y | A) is the wavefield excited by the virtual source Y received by the geophone A, but Y itself cannot excite seismic waves, so it is desirable to place the real source X as close as possible to the surface receiver Y, as shown in FIG. 7. FIG. 7 is a schematic of a down-going direct wave D (Y | A) in the reconstructed receive wavefield. The real source X is very close to the surface receiver Y, so D (X | A) ≈ D (Y | A). FIG. 8 is a single trace signal trace with the direct arrival D (X | A) from the source X to the subsurface geophone point A arranged in a line offset. It can be seen from fig. 8 that the waveform of the direct wave is very complex, representing the complexity of the near-surface structure. Since subsurface geophone point a is a single component geophone in fig. 8 and detects the vertical component of the seismic wave, small angle waves exhibit large amplitude variations, while large angle waves exhibit relatively small amplitude variations because this is a vertical component recording.

Based on the approximation in fig. 7, substituting equation (2) into equation (3) yields:

by substituting D (X | A; ω) into U (Y | X; ω), equation (4) can be simplified as:

V(A|B；ω)≈-4k²∫_surface∫_buriedD^*(X|A)U(Y|X；ω)D^*(B|X；ω)d²Xd²Y(5)

where V (A | B; ω) is the reconstructed wavefield at another subsurface receiver B when the subsurface receiver A is considered a virtual source; k is the wave number emitted by the surface seismic source; d^*(X | A) is the wavefield excited by the surface source X received by the geophone A; d^*(BxX; ω) is the reconstructed wavefield received at source X, with subsurface receiver B considered as the virtual source. In co-channel concentration, the surface data can be cross-correlated with the direct wave recorded by the geophone, and bidirectional wave field reconstruction is realized. Unlike classical seismic wavefield reconstruction methods, equation (5) is driven entirely by data, and is less affected by the surface layer and artifacts. FIG. 9(a) of FIG. 9 is a single shot signal record obtained after bi-directional wavefield reconstruction, where the signal-to-noise ratio of the single shot signal record of FIG. 9(a) is higher, and in particular, the reflection effect is better for shallow targets, than the single shot signal record of FIGS. 3 and 4, thereby reducing the effect of near-surface complexity on the seismic wavefield.

In the frequency domain, the spectrum of the bi-directionally reconstructed wavefield V (A | B) is less accurate than the original surface data. To obtain a more accurate comparison of the two-way reconstructed wavefield to the non-reconstructed wavefield, a zero-phase whitening filter function H needs to be applied to D (X | A) and D (B | X) to ensure that the spectrum of the two-way reconstructed wavefield is consistent with U (Y | X). The final reconstruction result can be obtained by adding the zero-phase whitening filter function H to equation (5):

V_H(A|B；ω)≈-4k²∫_surface∫_buried[HD^*(X|A；ω)]^*U(Y|X；ω)[HD^*(B|X；ω)]^*d²Xd²Y (6)

The receivers of the subsurface observation data are less than those of the surface observation data, and in order to compensate for the insufficient coverage times of the subsurface receivers, the seismic data signals need to be interpolated after being converted into the time domain in step S3.

V_H(A | B; omega) contains horizontal reflecting layers with small lateral variations, and is suitable for the dip-guided seismic trace interpolation method based on the plane wave decomposition filter.

A plane wave decomposition filter (PWD) is a method of estimating the dip of a seismic record by phase shifting the previous trace along the principal in-phase axis to predict the next trace with preserved amplitude. The estimated principal slope is then calculated by minimizing the prediction error. And the shape regularization ensures that the estimated inclination angle in the reconstructed wave field changes stably. Since PWDs are typically applied in the time domain, the reconstructed seismic data signal needs to be transformed back into the time domain via the following equation:

V_H(A|B；t)＝F_ω→t{V_H(A|B；ω)} (7)

Processing the seismic data signals by adopting a plane wave decomposition interpolation method, wherein the interpolation is written as a least square equation: v_HP ≈ δ, or

Wherein, V_nIs the reconstructed seismic data signal of the nth trace, Pn, n-1 is the dip domain plane wave decomposition operator between the nth trace and the n-2 th trace, by the formula:

and solving a least square equation.

The formula:

can be expressed as:

since the inverse of the diagonal matrix is very easy to solve, equation (9) can be solved directly to obtain the tilt angle field. FIG. 9(b) is a single shot signal record obtained after initial tilt estimation using a plane wave decomposition filter; smoothing the tilt field to form an extended tilt field

Extended tilt field

Length of (2)

Reconstruction of wavefields by interpolation with targets

Are the same length.

Another reconstruction

The method of (1) is to solve the linear system using the inverse PWD as a linear operator. The linear problem is to convert the matrix in equation (9) to

Or:

using structure constrained interpolation, the interpolated reconstructed wavefield is:

wherein,

is smoothed to form an expanded tilt angle field

The length of (a) of (b),

is that

And reconstructing data after the strip difference value. Equation (11) is a typical band diagonal matrix that can be solved using the conjugate gradient method. The trace continuity of the interpolated seismic signal record is enhanced while preserving slope information. FIG. 9(c) is a single shot signal record after interpolation. Comparing the single shot signal recordings in fig. 9(c) and fig. 9(a), it can be seen that the single shot signal recording in fig. 9(c) is clearer, the details are more obvious, the noise signal is weaker, and the signal-to-noise ratio is enhanced.

FIG. 10 is a schematic diagram of the principle of wavefield detection after bi-directional wavefield reconstruction and interpolation. Through the schematic diagram, the bidirectional wave field interference cross-correlation can be seen to distribute the seismic sources and the observation system underground, so that the complex near-surface interference is avoided.

Example two

The present example further illustrates and tests the method of the first example using time lapse seismic data from a desert environment.

In desert regions, 13 repeated two-dimensional surveys are conducted within 19 months, the first 6 surveys (S1-S6) are completed within 3 months, and 17 months are interrupted and 7 surveys (S7-S13) are completed within a week. The Mertz 26 controllable seismic source is used in each of the 13 surveys, and the repeated excitation precision of the seismic source reaches 1 m. The geophone layout is shown in FIG. 11, where the triangular pattern is 80 geophones and the dotted pattern is the seismic source of the earth's surface. All seismic sources make up 9 shot lines. And 300 surface detectors are further distributed near the middle shot line, the survey line containing 300 surface detectors is close to the fifth shot line, namely one shot line in the middle of the 9 shot lines, the track spacing is 7.5m, and the minimum offset distance is less than 5 m. The near-surface is a sand layer with the thickness ranging from several meters to dozens of meters, and a layered structure with the inclination angle smaller than 5 degrees is arranged below the near-surface. The imaging quality and reproducibility of the data is disturbed by near-surface structures.

In the embodiment, the field data is processed in three steps by acquiring seismic data signals: the first step involves nine shot line multimode orthogonal summations to remove the overly strong shot signal, random energy and scattered noise. And integrating the common shot channel level into a common detection wave channel level by using high-density shot points, and suppressing the surface wave noise by using FK filtering. The second step is a bi-directional wavefield reconstruction of the data. The direct waves in the 60ms window are picked up by an automatic algorithm, and reconstruction operators D (X | A) and D (B | X) are obtained according to the direct waves. And thirdly, performing static correction and dynamic correction, amplitude compensation, top cut and common center point superposition on the data. The amplitude compensation is one by one for 5 500ms windows per track. Common midpoint stacking includes arithmetic summation and energy averaging for gathers of different coverage times. Two non-reconstructed conventional common center point superposition methods for comparison use only the surface wavefield U (Y | X) and the subsurface wavefield V (Y | B), without a second step of reconstruction processing.

Fig. 12 is a schematic diagram of the three-dimensional data volume acquired in the present embodiment, in which X and Y axes respectively represent a surface shot point (2700) and a geophone point (300/80), an XZ plane represents a common gather, and a YZ plane represents a common shot gather. In the embodiment, a matrix transposition method is used for conversion between single-shot recording and single-track recording.

Fig. 13(a) is a cross-sectional view showing the common center of the earth surface observation system, and fig. 13(b) is a cross-sectional view showing the common center of the subsurface observation system, in which the layer indicated by the arrow is the target reflection layer. In fig. 13(a), the signal-to-noise ratio on the superimposed profile of the earth's surface survey line is the lowest, and particularly, the right side of the profile is affected by the desert environment of the earth's surface, and the imaging effect is very poor. As shown in FIG. 13(b), the cross-sectional view of the subsurface observation system of FIG. 13(b) superimposed at a common center point, even without reconstruction, has better signal-to-noise ratio and continuity than in FIG. 13(a) because the geophone portion avoids the influence of the subsurface on the seismic waves. Fig. 13(c) is a superimposed sectional view after the bidirectional wave field reconstruction, and the imaging quality of the sectional view in fig. 13(c) is significantly improved as compared with fig. 13(b) of fig. 13(a), and the wave field continuity is enhanced without being disturbed by the surface sand layer even in a shallow target layer. These results demonstrate the robustness of the method in this embodiment, and these imaging results can be further applied to formation velocity analysis and migration.

In addition to image quality, repeatability is also an important factor in measuring seismic monitoring quality. To compare the repeatability of the test system in this example, one was randomly drawn from two surveys (underground line-invariant) spaced one year apart. FIG. 14 is a comparison of the stacked profiles of a first survey (sets 1-6) and a second survey (sets 7-13), with the first survey profile on the left and the second survey profile on the right, FIG. 14(a) being a stacked profile of the common center of the earth's surface survey system, FIG. 14(b) being a stacked profile of the common center of the subsurface survey system, and FIG. 14(c) being a stacked profile of the common center of the wavefield via two-way reconstruction, as shown in FIG. 14. As can be seen in fig. 14(a), there is a significant difference between the left and right portions of the cross-sectional view, i.e. the results of the two surveys are poorly reproducible, particularly in the target zone indicated by the arrow in the figure. The section of fig. 14(b) also differs in the left and right parts, but the continuity of the wavefield is better in the section of fig. 14(b), especially for the target layer at the arrow mark, than in fig. 14 (a). Because the subsurface lines reduce part of the near surface interference at the receiving end. But there are still partial inconsistencies as the seismic source ends still pass near the surface. The difference between the left part and the right part of the cross-section of fig. 14(c) is small, and it can be seen that in the embodiment, the repeatability of detection can be obviously improved by re-referencing both the surface seismic source and the detector to be below the near-surface.

FIG. 15 is a schematic of root mean square mean values for a surface observation system, a subsurface observation system, and a bidirectionally reconstructed wavefield system. The reproducibility of the test can be quantitatively evaluated by a method of calculating Normalized Root Mean Square (NRMS) within the window of the target layer. NRMS is the percentage of the mean of the root mean square of the two seismic sections within a given window. NMRS is very sensitive to small changes in seismic data. NRMS is calculated for 13 exploration results respectively, and 78 NRMS values can be obtained for each common midpoint gather. NRMS histograms for the three methods are shown in fig. 15. The earth surface observation system has a bimodal distribution, and the peak values respectively appear at 35% and 100%. The NRMS values between the two surveys were large, with the two separate peaks having different shapes, indicating poor repeatability of the surface observation system. The same is true for subsurface observation systems, with peaks at about 20% and 70%. The distribution is more concentrated relative to the ground surface observation system, but the problem of inconsistency still exists between two groups of surveys. The data NRMS reconstructed by the bidirectional reconstructed wave field system presents unimodal distribution, and the peak value is 32%. Compared with the traditional virtual source method, the imaging quality and repeatability are remarkably improved through the bidirectional reconstruction wave field system.

EXAMPLE III

Based on the same inventive concept, the embodiment discloses a time-lapse seismic virtual source bidirectional wave field reconstruction system, which comprises:

the data acquisition module is used for acquiring seismic data signals;

the bidirectional wave field reconstruction module is used for converting the seismic data signals into a frequency domain, reconstructing wave fields of the seismic source end and the receiving end and reconstructing a seismic source and a receiver below the earth surface;

a time domain conversion module for converting the reconstructed seismic data signal back to the time domain;

and the common-center-point superposition module is used for carrying out common-center-point superposition on the seismic data signals to obtain a virtual seismic source excitation wave field.

The above description is only for the specific embodiments of the present application, but the scope of the present application is not limited thereto, and any person skilled in the art can easily conceive of the changes or substitutions within the technical scope of the present application, and shall be covered by the scope of the present application. Therefore, the protection scope of the present application shall be subject to the protection scope of the claims.

Claims

1. A time-lapse seismic virtual source bidirectional wave field reconstruction method is characterized by comprising the following steps:

s1, converting a pre-acquired seismic data signal into a frequency domain, reconstructing wave fields of a seismic source end and a receiving end, and reconstructing a seismic source and a receiver below the earth surface;

s2, transforming the reconstructed seismic data signal back to a time domain;

and S3, carrying out common center point superposition on the seismic data signals of the time domain to obtain a virtual seismic source excitation wave field.

2. The method for reconstructing the bidirectional wave field of the time-lapse seismic virtual source of claim 1, wherein the seismic data signals are divided into subsurface observation data and surface observation data, and the subsurface observation data comprise direct waves at a seismic source end and direct waves at a receiving end; and the earth surface observation data comprises reflected waves for removing surface waves, and direct waves at the seismic source end, direct waves at the receiving end and reflected waves are converted into a frequency domain.

3. The method for reconstructing a bidirectional wavefield of a time-lapse seismic virtual source as claimed in claim 2, wherein the wavefield reconstruction formula at the source end in step S1 is:

x, Y and B are the space coordinates of the seismic source, the earth surface receiver and the underground receiver respectively; v (BY; omega) is a wave field received by the earth surface receiver Y when the underground receiver B is used as a vertical virtual seismic source; v^*(BY; ω) is the complex conjugate of V (BY; ω); d^*(X | B; omega) is the direct wavefield excited by the seismic source received by the underground receiver B; u (X | Y; ω) is the reflected up-going wave from source X received by surface receiver Y; n is_xIs the normal vector of the integration surface.

4. The time-lapse seismic virtual source bidirectional wavefield reconstruction method of claim 3, wherein the receiving-end wavefield reconstruction formula in step S1 is:

where V (A | B; ω) is the reconstructed wavefield at another subsurface receiver B when the subsurface receiver A is considered a virtual source; k is the wave number emitted by the surface seismic source; d^*(X | A) is the wavefield excited by the surface source X received by the subsurface receiver A; d^*(BxX; ω) is the reconstructed wavefield received at source X, with subsurface receiver B considered as the virtual source.

5. The method of reconstructing a time-lapse seismic virtual source bi-directional wavefield of claim 4, wherein D is the number of seismic events in the seismic data set^*(X | A) and D^*And (B | X; omega) carrying out zero-phase whitening filtering, wherein the reconstructed wave field after filtering is as follows:

6. The method of reconstructing a time-lapse seismic virtual source bi-directional wavefield of claim 5, wherein the formula for transforming the reconstructed seismic data signals back into the time domain in step S2 is:

V_H(A|B；t)＝F_ω→t{V_H(A|B；ω)}

7. The method for reconstructing the two-way wavefield of a time-lapse seismic virtual source of any one of claims 2-5, wherein the number of receivers of the subsurface observation data is less than that of the surface observation data, and in order to compensate for the insufficient number of times of coverage of the subsurface receivers, the seismic data signals are interpolated after being converted into the time domain in step S2.

8. The method of reconstructing a two-way wavefield of a time-lapse seismic virtual source of claim 7, wherein the seismic data signals are processed using plane wave decomposition interpolation, the interpolation written as a least squares equation:

wherein, V_nIs the reconstructed seismic data signal of the nth trace, P_n,n-1Is the dip angle domain plane wave decomposition operator between the nth track and the (n-1) th track, and is determined by the formula:

and solving the least square equation.

9. The method for reconstructing a bidirectional wavefield of a time-lapse seismic virtual source as claimed in claim 8, wherein the plane wave decomposition interpolation method is used to obtain a reconstructed wavefield after interpolation, and the reconstructed wavefield is:

wherein,

is smoothed to form an expanded tilt angle field

The length of (a) of (b),

is that

And reconstructing data after the strip difference value.

10. A time-lapse seismic virtual source bi-directional wavefield reconstruction system, comprising:

the bidirectional wave field reconstruction module is used for converting the pre-acquired seismic data signals into a frequency domain, reconstructing wave fields of a seismic source end and a receiving end and reconstructing a seismic source and a receiver below the earth surface;

and the common-center-point superposition module is used for carrying out common-center-point superposition on the seismic data signals of the time domain to obtain a virtual seismic source excitation wave field.