Objective function

Seismic imaging can be formulated as an optimization problem (Nemeth et al., 1999; Ronen and Liner, 2000):

$$J(\mathbf{m}) = \big\vert\big\vert\mathbf{d} - \mathbf{L}\mathbf{m}\big\vert\big\vert _{2}^{2},$$ (5)

where $\mathbf{L}$ corresponds to the modeling operator, $\mathbf {d}$ corresponds to seismic data, $\mathbf{m}$ is an image of the subsurface containing both reflections and diffractions. In this work we use Kirchhoff modeling and migration as forward ( $\mathbf{L}$) and adjoint ( $\mathbf{L}^T$) operators correspondingly. To image reflections and diffractions as separate components of a wavefield we propose the following objective function:

$$J(\mathbf{m}_r,\mathbf{m}_d) = \Bigg\vert\Bigg\vert \mathbf{PD}_{... ...matrix}\mathbf{m}_r\\ \mathbf{m}_d \end{bmatrix} \Bigg\vert\Bigg\vert _{2}^{2}.$$ (6)

Here, $\mathbf{L}$ is forward Kirchhoff modeling operator, $\mathbf{m}_d$ and $\mathbf{m}_r$ are diffraction and reflection images (models) respectively, $\mathbf {P}$ stands for analytical path-summation migration operator computed from equation 2 ( $\mathbf {P}$ corresponds to $\hat{I}(\Omega,k,v_a,v_b,\beta)$ including forward Fourier, inverse Fourier and $\sigma=t^2$ transforms), $\mathbf {D}_{data}$ is a plane-wave destruction filter applied in the data domain. To penalize the objective function and perform reflection/diffraction separation the shaping regularization operator $\mathbf{R}_{\mathbf{S}_{\sigma},\mathbf{T}_{\lambda}}(\mathbf{m}_r,\mathbf{m}_d)$ is designed (see the “Regularization” subsection for details).

Rather than fitting seismic data $\mathbf {d}$ as in equation 5, we search for models minimizing a misfit weighted by path-summation migration operator. Path-summation operator $\mathbf {P}$ preceded by plane-wave destruction filter applied in the data domain ( $\mathbf {D}_{data}$) for reflection energy elimination can be treated as diffraction probability at a certain location. Therefore, model fitting becomes constrained to those samples, in which diffractions are most probable. The terms of the misfit in equation 6 can be rearranged to stress the misfit weighting by diffraction probability ( $\mathbf{P}\mathbf{D}_{data}$):

$$J(\mathbf{m}_r,\mathbf{m}_d) = \Bigg\vert\Bigg\vert \mathbf{PD}_{... ...athbf{R}_{\mathbf{S}_{\sigma},\mathbf{T}_{\lambda}}(\mathbf{m}_r,\mathbf{m}_d).$$ (7)

Incorporation of PWD ( $\mathbf {D}_{data}$) as a misfit weight into the inversion (equation 6) puts some of the reflected energy in the null space. PWD is a crucial part of the inversion. By weighting the residual with path-summation integral and PWD we guide the inversion to extract diffractions from the full wavefield, which is predominated by reflections. To restore reflections from the null space and account for their leakage to the diffraction image domain PWD can be disabled and another inversion cascade can be run:

$$J(\mathbf{m}_r,\mathbf{m}_d) = \Bigg\vert\Bigg\vert \mathbf{P} \m... ...ma},\mathbf{T}_{\lambda},\mathbf{O}_{\mathbf{m}_r}}(\mathbf{m}_r,\mathbf{m}_d),$$ (8)

where $\mathbf{R}_{\mathbf{S}_{\sigma},\mathbf{T}_{\lambda},\mathbf{O}_{\mathbf{m}_r}}(\mathbf{m}_r,\mathbf{m}_d)$ is the shaping regularization operator designed to perform reflection/diffraction separation after PWD operator ( $\mathbf {D}_{data}$) deactivation (see “Regularization” subsection for details).

Path-summation integral acts as a dip filter suppressing events with high wavenumbers (Merzlikin and Fomel, 2017a). High wavenumbers can be restored by running the third cascade of the inversion with path-summation integral filter disabled:

$$J(\mathbf{m}_r^w,\mathbf{m}_d^w) = \Bigg\vert\Bigg\vert \mathbf{d... ...\mathbf{T}_{\lambda},\mathbf{O}_{\mathbf{m}_r}}(\mathbf{m}_r^w,\mathbf{m}_d^w).$$ (9)

A practical alternative is to restore high wavenumbers separately for reflections and diffractions. Once reflection $\mathbf{m}_r$ and diffraction $\mathbf{m}_d$ components are restored (after the second inversion cascade 8), they can be used as starting models in the following inversions for wavenumber restoration:

$$J(\mathbf{m}_r^{w}) = \Vert \mathbf{d} - \mathbf{L} \mathbf{m}_r^{w} \Vert _{2}^{2}\ +\ \mathbf{R}_{\mathbf{S}_{\sigma}}(\mathbf{m}_r^{w}),$$ (10)

where $\mathbf{m}_r^{w}$ are reflections with high wavenumbers restored, which are initialized with $\mathbf{m}_r^{w}=\mathbf{m}_r$ at the iteration zero. High wavenumbers in diffractions $\mathbf{m}_d^{w}$ can be restored as

$$J(\mathbf{m}_d^{w}) = \Vert \mathbf{d} - \mathbf{m}_r^{w} - \math... ...}_d^{w} \Vert _{2}^{2}\ +\ \mathbf{R}_{\mathbf{T}_{\lambda}}(\mathbf{m}_d^{w}).$$ (11)

Reflections restored by equation 10 ( $\mathbf{m}_r^{w}$) are subtracted from the full wavefield data $\mathbf {d}$ to allow for diffraction-only restoration $\mathbf{m}_d^w$. No interference is expected between reflections $\mathbf{m}_r^{w}$ and diffractions $\mathbf{m}_d^w$ provided accurate starting models from the previous inversions (equations 6 and 8).

Instead of implementing $\mathbf{R}_{\mathbf{S}_{\sigma},\mathbf{T}_{\lambda}}(\mathbf{m}_r,\mathbf{m}_d)$, $\mathbf{R}_{\mathbf{S}_{\sigma},\mathbf{T}_{\lambda},\mathbf{O}_{\mathbf{m}_r}}(\mathbf{m}_r,\mathbf{m}_d)$, $\mathbf{R}_{\mathbf{S}_{\sigma}}(\mathbf{m}_r^{w})$ and $\mathbf{R}_{\mathbf{T}_{\lambda}}(\mathbf{m}_d^{w})$ as penalty terms for reflections and diffractions correspondingly, the problem is reformulated in the shaping framework (see “Regularization” subsection for details).