Methods and devices performing adaptive quadratic Wasserstein full-waveform inversion

Abstract

Methods and devices for seismic exploration of an underground structure apply W.sup.2-based full-wave inversion to transformed synthetic and seismic data. Data transformation ensures that the synthetic and seismic data are positive definite and have the same mass using an adaptive normalization. This approach yields superior results particularly when the underground structure includes salt bodies.

Claims

1. A method of seismic exploration of an underground structure using a Wasserstein metric-based full-wave inversion, FWI, the method comprising: obtaining seismic data acquired over the explored underground structure; generating synthetic data based on a velocity model; transforming the synthetic data and the seismic data by integrating a difference between the synthetic and the seismic data along each trace to generate transformed synthetic data and corresponding to null transformed seismic data; applying a linear normalization to traces of the transformed synthetic and the transformed seismic data using an adaptive normalization so that normalized corresponding traces to have same mass; updating the velocity model by applying the Wasserstein metric-based FWI to the normalized corresponding traces of the transformed synthetic and seismic data; and determining presence of natural resources in the explored underground structure based on the updated velocity model.

2. The method of claim 1, wherein the generating, the transforming, the applying of the linear normalization and the updating are performed iteratively until a loop-exit condition is met.

3. The method of claim 1, wherein the applying of the Wasserstein metric-based FWI includes computing a Frechet derivative of the synthetic data.

4. The method of claim 1, wherein a normalization constant used in the linear normalization is minimized while maintaining the normalized synthetic data and the normalized seismic data remaining positive definite.

5. The method of claim 1, wherein the underground formation includes a salt body.

6. The method of claim 1, further comprising generating a subsurface model based on the velocity model and other information, prior to generating the synthetic data.

7. The method of claim 1, wherein the transformed synthetic data {tilde over (f)} and the transformed seismic data {tilde over (g)} are
{tilde over (f)}(t)=∫.sub.0.sup.tf(s)ds−∫.sub.0.sup.tg(s)ds, {tilde over (g)}(t)=0 with f corresponding to the synthetic data and g to the seismic data, and the applying of the linear normalization is performed according to $\hat{f}\equiv \frac{\tilde{f}+c}{.Math.\tilde{f}+c.Math.},\hat{g}\equiv \frac{\tilde{g}+c}{.Math.\tilde{g}+c.Math.}$ with normalization constant c ensuring that {tilde over (f)}+c>0 and {tilde over (g)}+c>0.

8. A seismic data processing apparatus configured to perform seismic exploration of an underground structure using a Wasserstein metric-based full-wave inversion, FWI, the apparatus comprising: an interface configured to obtain seismic data acquired over the explored underground structure; and a data processing unit connecting to the interface and configured to generate synthetic data based on a velocity model; to transform the synthetic data and the seismic data by integrating a difference between the synthetic and the seismic data along each trace to generate transformed synthetic data corresponding to null transformed seismic data; to apply a linear normalization to traces of the transformed synthetic data and the transformed seismic data using an adaptive normalization so that normalized corresponding traces to have same mass; to update the velocity model by applying the Wasserstein metric-based FWI using the normalized corresponding traces of the transformed synthetic and of the seismic data; and to determine presence of natural resources in the explored underground structure based on the updated velocity model.

9. The seismic data processing apparatus of claim 8, wherein the data processing unit performs iteratively generating the synthetic data based on the velocity model, transforming the synthetic data and the seismic data, applying the linear normalization, and updating the velocity model by applying the W.sup.2-based FWI to the normalized corresponding traces of the synthetic and seismic data, until a loop-exit condition is met.

10. The seismic data processing apparatus of claim 8, wherein, when applying of the Wasserstein metric-based FWI , the data processing unit computes a Frechet derivative of the synthetic data.

11. The seismic data processing apparatus of claim 8, wherein a normalization constant used in the linear normalization is minimized while maintaining the normalized synthetic data and the normalized seismic data remain positive definite.

12. The seismic data processing apparatus of claim 8, wherein the underground formation includes a salt body.

13. The seismic data processing apparatus of claim 8, wherein the data processing unit generates a subsurface model based on the velocity model and other information, prior to generating the synthetic data.

14. The seismic data processing apparatus of claim 8, wherein the transformed synthetic data {tilde over (f)} and the transformed seismic data {tilde over (g)} are
{tilde over (f)}(t)=∫.sub.0.sup.tf(s)ds−∫.sub.0.sup.tg(s)ds, {tilde over (g)}(t)=0 with f corresponding to the synthetic data and g to the seismic data, and the applying of the linear normalization is performed according to $\hat{f}\equiv \frac{\tilde{f}+c}{.Math.\tilde{f}+c.Math.},\hat{g}\equiv \frac{\tilde{g}+c}{.Math.\tilde{g}+c.Math.}$ with normalization constant c ensuring that {tilde over (f)}+c>0 and {tilde over (g)}+c>0.

15. A computer readable storage medium storing executable codes, which, when executed by a processor, make the processor perform a method of seismic exploration of an underground structure using a Wasserstein metric-based full-wave inversion, FWI, the method comprising: obtaining seismic data acquired over the explored underground structure; generating synthetic data based on a velocity model; transforming the synthetic data and the seismic data by integrating a difference between the synthetic and the seismic data along each trace to generate transformed synthetic data corresponding to null transformed seismic data; applying a linear normalization to traces of the transformed synthetic and the transformed seismic data using an adaptive normalization so that normalized corresponding traces to have same mass; updating the velocity model by applying the Wasserstein metric-based FWI to the normalized corresponding traces of the transformed synthetic and seismic data; and determining presence of natural resources in the explored underground structure based on the updated velocity model.

16. The computer readable storage medium of claim 15, wherein the generating, the transforming, the applying of the linear normalization and the updating are performed iteratively until a loop-exit condition is met.

17. The computer readable storage medium of claim 15, wherein the applying of the Wasserstein metric-based FWI includes computing a Frechet derivative of the synthetic data.

18. The computer readable storage medium of claim 15, wherein a normalization constant used in the linear normalization is minimized while maintaining the normalized synthetic data and the normalized seismic data remain positive definite.

19. The computer readable storage medium of claim 15, wherein the underground formation includes a salt body.

20. The computer readable storage medium of claim 15, wherein the transformed synthetic data {tilde over (f)} and the transformed seismic data {tilde over (g)} are
{tilde over (f)}(t)=∫.sub.0.sup.tf(s)ds−∫.sub.0.sup.tg(s)ds, {tilde over (g)}(t)=0 with f corresponding to the synthetic data and g to the seismic data, and the applying of the linear normalization is performed according to $\hat{f}\equiv \frac{\tilde{f}+c}{.Math.\tilde{f}+c.Math.},\hat{g}\equiv \frac{\tilde{g}+c}{.Math.\tilde{g}+c.Math.}$ with normalization constant c ensuring that {tilde over (f)}+c>0 and {tilde over (g)}+c>0.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) For a more complete understanding of the present inventive concept, reference is now made to the following descriptions taken in conjunction with the accompanying drawings, in which:

(2) FIG. 1 is a graph illustrating a synthetic trace and a real data trace;

(3) FIG. 2 is a graph illustrating cost functions versus time shift;

(4) FIG. 3 is another graph representing cost functions versus time shift;

(5) FIG. 4 is a dataflow according to an embodiment;

(6) FIGS. 5A-5D are graphs illustrating the benefit of applying the W.sup.2-based FWI according to an embodiment to model data;

(7) FIGS. 6A-6D are graphs illustrating the benefit of applying the W.sup.2-based FWI according to an embodiment to real data;

(8) FIG. 7 is a flowchart of a method according to an embodiment; and

(9) FIG. 8 is a block diagram of a seismic data processing apparatus according to an embodiment.

DETAILED DESCRIPTION OF THE INVENTION

(10) The following description of the exemplary embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Instead, the scope of the invention is defined by the appended claims. The following embodiments are discussed using the terminology of seismic data processing for exploring underground structures, in particular, related to FWI methods based on the Wasserstein metric.

(11) Reference throughout the specification to “one embodiment” or “an embodiment” means that a particular feature, structure or characteristic described in connection with an embodiment is included in at least one embodiment of the subject matter disclosed. Thus, the appearance of the phrases “in one embodiment” or “in an embodiment” in various places throughout the specification is not necessarily referring to the same embodiment. Further, the particular features, structures or characteristics may be combined in any suitable manner in one or more embodiments.

(12) The quadratic Wasserstein metric (also known as W.sup.2 distance) is computed trace-by-trace for the synthetic data f and real data g, considered to be two density distributions defined on the same domain:
W.sup.2(f, g)=∫.sub.0.sup.T|t−G.sup.−1(F(t)|.sup.2f(t)dt   (1)
where F(t)=∫.sub.0.sup.tf(s)ds, G(t)=∫.sub.0.sup.tg(s)ds.

(13) Formula (1) can also be expressed in matrix form:
W.sup.2(f, g)=(t−G.sup.−1(F(t))).sup.Tdiag(f)(t−G.sup.−1(F(t))).   (2)

(14) A Frechet derivative with respect to the synthetic data f (i.e., the adjoint source) has to be computed in order to apply the W.sup.2-based cost function within the FWI framework. The adjoint-state method is a general framework for FWI, not specific to W2-based FWI. Different cost function metrics, such as, the conventional L2 or W2 may be used within this framework (as described, for example, in the 2006 article “A review of the adjoint-state method for computing the gradient of a function with geophysical application” by Plessix R.-E. published in Geophysical Journal international, 167(2), pp. 495-503).

(15) The adjoint source may be computed according to the following formula:

(16) $\begin{array}{cc}\nabla W^2=\left[-2\mathrm{diag}\left(\frac{d{G}^{-1}\left(y\right)}{dy}{|}_{F\left(t\right)}\right)U\mathrm{diag}\left(f\right)+\mathrm{diag}\left(t-G^{-1}\left(F\left(t\right)\right)\right)\right]\left(t-G^{-1}\left(F\left(t\right)\right)\right)& \left(3\right)\\ \mathrm{where}& \\ \left(\frac{d{G}^{-1}\left(y\right)}{dy}{|}_{F\left(t\right)}\right)=\frac{1}{g\left(G^{-1}\left(F\left(t\right)\right)\right)}& \left(4\right)\end{array}$
and U is the upper triangle matrix.

(17) The adjoint source is then back-propagated to compute the gradient according to the adjoint-state method.

(18) For the W.sup.2-based FWI to work properly, the density distributions f and g have to satisfy two conditions: positivity and mass conservation. Both f and g have to be positive everywhere (i.e., f>0, g>0), and their total masses must be equal (i.e., ∫.sub.0.sup.Tf=∫.sub.0.sup.Tg). Seismic data intrinsically has both positive and negative parts, and there is no guarantee that the total masses of synthetic and real data are equal. Therefore, normalization is necessary before applying the W.sup.2 metric. Various normalization techniques have been discussed in literature (for example, in the 2017 article, “Full-waveform inversion with an exponentially encoded optimal-transport norm,” by Qui et al., published in 87th Annual International Meeting, SEG, Expanded Abstracts, 1286-1290, or the 2018 article, “Application of optimal transport and the quadratic Wasserstein metric to full-waveform inversion,” by Yang et al., published in Geophysics, 83(1), R43-R62), with each technique having its advantages and drawbacks.

(19) Linear normalization seems most reliable when dealing with real data. The drawback of linear normalization is that it weakens the convexity of the W.sup.2 cost function. The following transformation may be applied to the seismic data before linear normalization to mitigate this problem:
{tilde over (f)}(t)=∫.sub.0.sup.tf(s)ds−∫.sub.0.sup.tg(s)ds, {tilde over (g)}(t)=0.   (5)

(20) Here, {tilde over (f)}(t) which is the integrals over the wavefields difference better promotes low frequency. The wavefield represents the sound wave amplitudes that travel through the subsurface earth after a shot is fired (i.e., a 4D volume of values in spatial (x, y, z) and time surveyed). A seismic trace is obtained by detecting (i.e., sampling amplitude of) these waves using receivers at a specific surface position (receiver position) for the recorded period of time. In other words, a trace is a projection of the wavefield to a particular position and along the time axis. Integrating the wavefield means integrating along the time dimension. This can be done on the trace or on the wavefield directly.

(21) An adaptive linear normalization is then applied to the transformed traces {tilde over (f)} and {tilde over (g)} to handle weak events and amplitude mismatch:

(22) $\begin{array}{cc}\hat{f}\equiv \frac{\tilde{f}+c}{<\tilde{f}+c>},\hat{g}\equiv \frac{\tilde{g}+c}{<\tilde{g}+c>},& \left(6\right)\end{array}$
where c is a normalization constant that ensures {tilde over (f)}+c>0 and {tilde over (g)}+c>0. Amplitude normalization yields a better handling of weak events and alleviates the negative effect of amplitude mismatches.

(23) The normalized signals {circumflex over (f)} and ĝ are used in equations (1)-(3) for FWI velocity updating.

(24) Thus, the data transformation performed prior to applying W.sup.2-based FWI includes two procedures: (A) integrating the wavefields, and (B) adaptive normalization. Integrating the wavefield boosts low frequency, thereby enhancing convexity. Adaptive normalization means subtracting the integrated signals and matching the subtracted result to a trace with an average amplitude of zero. FIGS. 1-3 illustrate these operations.

(25) The graph in FIG. 1 represents a synthetic trace 110 and a real data trace 120 (which are two identical Ricker wavelets) separated by a certain time shift. The graph in FIG. 2 represents the W.sup.2-based cost function plotted as a function of the time shift. Specifically, curve 210 is the L.sup.2-based cost function, curve 220 is the W.sup.2-based cost function with c=5, curve 230 is the W.sup.2-based cost function with c=1, and curve 240 is the W.sup.2-based cost function with c=0.5. The W.sup.2-based cost function with simple linear normalization loses global convexity, but remains better than the L.sup.2-based cost function. In addition, the convexity of the W.sup.2-based cost function weakens as the normalization constant c increases. Therefore, constant c is preferably as small as possible while still satisfying the positivity requirement.

(26) The adaptive normalization performed according to equation (5) allows using a smaller c by removing a large part of the energy from the traces. This benefit becomes clearer when the real seismic traces contain multiple events, with some significantly stronger than others. Adaptive normalization effectively removes the already matched events, which are usually strong, from the problem and gradually emphasizes the weaker events. This allows a smaller c value and boosted convexity.

(27) The graph in FIG. 3 illustrates the overall benefit brought by the adaptive normalization. Curve 310 represents the L.sup.2-based cost function without integral wavefield and curve 320 represents the L.sup.2-based cost function with integral wavefield. Curve 330 represents the W.sup.2-based cost function without integral wavefield and adaptive normalization, curve 340 improved represents the W.sup.2-based cost function with integral wavefield but without adaptive normalization, and curve 350 represents the W.sup.2-based cost function with integral wavefield and adaptive normalization. By comparing curves 310 with 320 and curve 330 with 340 and 350, one notes that integral wavefield enhances the convexity for both metrics.

(28) The adaptive quadratic Wasserstein FWI approach responds to the cycle skipping challenge by enhancing convexity via integral wavefield, and to the amplitude challenge by reducing amplitude sensitivity upon using adaptive normalization.

(29) FIG. 4 is a dataflow according to an embodiment. The initial velocity model 410 may be based on newly acquired seismic data and/or well data or on historic data from previous seismic exploration in the same area. Optional additional properties 420 may include other factors, such as anisotropic parameters, elastic parameters, etc. The recorded dataset 430 may include, for example, the geometry of the acquisition layout in addition to the seismic data. Inputs 410, 420 and 430 are used to generate a subsurface model of the underground structure at 440. The subsurface model specifies location of interfaces between volumes with different seismic wave propagation velocities, each interface being associated with a reflection coefficient that characterizes reflection and refraction occurring when a seismic wave is incident to an interface.

(30) Steps 440-480 codify the adaptive quadratic Wasserstein FWI and may then be performed (looped over) multiple times. The decision of whether or not to reiterate these steps is based on whether the W.sup.2-based cost function calculated at 480 has decreased below a predetermined threshold. At 490 a velocity model of the geophysical structure is output; if migration (e.g., reverse time migration of data using this velocity model) is applied to this velocity model then an image of the geophysical structure is obtained. Analysis of the image can suggest, to those trained in the field, the presence or absence of oil and/or gas among other natural resources inside the geophysical structure.

(31) FIGS. 5A-5D illustrate data processing of Marmousi model-based data (i.e., synthetic data) used to validate the above-described approach. FIGS. 5A-5D are four two-dimensional (2D) graphs, each having a horizontal dimension (of about 15 km) on x-axis, depth (down to 3 km) on y-axis and velocity in a range of 1,450 to 4,500 m/s illustrated as nuances of gray. A Ricker wavelet peaked at 4 Hz was used to generate the synthetic shots, with a fixed-spread acquisition geometry covering the model. FIG. 5A is the Marmousi model, FIG. 5B is the result of L.sup.2-based FWI, FIG. 5C is the result of a simple (i.e., no integral wavefield and adaptive normalization) W.sup.2-based FWI, and FIG. 5D is the result of W.sup.2-based FWI with adaptive normalization and integral wavefield. The FWIs were carried out at a single frequency of 6 Hz, with a 1D initial velocity model v(z). The L.sup.2-based FWI is severely cycle-skipped. Simple W.sup.2-based FWI has alleviated the cycle-skipping issue to some extent, but W.sup.2-based FWI with integral wavefield and adaptive normalization better recovers the model, thanks to its superior convexity.

(32) Similar to FIGS. 5A-5D, FIGS. 6A-6D are 2D graphs illustrating a test of W.sup.2-based FWI applied to a real data test. The data was acquired with a staggered-vessel configuration, which provided full-azimuth coverage for offsets up to 10 km, and long offsets up to 18 km for certain azimuths. The explored geophysical structure includes salt bodies 610 and 620, and a complex salt geometry 630 in-between, which is difficult to image using L.sup.2-based FWI. The initial velocity model illustrated in FIG. 6A was based on tomography and salt interpretation. FIG. 6B is the result of applying reverse time migration (RTM) to the initial velocity model in FIG. 6A.

(33) The FWI was carried out from 3 to 7 Hz. FIG. 6C illustrates the updated velocity model obtained by applying W.sup.2-based FWI, and FIG. 6D is the result of applying RTM to the updated velocity model in FIG. 6C. Applying W.sup.2-based FWI made significant corrections to the salt geometry and, as a result, the RTM image is much improved compared to the initial image, especially for the subsalt events.

(34) FIG. 7 is a flowchart of a method 700 according to an embodiment. The method includes obtaining seismic acquired over the explored underground formation at 710. Method 700 then includes generating synthetic data based on a velocity model at 720 and transforming the synthetic data and the seismic data to be positive definite and have the same mass using adaptive normalization at 730. Method 700 further includes updating the velocity model using W.sup.2-based FWI at 740. The updated velocity model is usable to determine the presence of natural resources in the explored underground formation.

(35) Steps 720 to 740 may be performed iteratively until a loop-exit condition is met. The loop-exit condition may be one of a predetermined number of iterations, with the W.sup.2-based cost function becoming lower than a predetermined threshold, a change of the W.sup.2-based cost function becoming less than a convergence threshold, etc., and a combination of such conditions.

(36) Step 730 (i.e., the transforming of the seismic data and the synthetic data) includes integrating wavefield and adaptive normalization. The adaptive normalization may be a linear normalization. The linear normalization may use a normalization constant (as in formula (6)) which is chosen to be as small as possible while the synthetic data and the seismic data remain positive definite. The underground formation may include a salt body.

(37) In one embodiment, a subsurface model is generated based on the velocity model and other information, prior to generating the synthetic data.

(38) FIG. 8 is a block diagram of a seismic data processing apparatus 800 according to an embodiment. Apparatus 800 includes an interface 810, a CPU 820, a user interface 830 and a data storage device 840. Interface 810 is configured to obtain seismic data acquired for the underground structure from other devices or the data acquisition equipment (collectively labeled 812). The user interface, which is optional, may include a display and may allow a user to vary the conditions of performing methods according to various embodiments (for example, may input an initial normalization constant that is then automatically optimized).

(39) CPU is configured to perform methods similar to method 700 described above. The CPU may include a software or hardware interface to a computer-readable recording medium such as the data storage device 840 or an external device. The computer-readable recording medium stores executable codes, which, when executed by the CPU, make the CPU perform methods similar to 700.

(40) The CPU may be configured to perform iteratively (1) generating synthetic data based on the velocity model, (2) transforming the synthetic data and the seismic data, and (3) updating the velocity model by applying the W.sup.2-based FWI to the transformed synthetic and seismic data, until a loop-exit condition (as already specified) is met.

(41) When transforming the synthetic and seismic data, the CPU may integrate the wavefield and perform an adaptive normalization. The adaptive normalization may be a linear normalization. The linear normalization may use a normalization constant chosen to be as small as possible if the synthetic data and the seismic data remain positive definite.

(42) In one embodiment the CPU may be configured to generate a subsurface model based on the velocity model and other information, prior to generating the synthetic data.

(43) The disclosed embodiments provide methods and systems for seismic exploration of an underground structure using adaptive quadratic Wasserstein (W.sup.2) FWI. It should be understood that this description is not intended to limit the invention. On the contrary, the embodiments are intended to cover alternatives, modifications and equivalents, which are included in the spirit and scope of the invention as defined by the appended claims. Further, in the detailed description of the embodiments, numerous specific details are set forth in order to provide a comprehensive understanding of the claimed invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

(44) Although the features and elements of the present embodiments are described in the embodiments in particular combinations, each feature or element can be used alone without the other features and elements of the embodiments or in various combinations with or without other features and elements disclosed herein.

(45) This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the subject matter is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims.