[2] The matrix expression of linearized migration/inversion

The linearized migration/inversion can be formulated from the integral expressed in equation (15). It can be regarded as an inverse generalized Radon transform. Equation (15) can be expressed as a matrix equation. The process for solving the equation set is just the migration/inversion imaging. Following Berkhout (1997), we first give a matrix expression of wave propagation from a source to a scatterer and then to a receiver:  

$$W^{U}R\left(\theta \right)W^{D} =\left[ \begin{array} {llcl}... ..._{Q1} & g^{D}_{Q2} & \cdots & g^{D}_{QN} \\ \end{array}\right],$$ (58)

where W^U is a discretized Green's function for upward wave propagation, W^D is a discretized Green's function for downward wave propagation, and $R\left(\theta \right)$ is a reflectivity matrix, which is related to the incident angle. If the variation of reflectivity with angle is neglected, $R\left(\theta \right)$ becomes a diagonal matrix. The reflectivity in this case is assumed to be the normal reflectivity. In practice, the reflectivity of a reflector varies with the incident angle. This is called an AVO/AVA phenomenon in seismology. The prestack migration/inversion aims at estimating the angle reflectivity to evaluate lithological variations. On the other hand, the residual moveout of the angle reflectivity indicates whether the macro migration/inversion velocity is reasonable or not. The synthetic wave field can be written as follows:  

$$P\left(\vec{x}_{r},\vec{x}_{s},\omega \right) =\sum\limits_{... ... \left( g^{U}rg^{D}\right) _{MN} \\ \end{array}\right]_{i_{z}}.$$ (59)

In equation (59), each element of the matrix P is a recorded seismic trace in the time domain and a recorded amplitude value for a shot-receiver pair in the frequency domain. Each column is a shot gather, and each row is a common receiver gather. Therefore, equation (59) can be regarded as the matrix expression of equation (15). The classical prestack migration can be formulated as the following:  

$$\left[ W^{U}\left( z_{0},z_{1}\right) \right] ^{H}P\left(\ve... ...}\left( z_{1},z_{0}\right) \right]^{H} = R\left( z_{1} \right).$$ (60)

The detailed matrix expression of equation (60) is

$$\begin{eqnarray} &&\left[ \begin{array} {llcl} \tilde{g}^{U}_{11} & \tilde{g}^{U... ...otfill \\ r_{P1} & r_{P2} & \cdots & r_{PQ} \\ \end{array}\right],\end{eqnarray}$$ (61)

where $R\left( z_{1} \right)$ is the image of the first layer. In $R\left( z_{1} \right)$, each row is an angle gather at an imaging point, and each column is a common angle gather. The multiplication of the pth row in the matrix $\left[ W^{U}\left( z_{0},z_{1}\right) \right] ^{H}$ by any column in the matrix $P\left(\vec{x}_{r},\vec{x}_{s},\omega \right)$ corresponds to a detection focusing of a shot gather; the multiplication of the qth column in the matrix $\left[ W^{D}\left( z_{1},z_{0}\right) \right]^{H}$ by any row in the matrix $P\left(\vec{x}_{r},\vec{x}_{s},\omega \right)$ corresponds to an emission focusing. Then, the image of the second layer can be obtained with  

$$\left[ W^{U}\left(z_{0},z_{2}\right) \right] ^{H}P\left(\vec... ...D}\left(z_{2},z_{0}\right) \right]^{H} = R\left( z_{2} \right).$$ (62)

Generally, the image of the z_ith layer is  

$$\left[ W^{U}\left( z_{0},z_{i}\right) \right]^{H} P\left(\ve... ...}\left( z_{i},z_{0}\right) \right]^{H} = R\left( z_{i} \right).$$ (63)

Here, the matrices $\left[ W^{U} \right]^{H}$ and $\left[ W^{D} \right]^{H}$ are non-recursive. Otherwise, equation (62) and (63) will be of the following forms:  

$$\left[ W^{U}\left( z_{1},z_{2}\right) \right] ^{H}\left[ W^{... ...}\left( z_{2},z_{1}\right) \right]^{H} = R\left( z_{2} \right),$$ (64)

or  

$$\left[ W^{U}\left( z_{1},z_{2}\right) \right] ^{H}P\left(\ve... ...}\left( z_{2},z_{1}\right) \right]^{H} = R\left( z_{2} \right),$$ (65)

and  

$$\left[ W^{U}\left( z_{i-1},z_{i}\right) \right] ^{H} \cdots ... ...left( z_{i},z_{i-1}\right) \right]^{H} = R\left( z_{i} \right),$$ (66)

or  

$$\left[ W^{U}\left( z_{i-1},z_{i}\right) \right] ^{H} P\left(... ...left( z_{i},z_{i-1}\right) \right]^{H} = R\left( z_{i} \right).$$ (67)

Defining the cost function as  

$$E\left(R\left( \theta \right) \right) =\Vert W^{U}R\left(\th... ...^{D}-P\left(\vec{x}_{r},\vec{x}_{s},\omega \right)\Vert^{2}_{2}$$ (68)

yields the formula of the linearized migration/inversion:  

$$R\left( \theta \right) =\frac{\left[ W^{U}\right] ^{H}P\left... ...ft[ W^{U}\right] \left[ W^{D}\right] \left[ W^{D}\right] ^{H}}.$$ (69)

The matrix expression of the migration/inversion in equation (69) is

$$\begin{eqnarray} &&\frac{\left[ \begin{array} {llcl} \tilde{g}^{U}_{11} & \tilde... ...otfill \\ r_{P1} & r_{P2} & \cdots & r_{PQ} \\ \end{array}\right],\end{eqnarray}$$ (70)

where the denominator term $\left[ M\right]$ is

$$\begin{eqnarray} \left[ M\right] &=&\left[\begin{array} {llcl} \tilde{g}^{D}_{11... ...D}_{Q1} & g^{D}_{Q2} & \cdots & g^{D}_{QN} \\ \end{array}\right] .\end{eqnarray}$$ (71)

From equation (61) and (70), the migration/inversion can be locally implemented, because all elements in the matrices W^U, $\left[ W^{U} \right]^{H}$, W^D and $\left[ W^{D} \right]^{H}$ relate only to a given layer. If the matrix $\left[ W^{U} \right]^{H}$ is the inverse of the matrix W^U, it can be expressed as  

$$\left[ W^{U}\right] ^{H}W^{U}=E ,$$ (72)

where the matrix E is an identity matrix. Similarly, if the matrix $\left[ W^{D} \right]^{H}$ is the inverse of the matrix W^D, we have  

$$W^{D}\left[ W^{D}\right] ^{H}=E.$$ (73)

In practice, $\left[ W^{U} \right]^{H}$ and $\left[ W^{D} \right]^{H}$ are the conjugates of W^Uand W^D respectively. Therefore, the matrix $\left[ M\right]$ is a band-width-limited diagonal matrix. The velocity structure and the acquisition geometry affect the inner structure of the matrix. In fact, M is a Hessian which will be discussed in detail later.