Wave Equation MVA

Introduction

Image gathers as a function of subsurface offset are an important tool for velocity analysis in areas of complex geology. The bottleneck in computing the image gathers force their computation for a limited number of subsufarce points and for a limited range of subsurface offsets. We offer a new persepective on image gathers where we organize the extended image as a function of all subsurface offsets for all subsurface points in a matrix whose \((i,j)^{\mathrm{th}}\) entry captures the interaction between gridpoints \(i\) and \(j\). We then propose an efficient way to probe the extended images for information for all subsurface offsets without explicitly calculating the source and receiver wavefields for all the sources.

Extended images

We arrange the time-harmonic source and receiver wavefields into the matrices \(U\) and \(V\), where each column corresponds to a source experiment. Given these matrices, the extended image at a single frequency, for all subsurface offsets and for all subsurface points can be written as the outer product of these two matrices, i.e., we have \[ \begin{equation*} E = VU^*, \end{equation*} \] where \(U,V\) are calculated using the two-way wave-equation. (Figure 1) illustrates how different subsurface-offset gathers are embedded in the 4-dimensional image volume.

In 2D, the extended image is a 5-dimensional function of all subsurface offsets and temporal shifts. So even in this case, it is prohibitively expensive to compute and store the extended image for all the subsurface points. To overcome this problem, we select \(l\) columns of \(E\) implicitly by multiplying this matrix with the tall matrix \(W = [{\bf{w}}_1, \ldots, {\bf{w}}_{l}]\) yielding, \[ \begin{equation*} \widetilde{E} = EW, \end{equation*} \] where \({\bf{w}}_i = [0, \ldots, 0, 1, 0,\ldots, 0]\) represents a single scattering point with the location of \(1\) corresponding to \(i^{th}\) grid location of the point scatterer. Now each column of \(\widetilde{E}\) represents a CIP gather at the locations represented by \({\bf{w}}_i\). We follow (Leeuwen and Herrmann, 2012) to efficiently compute \(\widetilde{E}\) by combining equations (1) and (2) as \[ \begin{equation*} \widetilde{E} = EW = H^{-*}P_r^TDQ^*P_sH^{-*}W, \end{equation*} \] where \(H\) is a discretization of the Helmholtz operator \((\omega^2 m + \nabla^2)\) with \(m\) is the squared slowness, The matrix \(Q\) represents the source function, \(*\) represents the conjugate-transpose, \(D\) is the reflected data matrix and \(P_s, P_r\) samples the wavefield at the source and receiver positions (and hence, their transpose injects the sources and receivers into the grid). We can compute this product efficiently as follows

• compute \(\widetilde{U} = H^{-*}W\) and sample this wavefield at the source locations \(\widetilde{D} = P_s\widetilde{U}\);
• correlate the result with source function \(\widetilde{Q} = Q^*\widetilde{D}\);
• use the result as data weights \(\widetilde{W} = D\widetilde{Q}\);
• inject the wavefield \(\widetilde{W}\) at the receiver locations and compute the extended image as \(\widetilde{E} = H^{-*}P_r^T\widetilde{W}\).

The computational cost of calculating \(\widetilde{E}\) is \(2l\) PDE solves plus the cost of correlating the source and data matrices. Thus, the cost of computing the CIPs does not depend on the number of sources or the number of subsurface offsets, as it does in the conventional methods for computing image gathers (Paul and Ivan, 2011). This is particularly beneficial when we are interested in computing only a few CIPs. An overview of the computational complexity is shown in Table 1.

AVA analysis and geological dip estimation

Because the extended images contain all possible subsurface offsets, we compute the angle-domain image gathers by selecting the subsurface offset that is aligned with the local dip. We do not need a priori information about the local geological dip to perform AVA analysis. (Kumar et al., 2013) showed that we can compute the local dip information on the fly from common-image-point gathers. To assess the quality of the angle-domain common-image-points gathers we compute the angle-dependent reflectivity coefficients and compare them with theoretical reflectivity coefficients yielded by the (linearized) Zoeppritz equations (Figure 2, 3).

Migration velocity analysis

Estimation of migration velocity is a crucial step in seismic imaging to get the accurate image of the subsurface. Wave-equation based migration velocity analysis (WEMVA) is one such tool for velocity estimation in complex geological settings. The objective of WEMVA is to measure the focusing of common-image gathers (CIG) by formulating an appropriate cost function (Biondi and Symes, 2004; Shen and Symes, 2008; Symes, 2008; Sava and Vasconcelos, 2011). In order to perform the automatic MVA, (Shen and Symes, 2008) consider penalizing the extended images as a function of subsurface offset with the lateral shift, which in our case is equivalent to demanding that our image matrix commutes with point-wise multiplication with the \(\bf{x}\)-position. For a focused image, we want to enforce that \[ \begin{equation*} E\mathrm{diag}({\bf{x}}) \approx \mathrm{diag}({\bf{x}})E. \end{equation*} \] The obvious way to achieve this is by formulating the optimization problem as \[ \begin{equation*} \min_{\bf{m}} \phi({\bf{m}}) = ||E({\bf{m}})\mathrm{diag}({\bf{x}}) - \mathrm{diag}({\bf{x}})E({\bf{m}})||^2, \end{equation*} \] where \(\| . \|\) is a matrix norm. If we choose the Frobenius norm as the matrix norm, then we can estimate this penalty via randomized trace estimation without explicitly forming the whole matrix as shown in (Kumar et al., 2014). Figure 4 shows the experimental results using proposed WEMVA formulation on Lens model.

References