LEAST-SQUARES IMAGING

The exploding reflector model provides the linear relation between the subsurface image and the wavefield on the surface. In forward modeling, the input is the subsurface image and the output is the wavefield on the surface. The migration is the conjugate operation to the forward modeling (Claerbout, 1992a). If the modeling operator is a unitary operator, the migration operator, which is the conjugate of the modeling operator, will backproject the wavefield to the subsurface image, which results in the original input of the modeling. However, in reality the forward operator is the continuous wave equation for a given earth model, and the backprojection operator, which is migration operator, is not the conjugate to the forward operator. Instead we use an operator that is conjugate to an approximated forward operator like the Kirchhoff, phase shift, and finite-difference method. This discrepancy along with the finite and the discrete limits produces many artifacts on the image after migration.

The forward modeling can be simply formulated as

$${\bf d = SW m }$$ (1)

where ${\bf W}$ is the wave equation operator that does the forward modeling, ${\bf m}$ is the continuous model space, ${\bf S}$ is the truncating and sampling operator that simulates the finite and discrete seismic survey, and ${\bf d}$ is the data we obtain in reality. In conventional migration, we approximate the forward operator ${\bf SW}$ as

$${\bf S\tilde W \simeq SW}$$ (2)

where ${\bf \tilde W}$ represents an approximated forward modeling operator. Therefore, the conventional migration actually solve the problem as follows

$${\bf \hat m \simeq \tilde W'S' d} {\bf = \tilde W'S'SW m}$$ (3)

Even though the sampling interval is small and the aperture is large, which makes the operator ${\bf S}$ identity matrix, some artifacts might be appeared if the operator ${\bf \tilde W}$differs from W. Thus, both the accurate approximation for the wave equation and the wide aperture and the dense sampling are important to get a good image.

The other approach to getting a good image for a given sampling interval and a given poor but cheap operator might be to use the least-squares optimization technique to find the best image. To solve formally for the unknown ${\bf m}$in the least-squares sense, we solve the normal equation as follows:

$${\bf \hat m } = {\bf (\tilde W'S'S\tilde W)}^{-1}{\bf \tilde W'S'd}$$ (4)

This may not recover the original image ${\bf m}$because ${\bf d}$ was generated by the operator ${\bf SW}$which differs from the operator ${\bf S\tilde W}$.If we assume that we have approximated the wave equation very accurately, that is ${\bf \tilde W = W}$,and the only problem that causes artifacts is due to the operator ${\bf S}$,we are solving correct normal equation as

$${\bf \hat m } = {\bf (W'S'SW)}^{-1}{\bf W'S'd}.$$ (5)

Usually the operator ${\bf S\tilde W}$ is too large in size to solve the inversion of ${\bf (\tilde W'S'S\tilde W)}$.An alternative might be an iterative inversion, with the hope of a fast convergence, such as the conjugate gradient method. Fortunately, the initial guess for the conjugate gradient is conventional migration, which is the image after applying the conjugate operator of the forward operator. It gets very close to the solution except for some artifacts. We can expect a fast convergence and a close solution obtained in a few iterations. The object function we want to minimize using the conjugate gradient can be formulated as follows :

$$\min_{\bf m} \vert\vert {\bf \hat d - S\tilde W m} \vert\vert$$ (6)

where ${\bf \hat d}$ represents the data generated by the forward operator ${\bf S\tilde W}$ and this operator could be any modeling operator like the Kirchhoff, phase-shift, or finite-difference method.