CN107272058B

CN107272058B - Imaging method, imaging apparatus, and computer storage medium

Info

Publication number: CN107272058B
Application number: CN201710542746.9A
Authority: CN
Inventors: 王之洋; 胡婷; 袁雨欣; 刘洪�
Original assignee: Institute of Geology and Geophysics of CAS
Current assignee: Institute of Geology and Geophysics of CAS
Priority date: 2017-07-05
Filing date: 2017-07-05
Publication date: 2020-02-11
Anticipated expiration: 2037-07-05
Also published as: CN107272058A

Abstract

The application discloses an imaging method, an imaging device and a computer storage medium. The method comprises the following steps: constructing a forward wave field of an underground shot point and a backward wave field of a wave detection point; judging the type of the medium for transmitting the forward transmission wave field and the backward transmission wave field; when the medium is an isotropic medium, performing longitudinal and transverse wave field separation on the forward wave field and/or the backward wave field according to a preset equation decoupling algorithm, and when the medium is non-uniform anisotropy, performing longitudinal and transverse wave field separation on the forward wave field and/or the backward wave field according to a preset variation function IDW interpolation algorithm; decomposing the separated longitudinal wave field and/or transverse wave field in four directions, namely up, down, left and right, to obtain a plurality of direction wave fields of the longitudinal wave field and/or transverse wave field; at least two directional wave fields are selected from the plurality of directional wave fields and subjected to normalized cross-correlation vector imaging. The imaging method can effectively suppress coupling noise, low-frequency noise and imaging false noise generated in the reverse time migration process of the elastic wave, and improve the imaging quality.

Description

Imaging method, imaging apparatus, and computer storage medium

Technical Field

The present application relates to the field of seismic exploration, and in particular, to an imaging method, an imaging apparatus, and a computer storage medium.

Background

The elastic wave seismic exploration technology has the advantages of high medium parameter obtaining precision, high certainty of medium parameter inversion, high imaging precision, high reservoir prediction precision, high fluid identification precision and the like. Since the elastic wave equation contains more information, it can describe the composition of the earth's medium more accurately than the acoustic wave equation. Therefore, the research on the elastic wave reverse time prestack depth migration technology is crucial to the research on improving the seismic imaging quality and solving some new problems of seismic exploration.

The problems of coupling noise, low-frequency noise and imaging false noise of the elastic wave imaging result become hot spots and difficulties concerned by the geophysical exploration society at home and abroad.

Elastic wave imaging noise mainly comes from coupling noise, low-frequency noise and imaging false noise, how to suppress the noise, for many years, scholars at home and abroad also propose various technical means to remove the noise, such as decoupling propagation, up-down direction wave decomposition and other technical methods, and although a certain effect is achieved, the improvement of the elastic wave imaging quality is not much.

Zhang et al (2007) deduces a decoupled elastic wave first-order velocity stress equation to realize separation of longitudinal and transverse waves. Longitudinal wave stress is introduced by Xiao and Leaney (2010), and a decoupling equation set based on an elastic wave first-order velocity stress equation is provided, is mathematically equivalent to the decoupling equation set provided by Zhang et al (2007), and realizes longitudinal and transverse wave decoupling propagation and imaging of VSP data. Zhou et al (2016) extend the set of decoupling equations from 2 dimensions to 3 dimensions based on the work of Zhang et al (2007). Fei et al (2010) propose that up-going and down-going wavefields can be separated and then imaging conditions applied separately to wavefields of different propagation directions, combined for imaging.

Fei et al (2015) proposed the concept of imaging artifacts, and they thought that Liu (2011) could suppress low frequency noise well, but neglected the imaging artifacts caused by backscattering when the velocity changes dramatically, i.e. the artifacts caused by cross-correlation imaging of the shot forward wavefield and the demodulator backward wavefield, called primary artifacts (primary RTM), and proposed de-primary RTM algorithm. Wangyibo et al (2016) derived the separation of the waves in the up, down, left, and right 4 quadrants of the scalar wave equation based on the Fei et al (2015) work.

The method has certain effect on suppressing low-frequency noise, coupling noise and imaging false noise, but the decoupling equation separation longitudinal and transverse waves can only be suitable for isotropic media and cannot be applied to anisotropic media, the up-down directional wave decomposition can well suppress the low-frequency noise generated by correlation of directional waves propagating in the same direction, and suppress the false noise caused by backward cross-correlation of a shot forward wave field and a wave detection point to a certain extent, but if the inclination angle of an underground reflecting layer is large and the target position is complex in distribution, more imaging noise from the left-right directional waves can still be generated, so that the existing method has certain technical limitation.

Disclosure of Invention

An object of the present invention is to provide an imaging method, an imaging apparatus, and a computer storage medium to solve the above-mentioned problem of poor imaging quality.

In order to achieve the above object, the present invention provides an imaging method, an imaging apparatus, and a computer storage medium, the imaging method including:

constructing a forward wave field of an underground shot point and a backward wave field of a wave detection point;

judging the type of a medium for transmitting the forward transmission wave field and the backward transmission wave field;

when the medium is an isotropic medium, performing vertical and horizontal wave field separation on the forward transmission wave field and/or the backward transmission wave field according to a preset equation decoupling algorithm, and when the medium is non-uniform anisotropy, performing vertical and horizontal wave field separation on the forward transmission wave field and/or the backward transmission wave field according to a preset mutation function IDW interpolation algorithm;

decomposing the separated longitudinal wave field and/or transverse wave field in four directions, namely up, down, left and right, to obtain a plurality of direction wave fields of the longitudinal wave field and/or the transverse wave field;

and selecting at least two directional wave fields from the plurality of directional wave fields, and performing normalized cross-correlation vector imaging on the selected directional wave fields.

In an alternative embodiment, the decomposing the separated compressional wave field and/or the separated shear wave field in four directions includes:

carrying out up-down directional wave separation on the separated longitudinal wave field and transverse wave field to obtain corresponding up-going directional wave field and down-going directional wave field;

and performing left-right directional wave decomposition on the uplink wave field and the downlink wave field to obtain corresponding left uplink wave field and left downlink wave field, and right uplink wave field and right downlink wave field.

In an optional embodiment, the performing up-down directional wave decomposition on the separated longitudinal wave field includes:

carrying out time complex field expansion on the separated longitudinal wave field;

carrying out Fourier transform on the longitudinal wave field after the time complex field expansion in the Z direction to obtain a wave number field longitudinal wave field;

selecting a part with wave number larger than zero from the longitudinal wave field of the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the Z direction to obtain a downlink direction wave field of the longitudinal wave field; and selecting a part with wave number less than zero from the longitudinal wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the Z direction to obtain an uplink wave field of the longitudinal wave field.

In an optional embodiment, the performing the directional wave decomposition in the left-right direction on the wave field of the upward direction includes:

performing time complex domain expansion on an uplink wave field of the longitudinal wave field;

performing Fourier transform on an uplink wave field of the longitudinal wave field after the time complex field expansion in the X direction to obtain a wave number field longitudinal wave field;

selecting a part with wave number larger than zero from an uplink wave field of the longitudinal wave field in the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a right uplink wave field of the longitudinal wave field; and selecting a part with wave number less than zero from the longitudinal wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a left uplink direction wave field of the longitudinal wave field.

In an optional embodiment, the performing the directional wave decomposition in the left-right direction on the wave field of the downlink direction includes:

performing time complex domain expansion on a wave field in the downstream direction of the longitudinal wave field;

performing Fourier transform on a downstream direction wave field of the longitudinal wave field after the time complex field expansion in the X direction to obtain a wave number field longitudinal wave field;

selecting a part with wave number larger than zero from a wave field in the downward direction of the longitudinal wave field in the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a wave field in the right downward direction of the longitudinal wave field; and selecting a part with wave number less than zero from the longitudinal wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a left downlink direction wave field of the longitudinal wave field.

In an alternative embodiment, the performing up-down directional wave decomposition on the separated shear wave field includes:

carrying out time complex field expansion on the separated shear wave field;

performing Fourier transform on the transverse wave field after the time complex field expansion in the Z direction to obtain a wave number field transverse wave field;

selecting a part with wave number larger than zero from the wave number domain transverse wave field, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the Z direction to obtain a downstream direction wave field of the transverse wave field; and selecting a part with wave number less than zero from the wave number domain transverse wave field, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the Z direction to obtain an uplink direction wave field of the transverse wave field.

performing time complex domain expansion on an uplink wave field of the transverse wave field;

performing Fourier transform on an uplink wave field of the transverse wave field after the time complex field expansion in the X direction to obtain a wave number field transverse wave field;

selecting a part with wave number larger than zero from an uplink wave field of the transverse wave field in the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a right uplink wave field of the transverse wave field; and selecting a part with wave number less than zero from the transverse wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a wave field of the left uplink direction of the transverse wave field.

performing time complex domain expansion on a wave field in the downlink direction of the transverse wave field;

performing Fourier transform on a wave field in a descending direction of the transverse wave field after the time complex field expansion in the X direction to obtain a wave number field transverse wave field;

selecting a part with wave number larger than zero from a wave field in the downward direction of the transverse wave field in the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a wave field in the right downward direction of the transverse wave field; and selecting a part with wave number less than zero from the transverse wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a left downlink wave field of the transverse wave field.

In an optional embodiment, when the medium is non-uniform and anisotropic, performing vertical and horizontal wavefield separation on the forward wavefield and/or the backward wavefield according to a preset variance function IDW interpolation algorithm, including:

when the wave field is in a non-uniform anisotropic medium, the forward transmission wave field and/or the backward transmission wave field at each moment are respectively transformed from a space domain to a wavenumber domain by utilizing Fourier transform;

selecting N reference models, calculating a separation operator according to the reference models, and performing truncation optimization on the separation operator by using a self-convolution combination window function in a wavenumber domain;

separating the forward wave field and/or the backward wave field under the reference model by using the separation operator in a wave number domain to obtain the separated longitudinal and transverse wave fields under the reference model;

transforming the separated longitudinal and transverse wave fields under the reference model from a wave number domain to a space domain, wherein N reference models correspondingly obtain N separated longitudinal and transverse wave fields under the reference model, and N is a positive integer;

in the spatial domain, an IDW interpolation algorithm is improved by using a variation function, the weight coefficients of N reference models are calculated according to the improved IDW interpolation algorithm, and the weighting interpolation processing is carried out on the separated longitudinal wave field and the separated transverse wave field under the N reference models to obtain the final separated longitudinal wave field and transverse wave field.

In an optional embodiment, the selecting N reference models includes:

traversing a plurality of anisotropic parameters of the initial model to be numerically simulated, constraining by using a critical value of a variation function, selecting a reference point with the maximum probability of occurrence in the critical value range, and searching according to the reference point to obtain the reference model.

In an alternative embodiment, the weight coefficients of the reference model are calculated by the variance function.

In an optional embodiment, the method for obtaining the variation function includes:

and fitting to obtain the variation function of the initial model by calculating the variation function value of the anisotropic parameter distribution of the initial model.

In an alternative embodiment, the vector imaging is elastic wave reverse time migration imaging.

The present application also provides an image forming apparatus including:

the construction unit is used for constructing an underground shot forward transmission wave field and a wave detection point backward transmission wave field;

the separation unit is used for judging the type of a medium for transmitting the forward transmission wave field and the backward transmission wave field; when the medium is an isotropic medium, performing vertical and horizontal wave field separation on the forward transmission wave field and/or the backward transmission wave field according to a preset equation decoupling algorithm, and when the medium is non-uniform anisotropy, performing vertical and horizontal wave field separation on the forward transmission wave field and/or the backward transmission wave field according to a preset mutation function IDW interpolation algorithm; decomposing the separated longitudinal wave field and/or transverse wave field in four directions, namely up, down, left and right, to obtain a plurality of direction wave fields of the longitudinal wave field and/or the transverse wave field;

and the imaging unit is used for selecting at least two directional wave fields from the plurality of directional wave fields and carrying out normalized cross-correlation vector imaging on the selected directional wave fields.

The present application further provides a computer storage medium having stored thereon a computer program which, when executed by a processor, performs the steps of:

The method utilizes the separation technology of an elastic wave vector wave field, namely a decoupling equation method is adopted in an isotropic medium, an improved IDW interpolation algorithm based on a variation function is adopted in a non-uniform anisotropic medium to realize the separation of longitudinal and transverse waves, then a Fourier transform method is applied to realize the decomposition of longitudinal and transverse wave omnibearing directional waves, at least two directional wave fields are selected from a plurality of directional wave fields on the basis of realizing the separation of the longitudinal and transverse waves and the decomposition of the omnibearing directional waves, and the selected directional wave fields are subjected to vector imaging, thereby solving the problem that the prior art can only decompose and image up, down, left and right directional waves on a standard wave and cannot be applied to vector waves (elastic waves). Meanwhile, the imaging method can effectively suppress coupling noise, low-frequency noise, imaging false noise and the like generated in the reverse time migration process of the elastic wave, and improves the imaging quality.

Drawings

In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings needed to be used in the description of the embodiments or the prior art will be briefly introduced below, it is obvious that the drawings in the following description are only some embodiments described in the present application, and for those skilled in the art, other drawings can be obtained according to the drawings without any creative effort.

Fig. 1 is a flow chart of an imaging method provided by an embodiment of the present application;

fig. 2 is a flowchart illustrating that a variogram IDW difference algorithm is adopted in an imaging method to separate the acquired elastic waves into longitudinal waves and transverse waves according to an embodiment of the present disclosure;

fig. 3 is a flowchart illustrating how to obtain a self-folding combination window function in an imaging method using a variogram IDW difference algorithm according to an embodiment of the present disclosure;

fig. 4 is a flowchart illustrating a longitudinal wave being separated into an upstream longitudinal wave and a downstream longitudinal wave in an imaging method according to an embodiment of the present application;

fig. 5 is a flowchart illustrating separation of the uplink shear wave into a left uplink shear wave and a right uplink shear wave in an imaging method according to an embodiment of the present application;

fig. 6 is a flowchart illustrating a separation of a downlink shear wave into a left downlink shear wave and a right downlink shear wave in an imaging method according to an embodiment of the present application;

fig. 7 is a flowchart illustrating a separation of shear waves into an uplink shear wave and a downlink shear wave in an imaging method according to an embodiment of the present application;

fig. 8 is a flowchart illustrating separation of the uplink shear wave into a left uplink shear wave and a right uplink shear wave in an imaging method according to an embodiment of the present application;

fig. 9 is a flowchart illustrating a separation of a downlink shear wave into a left downlink shear wave and a right downlink shear wave in an imaging method according to an embodiment of the present application;

FIG. 10(a) shows the result of Marmousi2 model elastic wave reverse time migration local imaging;

FIG. 10(b) shows the result of Marmousi2 model elastic wave reverse time migration local imaging;

FIG. 11(a) shows the result of PP imaging of the Marmousi2 model;

FIG. 11(b) shows the result of PP imaging of the salt dome model;

FIG. 11(c) is the PS imaging results of the salt dome model;

FIG. 12a is a Marmousi2 model anisotropic parameter distribution and a selected reference model;

FIG. 12b is the distribution of the anisotropic parameters of the Marmousi2 model and a selected reference model;

FIG. 12c is the Marmousi2 model anisotropic parameter distribution and the selected reference model;

FIG. 12d is the Marmousi2 model anisotropic parameter distribution and the selected reference model;

FIG. 13 is a diagram illustrating local extrema in a reference model selected during an interpolation operation of the IDW algorithm;

FIG. 14 is a diagram illustrating a reference model selected during an interpolation operation of an IDW algorithm according to an embodiment of the present disclosure;

fig. 15 is a schematic structural diagram of an elastic wave reverse time migration imaging device according to an embodiment of the present application.

Detailed Description

In order to make those skilled in the art better understand the technical solutions in the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application, and it is obvious that the described embodiments are only a part of the embodiments of the present application, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present application.

In the field of reverse time migration imaging, researchers have made many efforts to avoid interference from coupling noise, low frequency noise, imaging artifact noise, etc. during reverse migration imaging. For example, Fei et al (2015) proposes a de-primary RTM algorithm based on which an elastic wave is decomposed into up and down waves, and then noise in the decomposed waves is removed. However, the separation of up-going waves and down-going waves is hidden in imaging conditions, and although the separation of the up-going waves and the down-going waves can be realized, separated wave fields cannot be obtained.

Wang et al (2016) propose that combining longitudinal and transverse wave separation with up-down wave decomposition can effectively suppress coupling noise, low-frequency noise and imaging artifact noise, and explicitly separate up-down waves; however, the algorithm of Zhang and McMechan (2010) is too computationally intensive to use in anisotropic media. This method is not applicable to three-dimensional situations (in which anisotropy of the medium is to be considered), and therefore its application is limited to geological exploration, which requires consideration of three-dimensions. Moreover, Wang et al (2016) do not consider the decomposition of waves in the left-right direction, but only decompose waves in the up-down direction, and have a certain limitation in the presence of a vertical structure.

Wangyibo et al (2016) derived the separation of waves in the up, down, left, and right 4 quadrants of the scalar wave equation based on the Fei et al (2015) work, but this method is suitable for scalar waves (acoustic waves) and not for vector wave fields (elastic waves), and therefore cannot be applied to elastic wave reverse time migration imaging.

In view of the above, different researchers have tried to solve the noise problem in the imaging process from different aspects, but the various methods have not been effective or can not be used under specific conditions.

In order to solve the problem of noise in the imaging process, the application provides an imaging method. Fig. 1 is a schematic flow chart of an imaging method provided by an embodiment of the present invention, and as shown in fig. 1, the method may include the following steps:

s101: constructing a forward wave field of an underground shot point and a backward wave field of a wave detection point;

s102: judging the medium types of the forward transmission wave field and the backward transmission wave field;

s103: when the medium is an isotropic medium, performing longitudinal and transverse wave field separation on the forward wave field and/or the backward wave field according to a preset equation decoupling algorithm; when the medium is non-uniform anisotropy, carrying out longitudinal and transverse wave field separation on a forward wave field and/or a backward wave field according to a preset variable function IDW interpolation algorithm;

s104: decomposing the separated longitudinal wave field and the transverse wave field in four directions, namely, up, down, left and right, to obtain a plurality of direction wave fields of the longitudinal wave field and/or the transverse wave field;

s105: and selecting at least two directional wave fields from the plurality of directional wave fields, and carrying out vector imaging on the selected directional wave fields.

The method comprises the steps of applying a decoupling equation method in an isotropic medium and applying an IDW (inverse discrete wavelet) interpolation algorithm based on a variation function to separate a forward transmission wave field from a backward transmission wave field in a non-uniform anisotropic medium.

An embodiment of the present invention provides an imaging method. The method utilizes the separation technology of an elastic wave vector wave field, namely a decoupling equation method is adopted in an isotropic medium, and an improved IDW interpolation algorithm based on a variation function is adopted in a non-uniform anisotropic medium, so as to realize the separation of longitudinal waves and transverse waves; and then, a complex domain expansion method is applied to realize the decomposition of longitudinal and transverse waves in all directions, at least two direction wave fields are selected from a plurality of direction wave fields after the longitudinal and transverse wave separation and the all-direction wave decomposition are realized, and vector imaging is carried out on the selected direction wave fields, so that the problem that the longitudinal and transverse waves, the vertical and left and right direction waves can only be separated on a standard wave and cannot be applied to vector waves (elastic waves) in the prior art is solved. Meanwhile, the imaging method can effectively suppress coupling noise, low-frequency noise, imaging false noise and the like generated in the reverse time migration process of the elastic wave, and improves the imaging quality.

Specifically, IDW (inverse Distance weighted) is a common and simple spatial interpolation method. The weighted average is carried out by taking the distance between the interpolation point and the sample point as the weight, and the weight is given to the sample point which is closer to the interpolation point. For example, let a series of discrete points be distributed on the plane, and known as Xi, Yi, Zi (i ═ 1,2, …, n), the z-point value is determined by distance weighted values.

Specifically, taking the Marmousi2 model as an example, fig. 12a to 12d and fig. 13 are schematic diagrams of the distribution of anisotropic parameters of the Marmousi2 model and the selected reference model. From the above figure, the spatial distribution of the anisotropic parameters of the Marmousi2 model is not uniform. In the actual calculation process, the isotropic medium does not need to take into account the distance between two points. Because each point is the same for isotropic media, but different for anisotropic media. There is a certain distance between two points, which represents a similarity, assuming that closer points are more similar. Therefore, we need to interpolate a point, which can be obtained by weighted interpolation of the points near it, and the weights of the points near are larger. This method is an application of the inverse distance Interpolation (IDW) algorithm, but the algorithm does not take into account spatial variability. Therefore, the method also considers the space variability and adopts the IDW interpolation algorithm of the variation function based on the variability, thereby improving the separation quality when the vertical wavelength and the horizontal wavelength are separated, and further improving the imaging quality.

The variation function is a distance function between two points, not only describes the spatial distance, but also describes the spatial variability of the two points, and the interpolation calculation adopting the variation function can be more accurate. The specific solution process of the mutation function will be described below.

Referring to fig. 2, the step S103 further includes the following substeps:

s201: when the elastic vector wave field is in a non-uniform anisotropic medium, the elastic vector wave field at each moment is respectively transformed from a space domain to a wave number domain by utilizing Fourier transform;

s202: selecting N reference models, calculating a separation operator according to the reference models, and performing truncation optimization on the separation operator by using a self-convolution combination window function in a wavenumber domain;

s203: separating the elastic wave field under the reference model by using the separation operator in a wave number domain to obtain the separated longitudinal and transverse wave field under the reference model;

s204: transforming the separated longitudinal and transverse wave fields under the reference model from a wave number domain to a space domain, wherein N reference models correspondingly obtain N separated longitudinal and transverse wave fields under the reference model, and N is a positive integer;

s205: in the space domain, the weight coefficients of N reference models are calculated by using a variogram interpolation theory, and the separated longitudinal wave field and the separated transverse wave field under the N reference models are subjected to weighted interpolation processing to obtain the finally separated longitudinal wave field and the finally separated transverse wave field.

In the embodiment of the invention, the wave field is separated by using the separation operator in the wave number domain, and the finally separated wave field is obtained by applying IDW algorithm interpolation operation based on the mutation function in the space domain, so that the calculation complexity of the wave field separation can be effectively reduced. In the embodiment of the invention, the self-convolution combination window function is adopted in the wavenumber domain to cut off the separation operator and then separate the wave field, so that the efficiency of wave field separation can be improved, and the higher accuracy of the cut-off separation operator can be ensured.

Specifically, the reference model may be obtained by traversing a plurality of anisotropic parameters of the initial model to be numerically simulated, constraining by using the threshold R of the variation function, selecting the reference point with the highest occurrence probability within the threshold range, and searching according to the reference point (see fig. 13 and 14).

FIG. 3 is a flowchart illustrating a method for obtaining a self-convolution combining window function according to an embodiment of the present invention. As shown in fig. 3, the following steps may be included:

s301: selecting a window function with main lobe and side lobe performance higher than a set threshold value as an original window function;

s302: performing self-convolution operation on the original window function for L times to obtain a self-convolved window function, wherein L is a positive integer;

s303: and carrying out weighting operation on the self-convolution window function and the original window function to obtain the self-convolution combined window function.

In the above embodiment, the self-convolution operation is performed on the selected original window function, the self-convolution operation makes the main lobe performance of the original window function worse and the side lobe performance better, and then the weighting operation is performed on the original window function and the self-convolved window function according to the current requirement, so as to obtain a window function meeting the requirement. The precision and stability of the wave field separation can be improved by utilizing the separation operator of the self-convolution combination window function.

In step S302, the weighting coefficient for performing the weighting operation may have the following characteristics: under the condition that the main lobe performance is determined to be required to be prior, the weighting coefficient of the original window function is higher than that of the self-convolved window function; in case it is determined that side lobe performance precedence is required, the weighting coefficients of the self-convolved function are higher than the weighting coefficients of the original window function. Thus, the side lobe attenuation can be increased as much as possible while maintaining an appropriate main lobe width.

In order to make the technical solution of the present invention more easily understood by those skilled in the art, the variant function IDW difference algorithm and the sub-steps of the step will be described in a specific embodiment.

For non-uniform anisotropy, theoretically, each sampling point needs to calculate a polarization vector according to a corresponding anisotropy parameter to obtain a better separation result, but the processing calculation amount is too large, so that N reference models can be selected from an initial model, vector wave field (elastic wave) separation is respectively carried out in a wave number domain, then inverse transformation is carried out back to a space domain, and the separation result of the initial model is reconstructed according to a weight coefficient calculated in advance.

Therefore, the weight coefficient is calculated instead of the distance between two points in consideration of introducing a variation function describing the spatial structural variation and the stochastic variation of the regionalized variable. Specifically, the value of the variation function of the initial model anisotropic parameter distribution can be calculated according to the following formula:

wherein μ (h) represents a value of a variation function, Φ ═ f [ δ +2 (e- δ) sin ²(α-θ)]sin 2(α-θ)，

Epsilon and delta represent coefficients of anisotropy of the TI medium, theta represents a dip angle of a symmetric axis of the TTI medium,

the position increment representing the interpolated reference point, i-1,n, where n represents the number of interpolation reference points.

After a plurality of scattered variation function values are obtained through calculation, fitting is carried out according to the plurality of variation function values, and a theoretical variation function model is obtained.

Then, the weight coefficients of the respective reference models are calculated according to the following formula:

wherein, w ^kRepresents a weight coefficient of (ε) _k,δ _k,θ _k) Indicating the interpolated reference point.

As described above

Indicating position increments, for each

N may be calculated as 1,2 By different

The mode (a) is the abscissa and,

for the ordinate, a dot diagram is drawn, from which it can be found that: distance in space

The smaller, the smaller the variability,

the larger the size, the greater the variability when

When the value is increased to a critical value R, the two values basically have no great influence on each otherThis indicates that when selecting a reference point, or a point to be fitted, it should be selected within the R range, otherwise the selected point has not much practical significance. For example: the method comprises the steps of calculating to obtain a plurality of scattered variation function values according to an initial reference model, fitting the obtained plurality of scattered variation function values to obtain a theoretical variation function model, determining parameter values of all parameters in the variation function according to the fitted variation function model, traversing the probability of occurrence of epsilon, delta and theta of different values in the initial model, and selecting a point with the maximum probability of occurrence within a critical value of the fitted variation function as an interpolation reference point.

The following description is provided with reference to a specific embodiment, but it should be noted that the specific embodiment is only for better illustration of the present invention and should not be construed as a limitation of the present invention.

An inevitable problem in the separation of vector wavefields based on polarization characteristics is how to project wavefields more accurately in the corresponding polarization directions while ensuring the computational efficiency. In the wavenumber domain, this is called projection, and in the spatial domain, this is called spatial filtering.

In an isotropic medium, the separation effect of a vector wave field is only related to the precision of a differential operator, and the higher the precision is, the better the separation effect is, and meanwhile, the calculation amount is increased. However, in anisotropic media, the separation effect of the vector wave field is mainly related to two factors, namely anisotropic parameters and the precision of a quasi-differential operator, and considering that in TI media, the seismic wave dynamics can be expressed as the superposition of an isotropic part and an anisotropic part, a polarization vector can be obtained through calculation of a christofel equation and then is decomposed into the sum of a wave number vector and a deflection vector.

The deflection vector is a function with wave number as an independent variable and has space distribution characteristics, a conventional binomial window function and a Gauss window function are utilized in a wave number domain to cut off and approximate to a differential operator, and an IDW algorithm is utilized in a space domain to interpolate the anisotropic part of the polarization vector, so that the calculated amount of wave field separation is reduced, and the efficiency and the precision of the wave field separation are improved.

Specifically, there are two main factors affecting the technique of separating the vector wavefield based on polarization characteristics: one is to cut off the amplitude-frequency response characteristic of the window function, the main lobe shape of the amplitude-frequency response of the window function controls the range of the transition band, namely the spectrum coverage range, the shape of the side lobe determines the deviation degree of the differential operator approaching differential operator, and the performance of the main lobe and the side lobe directly influences the approaching precision; the other is an interpolation effect of an interpolation algorithm, and for an interpolation algorithm, not only higher interpolation accuracy is required, but also a smaller amount of calculation is required, and the execution efficiency of the IDW algorithm is higher, but only the calculation weight of the distance between two points is considered, so that the interpolation accuracy is not very high, and therefore, an interpolation algorithm capable of improving the interpolation accuracy while ensuring the execution efficiency is required.

In this example, a vector wave field separation method is proposed, which mainly solves the following two problems: the selection of the truncation window function is the selection of the interpolation algorithm.

Firstly, because the performance of the main lobe and the side lobe of the amplitude response of the window function directly affects the precision problem of the differential approximation pseudo-differential operator, if a proper window function is designed, firstly, how the performance of the main lobe and the side lobe affects the approximation precision is understood, and secondly, how to control the main lobe and the side lobe is researched. Then, the choice of interpolation algorithm introduces the optimization of spatial interpolation algorithm into the vector wave field separation technique, since for non-uniform anisotropic media it is better to separate the vector wave field if one quasi-differential operator is calculated for each point, but this will consume a huge amount of computation, assuming that the size of the quasi-differential operator is n ²The size of the model is N ²The calculated amount of vector wave field separation for the model is then 2n ²N ²This is much greater than 2mN for m order precision finite difference method ²Calculated amount of (1), n<And m is selected. Therefore, an interpolation algorithm is required to reconstruct the separation result of the initial model in a spatial domain by interpolation.

In this example, the basic principle underlying is: in an isotropic medium, the Helmholtz theorem is applied to respectively solve the divergence and the rotation of a wave field so as to separate longitudinal and transverse waves, and the formula is as follows:

P＝▽·U,

S＝▽×U。

in the wavenumber domain can be expressed as:

P＝▽·U＝ik _xU _x+ik _zU _z；

S＝▽×U＝ik _zU _x-ik _xU _z。

as can be seen from the above formula, in an isotropic medium, P-wave is the projection of the vector wave field in the wave number direction, S-wave is the projection of the vector wave field in the direction perpendicular to the wave number direction, and the polarization vector of P-wave is (k) _x,k _z) The polarization vector of the S wave is orthogonal thereto and thus is (k) _z,-k _x)。

For an anisotropic medium, the polarization vector of a longitudinal wave and a transverse wave in the anisotropic medium can also be obtained through a christofsel equation corresponding to the anisotropic medium, taking a two-dimensional TTI medium as an example, the christofsel equation corresponding to the TTI medium is as follows:

wherein:

wherein, c ₁₁,c ₁₅,c ₃₃,c ₃₅,c ₅₅Representing the tensor of the elastic coefficient, k _x,k _zExpressing the normalized wavenumber, Γ, the christoflel matrix, which is a typical eigenvalue and eigenvector problem, to have a non-zero solution to the formula, the coefficient determinant is zero,thus, the polarization vector P ═ (P) can be determined _x,P _z) This is k _x,k _zThe inverse transformation of iP back to the spatial domain can result in the pseudo-differential operator L, also called spatial filter operator.

Thomsen coefficients V characterizing the anisotropy of TI media (transversely isotropic media) are now introduced _P0,V _S0ε, δ, γ, wherein V _P0,V _S0ε, δ relate to the qP and qSV waves, V _S0γ is related to qSH waves because qSH waves are decoupled, so the separation of the vector wavefield only requires V _P0,V _S0Four parameters, ε and δ. Taking a two-dimensional TTI medium (a transversely isotropic medium with an inclined axis of symmetry) as an example, the tilt angle θ of the axis of symmetry of the TI medium is also introduced. In TI media, seismic wave dynamics are composed of two parts, the first part being an isotropic part and the second part being an anisotropic part, which can be expressed approximately as:

K≈K ^iso+L(ε,δ)+Q(ε,δ,V _S0) (4)

wherein K represents the dynamic characteristics of seismic waves, K ^isoRepresents an isotropic moiety, epsilon-delta-0; l (ε, δ) + Q (ε, δ, V) _S0) Denotes the anisotropic part, L (. epsilon., delta.) denotes the linear part, Q (. epsilon., delta., V) _S0) Representing the non-linear part. The same applies to the above approximation formula for the polarization vector, i.e., the polarization vectors of qP and qSV waves in the TI medium are also composed of two parts: an isotropic portion and an anisotropic portion.

For the isotropic part, a conventional optimization method of finite difference operator can be adopted, namely: the window function truncation method or optimization method is used, and the two methods are similar in nature and are expected to achieve a good approximation accuracy in a higher wave number range.

For the anisotropic part, theoretically, each sampling point needs to calculate a polarization vector according to a corresponding anisotropic parameter to obtain a better separation effect, but this inevitably brings a huge amount of calculation, and therefore, an interpolation mode is considered, that is: selecting N reference models from the initial model, respectively carrying out vector wave field separation in a wave number domain, then carrying out inverse transformation back to a space domain, and reconstructing a separation result of the initial model according to a weight coefficient calculated in advance. Since the interpolation algorithm, the interpolation precision and the execution efficiency of the interpolation algorithm are all factors to be considered, the IDW interpolation only considers the distance between two points to calculate the weight, the interpolation precision is poor, and the IDW interpolation only has good applicability to the uniform model, so that the introduction of a variation function describing the spatial structural change and the stochastic change of the regionalized variable can be considered to replace the distance between two points to calculate the weight coefficient.

The above is specifically described as follows:

1. isotropic part:

the derivative of a band-limited continuous signal f (x) at x-0 can be represented by a uniformly sampled signal f _nExpressed as:

there is a window function with length N +1 point, N is even number, to truncate the above two formulas, and through some simple processing, we can get finite difference operator:

wherein:

to optimize the precision of the finite difference operator approximation differential operator, the Chebyshev window function can be chosen to truncate:

where r represents the ripple rate, representing the degree of attenuation of the side lobes, N +1 represents the length of the window, and N is an even number.

The performance of the main lobe and the side lobe of the amplitude-frequency response of different window functions influences the approximation precision of the difference, and the specific influence mode comprises the following steps:

1) the main lobe size is related to the transition bandwidth: the main lobe is narrow, then the transition band is narrow, the spectrum coverage of the precision error of the finite difference operator using the window function to cut off approximation is large, the high-order precision can be achieved by using a low-order operator, the main lobe is wide, then the transition bandwidth is wide, and the spectrum coverage of the precision error of the finite difference operator using the window function to cut off approximation is large.

2) Relation of sidelobe attenuation and approximation precision stability: the attenuation of the side lobe directly influences the stability of the precision error of the finite difference operator of the window function truncation approximation, the larger the attenuation of the side lobe is, the smaller the fluctuation of the approximation precision error is, the high stability is realized, the smaller the attenuation of the side lobe is, the larger the fluctuation of the approximation precision error is, and the low stability is realized.

2. Anisotropic moiety:

in the TI anisotropic medium, the anisotropic part of the polarization vector of the seismic wave is composed of a linear part and a nonlinear part, and for the TTI medium, the nonlinear part is removed, that is, under the condition of weak anisotropy, the polarization angle of the qP wave can be approximately expressed as:

v _p＝α+f[δ+2(ε-δ)sin ²(α-θ)]sin 2(α-θ) (12)

where α denotes the phase angle, the angle of the propagation direction to the Z-axis, representing the isotropic part, theta denotes the tilt of the symmetry axis of the TTI medium,

removing the isotropic part in the above approximate formula, and only retaining the linear anisotropic part, we can get:

Φ＝f[δ+2(ε-δ)sin ²(α-θ)]sin 2(α-θ) (13)

wherein phi is approximately linear with epsilon and delta, and is approximately proportional with sin2(α -theta).

Therefore, it is possible to set the initial model to m ═ epsilon, δ, θ, and the anisotropy partial value of the polarization angle under the initial model condition to Φ, and further, it is possible to conditionally select N reference models m ^kAnd calculating to obtain phi corresponding to each reference model ^kAnd obtaining phi through interpolation according to the following formula:

wherein, w ^kRepresenting the weight coefficients.

The IDW interpolation algorithm only considers the spatial position relationship between the interpolation point and the reference point when calculating the weight, that is, when the spatial positions of different reference points and the interpolation point are consistent, the same weight coefficient exists, and although the execution efficiency of the IDW interpolation algorithm is high, the difference between the reference points is not considered.

In this example, the determination of the weight coefficient is improved, and specifically includes: assuming that phi satisfies a second-order stationary assumption or an intrinsic assumption, a variation function mu (h) is used to describe the characteristic that phi varies with different spatial positions, and the variation function value is calculated according to the following formula:

wherein the content of the first and second substances,

indicating position increments, for each N may be calculated as 1,2 By different

The mode (a) is the abscissa and,

The smaller, the smaller the variability,

the larger the size, the greater the variability when

Increasing to a critical value R, the two values will not substantially affect each other, which means that when selecting a reference point, or a point to be fitted, it should be selected within the range of R, otherwise the selected point has no much practical significance. It is described here that scattered points are obtained, and some theoretical variogram models can be used to fit these points to obtain the relevant parameters.

The weight coefficient w can be obtained after introducing the variation function R ^k：

Because μ (h) is a function of distance, more reliable interpolation results can be obtained by introducing μ (h) into the weight calculation of the IDW interpolation algorithm.

For the interpolated reference point (ε) _k,δ _k,θ _k) The method for selecting the interpolation reference point may include: calculating to obtain a plurality of scattered variation function values according to the initial model; fitting the obtained multiple scattered variation function values to obtain a theoretical variation function model; determining parameter values of all parameters in the variation function according to the variation function model; traversing the probability of the appearance of epsilon, delta and theta of different values in the numerical simulation initial model to be tested, and selecting the point with the maximum appearance probability within the fitted variation function critical value as an interpolation reference point.

Of course, other embodiments of the present invention provide other types of separation methods in addition to the above described longitudinal and transverse wave separation algorithm. For example, wavefields may be separated based on velocity differences or wave equation methods may be employed, etc. Therefore, the present application does not limit the algorithm for separating longitudinal and transverse waves, and the method for separating coupled longitudinal and transverse waves is all in accordance with the requirements of the present application.

Referring to fig. 1, step S104 of the method of the present application: in the vertical directional wave decomposition of the longitudinal and transverse wave fields after the longitudinal and transverse wave separation, referring to fig. 4, the vertical directional wave decomposition of the longitudinal wave field after the longitudinal and transverse wave separation includes the following sub-steps:

s401: carrying out time complex field expansion on the separated longitudinal wave field;

the complex field expansion strategy is a method for converting data into a complex form and performing various calculations, and the strategy is based on Hilbert transformation and can effectively improve the calculation efficiency. As can be seen from the spectral analysis, after Hilbert transform, the amplitude of the wave field in the frequency domain remains unchanged, and the phase is inverted. For complex domain extended wavefields, it can be found that in the frequency domain, its amplitude has a value only at positive frequencies, and its phase remains unchanged. The property can be utilized to apply the complex field expansion to the omnibearing directional wave decomposition, and the calculated amount is effectively low.

S402: carrying out Fourier transform on the longitudinal wave field after the time complex field expansion in the Z direction to obtain a wave number field longitudinal wave field;

where Kz represents the wavenumber of the longitudinal and transverse wavefield in the Z direction. When a seismic wave propagates in an underground medium, it is defined as positive downward, so that a downward propagating directional wave has positive Kz, and an upward propagating directional wave has negative Kz.

S403: selecting a part with wave number larger than zero from the longitudinal wave field of the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the Z direction to obtain a downlink direction wave field of the longitudinal wave field; and selecting a part with wave number less than zero from the longitudinal wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the Z direction to obtain an uplink wave field of the longitudinal wave field.

Referring to fig. 1, the imaging method of the present application further includes the steps of:

the upward wave field of the longitudinal and transverse waves is subjected to left-right directional wave decomposition to be decomposed into a right upward wave field and a left upward wave field;

in this step, referring to fig. 5, the decomposition of the upward wave field of the longitudinal wave into a right upward wave field and a left upward wave field by performing a left-right directional wave decomposition includes the following sub-steps:

s501: performing time complex field expansion on the separated uplink wave field of the longitudinal wave field;

wherein Kx represents the wave number of the longitudinal wave longitudinal and transverse wave fields in the X direction, and is defined as positive right, so that the Kx of the direction wave propagating to the right is positive, and the Kx of the direction wave propagating to the left is negative;

s502: performing Fourier transform on an uplink wave field of the longitudinal wave field after the time complex field expansion in the X direction to obtain a wave number field longitudinal wave field;

s503: selecting a part with wave number larger than zero from the upward wave field of the longitudinal wave field in the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a right upward wave field of the longitudinal wave field; and selecting a part with wave number less than zero from the longitudinal wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a left uplink direction wave field of the longitudinal wave field.

Referring to fig. 6, the down-going direction wave field of the longitudinal wave is decomposed into a right down-going direction wave field and a left down-going direction wave field by performing left-right direction wave decomposition, and includes the following sub-steps:

s601: performing time complex field expansion on the wave field in the downstream direction of the separated longitudinal wave field;

s602: performing Fourier transform on a downstream direction wave field of the longitudinal wave field after the time complex field expansion in the X direction to obtain a wave number field longitudinal wave field;

s603: selecting a part with wave number larger than zero from the wave field of the longitudinal wave field in the wave number domain in the downward direction, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a wave field of the longitudinal wave field in the right downward direction; and selecting a part with wave number less than zero from the longitudinal wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a left downlink direction wave field of the longitudinal wave field.

Referring to fig. 7, the separation of the shear wave signal into an uplink shear wave signal and a downlink shear wave signal in the up-down direction includes the following sub-steps:

s701: carrying out time complex field expansion on the separated shear wave field;

s702: performing Fourier transform on the transverse wave field after the time complex field expansion in the Z direction to obtain a wave number field transverse wave field;

s703: selecting a part with wave number larger than zero from the wave number domain transverse wave field, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the Z direction to obtain a downstream direction wave field of the transverse wave field; and selecting a part with wave number less than zero from the wave number domain transverse wave field, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the Z direction to obtain an uplink direction wave field of the transverse wave field.

The complex field, Kz, inverse fourier transform, etc. in these several steps have already been explained above and will not be described further here.

In step S104, referring to fig. 8, the separation of the uplink transverse wave signal into a left uplink transverse wave signal and a right uplink transverse wave signal in the left-right direction includes the following substeps:

s801: performing time complex field expansion on the separated uplink wave field of the transverse wave field;

s802: performing Fourier transform on an uplink wave field of the transverse wave field after the time complex field expansion in the X direction to obtain a wave number field transverse wave field;

s803: selecting a part with wave number larger than zero from the uplink wave field of the transverse wave field in the wave number domain, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a right uplink wave field of the transverse wave field; and selecting a part with wave number less than zero from the transverse wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a wave field of the left uplink direction of the transverse wave field.

The complex field, Kx, inverse fourier transform, etc. of these several steps have also been explained above and will not be described further here.

Referring to fig. 9, in step S104, the step of separating the downlink shear wave signal into a left downlink shear wave signal and a right downlink shear wave signal in the left-right direction includes the following sub-steps:

s901: performing time complex field expansion on the wave field in the downstream direction of the separated transverse wave field;

s1072: performing Fourier transform on a wave field in a descending direction of the transverse wave field after the time complex field expansion in the X direction to obtain a wave number field transverse wave field;

s903: selecting a part with wave number larger than zero from the wave field of the transverse wave field in the wave number domain in the downward direction, setting the part with wave number smaller than zero as zero, and performing inverse Fourier transform in the X direction to obtain a wave field of the transverse wave field in the right downward direction; and selecting a part with wave number less than zero from the transverse wave field of the wave number domain, setting the part with wave number greater than zero as zero, and performing inverse Fourier transform in the X direction to obtain a left downlink direction wave field of the transverse wave field.

Referring to fig. 1, step S105 is: selecting a direction wave field of a specific combination according to the actual geological condition, and imaging by applying a vector imaging condition of angle constraint to obtain an underground imaging result; in this step, the elastic wave can be divided into 8 directional wave fields according to the propagation direction of the longitudinal wave and the transverse wave and four directions, up, down, left and right. Taking longitudinal waves as an example, each wave field after separation is recorded as

S, R represents the elastic wave at the shot point and the elastic wave at the wave detection point, respectively. The superscript p represents the longitudinal wave. Subscripts l, r represent left and right, u, d represent up and down. The longitudinal waves at the shot point and the wave detection point are decomposed in four directions, namely, up, down, left and right, and then can be expressed as follows:

similarly, the transverse waves at the shot point and the wave detection point can be expressed as follows after being decomposed in four directions:

taking the longitudinal wave as an example, the directional waves of the upper, lower, left and right quadrants of the P wave at the shot point can be expressed by the following formula:

other locations and other similar wavefields may be expressed in the form of this equation and will not be described in further detail herein.

Meanwhile, in order to reduce the amount of calculation and the amount of storage, the data can be expanded to a complex field. For data expanded to a complex domain, the value of ω < 0 is zero, so that time Fourier transform on the data is not needed, the calculation amount is reduced, and the formula (19) -the formula (22) are simplified as follows:

after longitudinal and transverse wave decoupling, vertical and horizontal directional wave decomposition and simplification are carried out on the elastic wave fields of shot point forward transmission and wave detection point backward transmission, 16 separated directional wave fields can be obtained. The result of selecting the directional wave combination to perform the cross-correlation imaging is more than the result of performing the combined cross-correlation imaging by only performing the longitudinal and transverse wave decoupling and the decomposition of the upper and lower directional waves to obtain 8 directional wave fields in the prior art.

Therefore, the method can have more choices of imaging results. Meanwhile, when the imaging results are relatively more, noise results in the imaging results can be removed more easily, and some specific directional wave combinations can be removed according to specific geological structures of geological target areas. For example, I _ldld,I _lulu,I _rdrd,I _ruruThe combination of these same subscripts is the main source of low frequency noise for imaging; i is _luld,I _lurd,I _ruld,I _rurdThese shot and geophone field upgoing and downgoing waves, when there is a strong velocity change, can lead to primary artifacts due to the presence of backscatter. Wherein I represents the imaging result.

For the horizontal structure, theoretically, a good imaging result can be obtained through longitudinal and transverse wave separation and up-down azimuth wave decomposition, but the actual geological condition is complex, a complete horizontal or vertical structure is not provided, only up-down direction waves or left-right direction waves are not enough, and only the longitudinal and transverse wave separation and the decomposition of the omnibearing azimuth waves are combined, specific noise can be suppressed, and a good imaging result can be obtained.

When the imaging result in which noise exists is excluded, the remaining waveforms may be combined and imaged according to a preset directional wave. For example, the preset combination of directional waves may include: a left down longitudinal wave at the shot point, a right up longitudinal wave at the wave detection point and a right up transverse wave at the wave detection point. With this combination the following imaging conditions can be obtained:

where m represents spatial position, θ represents discrete angles in ADCIGs, θ represents a distance between the two _pRepresents the actual calculated angle, sigma represents the variance of the gaussian function,

representing the left down-going longitudinal wave wavefield at the shot, representing the up-right longitudinal wave field at the detection point, representing the up-right traveling shear wave field at the wave detection point.

The imaging effect of the directional wave combination is relatively excellent and can be used as a final imaging result.

Of course, other directional wave combinations than the above may be selected for reverse time shift imaging. The exploration personnel can select multiple groups of combined imaging according to the actual terrain, compare the imaging results of various directional wave combinations, and select the imaging with higher quality as the result.

While the process flows described above include operations that occur in a particular order, it should be appreciated that the processes may include more or fewer operations, which may be performed sequentially or in parallel.

Example 1:

referring to fig. 10(a) to 11(c), in order to verify the imaging quality of the imaging method of the present application, an experimenter used the method to image on a marmousi2 model as well as a salt dome model. Referring specifically to fig. 10(a) and 10(b), fig. 10(a) and 10(b) are partial results of PP imaging and PS imaging on a marmousi2 model. Where PP denotes a longitudinal wave at the shot point and a longitudinal wave at the wave detection point. PS represents the longitudinal wave at the shot point and the transverse wave at the wave detection point.

FIGS. 11(a) to (c) are respectively the result of PP imaging of the method on a Marmousi2 model; PP imaging results on the salt dome model and PS imaging results on the salt dome model. Where PP denotes the longitudinal wave at the shot point and the longitudinal wave at the pickup point. PS represents the longitudinal wave at the shot point and the transverse wave at the wave detection point.

It can be seen from the above figures that the elastic wave is subjected to separation in the vertical and horizontal directions on the basis of separation of the longitudinal and transverse waves, so that the elastic wave can be decomposed into a plurality of directional wave fields, specific directional waves are combined and imaged to form various low-frequency noises and false noises, the specific directional wave combinations are removed, and a proper directional wave combination is selected for cross-correlation imaging, so that the imaging quality of the imaging method provided by the invention is better.

Based on the same inventive concept as the imaging method shown in fig. 1, the present application also provides an elastic wave imaging apparatus. Referring to fig. 15, the image forming apparatus may include: the construction unit 401 is configured to construct an underground shot forward transmission wave field and a demodulator probe backward transmission wave field; a separation unit 402, configured to determine a type of a medium that propagates the forward wavefield and the backward wavefield; when the medium is an isotropic medium, performing vertical and horizontal wave field separation on the forward transmission wave field and/or the backward transmission wave field according to a preset equation decoupling algorithm, and when the medium is non-uniform anisotropy, performing vertical and horizontal wave field separation on the forward transmission wave field and/or the backward transmission wave field according to a preset mutation function IDW interpolation algorithm; decomposing the separated longitudinal wave field and/or transverse wave field in four directions, namely up, down, left and right, to obtain a plurality of direction wave fields of the longitudinal wave field and/or the transverse wave field; an imaging unit 403, configured to select at least two directional wave fields from the multiple directional wave fields, and perform normalized cross-correlation vector imaging on the selected directional wave fields.

Based on the same inventive concept as the imaging method shown in fig. 1, the present application also provides a computer storage medium having stored thereon a computer program which, when executed by a processor, performs the steps of: acquiring an underground shot forward wave field and a wave detection point backward wave field of a target stratum; according to a preset equation decoupling algorithm, carrying out longitudinal and transverse wave field separation on a forward wave field and/or a backward wave field propagated in an isotropic medium of the target stratum; according to a preset IDW (inverse discrete wavelet) interpolation algorithm, carrying out vertical and horizontal wave field separation on a forward wave field and/or a backward wave field propagated in a non-uniform anisotropic medium of the target stratum; decomposing the separated longitudinal wave field and the transverse wave field in four directions, namely, up, down, left and right, to obtain a plurality of direction wave fields of the longitudinal wave field and/or the transverse wave field; and selecting at least two directional wave fields from the plurality of directional wave fields, and carrying out vector imaging on the selected directional wave fields.

For convenience of description, the above devices are described as being divided into various units by function, and are described separately. Of course, the functionality of the units may be implemented in one or more software and/or hardware when implementing the present application.

As will be appreciated by one skilled in the art, embodiments of the present invention may be provided as a method, system, or computer program product. Accordingly, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.

The present invention is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each flow and/or block of the flow diagrams and/or block diagrams, and combinations of flows and/or blocks in the flow diagrams and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.

These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks.

These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.

In a typical configuration, a computing device includes one or more processors (CPUs), input/output interfaces, network interfaces, and memory.

The memory may include forms of volatile memory in a computer readable medium, Random Access Memory (RAM) and/or non-volatile memory, such as Read Only Memory (ROM) or flash memory (flash RAM). Memory is an example of a computer-readable medium.

Computer-readable media, including both non-transitory and non-transitory, removable and non-removable media, may implement information storage by any method or technology. The information may be computer readable instructions, data structures, modules of a program, or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), Static Random Access Memory (SRAM), Dynamic Random Access Memory (DRAM), other types of Random Access Memory (RAM), Read Only Memory (ROM), Electrically Erasable Programmable Read Only Memory (EEPROM), flash memory or other memory technology, compact disc read only memory (CD-ROM), Digital Versatile Discs (DVD) or other optical storage, magnetic cassettes, magnetic tape magnetic disk storage or other magnetic storage devices, or any other non-transmission medium that can be used to store information that can be accessed by a computing device. As defined herein, the computer readable medium does not include a transitory computer readable medium such as a modulated data wave field and a carrier wave.

It should also be noted that the terms "comprises," "comprising," or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising an … …" does not exclude the presence of other like elements in a process, method, article, or apparatus that comprises the element.

As will be appreciated by one skilled in the art, embodiments of the present application may be provided as a method, system, or computer program product. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present application may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.

The application may be described in the general context of computer-executable instructions, such as program modules, being executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. The application may also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media including memory storage devices.

The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. In particular, for the system embodiment, since it is substantially similar to the method embodiment, the description is simple, and for the relevant points, reference may be made to the partial description of the method embodiment.

The above description is only an example of the present application and is not intended to limit the present application. Various modifications and changes may occur to those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the scope of the claims of the present application.

Claims

1. An imaging method, comprising:

selecting at least two directional wave fields from the plurality of directional wave fields, and performing normalized cross-correlation vector imaging on the selected directional wave fields; wherein the imaging conditions of the normalized cross-correlation vector imaging include:

representing the left down-going longitudinal wave wavefield at the shot, representing the up-right longitudinal wave field at the detection point,

representing the up-right traveling shear wave field at the demodulation point,

representing the imaging results of the left downlink longitudinal wave at the shot point and the right uplink longitudinal wave at the wave detection point,

and imaging results of the left downlink longitudinal wave at the shot point and the right uplink transverse wave at the wave detection point are shown.

2. The imaging method of claim 1, wherein said decomposing the separated compressional and/or shear wave fields in four directions comprises:

3. The imaging method of claim 2, wherein said up-down directional wave decomposition of the separated longitudinal wave wavefield comprises:

4. The imaging method of claim 3, wherein performing left-right directional wave decomposition on the up-going wave wavefield comprises:

5. The imaging method of claim 3, wherein performing left-right directional wave decomposition on the down-going directional wave-field comprises:

6. The imaging method of claim 2, wherein said up-down directional wave decomposition of the separated shear wave field comprises:

carrying out time complex field expansion on the separated shear wave field;

7. The imaging method of claim 6, wherein performing left-right directional wave decomposition on the up-going wave wavefield comprises:

8. The imaging method of claim 6, wherein performing left-right directional wave decomposition on the down-going directional wave-field comprises:

9. The imaging method as claimed in claim 1, wherein the performing a vertical and horizontal wavefield separation on the forward and backward wavefields according to a preset variogram IDW interpolation algorithm when the medium is non-uniform anisotropy comprises:

10. The imaging method of claim 9, wherein said selecting N reference models comprises:

11. The imaging method of claim 10, wherein the weight coefficients of the reference model are calculated by the variance function.

12. The imaging method of claim 11, wherein the variogram acquisition method comprises:

13. The imaging method of claim 1, wherein the vector imaging is elastic wave reverse time migration imaging.

14. An image forming apparatus, comprising:

the imaging unit is used for selecting at least two directional wave fields from the plurality of directional wave fields and carrying out normalized cross-correlation vector imaging on the selected directional wave fields; wherein the imaging conditions of the normalized cross-correlation vector imaging include:

where m represents spatial position, θ represents discrete angles in ADCIGs, θ represents a distance between the two _pRepresents the actual calculated angle, sigma represents the variance of the gaussian function, representing the left down-going longitudinal wave wavefield at the shot,

representing the up-right longitudinal wave field at the detection point,

representing the up-right traveling shear wave field at the demodulation point, representing the imaging results of the left downlink longitudinal wave at the shot point and the right uplink longitudinal wave at the wave detection point,

15. A computer storage medium having a computer program stored thereon, the computer program, when executed by a processor, performing the steps of:

representing the left down-going longitudinal wave wavefield at the shot,

representing the up-right longitudinal wave field at the detection point,

representing the up-right traveling shear wave field at the demodulation point,