Resolving high frequency anomalies of gas cloud using full waveform inversion

High resolution models with structurally improved results significant to the physical properties of rocks in geologically complex areas require advance modeling methodologies. Low frequencies are required to understand the geological properties of the rocks while high frequencies is needed to address the structural challenges. Recent industry success in inversion have shown the accurate and robust results for the low frequencies. In this work, we provide a strategy to resolve geologically complex area such as gas cloud (at high frequencies) using full waveform inversion (FWI) based on 2D wave equation. Our contribution here, is the improvement in FWI imaging by: (i) solving the wave field equation to recover high resolution inversion results which is consistence to physical properties and shows structural enhancements; (ii) estimating the distribution of local minima which is largely affected by initial velocity model. To validate our approach, we demonstrate algorithms on synthetic gas


Introduction
To estimate the physical properties of the Earth subsurface through seismic-waves at multiscale frequencies for quantitative as well as qualitative information is a challenging task.In principle, full waveform inversion (FWI) is a suitable tool to obtain detailed subsurface properties of the Earth subsurface at high frequency.FWI is a data-strategy equation based on an ill-posed and highly non-linear approach.It has proved its effectiveness in significantly improving the spatial resolution of complex problem such as gas cloud [Prajapati and Ghosh 2015].Inversion at high frequencies is often challenging, though having the well log information helps in minimizing the problem to some extent as they contain 0-100 Hz frequencies or more.In reality, we often don't have these information for complex geological structure such as gas cloud.In that case inversion becomes much more challenging.
To exploit the full information contained in seismograms, FWI method was firstly introduced in time domain by Lailly [1983] and Tarantola [1984] and then later in frequency domain by Pratt et al. [1998], while industry success were shown around year 2009.
In complex geological environments, significant improvement in workflow can provide more accurate, robust and detailed inversion results which had been considered by the industry in deriving detailed velocity models.FWI has shown promising results at low frequencies (<12 Hz) and well resolved by Prajapati and Ghosh [2015], however, it is still unclear to what extent FWI can be applied efficiently at higher frequencies.
Some challenges have been recognized during inversion at high frequencies greater than 12 Hz.However, it has been seen through real and synthetic data that it improves the structural conformity by increasing the resolution [Prieux et al. 2013, Prajapati andGhosh 2015].These applications typically use only the low-frequencies (<12 Hz), but in practice, due to fact that we are not honoring the physics correctly and/or due to cycle-skipping or the local minima problem (prominent above 12 Hz) we generally avoid high frequencies in inversion results.
The presence of gas cloud has long been recognized as a significant problem in the seismic data processing around the world (North Sea, offshore South-East Asia).Due to complex wave propagation, anelastic losses, multiple scattering and strong contrast in these area inversion becomes more complicated and highly challenging.Non nonlinearity and non-uniqueness of the solution adds to this complexity.Seismic imaging below gas cloud is complex phenomenon and current approaches based on linear imaging do not provide satisfactory results.Many of the hydrocarbon bearing reservoirs in SE Asia have poor imaging [Ghosh et al. 2010], resulting from gas leakages and other complex geology.An example of such imaging problem (Figure 1) in SE Asia is shown, where two different gas cloud problem shown in the same field with similar structural play is a poor imaging due to "Seal" breach results in gas leakage.
The outline of the paper is summarized as follows.First, we review the methodology of full waveform inversion with different approach to solve the gas cloud challenges at high frequencies (Section 2).Next, we introduce the equation (Section 3) in detail to estimate the local minima effect on true model.Next, we solve the Gauss-Newton method to solve for high frequency (Section 4).Finally, we illustrate some insightful example (Section 5) using different approach of full waveform inversion as discussed above and conclude the paper with a discussion and conclusion (Section 6).To simplify the derivation, throughout the paper, we formulate the inversion in frequency-domain.

Full-waveform inversion
We consider the mathematical formulation for FWI in frequency domain [Prajapati and Ghosh 2015]: This equation defines a relation between the vector with gridded medium parameter m, A i (m) is the discretized Helmholtz operator, the wavefield u=[u 1 ;u 2 ;…], and the q i represent the discretized source vectors.

2D Inversion
We start with the initial model m 0 , assuming that this is closed to the true model with misfit of data where, p obs is the seismic data and p cal (m) is the model seismic data with a function of m.In the framework of misfit function C(m) with the starting velocity model m 0 plus a perturbation model m = m 0 +Dm is updated with every iteration with the l 2 norm given by where, † denote transposition, and Dp is the difference between the observed and those computed from for a model.If the problem is linearized, the minimum of misfit function of m 0 is approaches to zero when the derivative of the misfit function approaches to zero.

Formulation of simultaneous source encoding
It is well known that simultaneous source technique introduces random cross-talk that arises from the correlation between shots in the super-shots [Krebs et al. 2009, Prajapati andGhosh 2015].In order to mitigate cross-talk artifacts during the numerical inversion, a new encoding is generated at every iteration.High resolution result with high frequency can be obtained by re-sampling source positions and encoding functions at every iteration of the FWI algorithm.The simultaneous source encoding approach has been proposed to reduce the computational burden in time-domain [Guitton and Diaz 2012] and in the frequency-domain [Herrmann and Lin 2009].Another approach to create super-shots to reduce the cross-talk artifacts were introduced by Morton and Ober [1998] and Romero et al. [2000] using phase-encoding strategies of individual shots.
We formulated sophisticated algorithm of simultaneous inversion with a balance of low and high frequencies and uses the misfit function for full waveform inversion [Prajapati and Ghosh 2015] as where A(m) is the discretized Helmholtz operator assign for modeling parameter, m denote the model parameter, (5)

Formulation of local minima analysis
Full waveform inversion is highly sensitive to the starting velocity model.Most of FWI methods are unable to solve the convergence towards global local minima (Figure 2).If the starting model wavefield data within half wave cycle of the observed data (Figure 3), inversion solution will tend to converge local minima due to wrapping characteristics of phase in frequency domain [Bunks et al. 1995].Recent improvement based on Hessian estimation shows better convergence toward solution, because this search is quite computationally intensive [Virieux and Operto 2009], but difficult to estimate the distribution of local minima accurately in data.One way is to reduce the dimensionality model parameters; of course this is not an exact solution.
To estimate the distribution of local minima wavefield from the initial model, we referred the approach via the Lagrangian of Virieux and Operto [2009], and Huber et al. [2000].Which is typically, a forward modeling technique to recover a model and depends on the observed field u as, In practice, minimizing a l 2 -norm residual vector constrained optimization problem can be written as: (7) where m is vector model parameter; d i is the observed data; D i (m)=w i 2 diag(m)+L is the discretized Helmholtz operator with w i is the angular frequency and L is the discretized Laplacian operator; u i =[u 1 ;u 2 ;u 3 ;…] is the wavefield and q i is the source for all k experiments at each iteration of the optimization.
Due to constraint optimization, derivative-driven optimization is less affected by local minima [Huber et al. 2000].We can rewrite the above Equation (7) as: (8) And the gradient of the above objective function is given by ( 9) To compute the gradient, we solved both the forward D i (m i )u=q i and adjoint wave equations; however, in-situ of constrained method, the wavefield, gradient, and Hessian needs to be computed at each iteration  [ Pratt et al. 1998].Details of the problem reformulation and discretization can be found in Huber et al. [2000].At the end of each iteration, we obtained a forward and adjoint model with D, which contains a large, sparse matrix of inverse problem that can be solved for u i

Gauss-Newton method
The mathematical formulation based on Gauss-Newton method has been earlier explained by Virieux and Operto [2009] and Pratt et al. [1998].Compared to the conjugate-gradient or quasi-Newton method, Gauss-Newton method generally converges faster.We used l2 norm for FWI, which takes both amplitude and phase properties information into account in the inversion.Moreover, second-order information estimated by Hessian matrix is helpful in scaling and in the speed-up of the optimization.Before proceeding to higher frequencies with high resolution geological structures, we must recover the long wavelength or low-wavenumber component of the model, because theoretically at low frequency modeled and observed data match within half wave-cycle (Figure 3).Therefore phase matching is more important than amplitude; moreover, inversion of the frequencies recovers the long wavelength structures of the model and substantial help for cycle-skipping problem (Figure 3).The constraint of wavenumber content of the model is adapted based on the temporal frequencies included at each step of the FWI procedure [Wu and Toksöz 1987].

Inversion algorithm
An efficient alternative to solve FWI problem is by using the Gauss-Newton method.Instead of using firstorder steepest descent method (Equations 3 and 5), the second order Gauss-Newton method may speed up convergence rate during FWI.
The forward modeling function can be defined either in the time domain or in the frequency domain.
Here, we consider the mathematical formulation of forward modeling in the frequency domain.The objective function is defined as, (10) Model updating is carried out via the Gaussian-Newton method, (11) where Dp is the residual between the observed (p obs ) and the modeled (p k calc ) computed in the defined model space m k , J is the Fréchet derivative matrix obtained by perturbing each model parameter, J t J is the approximate Hessian is calculated implicitly and thus no need to store the Hessian matrix.

Example
The high frequency approach of FWI is tested on realistic synthetic gas cloud model which closely represent the shallow gas cloud environments of SE Asia (Figure 1) [Ghosh et al. 2010].The main aspects are (a) to enhance the structural properties of layer of media where the acoustic (Vp) velocity is significantly low as the reflection underneath the gas exhibit a 'push-down' effect due to localization of the anomaly, and (b) to image the gas cloud at higher frequencies ensuring structural conformity, while increasing resolution in the inversion results.

Example 1
To illustrate the structural information at high frequencies, we consider an acoustic velocity (Vp) model which is 5-km long offset and 0.7-km depth.Model is populated with Gaussian anomaly to make data more complex (Figure 4a).Note that due to Gaussian anomaly, the effect of attenuation, scattering phenomenon is very dominant; however, the adjoint velocity perturbation is less affected.An initial velocity model (Figure 4b) is used as background velocity for full waveform inversion.Inverting velocity range is varying from 800 m/s to 1600 m/s.In this experiment, first we generate synthetic data from 250 receivers at depth of 20-m and 125 sources at depth of 50-m that are equally spaced at the surface.For a source function, we used a full bandwidth source wavelet with a 55 Hz central frequency (Figure 5).
We use quasi-Newton technique such as Limitedmemory Broyden-Fletcher-Goldfarb-Shanno (l-BFGS) [Nocedal 1980] optimization method to invert the data in 12 discrete frequencies were selected between 3 Hz and 33 Hz and the inversion is carried out in a multiscale hierarchical strategy sequentially from low to high frequency data.We use the steepest decent method for optimization [Prajapati and Ghosh 2015], (12) where, a k is step length, Steepest-descent (i.e., the gradient ' −∇q [m] k ') method is used for the misfit function at point m 0 .The result of this exercise after five itera-tions is plotted in Figure 6a.Repeating this procedure leads to a better reconstruction of the wavefield through inversion at high frequencies, which is also observed in the corresponding model estimate as can be seen in Figure 6b-f.

Example 2
Although waveform inversion can provide accurate and high resolution geological models at both low and high frequencies (Figure 6); one of the main problems in FWI is the computational cost of the inversion for multiple sources and receivers.Using simultaneous source encoding, we solved the complex geological problem at low frequency [Prajapati and Ghosh 2015], but at high frequency, it is still a tough task.We briefly explain the formulation for simultaneous source encoding (Equations 4 and 5).For more details above equation, readers may follow Prajapati et al. [2014] and Prajapati and Ghosh [2015].
To illustrate the simultaneous FWI, we use same model (Figure 4) as used in example 1.The synthetic observed data (P obs ) is generated by solving the Helmholtz operator using an optimal 9-points finite difference method [Jo et al. 1996] in 12 equally spaced frequencies ranging from 3 Hz to 33 Hz.Both source and receiver are located at z = 20 m, with 40 m horizontal spacing interval.We use same wavelet as used in example 1.We used steepest-decent (i.e.gradient) method with l-BFGS (limited memory Broyden Fletcher Goldfrab Shanno) optimization method, which is a popular quasi-Newton method for solving non-linear optimization problem.Detail review and the application of the l-BFGS method on FWI can be found in Virieux and Operto [2009] and Nocedal [1980].The model parameters obtained from the inversion of the low frequency are used as starting velocity model for the higher frequency inversion.The inversion results at different frequency using simultaneous FWI, after 30 iterations are shown in Figure 7.These results indicate high resolution with structurally improved and significant to their properties.

Example 3
Objective function of waveform inversion is highly nonlinear, which can significantly affect the inversion results and solutions are often trapped in local minima.To avoid solutions that are trapped in local minima, one needs a good initial velocity model or data with low frequency components.However, getting good initial ve-locity model of complex geological structure is a challenging.Based on wave-Equations ( 6), ( 7) and ( 8), we estimate the distribution of the global minima effect of wavefield for a given medium in frequency domain.
To illustrate the variation of local minima, we consider a gas cloud model as depicted in Figure 1.The corresponding scatter wavefield at 12 Hz for a single source located at offset = 800 m, and depth = 50 m.We placed the receiver at 10-m depth with 100-m spacing interval.The result of this exercise is plotted in Figure 8a.Repeating this procedure with 100 shots and receivers at equal spacing at 10 m depth leads to a better reconstruction of the wavefield.We successfully estimated the distribution of local minima effect shown in Figure 8b with good initial velocity model.From result, it is clearly shown that local minima effect is significantly dominating by the boundary of sharp velocity contrast; however inside there is slight effect of local minima.

Example 4
To illustrate the FWI using Gauss-Newton method on high frequency, we consider a complex gas cloud model with gaussian anomaly (Figure 9a) and initial velocity model (Figure 9b).The original model had free-surface multiples turned on during the forward modeling.The data were regenerated using the pressure velocity field and the density with a 55 Hz maximum full bandwidth frequency (Figure 5).The data were acquired with 20 m depth for the source and receivers with 50 m for the shot intervals and 100 m for receiver interval.The inversion started with a smoothened original velocity field and 10 frequency bands between 3 Hz and 33 Hz.
Starting with the lower frequency band and marching towards the high frequency with minimizing the misfit function in each band throughout iterations.The inverted result after 30 iterations in each frequency band is shown in Figure 10a and the inversion result for selected frequency band is shown in Figure 10b-f.As the data residuals are directly proponational to the misfit of data during each iteration, so analyzing the residuals could be significant, and rate of convergence of the accuracy of the velocity field at a given iteration [Vigh et al. 2009].
In general, the data residual decreased substantially over the iteration (Figure 11).The residual energy RESOLVING GAS CLOUD AT HIGH FREQUENCY  can be explained by the small discrepancy between the true and final inverted velocity model, especially at the gas cloud boundary anomaly.

Discussion and conclusions
In this paper, we presented methodologies and strategies for FWI in frequency domain at high frequencies.The mathematical formulations include the conventional inversion, simultaneous source techniques, Gauss-Newton method.Source encoding relies on replacing all the sequential sources into a super-shot.
For each example, we selected discrete frequency band between 3 Hz and 33 Hz and inversion is carried out in a multiscale hierarchical strategy sequentially from low to high frequency.
To know the distribution of local minima effect on model data, we use Lagrange formulation to explore a larger search space using both the model and the wavefields.As the proposed method solves the governing wave equations independently, we avoid the storage of the gradient and Hessian matrix dataset.
As wave propagation though gas cloud lead to a significant deterioration of the seismic image below the gas cloud.Therefore, it is crucial to characterize the shallow gas cloud overburden, which has caused the wavefield to get trapped and scattered inside the gas cloud causing complex wave propagation.We envisage several possible inversions for complex and strong lateral variation of the overburden gas cloud problem.Since the problem is nonlinear, the waveform inversion is trivial as long as this method is efficient to solve the augmented acoustic wave equation.
As the minimization of the objective function is a nonlinear inverse problem and FWI is a local optimization method, the inversion results strongly depends on starting velocity model.Figure 2  The success of FWI depends on accurate initial velocity model to prevent cycle-skipping.Figure 3 illustrates the cycle-skipping problem of FWI.Incorrect starting velocity model i.e. more than half a wave-cycle in observed data, then inversion will be tapped in one of the many local minima and will be cycle-skipped (see the top schematic illustration of Figure 3).On the other hand if the starting velocity model is accurate enough i.e. within half a wave-cycle from the observed data and will end up at the local global minima (see the bottom schematic illustration of Figure 3).
Multilayer imaging below gas cloud is another chal- lenge to be solved using FWI technique.We use Gauss-Newton method to resolve this issue.We successfully imaged the layers below complex gas cloud overburden as shown in Figure 10.The presence of good signals at low frequencies retrieves the long wavelength structures of the model and this adds intermediate wavenumbers for higher frequency components, high resolution and the benefit of increased resolution in the inversion result, while ensuring structural conformity.

Figure 1 .
Figure 1.Gas masking effect: (a) dense shallow gas cloud masking with poor imaging; (b) poor imaging with gas leakage due to poor sealing of hydrocarbon reservoir.

Figure 2 .
Figure 2.An example of least-squares misfit function for nonlinear problem as a function of model parameter.Let consider X, Y and Z are three points.If the inversion is governed with starting model at point Y or Z then inversion will move to the strong or weak local minimum.On the other hand if the starting model at X, then inversion will be with global minimum.

Figure 3 .
Figure 3. Schematic illustration of the cycle-skipping problem in FWI.If the initial model and observed data have a lag is more than half a cycle, then the local minim moves in the wrong direction and if initial modeled and observed data have a lag is within half cycle, then the local minima moves in the optimum global direction.

Figure 4 .
Figure 4. (a) True velocity model and (b) corresponding initial velocity model for FWI.Note color bar is common for both figures.

Figure 5 .
Figure 5. (a) A full bandwidth source wavelet and (b) corresponding spectrum for forward modeling for FWI.

Figure 7 .
Figure 7. Same as Figure6, but for simultaneous source encoding FWI. 30 iterations for each frequency band were used.Frequency band is same as defined in Figure6caption.

Figure 8 .
Figure 8.(a) Scattering wavefield generated using single source at offset = 800 m and depth = 50 m; (b) scattering effect of local minima on whole data set with 100 sources and receivers.

Figure 9 .
Figure 9. (a) Complex gas cloud model and (b) corresponding background velocity model FWI.
is a schematic diagram which shows how the starting model affects the final inversion results.For example, if the starting model corresponds to position Y or Z, then we can observe strong or weak local minima, whereas if the starting model corresponds to X, then inversion ends up with local minima.

Figure 11 .
Figure 11.Misfit plot of selected frequency during FWI using Gauss-Newton method.It clearly shows that moving towards high frequency, misfit is decreasing drastically.