Consensus optimization of total variation–based reverse time migration

Abstract

A common method of imaging with seismic data is reverse time migration (RTM) in which recorded data are back-propagated to form an image of the subsurface. Least-squares RTM (LSRTM) extends this method to an iterative method that minimizes a data misfit term. An appropriate regularization term such as a sparsity-promoting functional (e.g., total variation (TV)) is required to stabilize the LSRTM solution. In this paper, in order to efficiently solve such regularized LSRTM via distributed optimization algorithms, we first reformulate the problem into a consensus form. The alternating direction method of multipliers (ADMM) allows us to split the problem into separate subproblems, resulting to a two-step iteration which efficiently solves the original problem. The most time-consuming step is due to LSRTM of the data associated with each source, separately, which is carried out in parallel via a set of workers. The resulting local images are sent to a master that is responsible for generating (via a proximal mapping) a unique global image that is close to the mean of local images while minimizing the regularizing function. This loop gives a general framework for the parallel computation of regularized LSRTM and allows different regularizing functions to be employed by using appropriate proximal operators. We demonstrate the performance of the proposed algorithm with a set of numerical examples. The results show that TV methods can generate more accurate images compared with the l1-norm sparse LSRTM, l2-norm LSRTM, and traditional RTM. Furthermore, the new algorithm shows improvements over the traditional formulations due to the consensus parallelization.

Appendix

Steps of computing the imaging condition

Computation of the imaging condition for each source includes performing the following three steps [30, 38]:

• 1. Computation of the forward wavefield, u = A^− 1b, which is obtained by forward time recursion:

  $$ u^{n} = 2u^{n-1} - u^{n-2} + (c_{0} T)^{2} \left( \nabla^{2} u^{n-1} - b^{n-1}\right) $$

  (27)

  for n = 0, 1,...,N − 1 where uⁿ ≡ u(x,nT).

• 2. Computation of the backward wavefield $v=A^{-T}P^{T}\ddot {d}$ which is the solution of the backward PDE:

  $$ A^{T}v=P^{T}\ddot{d} $$

  (28)

  where $\ddot {d}\equiv (D^{T} \otimes I)d$. Here, $P^{T}\ddot {d}$ serves as the source term in which P^T back-projects the second derivative of the data to the subsurface model [10]. One can think of $P^{T}\ddot {d}$ as a wavefield which is zero except at receiver positions where it takes a value equal to $\ddot {d}$. Since $A:=(\nabla ^{2} - \frac {1}{{c_{0}^{2}}(\bold {x})} \frac {\partial ^{2}}{\partial t^{2}})$, we get that $A^{T}:=(\nabla ^{2} - \frac {1}{{c_{0}^{2}}(\bold {x})} (\frac {\partial ^{2}}{\partial t^{2}})^{T})$ (note that ∇² is symmetric). Thus, the PDE form of Eq. 28 is

  $$ \nabla^{2} v(\bold{x},t) - \frac{1}{{c_{0}^{2}}(\bold{x})} (\frac{\partial^{2}}{\partial t^{2}})^{T} v(\bold{x},t) = \{P^{T}\ddot{d}\}(\bold{x},t). $$

  (29)

  Since ∂²/∂t² is approximated by the lower triangular matrix D as in Eq. 6, its transpose is an upper triangular matrix as follows:

  $$ (\frac{\partial^{2}}{\partial t^{2}})^{T} \approx D^{T}:=\frac{1}{T^{2}} \left[\begin{array}{cccccc} 1 & -2 &1& 0& {\cdots} & 0 \\ 0 & 1 &-2&1& {\cdots} & 0 \\ {\vdots} & & {\ddots} &{\ddots} & {\ddots} & {\vdots} \\ 0 & {\cdots} & 0 &1& -2 & 1 \\ 0 & {\cdots} & 0 & 0&1 & -2 \\ 0 & {\cdots} &0& 0 & 0 & 1 \end{array}\right]. $$

  (30)

  Accordingly, the backward wavefield defined by Eq. 29 can also be solved by the recursive (1), but, in this case, the recursion is performed backward in time N − 1,N − 2,..., 0 [30, 38].

• 3. Multiplication of the two wavefields obtained from steps 1 and 2 and summing over time (due to the summation operator e^T ⊗ I).