Method of high-resolution amplitude-preserving seismic imaging for subsurface reflectivity model

Abstract

The present disclosure provides a method of high-resolution amplitude-preserving seismic imaging for a subsurface reflectivity model, including: performing reverse time migration (RTM) to obtain an initial imaging result, performing Born forward modeling on the initial imaging result to obtain seismic simulation data, and performing RTM on the seismic simulation data to obtain a second imaging result; performing curvelet transformation on the two imaging results, performing pointwise estimation in a curvelet domain, and using a Wiener solution that matches two curvelet coefficients as a solution of a matched filter; and applying the estimated matched filter to the initial imaging result to obtain a high-resolution amplitude-preserving seismic imaging result.

Claims

1. A method of high-resolution amplitude-preserving seismic imaging for a subsurface reflectivity model, comprising: obtaining seismic observation data; performing reverse time migration (RTM) on the seismic observation data to obtain an initial imaging result, performing Born forward modeling on the initial imaging result to obtain seismic simulation data, and performing RTM on the seismic simulation data to obtain a second imaging result; performing curvelet transformation on the two imaging results, performing pointwise estimation in a curvelet domain, and using a Wiener solution that matches two curvelet coefficients as a solution of a matched filter; and applying the estimated matched filter to the initial imaging result and consequently obtaining a high-resolution amplitude-preserving seismic image; wherein: (1) performing RTM on the seismic observation data to obtain the initial imaging result m.sub.0(x), (2) substituting the initial imaging result m.sub.0(x) into a Born forward operator to obtain seismic simulation data; (3) performing RTM on the seismic simulation data to obtain the second imaging result m′(x); (4) performing curvelet transformation on the initial imaging result m.sub.0(x) and the second imaging result m′(x) to obtain transform coefficients C.sub.m.sub.0 and C.sub.m′ respectively; (5) using the Wiener solution $F\left(k,a,\theta \right)=\frac{C_{m_0}\left(k,a,\theta \right)C_{m^{\prime }}^{*}\left(k,a,\theta \right)}{C_{m^{\prime }}\left(k,a,\theta \right)C_{m^{\prime }}^{*}\left(k,a,\theta \right)+\mathrm{.Math.}},$  namely, an approximation of an inverse Hessian operator, wherein F represents the matched filter, K represents a wavenumber vector in a spatial frequency domain, a represents a scale, θ represents an angle, and ε is a regularization parameter; and (6) applying the inverse Hessian operator estimated in step (5) to C.sub.m0 and then performing inverse curvelet transformation to obtain an improved imaging result.

2. The method of high-resolution amplitude-preserving seismic imaging for a subsurface reflectivity model according to claim 1, wherein in step (1), RTM is performed on the seismic observation data by using a finite difference method to solve a wave equation, specifically comprising the following three steps: propagating forward a source wavefield, propagating backward a geophone wavefield, and applying a cross-correlation imaging condition.

3. The method of high-resolution amplitude-preserving seismic imaging for a subsurface reflectivity model according to claim 2, wherein the source wavefield is propagated forward in time, and the geophone wavefield is propagated backward in time from a recording moment; and the seismic observation data is obtained through finite difference forward modeling of the wave equation with a real velocity model as an input.

4. The method of high-resolution amplitude-preserving seismic imaging for a subsurface reflectivity model according to claim 1, wherein in step (2), Born forward modeling is performed on an obtained seismic imaging section to obtain the seismic simulation data, specifically comprising: (2.1) propagating forward a source wavefield to obtain a background wavefield p.sub.s(x), wherein the background wavefield is obtained through finite difference forward modeling of a wave equation with a migration velocity model as an input; (2.2) substituting $m\left(x\right)\frac{\partial }{\partial t}{p}_s\left(x\text{;}t\text{;}{x}_s\right)$  into the wave equation as a boundary condition and performing forward propagation to obtain a scattered wavefield p.sub.r(x); and (2.3) recording scattered wavefields at different moments at a receiver position to obtain the seismic simulation data.

5. The method of high-resolution amplitude-preserving seismic imaging for a subsurface reflectivity model according to claim 1, wherein in step (5), the matched filter is an operator that can map C.sub.m′ to C.sub.m.sub.0, pointwise estimation is performed, and the entire matching process is performed in a curvelet transform domain.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a real velocity model used in the present disclosure for testing;

(2) FIG. 2 is a migration velocity model used in the present disclosure for testing;

(3) FIG. 3 is a diagram that shows observed seismic data;

(4) FIG. 4 is a diagram that shows seismic simulation data obtained by performing Born forward modeling;

(5) FIG. 5 is a diagram that shows an initial imaging result obtained by performing RTM on observed seismic data;

(6) FIG. 6 is a diagram that shows a second imaging result obtained by performing RTM on seismic simulation data;

(7) FIG. 7 is a diagram that shows an improved imaging result obtained by the method of the present disclosure;

(8) FIG. 8 is a diagram that shows an imaging result obtained after 10 iterations of conventional LSRTM;

(9) FIG. 9 is a diagram that shows single-channel comparison of a real reflection coefficient, result of RTM, and result of the present disclosure; and

(10) FIG. 10 is a diagram that compares wavenumber spectra of results of RTM, LSRTM, and the present disclosure.

DETAILED DESCRIPTION OF THE EMBODIMENTS

(11) To make persons skilled in the art better understand the present disclosure, the following clearly and completely describes the technical solutions in the embodiments of the present disclosure 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 embodiments obtained by persons of ordinary skill in the art based on the embodiments of the present disclosure without creative efforts shall fall within the scope of protection of the present disclosure.

(12) It should be noted that the terms “first”, “second”, and so on in the specification and claims of the present disclosure and in the accompanying drawings are intended to distinguish between similar objects but do not necessarily indicate a specific order or sequence. It should be understood that the objects used in such a way may be exchanged under proper conditions to make it possible to implement the described embodiments of the present disclosure in other sequences apart from those illustrated or described here. Moreover, the terms “include”, “contain”, and any other variants mean to cover the non-exclusive inclusion, for example, a process, method, system, product, or device that includes a list of steps or units is not necessarily limited to those steps or units which are clearly listed, but may include other steps or units which are not expressly listed or inherent to such a process, method, system, product, or device.

(13) The present disclosure will be further explained in detail below with reference to the accompanying drawings.

(14) Step 1. Perform RTM on observation data to obtain an initial imaging result m.sub.0(x).

(15) A finite difference method is used to solve a wave equation to simulate the wave equation. During RTM imaging, a source wavefield is propagated forward, a geophone wavefield is propagated backward, and a cross-correlation imaging condition is applied. For an observation record D(x.sub.r;t;x.sub.s) a source is located at x.sub.s=(x.sub.s,y.sub.s,z.sub.s=0), and a geophone is located at x.sub.r=(x.sub.r,y.sub.r,z.sub.r=0). The forward propagation of the source wavefield can be described as follows:

(16) $\begin{array}{cc}\{\begin{array}{c}\left(\frac{1}{v^2\left(x\right)}\frac{{\partial }^2}{\partial t^2}-{\nabla }^2\right)p_s\left(x\text{;}t\text{;}{x}_s\right)=0\\ p_s\left(x\text{;}t\text{;}{x}_s\right)=\delta \left(x-x_s\right){\int }_0^{t}f\left(t^{\prime }\right){\mathrm{dt}}^{\prime }\end{array}& \left(1\right)\end{array}$

(17) where v(x) represents a velocity model, x represents coordinates of an underground position, p.sub.s represents the source wavefield, and f(t) represents a source wavelet. FIG. 1 and FIG. 2 show a real velocity model and a migration velocity model used in the present disclosure for testing respectively.

(18) The backward propagation of the geophone wavefield can be described as follows:

(19) $\begin{array}{cc}\{\begin{array}{c}\left(\frac{1}{v^2\left(x\right)}\frac{{\partial }^2}{\partial t^2}-{\nabla }^2\right)p_r\left(x\text{;}t\text{;}{x}_s\right)=0\\ p_r\left(x_r\text{;}t\text{;}{x}_s\right)=D\left(x_r,t\text{;}{x}_s\right)\end{array}& \left(2\right)\end{array}$

(20) where p.sub.r represents the geophone wavefield.

(21) The cross-correlation imaging condition is applied to obtain an image of underground structure.
m(x)=∫∫.sub.0.sup.T.sup.maxp.sub.s(x;t;x.sub.s)p.sub.r(x;t;x.sub.s)dtdx.sub.s  (3)

(22) where m(x) represents a stack imaging section, and T.sub.max represents total reception time of a seismic record. RTM is performed on the observation data shown in FIG. 3 to obtain the imaging result shown in FIG. 5.

(23) Step 2. Perform Born forward modeling to obtain seismic simulation data. First, propagate forward a source wavefield to obtain a background wavefield p.sub.s, described as follows:

(24) $\begin{array}{cc}\left(\frac{1}{v^2\left(x\right)}\frac{{\partial }^2}{\partial t^2}-{\nabla }^2\right)p_s\left(x\text{;}t\text{;}{x}_s\right)=\delta \left(x-x_s\right)f\left(t\right)& \left(4\right)\end{array}$

(25) Next, use

(26) $m\left(x\right)\frac{\partial }{\partial t}{p}_s\left(x\text{;}t\text{;}{x}_s\right)$
as a boundary condition to obtain a scattered wavefield p.sub.r, described as follows:

(27) $\begin{array}{cc}\left(\frac{1}{v^2\left(x\right)}\frac{{\partial }^2}{\partial t^2}-{\nabla }^2\right)p_r\left(x\text{;}t\text{;}{x}_s\right)=m_0\left(x\right)\frac{\partial }{\partial t}{p}_s\left(x\text{;}t\text{;}{x}_s\right)& \left(5\right)\end{array}$

(28) Finally, record scattered wavefields at a receiver position to obtain the seismic simulation data.
d(x.sub.r;t;x.sub.s)=p.sub.r(x.sub.r;t;x.sub.s)  (6)

(29) FIG. 4 shows the seismic simulation data.

(30) Step 3. Perform RTM on the seismic simulation data to obtain a second imaging result m′(x).

(31) Referring to FIG. 6, the specific implementation method is the same as that in step 1, and) only the observation data D(x.sub.r;t;x.sub.s); needs to be replaced by the simulation data d(x.sub.r;t;x.sub.s).

(32) Step 4. Perform curvelet transformation on the initial imaging result m.sub.0(x) and second imaging result m′(x) to obtain transform coefficients C.sub.m.sub.0 and C.sub.m′ respectively.
C.sub.m.sub.0(k,a,θ)=m.sub.0(x)φ*.sub.k,a,θ(x)dx  (7)
and
C.sub.m′(k,a,θ)=m′(x)φ*.sub.k,a,θ(x)dx  (8)

(33) where * represents complex conjugate.

(34) Step 5. Solve a matched filter in a curvelet domain as an approximation of an inverse Hessian operator.

(35) To match C.sub.m.sub.0 and C.sub.m′, this problem is described as minimizing the following least-squares objective function:
E(F)=½∥FC.sub.m′−C.sub.m.sub.0∥.sub.2.sup.2  (9)

(36) where F represents the matched filter. A Wiener solution for a minimum mean square error can be obtained by performing the following pointwise estimation:

(37) $\begin{array}{cc}F\left(k,a,\theta \right)=\frac{C_{m_0}\left(k,a,\theta \right)C_{m^{\prime }}^{*}\left(k,a,\theta \right)}{C_{m^{\prime }}\left(k,a,\theta \right)C_{m^{\prime }}^{*}\left(k,a,\theta \right)+\mathrm{.Math.}}& \left(10\right)\end{array}$

(38) where k represents a wavenumber vector in a spatial frequency domain, a represents a scale, θ represents an angle, and ε is a regularization parameter used to ensure stability of division.

(39) Step 6. Apply the estimated inverse Hessian operator to C.sub.m.sub.0 and then perform inverse curvelet transformation to obtain an improved imaging result. This process can be described as follows:
{circumflex over (m)}=Ψ.sup.−1FC.sub.m.sub.0  (11)

(40) where Ψ.sup.−1 represents the inverse curvelet transformation, F represents the matched filter estimated in step 5, C.sub.m.sub.0 represents a curvelet transform coefficient of the initial imaging result m.sub.0(x), and {circumflex over (m)} represents the high-resolution imaging result obtained by this method. The high-resolution imaging result shown in FIG. 7 is of comparable quality to that obtained after 10 iterations of conventional LSRTM in FIG. 8. The resolution is improved, and the amplitude is more balanced.

(41) The method of the present disclosure requires only two RTM imagings and one Born forward modeling, while 10 iterations of the conventional LSRTM require 10 RTM imagings and nine Born forward modelings. The method of the present disclosure greatly improves computation efficiency. Moreover, computation amount of curvelet transformation and inverse curvelet transformation is negligible compared with that of one wave propagation process. Single-channel comparison in FIG. 9 shows that the imaging result obtained by the method of the present disclosure is closer to a real reflection coefficient than the result obtained by the RTM, and maintains a relative amplitude relationship of reflection coefficients. Comparison of wavenumber spectra in FIG. 10 shows that both the method of the present disclosure and the LSRTM can recover high-wavenumber components to improve resolutions of seismic images.

(42) In conclusion, the method of the present disclosure resolves blurring and unbalanced amplitude, in imaging results, that are caused because an RTM imaging operator is an adjoint operator of a forward operator. RTM is performed on the initial imaging result. The matched filter of the two imaging result is found in the transform domain by performing curvelet transformation and used as an approximation of an inverse Hessian operator. The matched filter is applied to the initial imaging result to improve imaging quality.

(43) The above contents are merely used to illustrate the technical ideas of the present disclosure, rather than to limit the protection scope of the present disclosure. Any variations made based on the technical solutions according to the technical ideas proposed by the present disclosure shall fall within the protection scope as defined by the claims of the present disclosure.