Tomographic full waveform inversion: Practical and computationally feasible approach

Scale Separation

Although the next derivation is done in frequency domain, the actual implementation is done in time domain. First, we start with the tomographic full waveform objective function $J_{\rm TFWI}$ which can be written as:

$$J_{\rm TFWI}(\mathbf v(\mathbf h)) = \lVert \mathbf d \left( \mat... ...rm obs} \rVert^2_2 + \epsilon \lVert \mathbf A \mathbf v(\mathbf h) \rVert^2_2,$$ (1)

where $\mathbf h$ is the offset, $\mathbf v(\mathbf h)$ is the extended velocity model, $\mathbf d (\mathbf v(\mathbf h))$ is the modeled data, $\mathbf d_{\rm obs}$ is the observed surface data, $\epsilon$ is a scalar weight of the regularization term and $\mathbf A$ is a regularization operator. The modeled data $\mathbf d (\mathbf v(\mathbf h))$ is computed as:

$$d(\mathbf x_s, \mathbf x_r, \omega; \mathbf v(\mathbf h)) = f(\ma... ... x_s, \mathbf x, \omega; \mathbf v(\mathbf h)) \delta(\mathbf x_r - \mathbf x),$$ (2)

where $f(\omega)$ is the source function, $\omega$ is frequency, $\mathbf x_s$ and $\mathbf x_r$ are the source and receiver coordinates, and $\mathbf x$ is the model coordinate. In the acoustic, constant-density case the Green's function $G(\mathbf x_s, \mathbf x, \omega; \mathbf v(\mathbf h))$ satisfies:

$$\left( v^{-2}(\mathbf x, \mathbf h) * \omega^2 + \nabla^2 \right) G(\mathbf x_s, \mathbf x, \omega) = \delta(\mathbf x_s - \mathbf x),$$ (3)

where $*$ denotes a convolution operator over the subsurface offset axis (Symes, 2008; Biondi and Almomin, 2012). The first simplification of the extended velocity model is to separate it into a background and a perturbation as follows:

$$v^{-2}(\mathbf x, \mathbf h) = b(\mathbf x, \mathbf h) + r(\mathbf x, \mathbf h),$$ (4)

where $b(\mathbf x)$ is the background component, which is a smooth version of the slowness squared and $r(\mathbf x)$ is the perturbation component. This separation assumes that $b(\mathbf x)$ will contain the transmission effects and $r(\mathbf x)$ will contain the reflection effects. Depending on the error of the initial background velocity, the perturbation component can extend across several subsurface offsets so it is important to keep its offset axis. On the other hand, the background component is not expected to generate reflections that would be grossly time shifted with respect to the recorded data, and it thus safe to restrict it to zero offset. A physical interpretation is that the wave speed is not expected to vary much across different angles, at least in the isotropic case that we are analyzing. Therefore, the extent of background component across subsurface offsets can be reduced. If the velocity is expected to vary, the same reduction can be applied while keeping more than zero subsurface offset. In our derivation, we reduce the background to only the zero subsurface offset as follows

$$v^{-2}(\mathbf x, \mathbf h) \approx b(\mathbf x) + r(\mathbf x, \mathbf h).$$ (5)

After the separation and reduction, the Born approximation can be used to linearize wave equation where the data is assumed to contain primaries only. The linearized wave equation defines the data as follows:

$$d(\mathbf x_s,\mathbf x_r,\omega; \mathbf b, \mathbf r(\mathbf h)... ... r(\mathbf x, \mathbf h) G(\mathbf x + \mathbf h,\mathbf x_r,\omega;\mathbf b),$$ (6)

where the Green's functions now satisfy the conventional acoustic wave equation as follows:

$$\left( \omega^2 b(\mathbf x) + \nabla^2 \right) G(\mathbf x_s,\mathbf x,\omega) = \delta(\mathbf x_s-\mathbf x),$$ (7)

$$\left( \omega^2 b(\mathbf x) + \nabla^2 \right) G(\mathbf x,\mathbf x_r,\omega) = \delta(\mathbf x-\mathbf x_r)..$$ (8)

The forward modeling can be written in a compact notation as follows:

$$\mathbf d = \mathbf L(\mathbf b) \mathbf r(\mathbf h),$$ (9)

where $\mathbf L$ is the Born modeling operator. We now define a new objective function for the efficient TFWI inversion as:

$$J_{\rm ETFWI}(\mathbf b, \mathbf r(\mathbf h)) = \lVert \mathbf L... ...rm obs} \rVert^2_2 + \epsilon \lVert \mathbf A \mathbf r(\mathbf h) \rVert^2_2.$$ (10)

Notice that there are similarity with the regularized linearized inversion proposed by Clapp (2005), with the important difference that here we are including both the background and the perturbation components as variables in this objective function. The Born modeling operator is linear with respect to perturbation but nonlinear with respect to the background component. Therefore, another linearization around the ``background" background is required to compute the gradient. First, we rewrite the background as the sum of two components as follows:

$$b(\mathbf x) = b_0(\mathbf x) + \Delta b(\mathbf x),$$ (11)

where $b_0(\mathbf x)$ is the current background model and $\Delta b(\mathbf x)$ is the perturbation of the background. The Born approximation is used again to linearize the $\mathbf L$ operator with respect to the background resulting in a data-space tomographic operator. The data perturbation with respect to the background perturbation is now defined as:

$$\ \Delta d(\mathbf x_s,\mathbf x_r,\omega; \mathbf b_0, \mathbf r(\mathbf h)) = \sum_{\mathbf x,\mathbf y,\mathbf h}$$

$$\omega^4 f(\omega) G(\mathbf x_s,\mathbf y-\mathbf h,\omega;\math... ...;\mathbf b_0) \Delta b(\mathbf x) G(\mathbf x,\mathbf x_r,\omega;\mathbf b_0) +$$

$$\omega^4 f(\omega) G(\mathbf x_s,\mathbf x,\omega;\mathbf b_0) \D... ... r(\mathbf y, \mathbf h) G(\mathbf y+\mathbf h,\mathbf x_r,\omega;\mathbf b_0),$$ (12)

where $\mathbf y$ is the perturbation coordinates. As we can see in the previous equation, the tomographic operator correlates a background and a scattered wavefield from both source and receiver sides. The scattered wavefields are computed by correlating a background wavefield with the reflectivity model $\mathbf r$ and then propagating again to all model locations. The forward tomographic operator can be written in a compact notation as follows:

$$\Delta \mathbf d = \frac{\partial \mathbf L}{\partial \mathbf b} \mathbf r(\mathbf h) \Delta \mathbf b = \mathbf T \Delta \mathbf b,$$ (13)

Where $\mathbf T$ is the tomographic operator that relates changes in the background model to changes in the data. We can now compute the reflectivity gradient as follows:

$$g_{\mathbf r}(\mathbf x, \mathbf h) = \frac{\partial J}{\partial \mathbf r(\mathbf h)} = \mathbf L^* \Delta \mathbf d,$$ (14)

which is simply migration of the residuals. Then, we compute the background gradient as follows:

$$g_{\mathbf b}(\mathbf x) = \frac{\partial J}{\partial \mathbf b} = \mathbf T^* \Delta \mathbf d.$$ (15)

It is important to notice that reflectivity has two roles in computing the background model gradient. First, it is used to computed the data residuals. Second, it is used to scatter the background wavefields and compute the scattered wavefields of the tomographic operator.

Tomographic full waveform inversion: Practical and computationally feasible approach