Source Estimation by Full Wave Form Inversion

Abstract

We consider the inverse problem of estimating the parameters describing the source in a seismic event, using time-dependent ground motion recordings at a number of receiver stations. The inverse problem is defined in terms of a full waveform misfit functional, where the objective function is the integral over time of the weighted $L_2$ distance between observed and synthetic ground motions, summed over all receiver stations. The misfit functional is minimized under the constraint that the synthetic ground motion is governed by the elastic wave equation in a heterogeneous isotropic material. The seismic source is modeled as a point moment tensor forcing in the elastic wave equation. The source is described by 11 parameters: the six unique components of the symmetric moment tensor, the three components of the source location, the origin time, and a frequency parameter modeling the duration of the seismic event. The synthetic ground motions are obtained as the solution of a fourth order accurate finite difference approximation of the elastic wave equation in a heterogeneous isotropic material. The discretization satisfies a summation-by-parts (SBP) property that ensures stability of the explicit time-stepping scheme. We use the SBP property to derive the discrete adjoint of the finite difference method, which is used to efficiently compute the gradient of the misfit. A new moment tensor source discretization is derived that is twice continuously differentiable with respect to the source location. The differentiability makes the Hessian of the misfit a continuous function of all source parameters. We compare four different gradient-based approaches for solving the constrained minimization problem; two non-linear conjugate gradient methods (Fletcher–Reeves and Polak–Ribière), and two quasi-Newton methods (BFGS and L-BFGS). Because the Hessian of the misfit has a very large condition number, the parameters must be scaled before the minimization problem can be solved. Comparing several scaling approaches, we find that the diagonal of the Hessian provides the most reliable scaling alternative. Numerical experiments are presented for estimating the source parameters from synthetic ground motions in two different three-dimensional models; one in a simple layer over half-space, and one using a fully heterogeneous material. Good convergence properties are demonstrated in both cases.

Appendix 1: Proof of Theorem 1

We start by expanding the predictor into the corrector in Algorithm 1. Then rewrite the resulting expression as

$$\begin{aligned} \rho \frac{\mathbf{u}^{n+1}-2\mathbf{u}^n+\mathbf{u}^{n-1}}{\Delta _t^2} = \mathbf{L}_h(\mathbf{u}^n) + \mathbf{F}(t_n) + \frac{\Delta _t^2}{12}\bigl (\mathbf{L}_h(\mathbf{v}^n)+\mathbf{F}_{tt}(t_n)\bigr )+\mathbf{S}_G(\mathbf{u}^n-\mathbf{u}^{n-1}). \nonumber \\ \end{aligned}$$

(53)

Next, take the scalar product between (53) and ${\varvec{\kappa }}^n$, and sum over all time steps,

$$\begin{aligned}&\sum \limits _{n=0}^{M-1}\left( {\varvec{\kappa }}^n,\rho \frac{\mathbf{u}^{n+1}-2\mathbf{u}^n+\mathbf{u}^{n-1}}{\Delta _t^2}\right) _h = \sum \limits _{n=0}^{M-1}\left( {\varvec{\kappa }}^n,\mathbf{L}_h(\mathbf{u}^n)\right) _h + \sum \limits _{n=0}^{M-1} \left( {\varvec{\kappa }}^n,\mathbf{F}(t_n) + \frac{\Delta _t^2}{12}\mathbf{F}_{tt}(t_n)\right) _h \nonumber \\&\quad +\frac{\Delta _t^2}{12}\sum \limits _{n=0}^{M-1} \left( {\varvec{\kappa }}^n,\mathbf{L}_h(\mathbf{v}^n)\right) _h + \sum \limits _{n=0}^{M-1}\left( {\varvec{\kappa }}^n,\mathbf{S}_G(\mathbf{u}^n-\mathbf{u}^{n-1})\right) _h. \end{aligned}$$

(54)

The sum on the left hand side of (54) can be rewritten as

$$\begin{aligned}&\sum \limits _{n=0}^{M-1}\left( {\varvec{\kappa }}^n,\rho \frac{\mathbf{u}^{n+1}-2\mathbf{u}^n+\mathbf{u}^{n-1}}{\Delta _t^2}\right) _h =\sum \limits _{n=0}^{M-1}\left( \rho \frac{{\varvec{\kappa }}^{n+1}-2{\varvec{\kappa }}^n+{\varvec{\kappa }}^{n-1}}{\Delta _t^2}, \mathbf{u}^n \right) _h \nonumber \\&\quad +\frac{1}{\Delta _t^2}\left( \left( \rho \mathbf{u}^{-1},{\varvec{\kappa }}^0\right) _h - \left( \rho \mathbf{u}^0,{\varvec{\kappa }}^{-1}\right) _h +\left( \rho \mathbf{u}^M,{\varvec{\kappa }}^{M-1}\right) _h - \left( \rho \mathbf{u}^{M-1},{\varvec{\kappa }}^M\right) _h\right) , \end{aligned}$$

(55)

where the initial data $\mathbf{u}^0 = \mathbf{u}^{-1} = {\varvec{\kappa }}^M = {\varvec{\kappa }}^{M-1}=\mathbf{0}$ make

$$\begin{aligned} \sum \limits _{n=0}^{M-1}\left( {\varvec{\kappa }}^n,\rho \frac{\mathbf{u}^{n+1}-2\mathbf{u}^n+\mathbf{u}^{n-1}}{\Delta _t^2}\right) _h =\sum \limits _{n=0}^{M-1}\left( \rho \frac{{\varvec{\kappa }}^{n+1}-2{\varvec{\kappa }}^n+{\varvec{\kappa }}^{n-1}}{\Delta _t^2}, \mathbf{u}^n \right) _h. \end{aligned}$$

The first sum on the right hand side of (54) is treated by the self-adjoint property,

$$\begin{aligned} \sum \limits _{n=0}^{M-1}({\varvec{\kappa }}^n,\mathbf{L}_h(\mathbf{u}^n))_h = \sum \limits _{n=0}^{M-1}(\mathbf{L}_h({\varvec{\kappa }}^n),\mathbf{u}^n)_h. \end{aligned}$$

The second last sum of (54) can be rewritten

$$\begin{aligned} ( {\varvec{\kappa }}^n, \mathbf{L}_h(\mathbf{v}^n) )_h&= ( \mathbf{L}_h({\varvec{\kappa }}^n), \mathbf{v}^n )_h = \left( \mathbf{L}_h({\varvec{\kappa }}^n), \frac{1}{\rho }(\mathbf{L}_h(\mathbf{u}^n) + \mathbf{F}(t_n)) \right) _h \nonumber \\&= \left( \mathbf{L}_h({\varvec{\kappa }}^n), \frac{1}{\rho }\mathbf{L}_h(\mathbf{u}^n) \right) _h + \left( \mathbf{L}_h({\varvec{\kappa }}^n), \frac{1}{\rho }\mathbf{F}(t_n) \right) _h \nonumber \\&= \left( \mathbf{L}_h\left( \frac{1}{\rho }\mathbf{L}_h({\varvec{\kappa }}^n)\right) , \mathbf{u}^n \right) _h + \left( \frac{1}{\rho }\mathbf{L}_h({\varvec{\kappa }}^n), \mathbf{F}(t_n) \right) _h\nonumber \\&= (\mathbf{L}_h({\varvec{\zeta }}^n), \mathbf{u}^n )_h + ({\varvec{\zeta }}^n, \mathbf{F}(t_n) )_h. \end{aligned}$$

(56)

The super-grid damping term can be written

$$\begin{aligned}&\sum \limits _{n=0}^{M-1}({\varvec{\kappa }}^n,\mathbf{S}_G(\mathbf{u}^n-\mathbf{u}^{n-1}))_h = \sum \limits _{n=0}^{M-1}({\varvec{\kappa }}^n,\mathbf{S}_G(\mathbf{u}^n))_h - \sum \limits _{n=0}^{M-1}({\varvec{\kappa }}^{n+1},\mathbf{S}_G(\mathbf{u}^n))_h\nonumber \\&\quad -({\varvec{\kappa }}^0,\mathbf{S}_G(\mathbf{u}^{-1}))_h + ({\varvec{\kappa }}^M,\mathbf{S}_G(\mathbf{u}^{M-1}))_h. \end{aligned}$$

(57)

The boundary terms are zero because of the initial data $\mathbf{u}^{-1}={\varvec{\kappa }}^{M} = 0$, and we use the symmetry (9) to obtain

$$\begin{aligned} \sum \limits _{n=0}^{M-1}({\varvec{\kappa }}^n,\mathbf{S}_G(\mathbf{u}^n-\mathbf{u}^{n-1}))_h = \sum \limits _{n=0}^{M-1}({\varvec{\kappa }}^n-{\varvec{\kappa }}^{n+1},\mathbf{S}_G(\mathbf{u}^n))_h = -\sum \limits _{n=0}^{M-1}(\mathbf{S}_G({\varvec{\kappa }}^{n+1}-{\varvec{\kappa }}^{n}),\mathbf{u}^n)_h \end{aligned}$$

Collecting terms gives that (54) is equivalent to

$$\begin{aligned}&\sum \limits _{n=0}^{M-1}\left( \rho \frac{{\varvec{\kappa }}^{n+1}-2{\varvec{\kappa }}^n+{\varvec{\kappa }}^{n-1}}{\Delta _t^2}, \mathbf{u}^n \right) _h = \sum \limits _{n=0}^{M-1}(\mathbf{L}_h({\varvec{\kappa }}^n),\mathbf{u}^n)_h + \sum \limits _{n=0}^{M-1} \left( {\varvec{\kappa }}^n,\mathbf{F}(t_n) + \frac{\Delta _t^2}{12}\mathbf{F}_{tt}(t_n)\right) _h \nonumber \\&\quad +\frac{\Delta _t^2}{12}\sum \limits _{n=0}^{M-1} (\mathbf{L}_h({\varvec{\zeta }}^n), \mathbf{u}^n )_h + \frac{\Delta _t^2}{12}\sum \limits _{n=0}^{M-1}({\varvec{\zeta }}^n,\mathbf{F}(t_n))_h -\sum \limits _{n=0}^{M-1}(\mathbf{S}_G({\varvec{\kappa }}^{n+1}-{\varvec{\kappa }}^{n}),\mathbf{u}^n)_h. \end{aligned}$$

(58)

Expanding the predictor (11) into the corrector (12) gives

$$\begin{aligned} \rho \frac{{\varvec{\kappa }}^{n+1}-2{\varvec{\kappa }}^n+{\varvec{\kappa }}^{n-1}}{\Delta _t^2} = \mathbf{L}_h({\varvec{\kappa }}^n) + \mathbf{G}(t_n) + \frac{\Delta _t^2}{12}\mathbf{L}_h({\varvec{\zeta }}^n) - \mathbf{S}_G({\varvec{\kappa }}^{n+1}-{\varvec{\kappa }}^n). \end{aligned}$$

(59)

Identity (13) of Theorem 1 is obtained by inserting (59) into the left hand side of (58).