A graphics processing unit implementation of time-domain full-waveform inversion

Nonlinear conjugate gradient method

The conjugate gradient (CG) algorithm decreases the misfit function along the conjugate gradient direction:

$$\textbf{d}_k= \begin{cases}-\nabla E(\textbf{m}_0), & k=0\\ -\nabla E(\textbf{m}_k)+\beta_k \textbf{d}_{k-1}, & k\geq 1 \end{cases}$$ (6)

There are a number of ways to compute $\beta_k$ . We use a hybrid a hybrid scheme combing Hestenes-Stiefel method and Dai-Yuan method (Hager and Zhang, 2006)

$$\beta_k=\max(0, \min(\beta_k^{HS},\beta_k^{DY})).$$ (7)

in which

$$\left\{ \begin{split}\beta_k^{HS}&=\frac{\langle\nabla E(\textb... ...xtbf{m}_k)-\nabla E(\textbf{m}_{k-1})\rangle} \end{split} \right.$$ (8)

This provides an automatic direction reset while avoiding over-correction of $\beta_k$ in conjugate gradient iteration. It reduces to steepest descent method when the subsequent search directions lose conjugacy. The gradient of the misfit function w.r.t. the model is given by (Bunks et al., 1995)

$$\nabla E_{\mathrm{m}} =\frac{2}{v^3(\textbf{x})}\sum_{r=1}^{ng}\s... ...;\textbf{x}_s)}{\partial t^2}p_{res}(\textbf{x}_r,t;\textbf{x}_s)\mathrm{d}t\\$$ (9)

where $p_{res}(\textbf{x},t;\textbf{x}_s)$ is the back propagated residual wavefield, see the Appendix B and C for more details. A Gaussian smoothing operation plays an important role in removing the noise in the computed gradient. A precondition is possible by normalizing the gradient by the source illumination which is the energy of forward wavefield accounting for geometrical divergence (Bai et al., 2014; Gauthier et al., 1986):

$$\nabla E(\mathrm{m}_k)=\frac{\nabla E_{\mathrm{m}}}{\sqrt{\sum_{s=1}^{ns}\int_{0}^{t_{\max}} p_{cal}^2(x,t;x_s)\mathrm{d}t+\gamma^2}}$$ (10)

where $\gamma$ is a stability factor to avoid division by zero. To obtain a reasonable step size $\alpha_k$ in equation (5), we estimate a small step length $\epsilon$ proposed by Pica et al. (1990):

$$\max(\epsilon \vert\textbf{d}_k\vert)\leqslant \frac{\max(\vert\textbf{m}_k\vert)}{100}.$$ (11)

and the Taylor approximation

$$\textbf{J}_k\textbf{d}_k\approx\frac{\textbf{f}(\textbf{m}_k+\epsilon \textbf{d}_k)-\textbf{f}(\textbf{m}_k)}{\epsilon}$$ (12)

We summarize the FWI flowchart in Figure 1.

flowchart
Figure 1. Backward reconstruction can be realized using the saved boundaries. Note that no absorbing boundary condition is applied on the top boundary of the model in the forward modeling.

A graphics processing unit implementation of time-domain full-waveform inversion