Wave Equation MVA


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.

Figure1Different slices through the 4-dimensional image volume \(e(z,z',x,x')\) around \(z=z_0\) and \(x=x_0\). (a) Conventional image \(e(z,z,x,x)\), (b) Image gather for horizontal and vertical offset \(e(z,z',x_0,x')\), (c) Image gather for horizontal offset \(e(z,z,x_0,x')\) and (d) Image gather for a single scattering point \(e(z_0,z',x_0,x')\). (e-g) shows how these slices are organized in the matrix representation of \(e\).}

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.

# of PDE solves flops
conventional \(2 N_s\) \(N_s N_{h_x} N_{h_y} N_{h_z}\)
this work \(2 N_x\) \(N_r N_s\)
Table1Computational complexity of the two schemes in terms of the number of sources \(N_{\mathrm{s}}\), receivers \(N_{\mathrm{r}}\), sample points \(N_{x}\) and desired number of subsurface offsets in each direction \(N_{h_{\{x,y,z\}}}\).

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).

Figure2Estimation of local geological dip. (a) Dipping one-layer model overlay CIP location at \(x = 1750\) m and \(z = 460\) m indicated by star. (b) CIP gather overlay on dipping model. (c) Stack-power versus dip-angle. We can see the maximum stack-power corresponds to the dip value of \(12^\circ\), which is the true dip value.
Figure3Modulus of angle-dependent reflectivity coefficients in dipping one-layer model at \(z = 460\) m and \(x = 1750\) m. Reflectivity coefficients (a) with no dip \(\theta=0^\circ\) and (b) with the dip obtained via the method described above (\(\theta = 12^\circ\)).

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.

Figure4Lens model.(a,b,c) True model, initial model and inversion results after 30 iterations. (d,e,f) CIG along x=1500 and z for true model, initial model and inversion results. (g) Vertical trace profile extracted along x=1500m.


Biondi, B., and Symes, W. W., 2004, Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging: Geophysics, 69, 1283. Retrieved from http://sepwww.stanford.edu/sep/biondo/PDF/Journals/Geophysics/adcig-geo-online.pdf

Kumar, R., Leeuwen, T. van, and Herrmann, F. J., 2013, AVA analysis and geological dip estimation via two-way wave-equation based extended images: SEG technical program expanded abstracts. SEG. Retrieved from https://slim.gatech.edu/Publications/Public/Conferences/SEG/2013/kumar2013SEGAVA/kumar2013SEGAVA.pdf

Kumar, R., Leeuwen, T. van, and Herrmann, F. J., 2014, Extended images in action: Efficient WEMVA via randomized probing: EAGE. Retrieved from http://www.earthdoc.org/publication/publicationdetails/?publication=76362

Leeuwen, T. van, and Herrmann, F. J., 2012, Wave-equation extended images: Computation and velocity continuation: EAGE technical program. EAGE. Retrieved from https://slim.gatech.edu/Publications/Public/Conferences/EAGE/2012/vanleeuwen2012EAGEext/vanleeuwen2012EAGEext.pdf

Paul, S., and Ivan, V., 2011, Extended imaging conditions for wave-equation migration: Geophysical Prospecting, 59, 35–55. Retrieved from http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2478.2010.00888.x/abstract

Sava, P., and Vasconcelos, I., 2011, Extended imaging conditions for wave-equation migration: Geophysical Prospecting, 59, 35–55. Retrieved from http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2478.2010.00888.x/abstract

Shen, P., and Symes, W. W., 2008, Automatic velocity analysis via shot profile migration: Geophysics, 73, VE49–VE59. Retrieved from http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2478.2010.00888.x/abstract

Symes, W. W., 2008, Migration velocity analysis and waveform inversion: Geophysical Prospecting, 56, 765–790. Retrieved from http://www.trip.caam.rice.edu/reports/2007/amva2.pdf