Target-oriented least-squares migration/inversion with sparseness constraints

introduction

Migration is an important and robust tool for imaging subsurface structures using reflection seismic data. However, migration operator is only the adjoint of the forward Born modeling operator (Lailly, 1983), which produces reliable structural information of the subsurface (assuming an accurate background velocity is known), but blurs the image because of the non-unitary nature of the Born modeling operator. To deblur the migrated image and correct the effects of limited acquisition geometry, complex overburden and bandlimited wavefields, the imaging problem can be treated as an inverse problem, which, instead of using the adjoint operator, uses the pseudo-inverse of the Born modeling operator to optimally reconstruct the reflectivity. This inversion-based imaging mehtod is also widely known as least-squares migration (Clapp, 2005; Nemeth et al., 1999; Valenciano, 2008; Kuhl and Sacchi, 2003).

The standard least-squares migration/inversion (LSI) approach tries to minimize an objective function defined in the data space, which compares the mismatch between the modeled and the observed primaries (Clapp, 2005; Nemeth et al., 1999; Kuhl and Sacchi, 2003). The objective function is then minimized with a gradient-based optimization solver, which iterates until an acceptable image is obtained. However, the data-space inversion scheme lacks flexibility and cannot be implemented in a target-oriented fashion. Full-domain migration/demigration has to be carried out within each iteration; and the optimization converges slowly without a proper preconditioner. Therefore, the data-space inversion scheme is computationally challenging for large-scale applications.

One way to reduce the computational cost is by solving the LSI problem in a target-oriented fashion (Yu et al., 2006; Valenciano, 2008). This can be achieved by minimizing an objective function defined in the model space, instead of the data space. The target-oriented model-space formulation allows us to invert only areas of particular interest, such as subsalt regions, where potential reservoirs are located and migration often fails to provide reliable images. Solving the LSI in the model space requires explicitly computing the Hessian, the normal operator of the forward Born modeling operator. The full Hessian, however, is expensive to compute without certain approximations. Fortunately, as demonstrated by Valenciano (2008) and Tang and Biondi (2009), for a typical conventional acquisition geometry (shot records do not interfer), the Hessian matrix is sparse and diagonally dominant for most areas. Thus a truncated Hessian with a limited number of off-diagonal elements (the number is usually very small) can be used to approximate the exact Hessian for inverse filtering.

The truncated Hessian can be computed by storing the Green's functions (Valenciano, 2008), which, however, may bring considerable computational issues (e.g. disk storage, I/O and etc.), because the Green's functions can be huge for practical applications, especially in 3-D. To reduce the computational overburden, this paper computes the Hessian using the phase-encoding method (Tang, 2008b). As demonstrated by Tang (2008b) and Tang (2008a), computing the phase-encoded Hessian does not require storing any Green's functions and it is also more efficient: the cost for computing the receiver-side randomly phase-encoded Hessian is about one shot-profile migration, and if a mixed simultaneous phase-encoding strategy is used, the cost is about one plane-wave source migration.

Besides the computational cost, two main issues, i.e., the operator mismatch and the underdetermined nature of the seismic inverse problem, make the practical application of LSI less effective. The first issue often arises when our modeling operator is not sufficient to predict the physics of the data, for example, anisotropy or elasticity presents in the data but is not accurately modeled by our numerical operators. This can cause data-inconsistency problems. The second issue is due to the limited surface seismic acquisition geometry, which makes the inversion have an infinite number of solutions that fit the observed data equally well. Regularization is therefore important to stabilize the inversion and make it converge to geologically reasonable solutions. In this paper, I exploit the application of a non-quadratic regularization operator that imposes sparsness to the model space (Sacchi and Ulrych, 1995; Ulrych et al., 2001). The model-space sparsity is achieved by minimizing the model residual in the $\ell_1$ or Cauchy norm, whose distribution is longer-tailed than the Gaussian distribution (the $\ell _2$ norm), hence it penalizes weak energy and leads to spiky solutions (Amundsen, 1991). The application of the sparseness constraint to seismic imaging has also been reported by Tang (2006) and Wang and Sacchi (2007), who use it to regularize prestack image gathers. In this paper, however, I use it to regularize the prestack image (zero subsurface offset) to enhance the resolution of the inverted reflectivity. I compare the one-way wave-equation inversion results on the Marmousi model regularized using the sparseness constraint with those regularized using a standard $\ell _2$ norm damping. The experiments have been carried out on data sets synthesized using both one-way wave-equation Born modeling and two-way acoustic wave-equation finite-difference modeling. The first case represent the ideal scenario, where our modeling operator (one-way wave-equation propagator for this case) matches all the physics in the data. I show that both inversion schemes work well under this situation, and sparseness constrained inversion can offer slightly higher resolution. The second case is much more challenging for both schemes, because our modeling operator can not model all the complexities present in the data (e.g., amplitudes, multiples and etc.). My experiments show that under this difficult situation, the sparseness-constrained approach provides us a better inversion result for the Marmousi model, which suggests the importance of accurate model covariance (or the a priori information) for the LSI problem.

This paper is organized as follows: I first briefly review the theory of target-oriented LSI and phase-encoded Hessian, then I discuss the sparseness constraint which minimizes the model residual in the Cauchy norm. Finally I apply the regularized inversion scheme to the Marmousi model.

Target-oriented least-squares migration/inversion with sparseness constraints