Least-squares time-reversal imaging

Since backward wave propagation is not the exact inverse of forward wave propagation, the weighted time-reversal wavefield does not contain the correct phase and amplitude information of the true source. To solve this issue, we reformulate the problem of passive source imaging from the perspective of full waveform inversion (FWI) using the adjoint state method (Plessix, 2006; Virieux and Operto, 2009). In acoustic media, we can define the objection function as the data misfit measured by the $L_2$ norm:

$$J = \frac{1}{2}\vert\vert \mathbf{S}u(\mathbf{x},t) - d(\mathbf{x},t) \vert\vert _2^2 \;,$$ (7)

where $u(\mathbf{x})$ is the wavefield, $d(\mathbf{x},t)$ is the observed seismic data and $\mathbf{S}$ is the matrix that defines aquisition geometry. The corresponding state and adjoint state equations can be expressed as

$$(\frac{1}{c^2(\mathbf{x})}\,\frac{\partial^2}{\partial t^2} - \nabla^2) u(\mathbf{x},t) = f_s(\mathbf{x},t)\;,$$ (8)

$$(\frac{1}{c^2(\mathbf{x})}\,\frac{\partial^2}{\partial t^2} - \na... ...{x},t) = \mathbf{S}^{\intercal}(\mathbf{S}u(\mathbf{x},t) - d(\mathbf{x},t))\;.$$ (9)

where $\lambda(\mathbf{x})$ is the adjoint wavefield, $c(\mathbf{x})$ is the velocity and $f_s(\mathbf{x},t)$ is the source function . The gradients for source and velocity are respectively:

$$g_{f_s}(\mathbf{x},t) = \lambda(\mathbf{x},t) \;,$$ (10)

$$g_c(\mathbf{x}) = \frac{2}{c^3(\mathbf{x})} \int_0^T \lambda(\mathbf{x},t) \frac{\partial^2u(\mathbf{x},t)}{\partial t^2} dt \;.$$ (11)

Note that the gradient with respect to the seismic source is exactly the adjoint wavefield, i.e., the time-reversal wavefield modeled by the adjoint state equation. This explains why the conventional time-reversal imaging methods are appropriate for locating seismic sources. However, the difference between the true seismic source and a time-reversal wavefield is quite significant. First, the true seismic source $f_s(\mathbf{x},t)$ is sparse, while a time-reversal wavefield is not. Second, even at locations and times that correspond to true seismic sources, the time-reversal wavefield does not provide accurate amplitude and phase information.

To address this issue, we note that the wavefield is linearly dependent on the source term (Rickett, 2013). Rewriting equation 8 as a linear operator:

$$\mathbf{d} = \mathbf{A} \mathbf{f}_s \;,$$ (12)

where linear operator $A$ represents the combination of acoustic wave equation and acquisition geometry matrix $\mathbf{S}$. Given a starting velocity model, one can obtain accurate information about the source by solving the least-squares problem

$$\min\vert\vert \mathbf{d}-\mathbf{A}\mathbf{f}_s \vert\vert _2^2$$ (13)

using the pseudo-inverse

$$\hat{\mathbf{f}}_s = \mathbf{A}^{\dagger} \mathbf{d} = (\mathbf{A}^{\intercal}\mathbf{A})^{-1}\mathbf{A}^{\intercal} \mathbf{d} \;.$$ (14)

In the rest part of the paper, we will refer to the method defined by equation 14 as least-squares time-reversal imaging (LSTRI). Although LSTRI is a linear inversion, we recognize that the model space is actually four-dimensional, and includes the three spatial axes and one time axis. Without proper preconditioning, the computational cost of such an inversion can be impractical. We propose to use the weighting function derived in equation 5 as a diagonal model weighting matrix to precondition the conjugate gradient (CG) iterations. Since $W(\mathbf{x},t)$ effectively restricts CG to only update the locations that correspond to possible source locations, the eigenvalues of the matrix being inverted become better clustered and CG is able to achieve a faster convergence rate (Shewchuk, 1994).

The starting velocity model used to invert for the passive sources might be inaccurate. Most often, they are derived from well logs or travel-time tomography. FWI is capable of providing more detailed velocity structures. After initial source inversion, the inverted source can be used in FWI to update the velocity. Furthermore, one can jointly update the source and velcoty under the framework of the variable projection method (Aravkin et al., 2012; Golub and Pereyra, 1973; Rickett, 2013). In practice, we perform sufficient linear iterations to obtain accurate source functions and only a few nonlinear iterations to update the velocity model in an alternating fashion.