Preconditioned least-squares reverse-time migration using random phase encoding

Hessian Estimation

For any linear operator $\mathbf F$ , the Hessian matrix is computed as follows:

$$H(\mathbf x, \mathbf y) = \mathbf F^*(\mathbf y) \mathbf F(\mathbf x),$$ (13)

where $\mathbf x$ and $\mathbf y$ are model coordinates. There are several ways to utilize the Hessian matrix in the inversion process, but I will focus on using its diagonal as a preconditioner to the gradient:

$$\mathbf s_k = \frac{\mathbf g_k}{{\rm diag} \lbrace \mathbf H \rbrace+\epsilon \mathbf I},$$ (14)

where $\mathbf g_k$ is gradient at the $k^{\rm th}$ iteration, $\mathbf s_k$ is the preconditioned search direction, $\mathbf I$ is an identity operator, and $\epsilon$ is a constant used to stabilize the division. For the modeling operator, the diagonal of the Hessian matrix can be written as follows:

$${\rm diag} \lbrace \mathbf H \rbrace = D(\mathbf x) = \sum_{\omeg... ...\rvert^2 \sum_{\mathbf x_r} \lvert G_0(\mathbf x ,\mathbf x_r,\omega) \rvert^2.$$ (15)

Unlike with forward and adjoint modeling operators, the computations must be done on each source-receiver pair separately. As with LSRTM, the expense of this operation can be reduced by encoding the source or receiver side, or both sides. I first define a receiver-side encoding function $\beta$ as follows:

$$\beta(\mathbf x_r, p_r) = \frac{1}{\sqrt{N_{\rm realize}}} \gamma(\mathbf x_r, p_r),$$ (16)

where $p_r$ is the realization index, and the other variables are the same as in the encoding function $\alpha$ . I now define an encoded receiver wavefield:

$$R(\mathbf x, p_r,\omega)= \sum_{\mathbf x_r} \beta(\mathbf x_r, p_r) G_0(\mathbf x,\mathbf x_r,\omega).$$ (17)

This encoding results in the receiver-side blended function $\widetilde{\mathbf D}$ , which can be written as follows:

$$\widetilde{D}(\mathbf x, p_r) = \sum_{\omega} \omega^4 \lvert f(\... ...s} \lvert G_0(\mathbf x,\mathbf x_s,\omega) R(\mathbf x , p_r,\omega) \rvert^2.$$ (18)

The cost of one realization of the function $\widetilde{\mathbf D}$ is equivalent to an unblended migration of all the shots. Additionally, the source side can also be blended:

$$\widetilde{\widetilde{D}}(\mathbf x, p_s, p_r) = \sum_{\omega} \omega^4 \lvert S(\mathbf x, p_s,\omega) R(\mathbf x , p_r,\omega) \rvert^2.$$ (19)

The cost of one realization of the function $\widetilde{\widetilde{\mathbf D}}$ is equivalent to migrating one shot only. However, the additional blending results in more crosstalk. Hence, the function $\widetilde{\widetilde{\mathbf D}}$ requires more realizations to reduce the crosstalk artifacts than does the function $\widetilde{\mathbf D}$ .

Although the cost of computing the function $\mathbf D$ can be reduced with blending, additional propagation(s) are still required to compute the receiver side. Therefore, preconditioning with the source intensity function $\mathbf D_S$ can be done by ignoring the receiver side of the Hessian matrix. The source intensity function can be written as follows:

$$D_S(\mathbf x) = \sum_{\omega} \omega^4 \lvert f(\omega) \rvert^2 \sum_{\mathbf x_s} \lvert G_0(\mathbf x,\mathbf x_s,\omega) \rvert^2.$$ (20)

The previous formulation computes the exact source function in one iteration of LSRTM. However, if the inversion is done with the blended operator, the source intensity function can be computed using the blended source wavefield as follows:

$$\widetilde{D}_S(\mathbf x, p_s) = \sum_{\omega} \omega^4 \lvert S(\mathbf x,\mathbf p_s,\omega) \rvert^2.$$ (21)

By comparing equation (21) to equation (19), we can see that ignoring the receiver side in the Hessian matrix can be physically interpreted as having receivers everywhere in the subsurface. As a result, the source intensity function overestimates the Hessian matrix. Therefore, I propose a better approximation to the Hessian matrix by using the blended source wavefiled to approximate the receiver-side. This source-based Hessian can be written as follows:

$$\widetilde{D}_{SS}(\mathbf x, p_s) = \sum_{\omega} \omega^4 \lvert S(\mathbf x,\mathbf p_s,\omega) \rvert^4.$$ (22)

The function $\widetilde{\mathbf D}_{SS}$ approximates the receiver wavefield by the source wavefield. Physically, this assumes that the receivers exist at the same location as the sources. This is a better approximation than the original source intensity function, especially for the fixed-spread geometry. This formulation requires no additional propagation if the source side is blended. However, there are two sources of error in equation (22). First, the receiver spacing could be different than the source spacing, even if their spreads cover the same area. Second, the receiver side should have a different encoding function than the source side. These errors will prevent the source-based Hessian from approaching the true Hessian matrix, regardless of the number of realizations.

Preconditioned least-squares reverse-time migration using random phase encoding