Method and system for stabilizing Poynting vector of seismic wavefield

Abstract

The present disclosure provides a method and system for stabilizing a Poynting vector of a seismic wavefield. The method includes: adjusting an amplitude of a time derivative of the seismic wavefield, and computing a Poynting vector of the adjusted time derivative of the seismic wavefield to obtain a first Poynting vector, where a difference between the amplitude of the first Poynting vector and the amplitude of a second Poynting vector is within a set range, and the second Poynting vector belongs to the seismic wavefield; and conducting operation on the second Poynting vector and the first Poynting vector to obtain a final Poynting vector of the seismic wavefield. The present disclosure addresses instability of Poynting vectors computation.

Claims

1. A method for stabilizing a Poynting vector of a seismic wavefield, comprising: adjusting an amplitude of a time derivative of the seismic wavefield, and computing a Poynting vector of the adjusted time derivative of the seismic wavefield to obtain a first Poynting vector, wherein a difference between the amplitude of the first Poynting vector and the amplitude of a second Poynting vector is within a set range, and the second Poynting vector belongs to the seismic wavefield; and conducting operation on the second Poynting vector and the first Poynting vector to obtain a final Poynting vector of the seismic wavefield; wherein said adjusting an amplitude of a time derivative of the seismic wavefield, and computing a Poynting vector of the adjusted time derivative of the seismic wavefield to obtain a first Poynting vector specifically comprises:. computing a first-order time derivative of the seismic wavefield; computing an amplitude adjustment factor of the first-order time derivative of the seismic wavefield; adjusting the amplitude of the first-order time derivative of the seismic wavefield by adopting the amplitude adjustment factor to obtain an adjusted time derivative; and computing the Poynting vector of the adjusted time derivative to obtain the first Poynting vector.

2. (canceled)

3. The method for stabilizing a Poynting vector of a seismic wavefield according to claim 1, wherein said conducting operation on the second Poynting vector and the first Poynting vector to obtain a final Poynting vector of the seismic wavefield specifically comprises: adding the second Poynting vector and the first Poynting vector together to obtain the final Poynting vector.

4. The method for stabilizing a Poynting vector of a seismic wavefield according to claim 1, wherein said computing a first-order time derivative of the seismic wavefield specifically comprises: computing the first-order time derivative of the seismic wavefield by adopting a finite difference method.

5. The method for stabilizing a Poynting vector of a seismic wavefield according to claim 4, wherein a computational formula of the amplitude adjustment factor is as follows: $a=-\frac{\Delta t}{\sin \left(2\pi f_d\Delta t\right)};$ wherein α denotes an amplitude adjustment factor; f.sub.d denotes a basic frequency of a seismic wavefield; and Δt denotes a time interval of finite difference.

6. The method for stabilizing a Poynting vector of a seismic wavefield according to claim 1, wherein a computational formula of the adjusted time derivative is as follows: ${\stackrel{\_}{u}}^t=a\frac{\partial u^t}{\partial t};$ wherein ū.sup.t denotes an adjusted time derivative at moment t; α denotes an amplitude adjustment factor; t denotes time; u.sup.t denotes a seismic wavefield at moment t; and $\frac{\partial u^t}{\partial t}$ denotes a first-order time derivative obtained by finding a partial time derivative of u.sup.t.

7. The method for stabilizing a Poynting vector of a seismic wavefield according to claim 1, wherein a computational formula of the first Poynting vector is as follows: $P^{\stackrel{\_}{u}}=-\frac{\partial {\stackrel{\_}{u}}^t}{\partial t}\nabla {\stackrel{\_}{u}}^t;$ wherein P.sup.ū denotes a first Poynting vector; ū.sup.t denotes an adjusted time derivative at moment t; t denotes time; and ∇ū.sup.t denotes a gradient of ū.sup.t.

8. The method for stabilizing a Poynting vector of a seismic wavefield according to claim 1, wherein a computational formula of the first Poynting vector is as follows: $P^{\overline{u}}=-\frac{2u^t\sin \left(\pi f_d\Delta t\right)}{\Delta t}\nabla {\stackrel{\_}{u}}^t;$ wherein P.sup.ū denotes a first Poynting vector; ū.sup.t denotes an adjusted time derivative at moment t; t denotes time; ∇ū.sup.t denotes the process of finding the gradient of ū.sup.t; u.sup.t denotes a seismic wavefield at moment t; f.sub.d denotes a basic frequency of the seismic wavefield; and Δt denotes a time interval of finite difference.

9. The method for stabilizing a Poynting vector of a seismic wavefield according to claim 4, wherein a computational formula of the first-order time derivative is as follows: $\frac{\partial u^t}{\partial t}=\frac{u^{t+\Delta t}-u^{t-\Delta t}}{2\Delta t};$ wherein u.sup.t denotes a seismic wavefield at moment t; t denotes time; Δt denotes a time interval of finite difference; u.sup.t+Δt denotes a seismic wavefield at moment t+Δt; and u.sup.t−Δt denotes a seismic wavefield at moment t−Δt.

10. A system for stabilizing a Poynting vector of a seismic wavefield, comprising: an amplitude adjustment module, which is configured to adjust an amplitude of a time derivative of the seismic wavefield, and compute a Poynting vector of the adjusted time derivative of the seismic wavefield to obtain a first Poynting vector, wherein a difference between the amplitude of the first Poynting vector and the amplitude of a second Poynting vector is within a set range, and the second Poynting vector belongs to the seismic wavefield; and a Poynting vector operation module, which is configured to conduct operation on the second Poynting vector and the first Poynting vector to obtain a final Poynting vector of the seismic wavefield; wherein the amplitude adjustment module is specifically configured to: compute a first-order time derivative of the seismic wavefield; compute an amplitude adjustment factor of the first-order time derivative of the seismic wavefield; adjust the amplitude of the first-order time derivative of the seismic wavefield by adopting the amplitude adjustment factor to obtain an adjusted time derivative; and compute the Poynting vector of the adjusted time derivative to obtain the first Poynting vector.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0034] To describe the embodiments of the present disclosure or the technical solutions in the related art more clearly, the accompanying drawings required in the embodiments are briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present disclosure. A person of ordinary skill in the art may further obtain other accompanying drawings based on these accompanying drawings without creative labor.

[0035] FIG. 1 is a schematic diagram illustrating time-amplitude of a seismic wavelet and a first-order time derivative thereof;

[0036] FIG. 2 is a flowchart of a method for stabilizing a Poynting vector of a seismic wavefield according to an embodiment of the present disclosure;

[0037] FIG. 3 is a schematic diagram illustrating time-amplitude of a Poynting vector of a seismic wavefield and a Poynting vector of a first-order time derivative of the seismic wavefield according to an embodiment of the present disclosure;

[0038] FIG. 4 is a schematic diagram illustrating a seismic wavefield in a homogeneous medium at moment 0.4 s;

[0039] FIG. 5 is a schematic diagram illustrating Poynting vectors computed by an original Poynting vector equation, a spatial smoothing method, an optical flow method, a time shifting method and the method in the present disclosure for the seismic wavefield as shown in FIG. 4; where FIG. 5(a) is a schematic diagram illustrating an x component of the Poynting vector obtained by the original Poynting vector computation equation, FIG. 5(b) is a schematic diagram illustrating a z component of the Poynting vector obtained by the original Poynting vector computation equation, FIG. 5(c) is a schematic diagram illustrating an x component of the Poynting vector computed by the smoothing method, FIG. 5(d) is a schematic diagram illustrating a z component of the Poynting vector computed by the smoothing method, FIG. 5(e) is a schematic diagram illustrating an x component of the Poynting vector obtained after conducting iteration five times by the optical flow method, FIG. 5(f) is a schematic diagram illustrating a z component of the Poynting vector obtained after conducting iteration by the optical flow method, FIG. 5(g) is a schematic diagram illustrating an x component of the Poynting vector computed by the time shifting method, FIG. 5(h) is a schematic diagram illustrating a z component of the Poynting vector computed by the time shifting method, FIG. 5(i) is a schematic diagram illustrating an x component of the Poynting vector computed by the method in the present disclosure, and FIG. 5(j) is a schematic diagram illustrating a z component of the Poynting vector computed by the method in the present disclosure;

[0040] FIG. 6 is a schematic diagram illustrating the error between a propagation angle computed based on a Poynting vector and an actual propagation angle, where FIG. 6(a) is a schematic diagram illustrating the error between a propagation angle computed based on an original Poynting vector and an actual propagation angle, FIG. 6(b) is a schematic diagram illustrating the error between a propagation angle computed by the smoothing method and an actual propagation angle, FIG. 6(c) is a schematic diagram illustrating the error between a propagation angle computed by the optical flow method and an actual propagation angle, FIG. 6(d) is a schematic diagram illustrating the error between a propagation angle computed by the time shifting method and an actual propagation angle, and FIG. 6(e) is a schematic diagram illustrating the error between a propagation angle computed by the method in the present disclosure and an actual propagation angle;

[0041] FIG. 7 is a schematic diagram of a Marmousi model;

[0042] FIG. 8 is a schematic diagram illustrating a seismic wavefield of a Marmousi model at moment 1.0 s;

[0043] FIG. 9 is a schematic diagram illustrating a Poynting vector computed by using an original Poynting equation for the seismic wavefield as shown in FIG. 8; where FIG. 9(a) is a schematic diagram illustrating an x component of a Poynting vector computed by using an original Poynting equation for the seismic wavefield as shown in FIG. 8, and FIG. 9(b) is a schematic diagram illustrating a z component of a Poynting vector computed by using an original Poynting equation for the seismic wavefield as shown in FIG. 8;

[0044] FIG. 10 is a schematic diagram illustrating a Poynting vector obtained by using the method in the present disclosure for the seismic wavefield as shown in FIG. 8; where FIG. 10(a) is a schematic diagram illustrating an x component of a Poynting vector obtained by using the method in the present disclosure for the seismic wavefield as shown in FIG. 8, and FIG. 10(b) is a schematic diagram illustrating a z component of a Poynting vector obtained by using the method in the present disclosure for the seismic wavefield as shown in FIG. 8; and

[0045] FIG. 11 is a structure diagram of a system for stabilizing a Poynting vector of a seismic wavefield according to an embodiment of the present disclosure.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0046] The technical solutions in the embodiments of the present disclosure will be described below clearly and completely with reference to the accompanying drawings in the embodiments of the present disclosure. Apparently, the described embodiments are merely some rather than all of the embodiments of the present disclosure. All other examples obtained by a person of ordinary skill in the art based on the examples of the present disclosure without creative efforts shall fall within the protection scope of the present disclosure.

[0047] To make the above-mentioned objective, features, and advantages of the present disclosure clearer and more comprehensible, the present disclosure will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0048] Propagation vectors of a seismic wavefield play a very important role in the field of seismic exploration. As a key parameter for generation of angle-domain common-image gathers, propagation vectors are widely used in reverse time migration for suppression of low frequency noise caused by backscattering, illuminance compensation, correction of polarity reversal of converted waves, etc. Propagation vectors can not only be obtained in a spatial-temporal domain, but also be obtained in a frequency-wavenumber domain. Frequency-wavenumber domain methods require a lot of Fourier transform, which leads to low computational efficiency. Since the direction of the Poynting vectors is consistent with that of the propagation vectors, and it only requires the computation of space and time derivatives of a seismic wavefield, Poynting vector computation achieves much higher computational efficiency than that of the frequency-wavenumber domain method, and is thus applied to a wide range of areas. However, Poynting vector computation has the defect of instability.

[0049] In order to solve the above problem, embodiments provide a method for stabilizing a Poynting vector of a seismic wavefield. FIG. 2 is a flowchart of a method for stabilizing a Poynting vector of a seismic wavefield according to an embodiment of the present disclosure. Referring to FIG. 2, the method includes the following steps:

[0050] Step 101: adjust an amplitude of a time derivative of the seismic wavefield, and compute a Poynting vector of the adjusted time derivative of the seismic wavefield to obtain a first Poynting vector, where the difference between the amplitude of the first Poynting vector and the amplitude of a second Poynting vector is within a set range, and the second Poynting vector belongs to the seismic wavefield.

[0051] Step 102: conduct operation on the second Poynting vector and the first Poynting vector to obtain a final Poynting vector of the seismic wavefield.

[0052] In an optional implementation, the set range in step 101 can be an order of magnitudes, that is, the amplitude of the time derivative of the seismic wavefield may be adjusted in this implementation, such that the amplitude of the Poynting vector (the first Poynting vector) of the adjusted time derivative of the seismic wavefield is within the same order of magnitudes as that of the second Poynting vector. The purpose of adjusting the two to the same order of magnitudes described herein is to improve computational accuracy.

[0053] In an optional implementation, step 101 specifically includes:

[0054] 1) Compute a first-order time derivative of the seismic wavefield. Specifically, the first-order time derivative of the seismic wavefield may be computed by adopting a finite difference method. A computational formula of the first-order time derivative is as follows:

[00007]$\frac{\partial u^t}{\partial t}=\frac{u^{t+\Delta t}-u^{t-\Delta t}}{2\Delta t};$

[0055] where

[00008]$\frac{\partial u^t}{\partial t}$

denotes a first-order time derivative obtained by finding a partial time derivative of u.sup.t; u.sup.t denotes a seismic wavefield at moment t; t denotes time; Δt denotes a time interval of finite difference; u.sup.t+Δt denotes a seismic wavefield at moment t+Δt denotes a seismic wavefield at moment t−Δt.

[0056] 2) Compute a Poynting vector of the seismic wavefield to obtain a second Poynting vector. A computational formula is specifically as follows:

[00009]$P^u=-\frac{\partial u^t}{\partial t}\nabla u^t;$

[0057] where P.sup.u denotes a second Poynting vector; and ∇ū.sup.t denotes a gradient of u.sup.t.

[0058] 3) Compute an amplitude adjustment factor of the first-order time derivative of the seismic wavefield. Specifically,

u.sup.t+Δt and u.sup.t−Δt may be expressed as:

u.sup.t+Δt=u.sup.texp(−iωΔt);

u.sup.t−Δt=u.sub.texp(iωΔt);

[0059] Temporal difference is expressed as:

[00010]$\frac{u^{t+\Delta t}-u^{t-\mathrm{\Delta t}}}{2\Delta t}=\frac{u^t}{2\Delta t}\left[\exp \left(-i\mathrm{\omega \Delta }t\right)-\exp \left(i\mathrm{\omega \Delta }t\right)\right]=\frac{u^t}{2\Delta t}\left[\cos \left(\mathrm{\omega \Delta }t\right)-i\sin \left([\mathrm{\omega \Delta }t\right)\right]-\frac{u^t}{2\Delta t}\left[\cos \left(\mathrm{\omega \Delta }t\right)+i\sin \left([\mathrm{\omega \Delta }t\right)\right]=-\frac{\sin \left(\mathrm{\omega \Delta }t\right)}{\Delta t}i{u}^t;$

[0060] iu.sup.t and u.sup.t have the same amplitude, but have a phase difference of 90 degrees therebetween. The advantage of phase difference of 90 degrees is that: the amplitudes of the two Poynting vectors are complementary, with the maximum value of the one corresponding to the minimum value of the other. The above formula is followed by:

[00011]$i{u}^t=-\frac{\Delta t}{\sin \left(\omega \Delta t\right)}\frac{u^{t+\Delta t}-u^{t-\Delta t}}{2\Delta t}=-\frac{\Delta t}{\sin \left(\omega \Delta t\right)}\frac{\partial u^t}{\partial t};$

[0061] where ω=2πf, ω denotes an angular frequency of a seismic wavefield, f denotes a frequency of a seismic wavefield, from which the following can be derived:

[00012]$i{u}^t=-\frac{\Delta t}{\sin \left(2\pi f\Delta t\right)}\frac{\partial u^t}{\partial t};$

[0062] For the sake of simplicity, f is made equivalent to f.sub.d, and then

[00013]$i{u}^t=-\frac{\Delta t}{\sin \left(2\pi f_d\Delta t\right)}\frac{\partial u^t}{\partial t};$

[0063] a computational formula of the amplitude adjustment factor is as follows:

[00014]$a=-\frac{\Delta t}{\sin \left(2\pi f_d\Delta t\right)};$

[0064] where α denotes an amplitude adjustment factor; and f.sub.d denotes a basic frequency of a seismic wavefield.

[0065] 4) Adjust the amplitude of the first-order time derivative of the seismic wavefield by adopting the amplitude adjustment factor to obtain an adjusted time derivative. Specifically, a computational formula of the adjusted time derivative is as follows:

[00015]${\stackrel{\_}{u}}^t=a\frac{\partial u^t}{\partial t};$

[0066] where ū.sup.t denotes an adjusted time derivative at moment t.

[0067] 5) Compute the Poynting vector of the adjusted time derivative to obtain the first Poynting vector. The first Poynting vector may be computed using either of the following two methods.

[0068] Method I: a computational formula of the first Poynting vector is as follows:

[00016]$P^{\overline{u}}=-\frac{\partial {\stackrel{\_}{u}}^t}{\partial t}\nabla {\stackrel{\_}{u}}^t;$

[0069] where P.sup.ū denotes a first Poynting vector; and ∇ū.sup.t denotes a gradient of ū.sup.t.

[0070] Method II: a computational formula of the first Poynting vector is as follows:

[00017]$P^{\overline{u}}=-\frac{2u^t\sin \left(\pi f_d\Delta t\right)}{\Delta t}\nabla {\overline{u}}^t.$

[0071] Method II has the advantages that in the computational process, the procedure of solving partial time derivatives of ū.sup.t is omitted, which reduces the amounts of computation, and allows for higher computational efficiency.

[0072] In an optional implementation, step 102 specifically includes:

[0073] adding the second Poynting vector and the first Poynting vector together to obtain the final Poynting vector. A specific computational formula is as follows:

P=P.sup.u+P.sup.ū;

[0074] where P denotes a final Poynting vector.

[0075] By means of the method for stabilizing a Poynting vector of a seismic wavefield according to the foregoing embodiment, given the consistency between the direction of a Poynting vector of the seismic wavefield and that of a Poynting vector of its time derivative, although the amplitude of the Poynting vector of the seismic wavefield is much smaller than that of the Poynting vector of its time derivative, the variations in the amplitudes of the two complement each other, that is, as shown in FIG. 3, the minimum value of the Poynting vector of the seismic wavefield corresponds to the maximum value of the Poynting vector of the time derivative. Therefore, the instability in the zero value or minimum value can be avoided by adjusting amplitudes of the two within a set range (e.g., the same order of magnitudes) and then conducting operation (e.g., summation).

[0076] For the purpose of illustrating the advantages of the present disclosure, the method in the present disclosure is compared with other methods proposed to solve the instability of Poynting vector computation.

[0077] Existing methods for addressing instability of Poynting vector computation are as follows:

[0078] (1) Spatial smoothing method. After a Poynting vector is obtained, it is subject to smoothing using Gaussian function.

[0079] (2) Multiplex one-dimensional smoothing method. For two-dimensional data and three-dimensional data, Gaussian smoothing function is used directly, which leads to a larger amount of computation. In order to lower the amount of computation, the multiplex one-dimensional smoothing method is adopted. For example, regarding an x component of a Poynting vector, one-dimensional smoothing is first conducted in the x direction, and then conducted in the z direction.

[0080] (3) Temporal smoothing method. Firstly, a Poynting vector at each position in space is computed within a period of time, and then subject to temporal smoothing.

[0081] (4) Time integration method. Firstly, a Poynting vector at each position in space is computed within a period of time, and then subject to time integration. The integral value is regarded as a Poynting vector at a precomputed position.

[0082] (5) Time shifting method. Firstly, a Poynting vector at each position in space is computed within a period of time, and then the maximum value of an absolute value of the Poynting vector within this period of time is searched. The Poynting vector with maximum value is regarded as a Poynting vector at a computed position.

[0083] (6) Optical flow method. Developed from computer vision, this method helps to obtain the accurate propagation direction step by step through iteration. An iterative formula of the method is as follows:

[00018]$\{\begin{array}{c}{p}_x^{n+1}={\stackrel{\_}{p}}_x^n-\frac{{\phi }_x\left({\phi }_x{\stackrel{\_}{p}}_x^n+{\phi }_y{\stackrel{\_}{p}}_y^n+{\phi }_x{\stackrel{\_}{p}}_y^n+{\phi }_t\right)}{{\alpha }^2+{.Math.{\phi }_x.Math.}^2+{.Math.{\phi }_y.Math.}^2+{.Math.{\phi }_x.Math.}^2}\\ p_y^{n+1}={\stackrel{\_}{p}}_y^n-\frac{{\phi }_y\left({\phi }_x{\stackrel{\_}{p}}_x^n+{\phi }_y{\stackrel{\_}{p}}_y^n+{\phi }_x{\stackrel{\_}{p}}_y^n+{\phi }_t\right)}{{\alpha }^2+{.Math.{\phi }_x.Math.}^2+{.Math.{\phi }_y.Math.}^2+{.Math.{\phi }_x.Math.}^2}\\ p_z^{n+1}={\stackrel{\_}{p}}_z^n-\frac{{\phi }_z\left({\phi }_x{\stackrel{\_}{p}}_x^n+{\phi }_y{\stackrel{\_}{p}}_y^n+{\phi }_x{\stackrel{\_}{p}}_y^n+{\phi }_t\right)}{{\alpha }^2+{.Math.{\phi }_x.Math.}^2+{.Math.{\phi }_y.Math.}^2+{.Math.{\phi }_x.Math.}^2}\end{array};$

[0084] where p.sub.x, p.sub.y and p.sub.z denote components of a Poynting vector P in three directions of x, y and z under Cartesian coordinates; the superscript n denotes a number of iterations; φ.sub.x=∂u/∂x, φ.sub.y=∂u/∂y, φ.sub.z=∂u/∂z; and α denotes a weighting coefficient.

[0085] Among the above methods, Method (1) to Method (4) increase the stability of Poynting vector, but reduce the accuracy of Poynting vector computation at the same time. Method (5) and Method (6), despite their computational accuracy and stability, require a large number of search operations (as required in time shifting method) or iterative operation (as required in optical flow method), which leads to low computational efficiency. By contrast, the method for stabilizing a Poynting vector based on a seismic wavefield and a first-order time derivative thereof according to the present disclosure has high computational efficiency while guaranteeing the computational accuracy and stability.

[0086] To further illustrate the advantages of the present disclosure, the following two specific examples are adopted to compare the method proposed in the present disclosure with an existing method.

[0087] In this example, a simple model is adopted. FIG. 4 is a schematic diagram illustrating a seismic wavefield in a homogeneous medium at moment 0.4 s, with a velocity of 2.0 km/s and a basic frequency of 30 Hz. FIG. 5 is a schematic diagram illustrating Poynting vectors computed by an original Poynting vector equation, a spatial smoothing method, an optical flow method, a time shifting method and the method in the present disclosure, the original Poynting vector equation being

[00019]$P=-\frac{\partial u}{\partial t}\nabla u;$

where P denotes a Poynting vector, and u denotes a seismic wavefield.

[0088] As can be seen from FIG. 4 and FIG. 5, the original Poynting vector has distinct stripes in white bands and black bands, and values on these stripes are close to or equal to zero. A point where the value is equal to zero is regarded an instable point. In other methods, there are no stripes or stripes are less distinct. FIG. 6 shows the error between a propagation angle of seismic wavefield computed by using the Poynting vector shown in FIG. 5 and an actual propagation angle. The point with relatively high amplitude in FIG. 6a is the point where Poynting vector computation is instable. It can be seen from FIG. 6a-FIG. 6e that smoothing method, optical flow method, time shifting method and the method in the present disclosure all improve the stability of Poynting vector computation. Among all the foregoing methods, the method in the present disclosure and the time shifting method produce the smallest error, and achieve the highest computational efficiency.

[0089] Compared with the prior art, the present disclosure achieves higher accuracy as well as higher computational efficiency. For the smoothing algorithm, smoothing plays the role of filtering out high frequency and retaining low frequency, which, despite the capacity to avoid points with computational instability, cannot significantly improve the accuracy. Given the consistency between the propagation direction of the Poynting vector of the seismic wavefield and that of the Poynting vector of its first-order time derivative, amplitude adjustment factors will not change the direction of the Poynting vector, which enables the higher computational accuracy of the method described herein. Additionally, as the variations in the amplitude of the Poynting vector of the seismic wavefield and that of the Poynting vector of its first-order time derivative complement each other, that is, the maximum value of one corresponds to the minimum value of the other, the computational instability caused by the minimum value can be offset by adding the two together.

[0090] Compared with the time shifting method and the optical flow method, the method in the present disclosure achieves computational accuracy which is similar to that of the time shifting method and higher than that of the optical flow method; besides, as compared with the time shifting method which requires a large number of search operations, and the optical flow method which requires a large number of iterative operations, the method in present disclosure only needs to compute the Poynting vector for one more time, thus enabling high computational efficiency.

[0091] In this example, a complex model is adopted. FIG. 7 shows a Marmousi model, FIG. 8 shows a seismic wavefield at moment 1.0 s, FIG. 9 shows a Poynting vector computed by using an original Poynting vector equation, and FIG. 10 shows a Poynting vector computed by the method in the present disclosure. Comparison between FIG. 10 and FIG. 9 shows that the stripes in the black and white bands in FIG. 9 are very distinct and there are many instable points, while the stripes in FIG. 10 are less distinct and the stability has been greatly improved.

[0092] Poynting vector obtained by using the method in the present disclosure plays a very important role in seismic exploration, and can be applied to a variety of scenarios. For example, Poynting vector can be widely used in reverse time migration for suppression of low frequency noise caused by backscattering, illuminance compensation, correction of polarity reversal of converted waves, etc.

[0093] The low frequency noise caused by backscattering is suppressed in reverse time migration, and a Poynting vector P.sub.src of a seismic source wave field and a Poynting vector P.sub.rec of a rectification point wave field are solved respectively by using the method provided in this embodiment, where subscript src denotes the source wave field, rec denotes the rectification point wave field, and the reflection angle is obtained by using the following formula:

[00020]$\theta =0.5\times \arccos \left(\frac{P_{\mathrm{src}}.Math.P_{\mathrm{rec}}}{.Math.P_{\mathrm{src}}.Math..Math.P_{\mathrm{rec}}.Math.}\right);$

[0094] Low frequency noise can then be suppressed by using imaging in the following formula:

I(x)=∫.sub.0.sup.T.sup.maxcos.sup.nθ.Math.u.sub.src(x,t)u.sub.rec(x,t)dt;

[0095] where I denotes an imaging profile under reverse time migration, x denotes a spatial position, T.sub.max denotes a maximum imaging time, n is a positive integer greater than 1, which generally delivers a better effect when set to be greater than or equal to 3.

[0096] The present disclosure further provides a system for stabilizing a Poynting vector of a seismic wavefield. FIG. 11 is a structure diagram of a system for stabilizing a Poynting vector of a seismic wavefield according to an embodiment of the present disclosure. Referring to FIG. 11, the system includes:

[0097] an amplitude adjustment module 201, which is configured an amplitude of a time derivative of the seismic wavefield, and compute a Poynting vector of the adjusted time derivative of the seismic wavefield to obtain a first Poynting vector, where a difference between the amplitude of the first Poynting vector and the amplitude of a second Poynting vector is within a set range, and the second Poynting vector belongs to the seismic wavefield; and

[0098] a Poynting vector operation module 202, which is configured to conduct operation on the second Poynting vector and the first Poynting vector to obtain a final Poynting vector of the seismic wavefield.

[0099] By means of the system for stabilizing a Poynting vector of a seismic wavefield according to an embodiment of the present disclosure, given the consistency between the direction of a Poynting vector of the seismic wavefield and that of a Poynting vector of its time derivative, although the amplitude of the Poynting vector of the seismic wavefield is much smaller than that of the Poynting vector of its time derivative, the variations in the amplitudes of the two complement each other, that is, the minimum value of the Poynting vector of the seismic wavefield corresponds to the maximum value of the Poynting vector of the time derivative. Therefore, the instability in the zero value or minimum value can be avoided by adjusting amplitudes of the two within a set range and then conducting operation.

[0100] Each example of the present specification is described in a progressive manner, each example focuses on the difference from other examples, and the same and similar parts between the examples may refer to each other. Since the system disclosed in an embodiment corresponds to the method disclosed in another embodiment, the description is relatively simple, and reference can be made to the method description.

[0101] Specific examples are used herein to explain the principles and embodiments of the present disclosure. The foregoing description of the embodiments is merely intended to help understand the method of the present disclosure and its core ideas; besides, various modifications may be made by a person of ordinary skill in the art to specific embodiments and the scope of application in accordance with the ideas of the present disclosure. In conclusion, the content of the present description shall not be construed as limitations to the present disclosure.