Acoustic 2 D full waveform inversion to solve gas cloud challenges

The existing conventional inversion algorithm does not provide satisfactory results due to the complexity of propagated wavefield though the gas cloud. Acoustic full waveform inversion has been developed and applied to a realistic synthetic offshore shallow gas cloud feature with Student-t approach, with and without simultaneous sources encoding. As a modeling operator, we implemented the grid based finite-difference method in frequency domain using second order elastic wave equation. Jacobin operator and its adjoint provide a necessary platform for solving full waveform inversion problem in a reduced Hessian matrix. We invert gas cloud model in 5 frequency band selected from 1 to 12 Hz, each band contains 3 frequencies. The inversion results are highly sensitive to the misfit. The model allows better convergence and recovery of amplitude losses. This approach gives better resolution then the existing least-squares approach. In this paper, we implement the full waveform inversion for low frequency model with minimum number of iteration providing a better resolution of in-


Introduction
Generally a two-step approach is applied in seismic exploration: (1) firstly, the acquired data is processed for primaries only for structured imaging; (2) secondly, to derive the lithological and fluid properties derived through AVO (amplitude versus offset) inversion.However, the above process is somewhat "flawed" as the imaging step is linear whereas the problem at hand is highly non-linear.
Hence, a more theoretically correct approach would be to use the full waveform non-linear inversion.This inversion process should contain not only the primaries but all the other arrivals (direct, refracted and reflected), including multiples in its solution for a unified approach.
Many of the hydrocarbon bearing reservoir in the SE Asia have poor imaging [Ghosh et al. 2010, Ghazali 2011], resulting from gas leakages and other complex geology.An example of such imaging problem (Figure 1) in SE Asia production field, where in one part there is no gas leakage hence, we have prefect imaging, whereas in another part of the same field with similar structural play is a poorly imaged due to "Seal" breach resulting from gas leakage.
Note that, when the wave propagates through a shallow anomaly such as gas cloud, the wave undergoes the following changes: (1) delay in arrival travel time due to significant decreased velocity in gas cloud area; (2) reduction in amplitude of due to attenuation; (3) phase distortion; (4) pre-secure of internal peg-leg multiples often refereed as "Coda" [Ghazali 2011].
The obvious approach is through Q-compensation [Reilly et al. 2008], which uses a tandem velocity Q model in pre-stack depth imaging to compensate for travel time delay, amplitude loss and phase distortion.However, Q-application in gas wipeout has limitation as vertical seismic profiling (VSP) measurement shows, Q to be as low as 40, and wherever Q model is correct the condition of a constant frequency Q is violated.
A second approach is from [Ghazali 2011].where he developed an "equivalent media" model to account for time delays and amplitude losses using inverse scattering concept.In first step he solves the kinematics to account for travel time and then uses an inverse approach to account for amplitude losses.Thus, a key contribution of his work is to incorporate transmission losses and hence geologically constraint reservoir model with appropriate acoustic properties.
The full waveform inversion (FWI) technique based on least-squares in time domain was first introduced by Tarantola [1984] and later by Pratt et al. [1998] in frequency domain.The basic idea of FWI is to minimize of the misfit between observed and calculated result after each optimization during inversion.

Subject classification:
Exploration geophysics, Seismic methods, Waves and wave analysis, Full waveform inversion, Gas cloud.
The calculation of misfit and gradient depends on the acquisition geometry (or the quality) of acquired data.Previously, some studies have been performed on complex gas cloud by: (1) using finite difference elastic model, (2) first arrival tomography [Virieux and Operto 2009], (3) using l 2 norm [Brossier et al. 2009].
In this study, we approached the problem in two dimensional acoustic-frequency domains, coupled with both simultaneous and non-simultaneous techniques.Results clearly demonstrate the effectiveness of proposed algorithm.where J is known as Jocabian matrix.

Student-t formulation
Due to convex or hyperbolic nature of data, the presence of large outliers and unexplained events in the data is a challenging task.Optimization algorithm like least-square cannot explain or recover the outliers.Disadvantage with the least-squares regression in seismic inversion has been discussed in Virieux and Operto [2009].In this paper, we use a different approach to recover the outliers and unexplained events by "Student-t" approach [Aravkin et al. 2011].
The modified Student-t penalty is written as, where, v− scale parameter, k− the degree of freedom, low value of k (10 -2 ) provide good reconstruction of perturbation.

Formulation of simultaneous source encoding
Simultaneous source encoding is a summing encoded source into super-shots [Romero et al. 2000, Schuster et al. 2011] conclude that the number of migration iterations is synchronous with the number of shot gathers in a super-shots gather to gain computational sped up.where, d ~denote the corresponding data, the random weighting vector ~=(~1, ~2, ~3... ~k) t ; such that M(~~t)=1/NR [Leeuwen and Herrmann 2013].
The corresponding supershot l ~and adjoint wavefield m ~is given by and the encoded misfit is given by .

Optimization method
Quality and rate of convergence of inversion depend significantly on the gradient method used [Pratt et al. 1998].The problem lies in the inability of gradient method to determine a reliable step length.We use the steepest decent method for optimization, where, a k is step length.Steepest-descent (i.e., the gradient '−Dq[m] k ') method is used for the misfit function at point m 0 .The curvature (second derivative) of the misfit function is the known as 'Hessian 'function.The perturbation model of Equation ( 4) gives the minimum misfit function in each iteration.

Numerical example: offshore Malay Basin model
We use a realistic synthetic gas clouds model which closely represent the shallow gas clouds environments of SE Asia (Figure 1) [Ghosh et al. 2010, Ghazali 2011].The main target is to enhance the properties of media where the acoustic (V p ) velocity is significantly low.In the gas clouds region, the reflections underneath the gas exhibit a 'push-down' effect due localization of the anomaly (Figure 1).Hence, full waveform inversion is right approach; however, finding a good velocity model is the major hurdle.
The work flow we adopted for this is as follows.
(3) We populate the model with Gaussian anomalies to add complexity to the data, note that in Gaussian anomaly, the effect of attenuation, scattering phenomenon is very dominant; however, the adjoint velocity perturbation is less affected.
(4) Finally, we obtain our initial velocity model (Figure 2B).A number of initial realizations are constructed by adding random variation around the model.
We use 12 iterations per frequency band, as larger iterations may lead to erroneous results and over data fitting [Li et al. 2013].
We consider three scenarios to illustrate the FWI results by (a) least-squares, (b) Student-t misfit and by (c) simultaneous source encoding on same model.Figure 2C shows the data at 3 Hz, 6 Hz, and 9 Hz in the source-receiver domain.

Discussion of results
The magnitude of the gradient with respect to model parameter is governed by the location of source and receivers array.The use of least-squares, Student's and simultaneous source encoding in frequency domain FWI highlight the sensitivity of the optimization convergence toward the true model.Hereinafter is our finding.
(a) In the first scenario, the application with leastsquares technique in a gas cloud model intrinsically amplifies the weight of the high residuals during inversion, hence causing divergence during the optimization.This  is often the case if the residuals do not correspond to signals (Figure 5a), causing the large outliers in the inversion result.
(b) On the other hand, Student-t method diminishes the influence of large outliers, which suggest a better choice when the quality of data is structurally complex (Figure 5b).
(c) In the third scenario, we implemented simultaneous source encoding FWI methods [Krebs et al. 2009].The advantage with simultaneous source encoding is an efficient computation at efficiency.Moreover, computational costs of full waveform inversion are a major issue to estimate the inversion result for large scale problem such as gas cloud.The basic idea of simultaneous source encoding is based on number of migration iterations should be equal to number of shot gathers in a super-shot gather.Result using simultaneous source encoding is shown in Figure 6, which is more efficient and effective in comparison to the conventional inversion.
Further, the results with Student-t penalty function and simultaneous source encoding are more informative than others (Figures 5c and 6b).We observe that the results of least-squares are apparently less sensitive to the data.However, using simultaneous source encoding, some artifacts are generated at top of the data (Figure 6a).
Further the resolution of the inversion results is still lower than the Student-t penalty but it's also well explain the subsurface structure of gas cloud model very well.Simultaneous source encoding inversion consists of sources that can reduce the computational cost by a factor proponational to the number sources stacked within a super-shots for forward modeling performed in the inversion procedure.On the other hand, simultaneous sources introduce the artifact (also known as crosstalk) noises in the final inversion result due to interference among the individual sources within a super-shot or due to less constraint of top of boundary during inversion.

Conclusion
(1) Least-squares algorithm is a fast and widely used algorithm, in seismic inversion; however it fails in converse hyperbolic behavior of the data from complex feature like gas cloud model (Figure 5a).( 2) Student-t approach provides a significant improvement in final inverted model (Figure 5b).(3) Source encoding waveform inversion, in which the sources are modeled simulta-   neously, is considerably faster and more accurate than the conventional inversion (Figure 6a,b).Combined misfit plot with true model (Figure 5c) showed the confidence and effectiveness of the proposed algorithm.(4) Our method recovers the amplitude losses caused by gas cloud masking with reduced the computational cost to speed up.(5) Theoretically, the error in the misfit and gradient decays much rapidly as the number of sources increases.( 6) We used l-BFGS optimization strategies for dealing with complex model to approximate the gradient.
While the results of the two different methods are very encouraging, there is still scope for improvement.An inversion result on a small complex synthetic model suggest that using randomly chosen sequential sources can be as efficient as simultaneous source encoding.A notable advantage is that it applies on marine seismic data, which allows us in practice as source estimation and source-dependent data-weighting with control over the error between the approximated and full misfit gradient.
Due to ill posed and nonlinear relationship between data and model, inversion needs to be iterated several times to converge towards the local minima of misfit function.This is particularly true, when working at low frequencies because theoretically at low frequency, there is a more probability of chance that modeled and observed data match within half a wave-cycle.
The basic formulation is as follows where, Dd = d obs − d k calc is the data misfit, the residual between the observed data d obs and the modeled data d k calc computed in the defined model space m k .The misfit function C(m) is updated with every iteration with initial velocity model m 0 and perturbation model Dm is given by where, t denote the transpose.The updated perturbation model is given by using , in Equation (3), we get This is Equation (3) of Warner et al. [2013], { is the Hessian matrix of second order differential equation of the gradient function C(m 0 ).The l 2 norm is usually written as: Pratt et al. [1998] defined as ∝, as step length,

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

Figure 4 .
Figure 4. (a) Shot gather time-domain seismic data; (b) model descriptions used for modeling and inversion.

Figure 5 .
Figure 5. (a) Full waveform inversion result using conventional approach (least-squares technique); (b) full waveform inversion result using the proposed approach; (c) misfit plot with the proposed approach (red color).The inversion results are well fit with the true model.