Anisotropic tomography with rock physics constraints

Theory

Here, we define a VTI tomography in using event traveltime as data. To the first order, when we assume the ray path does not change with perturbation in the background velocity model, the traveltime misfit is related to the velocity perturbation as follows:

$$\Delta t = \int_{\rm raypath} ~ \frac{1}{v^2} \Delta v \operatorname{d}l.$$ (1)

For weak anisotropy, the angle dependent velocity in VTI media can be described by vertical velocity $v_0$ and the Thomson anisotropic parameters $\epsilon$ and $\delta$ as follows (Tsvankin and Thomsen, 1994):

$$v(\theta) = v_0(1 + 2\delta \sin^2\theta \cos^2\theta + 2\epsilon \sin^4\theta),$$ (2)

where $\theta$ is the angle between propagation direction and the vertical direction. The angles are defined by the raypath in equation 1. Consequently, the velocity perturbation can be expressed as follows:

$$\Delta v(\theta) = \frac{\partial v(\theta)}{\partial v_0} \Delta... ...on} \Delta \epsilon + \frac{\partial v(\theta)}{\partial \delta} \Delta \delta.$$ (3)

Now we can see explicitly that the traveltime misfit has three contributions from $v_0$ , $\epsilon$ and $\delta$ , respectively:

$$\Delta t = \int \frac{1}{v^2(\theta)} \frac{\partial v(\theta)}{\... ...a)} \frac{\partial v(\theta)}{\partial \delta} \Delta \delta \operatorname{d}l.$$ (4)

In practice, we usually formulate the seismic tomography problem as an inversion problem with another model space regularization term. Then the objective function reads:

$$S(\Delta m) = \vert\vert \mathbf{T} \Delta m - \Delta t \vert\vert _D^2 + \vert\vert m_0 + \Delta m - m_{\rm prior} \vert\vert _M^2,$$ (5)

where $D$ denotes the data space whose covariance is $\mathbf{C_D}$ and $M$ denotes the model space whose covariance is $\mathbf{C_M}$ . In the data space fitting goal, operator $\mathbf{T}$ is a row operator vector $[\mathbf{T}_{v_0} ~ \mathbf{T}_{\epsilon} ~ \mathbf{T}_{\delta}]$ , each element corresponding to one of the three terms in equation 4 respectively; $\Delta m$ is a column model vector $[\Delta v_0 ~ \Delta \epsilon ~ \Delta \delta]'$ ; $\Delta t$ is the traveltime misfit. In the model space fitting goal, $m_0$ is the background (initial) model $[v_{0_{0}},\epsilon_0,\delta_0]'$ ; $m_{\rm prior}$ is the mean of the prior distribution in the model space $[v_{0p},\epsilon_p,\delta_p]'$ . We evaluate both fitting goals in L2 using an LSQR solver.

We can also obtain the normal equation representation of the objective function by taking the derivative of equation 5 with respect to $\Delta m$ :

$$\begin{vmatrix}\mathbf{C_D^{-1/2}} \mathbf{T} \\ \lambda ~ \mathb... ... \Delta t \\ \lambda ~ \mathbf{C_M^{-1/2}} (m_{\rm prior} - m_0), \end{vmatrix}$$ (6)

where $\lambda$ is the balancing factor between the data fitting equation and the model fitting equation. For the synthetic study here, we may assume $\mathbf{C_D} = \mathbf{I}$ , since there is no noise in these data.

To speed up the convergence, we turn this regularized problem into a preconditioned problem by doing a variable substitution:

$$\begin{vmatrix}\mathbf{T} \mathbf{C_M^{1/2}} \\ \lambda ~ \mathbf... ...}\Delta t \\ \lambda ~ \mathbf{C_M^{-1/2}} (m_{\rm prior} - m_0) \end{vmatrix}.$$ (7)

Now LSQR inverts the preconditioning variable $\Delta m'$ , and the real model updates are obtained by $\Delta m = \mathbf{C_M^{1/2}} \Delta m'$ after inversion.

For a multi-parameter model inversion, we assume the model covariance can be fully decoupled into two different parts: spatial covariance for each parameter, and the local cross-parameter covariance. Mathematically, that translates into:

$$\mathbf{C_M^{1/2}} = \mathbf{S} \mathbf{\sigma} = \left\vert \beg... ...} & \sigma_{\delta \epsilon} & \sigma_{\delta \delta} \end{array} \right\vert ,$$ (8)

where $\mathbf{S}$ is a spatial smoothing matrix, and $\mathbf{\sigma}$ is the square-root of the point-by-point cross-parameter covariance matrix. Each element of $\mathbf{S}$ can be estimated according to the user's assumption for the smoothness of different parameters. For example, the steering filters (Clapp et al., 2004; Fomel, 1994) provide a good choice to incorporate the structural information of the subsurface. The local cross-parameter covariance can be estimated from the rock physics modeling (Bachrach, 2010a).

covariance
Figure 1. Left: Stochastic realizations of a rock physics modeling for shale and sandy shale. The ellipsoids are fitted locally to describe the multi-Gaussian relationship among the parameters. Right: Small dots in the background are the estimated prior distribution by the operator using the ellipsoids on the left. Larger dots are the null-space projection of the operator.

The left panel on Figure 1 shows an example of the stochastic realizations of a rock physics model for shale and sandy shale. The dots are scattered $v_0$ , $\epsilon$ and $\delta$ results of hundreds of realizations, while the nine ellipsoids are fitted locally with $v_0$ the controlling variable to describe the multi-Gaussian relationship among the parameters. With a statistical study of the laboratory measurements on the cores, the range of the required input parameters for the rock physics model has been limited to a relatively narrow range, which leads to tight ellipsoids in Figure 1. We refer interested readers to details in Bachrach (2010a).

In this paper, we build the cross-parameter covariance matrix for each grid point in the subsurface according to the velocity, and linearly interpolate between the ellipsoids to which the velocity value belongs. Notice that additional non-linearity has been introduced during this process. Better methods to utilize the stochastic realizations could fit ellipsoids centered at every velocity value, which can be more precise but outside the scope of this study.

Anisotropic tomography with rock physics constraints