Extended imaging

Introduction

The formation of subsurface-offset gathers, such as common-image gathers (CIGs) (Rickett and Sava, 2002; W. W. Symes, 2008; Stolk et al., 2009; Kalita and Alkhalifah, 2016), angle-domain common-image gathers (ADCIGs) (de Bruin et al., 1990; Kumar et al., 2018), and common-image point gathers (CIPs) (P. Sava and Vasconcelos, 2011; Leeuwen and Herrmann, 2012) has become an essential component of modern seismic imaging workflows. Each gather provides crucial information on the velocity model and the scattering (amplitude versus offset). Contrary to CIGs, CIPs provide information on the complete scattering mechanism because they are functions of the full omni-directional subsurface offset (see Figure 1). Because of the high computational costs, subsurface-offset CIP gathers are often only evaluated at limited numbers of subsurface points and limited ranges of subsurface offsets. At SLIM, we overcome these limitations by offering a new perspective on image gathers by forming extended image volumes (EIVs, (Leeuwen et al., 2017a)), which are quadratic in the image size, as a function of all subsurface offsets for all subsurface points. We organize these EIVs in a matrix whose \((i,j)^{\mathrm{th}}\) entry captures the interaction between gridpoints \(i\) and \(j\). We then propose efficient ways to probe and factorize these EIVs providing us with information for all subsurface midpoints (\(\vec{m}\)) and offsets (\(\vec{h}\)) \[ \begin{equation} \vec{m}=\frac{\vec{x}+\vec{x}^\prime}{2}\quad \text{and}\quad \vec{h}=\frac{\vec{x}-\vec{x}^\prime}{2}, \label{eq:offset} \end{equation} \] without explicitly calculating the source and receiver wavefields for all sources.

Full subsurface extended image volumes

To arrive at our formulation, 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 EIV at a single frequency, for all subsurface offsets and for all subsurface points can be written as the outer product between these two matrices, i.e., we have \[ \begin{equation} E = VU^*, \label{eq:E} \end{equation} \] where \(U,V\) are calculated using the two-way wave-equation. Figure 2 illustrates how the different subsurface-offset gathers are embedded in the 4-dimensional EIV.

In 2D, the EIVs are 5-dimensional functions 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, \label{eq: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 Common Image Point (CIP) gather at the locations represented by \({\bf{w}}_i\). We follow (Leeuwen and Herrmann, 2012) to efficiently compute \(\widetilde{E}\) by combining equations \(\ref{eq:E}\) and \(\ref{eq:EW}\) 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 \mathrm{diag}(\mathbf{m}) + \nabla^2)\) with \(\mathbf{m}\) the squared slowness and \(\nabla^2\) the Laplacian. The matrix \(Q\) represents the source function, \(*\) the conjugate-transpose, \(D\) is the data matrix and \(P_s, P_r\) sample 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 for conventional methods used to form CIG’s (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.

Below, we list a number of applications of EIVs.

AVA analysis and geological dip estimation

Because EIVs contain all possible subsurface offsets, its possible to correct angle-domain CIGs ofr the local geologic dip by selecting the appropriate subsurface offsets. To accomplish this, we do not need a priori information on local geological dip to perform AVA analysis. (Kumar et al., 2013; Leeuwen et al., 2017b) showed that local dip information can be computed on the fly from CIG’s. 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 3, 4). These Figures demonstrate that these dip corrections are feasible and can lead to significant improvements in the CIG’s as shown in Figure 5

Migration velocity analysis

To achieve accurate high-fidelity images of the subsurface, estimation of the migration velocity model is crucial. 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 (CIGs) by formulating an appropriate cost function (Biondi and Symes, 2004; Shen and Symes, 2008; W. W. Symes, 2008; P. Sava and Vasconcelos, 2011). In order to perform the automatic MVA, (Shen and Symes, 2008) considered penalizing the extended images as a function of subsurface offset with lateral shifts, which in our formulation is equivalent to demanding that our EIV matrix commutes with point-wise multiplication with the \(\bf{x}\)-position. For a focused image (Leeuwen et al., 2017b), 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; Leeuwen et al., 2017b). Figure 6 shows the experimental results using proposed WEMVA formulation on Lens model.

Low-rank factorization

While the formation of EIVs allowed for new and enabling perspectives on seismic imaging, migration velocity analyses, and redatuming (Kumar et al., 2019) computational challenges remained. To address these, Yang et al. (2021) proposed a low-rank factorizing for EIVs based on randomized SVDs. This LR factorization provides access to the migrated image, CIPs, and CIG all from these factors alone. Given the number of probing vectors \(n_p\) and \(q\), this randomized SVD involves the steps outlined in Algorithm 1.

Velocity continuation

Given the L, R factorization, invariance-properties of EIVs (Kumar et al., 2019) can be used to map EIVs, or factors thereof, for one velocity model to another without the need to refactor—i.e., the factors for velocity model \(2\) can be related to these of velocity model \(1\) via \[ \begin{equation} \begin{aligned} \mathbf{L}_{2} &= \mathbf{H}_{2}^{-*} \cdot \mathbf{H}_{1}^{*} \cdot \mathbf{L}_{1}, \\ \\ \mathbf{R}_{2} &= \mathbf{H}_{2}^{-1} \cdot \mathbf{H}_{1} \cdot \mathbf{R}_{1}. \end{aligned} \label{eq:factors} \end{equation} \] involving fewer applications of the wave equation and its inverse (Yang et al., 2021). Based on these expressions, we are able to map the LR factors between the velocity models plotted in Figure 7. Given the new factors,

References