Plane-wave Sobel attribute for discontinuity enhancement in seismic images

Theory

The traditional Sobel operator approximates a smoothed gradient of the image intensity function. It is defined as the convolution of an image with two 3$\times$3 filters. The first of these filters ($S_i$) differentiates in the inline direction and averages in the crossline direction. The second filter ($S_x$) differentiates in the crossline direction and averages in the inline direction.

$$S_i=\left[ \begin{array}{rrr} -1 & 0 & 1 \\ -2 & 0 & 2 \\ ... ...\right]\left[ \begin{array}{ccc} -1 & 0 & 1 \end{array}\right]$$ (1)

$$S_x=S_i^{\mathrm{T}}=\left[ \begin{array}{rrr} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{array}\right]$$ (2)

In the $Z$-transform notation, filters (1) and (2) can be expressed as

$$\begin{array}{l} S_i(Z_i,Z_x)=(Z_x+2+Z_x^{-1})(Z_i-Z_i^{-1}) \\ S_x(Z_i,Z_x)=(Z_x-Z_x^{-1})(Z_i+2+Z_i^{-1}) \end{array}$$ (3)

where $Z_j$ is a phase shift in the $j$ direction.

The inline and crossline images are combined to approximate the magnitude of the image gradient (Chopra and Marfurt, 2007) where $\mathbf{d}$ is the data and $\mathbf{S}_i$ and $\mathbf{S}_x$ are convolution operators with the filters $S_i$ and $S_x$, respectively.

$$\Vert\mathbf{\nabla}\mathbf{d}\Vert\approx \sqrt{(\mathbf{S}_i\mathbf{d})^2+(\mathbf{S}_x\mathbf{d})^2}$$ (4)

We propose to modify the filter for application to 3D seismic images by orienting the filter along the structure of seismic reflectors. Local slopes of seismic reflections are estimated in the inline and crossline directions using accelerated plane-wave destruction filters (Chen et al., 2013). The local-plane wave model assumes seismic traces can be effectively predicted by dynamically shifting adjacent seismic traces. This dynamic shift corresponds to the dip of seismic reflections. This model is useful for seismic data characterization and is the basis for plane-wave destruction filters. The local plane-wave differential equation is defined by Claerbout (1992) as

$$\frac{\partial u}{\partial x}+p\frac{\partial u}{\partial t} = 0 \ ,$$ (5)

where $u$ is the seismic wavefield and $p$ is the temporally and spatially variable local slope. The optimal local slopes in the inline and crossline directions can be determined by minimizing the regularized plane-wave residual (Fomel, 2002; Chen et al., 2013).

We further modify the Sobel filter by replacing the derivative operation with plane-wave destruction (Fomel, 2002) and the smoothing operation with plane-wave shaping (Fomel, 2007; Swindeman, 2015). High order plane-wave destruction filters are described in the $Z$-transform notation as

$$\begin{array}{l} C_i(p_i)=B(p_i,Z_t^{-1})-Z_iB(p_i,Z_t)\\ C_x(p_x)=B(p_x,Z_t^{-1})-Z_xB(p_x,Z_t) \end{array}$$ (6)

where $C$ is the plane-wave destruction filter, $B$ is an all-pass filter, and $p_i$ and $p_x$ are the local slopes in the inline and crossline directions, respectively.

Inline and crossline shaping filters are applied to the crossline and inline plane-wave destruction images, respectively. Thus, the plane-wave Sobel filter modifies equation (3) by effectively replacing $Z_i$ with

$$Z_i\frac{B(p_i,Z_t)}{B(p_i,Z_t^{-1})}$$ (7)

and $Z_x$ with

$$Z_x\frac{B(p_x,Z_t)}{B(p_x,Z_t^{-1})} \ .$$ (8)

This modification effectively orients the plane-wave Sobel filter along seismic reflection structures.

In the conventional implementation, the inline and crossline images are combined to produce the final image (equation 4). We propose an alternative approach based on the efficient azimuth scanning workflow (Merzlikin et al., 2017b). We scan through a window of azimuths and compute the $L_q$ norm of linear combinations of the inline and crossline images weighted by the sine and cosine of the azimuth $\alpha$

$$\hat{\mathbf{S}}(\alpha,p_i,p_x)=\Vert\mathbf{S}_i(p_i,p_x)\... ...d}\cos\alpha+\mathbf{S}_x(p_i,p_x)\mathbf{d}\sin\alpha\Vert _q$$ (9)

where $\mathbf{S}_i(p_i,p_x)$ and $\mathbf{S}_x(p_i,p_x)$ correspond to convolution operators with the filters $S_i(p_i,p_x)$ and $S_x(p_i,p_x)$, respectively. The azimuth which corresponds with the discontinuous features produces the best image at each point is picked on a semblance-like panel using a regularized automatic picking algorithm (Fomel, 2009). In images with cross-cutting relationships between discontinuous geologic features, the azimuth which corresponds to the more prominent discontinuity will be preferentially selected. The ensemble of images is then sliced using the pick to generate the optimal image. This improves the resolution of discontinuous features by effectively orienting the plane-wave Sobel filter perpendicular to edges in the seismic image.

As suggested by equation (9), this workflow can be applied the seismic image $\mathbf{d}$; however, we find that the azimuth of the seismic expression of discontinuous geologic features can be estimated more accurately by applying the workflow outlined above to a coherence, semblance, or Sobel filtered image. In other words, we suggest to replace $\mathbf{d}$ in equation (9) with the plane-wave Sobel filter image ( $(\mathbf{S}_i(p_i,p_x)\mathbf{d})^2+(\mathbf{S}_x(p_i,p_x)\mathbf{d})^2$). By cascading multiple iterations of the plane-wave Sobel filter, we generate a more segmented image with less noise and improved estimates of the azimuth of the discontinuous features as expressed in the seismic image. Instead of replacing $\mathbf{d}$ with the plane-wave Sobel filter image, we invite users to cascade their favorite coherence or semblance attributes with the plane-wave Sobel filter to further enhance discontinuous features in seismic images.

Plane-wave Sobel attribute for discontinuity enhancement in seismic images