Migration/inversion & preconditioning

Wave equation modeling can be regarded as a linear operator, the domain of the operator is a the scattering field, and the range is a seismic dataset. When the scattering field and seismic data are discretized they can be represented as finite dimensional vectors. The modeling operator, A, maps from a scatterer vector, s to a data vector d.

As = d

When imaging we wish to obtain an estimate of s from the measured data d. One possible operation is to apply the adjoint operator,

$$\hat{s} = A'd$$

If the operator is nearly unitary this will be a good estimate of the model. When the operator is not unitary, as when our modeling operator is truncated, this will not be a good estimate. The estimate can be improved by solving the least squares inverse problem

$$\min_{\hat{s}} \vert\vert A \hat{s} - d \vert\vert^2$$

If a conjugate gradient algorithm is used the application of the adjoint can be regarded as the first iteration of that algorithm.

Another method for obtaining an estimate is to note that the solution of the least squares problem is given by

$$\hat{s} = (A'A)^{-1} A' d$$

Calculation of (A'A)^-1 is very expensive but we might obtain an acceptable result by approximating this inverse. If we assume A'A is diagonally dominant we can approximate the inverse by ( Diag(A'A) )^-1, where Diag extracts the diagonal of A'A and sets all other elements to zero. So the estimate of the solution is given by,

 

$$\hat{s} = ( Diag(A'A) )^{-1} A' d.$$ (8)

This procedure can also be interpreted as the first step of an iterative algorithm. In this case the problem being solved is a preconditioned version of the original problem. M = ( Diag(A'A) )^-1/2 is a diagonal preconditioner that equalizes the norm of the columns of the operator A. The preconditioned least squares problem is,

$$\begin{eqnarray} \min_{\hat{t}} \vert\vert AM \hat{t} - d \vert\vert^2 \\ \hat{s} = M t.\end{eqnarray}$$ (9) (10)

One step of the preconditioned iterative algorithm followed by application of the preconditioner to solve for $\hat{s}$ is exactly the same as equation (8).

In many applications it is relatively simple to calculate ( Diag(A'A) ). For the modeling algorithm used here it is,

$$(A'A)(x,z;h) = \int \int G^*( y,z,\omega; h ) G(x-y,z,\omega; h ) \ dy \ d\omega.$$

If the modeling is not limited in spatial aperture, the weight is just a function of depth and offset panel. It is a depth and offset dependent scale applied after adjoint modeling rather than the time dependent scale usually applied before migration in seismic processing.

Figure  shows the constant offset migrations calculated using this weight function. The amplitudes of the two scatterers have been correctly recovered. The diagonal weight only corrects the relative amplitudes of the result, it does not correct the spectrum of the result, the full inverse operator (A'A)^-1A' would correct the spectrum as well.

 

ata-mig
Figure 7 Constant offset migration/inversion in a V(z) medium.