Noise-robust detection and tracking of salt domes in postmigrated volumes using texture, tensors, and subspace learning

TL;DR

This paper presents noise-robust algorithms for detecting and tracking salt dome boundaries in seismic data, combining texture analysis, tensor-based subspace learning, and noise-adjusted PCA to improve accuracy and efficiency in noisy environments.

Contribution

It introduces a novel combination of texture-based detection, tensor-based tracking, and noise-adjusted PCA for robust salt dome boundary identification in seismic volumes.

Findings

Algorithms are effective in noisy seismic data

Boundary tracking reduces parameter tuning

Validated on offshore seismic datasets with encouraging results

Abstract

The identification of salt dome boundaries in migrated seismic data volumes is important for locating petroleum reservoirs. The presence of noise in the data makes computer-aided salt dome interpretation even more challenging. In this paper, we develop noise-robust algorithms that can label boundaries of salt domes both effectively and efficiently. Our research is twofold. First, we utilize a texture-based gradient to accomplish salt dome detection. We show that by employing a dissimilarity measure based on two-dimensional (2D) discrete Fourier transform (DFT), the algorithm is capable of efficiently detecting salt dome boundaries with accuracy. At the same time, our analysis shows that the proposed algorithm is robust to noise. Once the detection is performed for an initial 2D seismic section, we propose to track the initial boundaries through the data volume to accomplish an efficient labeling process by avoiding parameters tuning that would have been necessary if detection had been performed for every seismic section. The tracking process involves a tensor-based subspace learning process, in which we build texture tensors using patches from different seismic sections. To accommodate noise components with various levels in a texture tensor, we employ noise-adjusted principal component analysis (NA-PCA), so that principal components corresponding to greater signal-to-noise ratio values may be selected for tracking. We validate our detection and tracking algorithms through experiments using seismic datasets acquired from Netherland offshore F3 block in the North Sea with very encouraging results.

Figures

Tables

Table 1. Table 1: The statistical measures of SalSIM indices in Figure  LABEL:fig:detection_cmpr obtained from various detection methods.

Detection Methods	Mean	Standard	AMD (pixels)
		Deviation
Detection method	0.9348	0.0104	11.64
based on GoT Maps
Detection method based	0.9290	0.0104	13.92
on GLCM Contrast Maps Berthelot et al., (2013)
Detection method based on	0.9054	0.0159	20.73
Gradient Maps Aqrawi et al., (2011)

Table 2. Table 2: The statistical measures of the SalSIM indices of boundaries delineated by various detection methods in the crossline dataset.

Detection Methods	Mean	Standard	AMD (pixels)
		Deviation
Detection method	0.9678	0.00211	4.16
based on GoT Maps
Detection method based	0.9596	0.00326	6.57
on GLCM Contrast Maps Berthelot et al., (2013)
Detection method based on	0.9410	0.00822	8.74
Gradient Maps Aqrawi et al., (2011)

Table 3. Table 3: The statistical measures of SalSIM indices in Figure  LABEL:fig:tracking_cmpr obtained from various tracking strategies.

Tracking Methods	Mean	Standard	AMD
		Deviation	(pixels)
Proposed tracking method	0.9531	0.0111	8.32
Tracking method based	0.9496	0.0148	9.30
on vectorization
Tensor-based tracking method	0.9364	0.0295	11.26
without GoT maps
Tensor-based tracking method	0.9422	0.0189	10.53
with GLCM contrast maps

Equations

$${\mathbf{G}}_{\mathbf{x}}[i,j]=d\left(W_{x-}^{i,j},W_{x+}^{i,j}\right),$$

$${\mathbf{G}}_{\mathbf{y}}[i,j]=d\left(W_{y-}^{i,j},W_{y+}^{i,j}\right),$$

$$\mathbf{G}[i,j]={\left({{\mathbf{G}}_{\mathbf{x}}}^2[i,j]+{{\mathbf{G}}_{\mathbf{y}}}^2[i,j]\right)}^{\frac{1}{2}},$$

$${\mathbf{G}}_{\mathbf{x}}[i,j]=\sum _{n=1}^N{w}_n.d\left(W_{n,x-}^{i,j},W_{n,x+}^{i,j}\right),$$

$${\mathbf{G}}_{\mathbf{y}}[i,j]=\sum _{n=1}^N{w}_n.d\left(W_{n,y-}^{i,j},W_{n,y+}^{i,j}\right),$$

$$\mathbf{G}[i,j]={\left({{\mathbf{G}}_{\mathbf{x}}}^2[i,j]+{{\mathbf{G}}_{\mathbf{y}}}^2[i,j]\right)}^{\frac{1}{2}},$$

$$d(W_{-},W_{+})=E\left\{\left|\mathcal{F}\left\{\left|\mathcal{F}\left\{\left|W_{-}-W_{+}\right|\right\}\right|\right\}\right|\right\},$$

$$\mathbf{I}=\left[\begin{array}{cc}{I}_{xx}& I_{xy}\\ I_{yx}& I_{yy}\end{array}\right],$$

$$I_{xx}=\sum _{[k,l]\in W^{i,j}}({\mathbf{\Delta }}_{\mathbf{y}}[k,l]-{\bar{\Delta }}_y{)}^2,$$

$$I_{xy}=I_{yx}=\sum _{[k,l]\in W^{i,j}}-({\mathbf{\Delta }}_{\mathbf{x}}[k,l]-{\bar{\Delta }}_x)({\mathbf{\Delta }}_{\mathbf{y}}[k,l]-{\bar{\Delta }}_y),$$

$$I_{yy}=\sum _{[k,l]\in W^{i,j}}({\mathbf{\Delta }}_{\mathbf{x}}[k,l]-{\bar{\Delta }}_x{)}^2,$$

$$\mathbf{D}[i,j]=\sum _{n=1}^{N}1-\frac{\min ({\Lambda }_{1,n}[i,j],{\Lambda }_{2,n}[i,j])}{\max ({\Lambda }_{1,n}[i,j],{\Lambda }_{2,n}[i,j])}.$$

$$\{\begin{array}{rl}& \arg \underset{[i,j]}{\min }\mathbf{D}\ast \mathbf{\Omega }\\ & \mathbf{\Omega }[u,v]=e^{-\frac{d_{\mathbf{\Omega }}^2(u,v)}{2{\sigma }_{\mathbf{\Omega }}^2}},u,v\in \left\{1,2,\cdots ,2N+1\right\}\end{array}$$

$$\arg \underset{T}{\min }\left\{{\sigma }_1^2(T)\sum _{i=0}^{T-1}p(i)+{\sigma }_2^2(T)\sum _{i=T}^{N-1}p(i)\right\},$$

$$\arg \underset{T}{\max }\left\{\left(\sum _{i=0}^{T-1}p(i)\right)\left(\sum _{i=0}^{T-1}p(i)\right)\left({\mu }_1(T)-{\mu }_2(T)\right)\right\},$$

$${\hat{W}}_1=W_1+n_1,$$

$${\hat{W}}_2=W_2+n_2,$$

$$d({\hat{W}}_1,{\hat{W}}_2)=E\left\{\left|\mathcal{F}\left\{\left|\mathcal{F}\left\{{\Delta }_W+{\Delta }_n\right\}\right|\right\}\right|\right\}.$$

$$\mathcal{F}\left\{{\Delta }_W+{\Delta }_n\right\}=\mathcal{F}\left\{{\Delta }_W\right\}+\mathcal{F}\left\{{\Delta }_n\right\}$$

$$\mathcal{F}\left\{{\Delta }_n\right\}=\sum _{k=0}^{2n}\sum _{l=0}^{2n}a[k,l]+j\sum _{k=0}^{2n}\sum _{l=0}^{2n}b[k,l],$$

$$n_m=|\mathcal{F}\left\{{\Delta }_n\right\}|\cos \varphi ,$$

$$\scalebox{0.96}{$\mathcal{B}(i_{1}\cdots i_{n-1}j_{n}i_{n+1}\cdots i_{N})=\sum\limits_{i_{n}}\mathcal{A}(i_{1}\cdots i_{N})\cdot\mathbf{U}(j_{n},i_{n})$}.$$

$$\mathcal{C}=\mathcal{A}{\times }_1{\tilde{\mathbf{U}}}^{(1{)}^T}{\times }_2{\tilde{\mathbf{U}}}^{(2{)}^T}\cdots {\times }_N{\tilde{\mathbf{U}}}^{(N{)}^T},$$

$$\scalebox{0.85}{$\begin{aligned} e_{c}=&\sum\limits_{n=\{1,2\},\mathcal{M}=\{\mathcal{S},\mathcal{G}\}}\left\|\mathcal{M}_{p_{i}}-\mathcal{M}_{p_{i}}\times_{n}\left(\widetilde{\mathbf{U}}^{(n)}_{\widetilde{\mathcal{M}}_{k}}\cdot\widetilde{\mathbf{U}}^{(n)^{T}}_{\widetilde{\mathcal{M}}_{k}}\right)\right\|_{F}^{2}+\\ &\sum\limits_{\mathcal{M}=\{\mathcal{S},\mathcal{G}\}}\left\|\mathbf{M}_{p_{i}}^{(3)}-\mathbf{M}_{p_{i}}^{(3)}\cdot\widetilde{\mathbf{U}}_{\widetilde{\mathcal{M}}_{k}}^{(3)}\cdot\widetilde{\mathbf{U}}_{\widetilde{\mathcal{M}}_{k}}^{(3)^{T}}\right\|_{F}^{2}\end{aligned}$},$$

$$\begin{array}{rl}{e}_t=& \sum _{n=\{1,2\},\mathcal{M}=\{\mathcal{S},\mathcal{G}\}}{\lambda }_{\mathcal{M}}\cdot {\Vert {\mathcal{M}}_{p,j}-{\mathcal{M}}_{p,j}{\times }_n\left({\tilde{\mathbf{U}}}_{\tilde{\mathcal{M}}}^{(n)}\cdot {\tilde{\mathbf{U}}}_{\tilde{\mathcal{M}}}^{(n{)}^T}\right)\Vert }_F^2+\\ & \sum _{\mathcal{M}=\{\mathcal{S},\mathcal{G}\}}{\lambda }_{\mathcal{M}}\cdot {\Vert {\mathbf{M}}_{p,j}^{(3)}-{\mathbf{M}}_{p,j}^{(3)}\cdot {\tilde{\mathbf{U}}}_{\tilde{\mathcal{M}}}^{(3)}\cdot {\tilde{\mathbf{U}}}_{\tilde{\mathcal{M}}}^{(3{)}^T}\Vert }_F^2\end{array},$$

$$\scalebox{0.79}{$\mathbf{\widehat{B}}[i,j]=\left\{\begin{aligned} &0\qquad\qquad\qquad\quad\ \mathbf{B}[i,j]-\underset{R_{m}\times R_{m}}{\mbox{median}}(\mathbf{B}[i,j])=-1\\ &\underset{R_{m}\times R_{m}}{\mbox{median}}(\mathbf{B}[i,j])\quad\mbox{Otherwise}\end{aligned}\right.$},$$

$$\mathbf{L}\oplus \mathbf{H}=\bigcup _{z\in \mathbf{L}}{\mathbf{H}}_z,$$

$${\mathbf{\Sigma }}_{\hat{\mathbf{N}}}={\mathbf{P}}^T{\mathbf{\Sigma }}_{\mathbf{N}}\mathbf{P}={\mathbf{P}}^T\mathbf{V}{\mathbf{\Delta }}_{\mathbf{N}}{\mathbf{V}}^T\mathbf{P}=\mathbf{I},$$

$$\text{mbox}SalSIM=\underset{\text{Local Item}}{\underbrace{e^{-\alpha \cdot ({\mu }_d+{\sigma }_d)}}}\cdot \underset{\text{Global Item}}{\underbrace{e^{-\beta \cdot d_{max}}}},$$

$$\scalebox{0.97}{$\overline{\mathbf{C}}[i,j]=\sum\limits_{[\Delta i,\Delta j]}\left(\frac{1}{N_{g}}\sum\limits_{u}\sum\limits_{v}(u-v)^{2}\mathbf{P_{[\Delta i,\Delta j]}}[u,v]\right)$},$$

$$F=\sqrt{F_x^2+F_y^2+F_z^2},$$

Full text

Citation

Z. Wang, T. Hegazy, Z. Long and G. AlRegib, “Noise-robust detection and tracking of salt domes in postmigrated volumes using texture, tensors, and subspace learning,” Geophysics, vol. 80, no. 6, pp. WD101-WD116, Sep. 2015.

DOI

https://doi.org/10.1190/geo2015-0116.1

Review

Date of publication: 23 September 2015

Bib

@article{wang2015noise,

title={Noise-robust detection and tracking of salt domes in postmigrated volumes using texture, tensors, and subspace learning},

author={Wang, Z. and Hegazy, T. and Long, Z. and AlRegib, G.},

journal={Geophysics},

volume={80},

number={6},

pages={WD101–WD116},

year={2015},

publisher={Society of Exploration Geophysicists}}

Copyright

©2015 SEG. Personal use of this material is permitted. Permission from SEG must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

Contact

zhiling.long@gatech.edu OR alregib@gatech.edu

https://ghassanalregib.com/

\lefthead

Wang et al. \rightheadSalt Dome Detection and Tracking

1 Introduction

The deposition of salt in marine basins commonly intrudes into surrounding rock strata and forms an important geological structure, salt domes. Because of the impermeability of salt, salt domes may seal porous reservoir rocks and lead to the formation of petroleum reservoirs. To localize the possible positions of reservoir regions, experienced interpreters need to accurately label the boundaries of salt domes in migrated seismic data. With the dramatically increasing amount of acquired seismic data, however, manual interpretation is becoming time consuming and labor intensive.

To speed up interpretation, in recent years, researchers have proposed various interactive methods that assist with the detection of salt dome boundaries using graph theory and image processing techniques. Lomask et al., (2004) defined seismic sections as weighted undirected graphs and applied the normalized cut image segmentation (NCIS) Shi and Malik, (2000) to globally optimize the delineation of salt dome boundaries. As an extension of this method, Lomask et al., (2007) employed local dips in the weight matrix and utilized bound constraints to remove boundary artifacts. Similarly, Halpert et al., (2009) introduced the modified NCIS by combining multiple seismic attributes with adaptive weights. Although these NCIS-based methods can be implemented in parallel, the high computational cost limits their future application on high-resolution or three-dimensional (3D) seismic data. To improve the efficiency of global segmentation, Halpert et al., (2010) proposed to detect the boundaries of salt domes using a pairwise region comparison based on the minimum spanning tree Felzenszwalb and Huttenlocher, (2004), which reduces the algorithm complexity from $\mathcal{O}\left(n^{2}\right)$ to $\mathcal{O}\left(n\log n\right)$.

In addition to the segmentation-based methods, edge detection has also been explored. Zhou et al., (2007) and Aqrawi et al., (2011) applied dip-guided 2D and 3D Sobel filters to detect salt dome boundaries in time slices. In addition, Berthelot et al., (2013) proposed to delineate salt domes by testing a combination of multiple attributes (i.e., texture, frequency, and dip) in a supervised Bayesian classification model. To characterize the homogeneous texture of salt bodies, most recently, Hegazy and AlRegib, 2014b derived a directionality attribute from the moment of inertia tensor defined by gradient components. Although the methods mentioned above may be able to achieve various levels of improvement in the detection performance, normally they involve many parameters to be manually tuned, which unavoidably compromises the overall efficiency. In addition to the effectiveness and efficiency issues, salt dome interpretation may be further complicated by the presence of noise. Typically, noise from different sources is present, which may vary across a seismic volume, thus being space-time dependent. Accommodating such noise effect is a challenging task for computer-aided interpretation.

In our work, we view seismic volumes as big data volumes containing certain information that needs to be extracted manually or automatically. Hence, we import our long experience in digital signal processing and visual data analysis to seismic interpretation. In our view, a large seismic volume contains information about salt dome structures. Such structure is characterized by certain attributes. However, calculating these attributes accurately and efficiently faces many challenges. First, the data are subjected to a long processing pipeline of stacking and migration. Not only do these steps introduce artifacts, but they may magnify noise in the raw data. Second, the raw data are an undersampled depiction of the subsurface structure. This is analogous to detecting detailed activities from low-quality videos captured with a low-resolution camera. Third, noise in the data volume can be a function of the subsurface structure, and it may contain information about the subsurface structure. Instead of denoising, our long-term plan is to take advantage of such noise.

In this paper, we focus on developing noise-robust algorithms that can label boundaries of salt domes both effectively and efficiently. Our research is twofold. First, we define a texture-based attribute, namely gradient of texture (GoT), and combine it with a dissimilarity measure based on 2D discrete Fourier transform (DFT) to implement an algorithm capable of detecting salt dome boundaries with high accuracy and efficiency. We further analyze and show that the proposed algorithm is robust to noise. Once salt dome detection is performed for an initial 2D seismic section, we propose to track the initial boundaries through the data volume to accomplish the labeling process. This way, we can avoid tweaking parameters for boundary detection in every seismic section, thus being more efficient. The tracking involves a tensor-based subspace learning process. It starts by forming texture tensors according to the initially labeled boundaries, to better capture the strong correlation between patches of similar texture patterns. It then projects the initial boundaries onto the neighboring seismic section, and determines the tracked boundaries by evaluating texture similarity in the subspace of the texture tensors. To maintain consistency, we rename initial and neighboring seismic sections as reference and predicted sections, respectively, in the rest of this paper. To apply the tensor-based subspace learning, we also take noise into consideration. Unlike traditional applications, with seismic data it is likely that noise of different levels may be included in the same texture tensor. To handle such situations, we employ noise-adjusted subspace learning techniques such as those presented in Roger, (1994); Chang and Du, (1999).

An overall block diagram for our salt dome boundary labeling system is depicted in Figure LABEL:fig:diagram. It includes both the detection of salt dome boundaries for the reference sections, and the tracking of the boundaries through the predicted sections. This figure will be used throughout the paper. The rest of the paper is organized as follows. First, we discuss our salt dome detection method. Then, we present our salt dome tracking algorithm. Thereafter, we conduct experiments to evaluate the performance of the proposed methods and the noise effect. Based on the results, we draw conclusions at the end of the paper.

2 Theory

2.1 Salt Dome Detection

The top boundary of a salt dome usually has a cap rock, which appears in migrated seismic surveys as areas with high contrast and strong edges (Figure LABEL:fig:SaltDome). However, the lateral boundaries usually lack such high contrast and are more difficult to detect. Figure LABEL:fig:SaltDome shows that the lateral boundaries of the salt dome have no explicit edge existing between the salt region and the non-salt region containing horizons. However, seismic interpretation experts would still be able to visually delineate such boundaries using the change of texture between salt and non-salt regions shown by other attributes. As a matter of fact, even non-experts can still perceptually differentiate between the two regions and possibly draw a boundary line. Therefore, we believe that an automated method should take human visual perception into account. Using traditional edge detection methods (e.g., Sobel filter) to delineate such boundary would fail because the boundary exhibits a change in texture, rather than an explicit edge.

Figure LABEL:fig:diagram shows an overview of the proposed detection method. To capture the change in texture, we propose a new attribute, which we call Gradient of Texture (GoT). Then, a threshold is applied to the GoT map $\mathbf{G}$ to obtain a binary image $\mathbf{B_{t}}$. From an initialization point inside the salt dome, we perform region growing such that the grown region boundaries stop at the threshold value of GoT. We perform binary morphological operations on the grown region to enhance it. The boundary of the salt region is determined based on enhanced grown region. In the following, we explain each of the modules in detail.

In relation to existing methods from the image processing literature, we show in the following subsections that our proposed GoT is faster and more general than the “texture gradient” introduced by Martin et al., (2004). First, Martin et al., (2004) computes features along 13 different directions, while we show that only 2 directions (horizontal and vertical) are sufficient. Second, Martin et al., (2004) computes the gradient as the norm of the difference between 13-dimensional feature vectors while we directly assess the dissimilarity between the neighborhood windows without computing feature vectors. Finally, Martin et al., (2004) compares texture in two disk halves centered around a center point, while we use square neighborhood windows sharing a side centered around the center point. Using square (or even rectangular) windows makes it easier to take advantages of existing efficient algorithms such as Fast Fourier Transform (FFT) to assess dissimilarity. It is noteworthy that comparing texture windows is also used in several other applications. For example, Akl et al., (2014) developed a method for a texture synthesis application where it is required to measure the dissimilarity between synthesized texture patches and the input texture to be imitated. First, features such as coherency and orientation are computed. Then, the comparison performed by Akl et al., (2014) measures the dissimilarity using various distance measures. In contrast, our method compares the texture windows directly using a very efficient dissimilarity measure without having to compute feature vectors, which can be very computationally demanding. Further, our method takes into account perceptual similarity not just objective similarity as we believe this is one of the important factors in seismic interpretation.

2.2 Gradient of Texture

We define the gradient of texture at a given point as the perceptual dissimilarity between two square neighborhood windows that share a side centered around the given point. The shared side is perpendicular to the direction of interest Gonzalez and Woods, (2006); Petrou and Petrou, (2010). Figure LABEL:fig:GoT illustrates an example with two texture regions separated by a purely vertical boundary. For a given row, if we slide the center point (and hence the two neighborhood windows in the figure) along the horizontal direction, we expect the GoT profile to follow the curve shown at the bottom of the figure. For other rows, we expect similar curves because the boundary is purely vertical in this example. Theoretically, the highest dissimilarity, and hence the highest GoT value, is obtained when the center point falls exactly on the texture boundary because the left neighborhood window $W_{x-}$ will have completely different texture content from that of the right neighborhood window $W_{x+}$. When we move the center point away from the texture boundary to the position shown in Figure LABEL:fig:GoT, the content of $W_{x+}$ is purely from the right texture region, while $W_{x-}$ contains both textures. The partially shared content will cause the dissimilarity value, and hence the GoT value, to drop. If we keep moving the center point along with the neighborhood windows away from the boundary, we reach a point where the contents of $W_{x-}$ and $W_{x+}$ are from the same texture region; i.e., the texture contents will be exactly similar. This exact similarity means zero dissimilarity, and hence, zero GoT.

For the texture boundary in Figure LABEL:fig:GoT, computing the gradient of texture in the horizontal direction is sufficient because the texture boundary is purely vertical. For the general case, we compute the gradient of texture as a vector with components in both the horizontal and vertical directions. Figure LABEL:fig:GoT_xycomp shows how the gradient of texture is computed for the general case. The solid neighborhood windows are used to compute the horizontal component, while the dashed neighborhood windows are used to compute the vertical component. Mathematically,

[TABLE]

where $\mathbf{G_{x}}[i,j]$ and $\mathbf{G_{y}}[i,j]$ denote the GoT values at point $[i,j]$ along the horizontal and vertical directions respectively, $\mathbf{G}[i,j]$ represents the GoT value combining the horizontal and vertical components, function $d(\cdot)$ defines a dissimilarity measure, and the four $W^{i,j}$ expressions with four different subscripts refer to the four neighborhood windows shown in Figure LABEL:fig:GoT_xycomp.

With the large scale variations in texture content, the neighborhood window size used to compute the GoT attribute needs to be carefully chosen. To improve robustness, we use different window sizes (5 in our case) and we choose the weighted average of the quantities as the attribute value. The neighborhood window sizes range from $3\times 3$ to $11\times 11$. Figure LABEL:fig:neighb_wins shows neighborhood windows of minimum and maximum sizes around a labeled salt dome boundary. Note that different window sizes enable capturing texture dissimilarity across different scales. The multi-scale version of GoT is mathematically expressed as follows:

[TABLE]

where, in our case, $N=5$. Also, for a given $n$, the neighborhood window size is $(2n+1)\times(2n+1)$, and $W^{i,j}_{n,x-}$, $W^{i,j}_{n,x+}$, $W^{i,j}_{n,y-}$, $W^{i,j}_{n,y+}$ respectively denote the neighborhood windows of size $(2n+1)\times(2n+1)$ that are to the left, right, bottom and top of the point $[i,j]$ where the attribute is being computed. Finally, $w_{n}$ denotes the weight of the GoT value computed at the window whose size is $(2n+1)\times(2n+1)$. The GoT value increases as we increase the neighborhood window size because of the dissimilarity measure we are using (as we show in the next section). Thus, with equal weights, the overall attribute would be more biased to the larger window sizes. As a result, to partially compensate for this bias, we choose $w_{n}$ to be inversely proportional to $n$ as $w_{n}=1/n$.

2.3 Dissimilarity Measure

To measure the perceptual dissimilarity between neighboring windows, we use the following measure

[TABLE]

where $\mathscr{F}\{\cdot\}$ denotes the 2D DFT, $E$ is the expectation operator, and $W_{-}$, $W_{+}$ denote the two neighboring windows whose dissimilarity is to be evaluated. This measure is a close variation of the one introduced in Hegazy and AlRegib, 2014a , which was shown to be consistent with human perception while being computationally efficient. The dissimilarity measure characterizes the amount of variation, or chaos, in the magnitude spectrum of the absolute difference between the two windows. The outer Fourier transform is considered the magnitude spectrum of the magnitude spectrum. Since the proposed dissimilarity measure is highly consistent with human perception Hegazy and AlRegib, 2014a , we believe it partially imitates how the human interpreter delineates the salt boundaries according to the texture variations inside and outside the salt region.

2.4 Initialization and Region Growing

GoT detects several texture boundaries including ones outlining salt and others that are irrelevant to salt, e.g., from strong reflectors that are away from the salt dome. In order for the detection method to focus on the salt boundaries only, we identify an initialization point inside the salt dome and then use region growing until a predefined threshold for GoT magnitude is met. The initialization point could be selected automatically or interactively with human interpreter involved. Even with interactive selection of initialization point, the time spent by human interpreter will be considerably reduced since the interpreter would have to quickly click on any arbitrary point within the salt region, as opposed to carefully traversing a long, tortuous salt dome boundary. Note that the proposed method is not sensitive to the initialization point selection as long as the selected point falls anywhere within the salt region. In other words, the initialization point needs to fall within the salt region, but is otherwise arbitrary.

For automated initialization, we propose a selection method that is based on the directionality attribute Hegazy and AlRegib, 2014b . The rationale is that salt regions would have the lowest directionality in the whole seismic section under consideration. To compute the directionality map $\mathbf{D}$, we first compute the moment of inertia tensor $\mathbf{I}$ of the scattered plot of the horizontal and vertical components of the intensity gradient $\mathbf{\Delta}$ computed for a small neighborhood window $W^{i,j}$ centered around the point $[i,j]$. Note that the gradient of intensity $\mathbf{\Delta}$ used here is different from the gradient of texture (GoT) $\mathbf{G}$ introduced in the previous subsections. Mathematically, we compute the moment of inertia tensor $\mathbf{I}$ around a given point [i,j] as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where:

• •

$\mathbf{I}$ is the moment of inertia tensor for the scattered plot of the intensity gradient components, which is obtained for a neighborhood window $W^{i,j}$ centered around point $[i,j]$.

• •

$W^{i,j}$ is a neighborhood window centered around point $[i,j]$. The summation is performed over all the points $[k,l]$ in $W^{i,j}$.

• •

$\mathbf{\Delta_{x}}[k,l]$ and $\mathbf{\Delta_{y}}[k,l]$ respectively denote the horizontal and vertical components of the intensity gradient at a point $[k,l]$ inside $W^{i,j}$ .

• •

$\bar{\Delta}_{x}$ and $\bar{\Delta}_{y}$ respectively denote the mean of the horizontal and vertical components of the gradient computed over the neighborhood window $W^{i,j}$.

The eigenvalues of the moment of inertia tensor correspond to the axes of minimum and maximum moments of inertia. Thus, we can use the eigenvalues of the moment of inertia tensor to capture the directionality of texture as we show below. To make the method more robust, we average the directionality maps obtained using multiple window sizes as outlined in Hegazy and AlRegib, 2014b . Let $\Lambda_{1,n}[i,j]$, $\Lambda_{2,n}[i,j]$ denote the eigenvalues of the moment of inertia tensor computed around point $[i,j]$ using a neighborhood window whose size is $(2n+1)\times(2n+1)$ and centered around [i,j]. We can mathematically define the multi-scale directionality attribute at point $[i,j]$ as:

[TABLE]

For non-salt regions, one eigenvalue will be much larger than the other eigenvalue. Thus, the directionality attribute will be close to 1. On the other hand, for salt regions, the eigenvalues will be almost equal, leading to directionality attributes close to zero.

For added reliability and to protect against selecting a noisy point that happened to have very low directionality, we smooth out the directionality map with a Gaussian filter. We obtain the initialization point by minimizing the smoothed multi-scale directionality map. Therefore, the initialization point is selected according to:

[TABLE]

where $\mathbf{\Omega}$ is a Gaussian kernel with a size of $(2N+1)\times(2N+1)$ used to smooth out the multi-scale directionality attribute map and $\ast$ is the 2D convolution operator. The elements in $\mathbf{\Omega}$ indexed by $u$ and $v$ are determined by function $d_{\mathbf{\Omega}}$, which calculates the distance from $[u,v]$ to the center of the Gaussian kernel. In addition, $\sigma_{\mathbf{\Omega}}$ denotes the standard deviation of the Gaussian kernel.

We grow the initialization point into a region. Since the initialization point is inside the salt dome, i.e., inside a texture region, we expect the texture gradient around the initialization point to be very low (theoretically zero). Thus, the initialization point can be grown since GoT in the neighboring region is less than a certain threshold value (determined according to the next subsection). As we grow the region, we start meeting points that are close to the boundary of the dome where the GoT value is high. We stop growing along those points that reached the threshold of the GoT map. In the following subsection, we discuss how the threshold value is determined.

2.5 Thresholding

The threshold value can be determined interactively with a human interpreter or automatically. For automated threshold value selection, we propose using Otsu’s method Otsu, (1979). The algorithm assumes that there are two classes of pixels following bimodal histogram. The optimal threshold is therefore calculated as the one that minimizes the intra-class variance, which is the weighted sum of the variances of the two classes. Thus, the optimal threshold is given by:

[TABLE]

where $T$ is the threshold value, $p(i)$ is the probability of the intensity value $i$, and $\sigma_{1}^{2}(T)$, $\sigma_{2}^{2}(T)$ are the variances of the first and second classes, respectively. The first and second classes are determined given the threshold $T$. The optimal threshold value is found by exhaustively searching for $T$ between [math] and $N-1$.

For more efficient implementation, Otsu shows that minimizing the intra-class variance is equivalent to maximizing the inter-class variance, which can be simplified as follows:

[TABLE]

where $\mu_{1}(T)$ and $\mu_{2}(T)$ are the mean values of the first and second classes, respectively given a threshold value $T$ Otsu, (1979).

2.5.1 Noise Impact

In this subsection, we discuss the impact of acquisition noise on the salt boundaries detected on the migrated seismic data. Behera and Tsvankin, (2008) study the impact of shot gathers contaminated with Gaussian noise on the migrated seismic data. Under Gaussian noise, the reflectors in the migrated image are correctly positioned but are affected by noise Behera and Tsvankin, (2008). In the following analysis, we assume that the random acquisition noise will result in a migrated image that is contaminated with additive random noise.

Let $\hat{W}_{1}$, $\hat{W}_{2}$ be two noise-contaminated adjacent neighborhood windows whose dissimilarity is to be evaluated to compute the GoT attribute at their center point, such that

[TABLE]

where $W_{1}$, $W_{2}$ denote the noise-free windows, and $n_{1}$, $n_{2}$ denote random noise added to $W_{1}$ and $W_{2}$, respectively. We re-write the dissimilarity measure in Eq.(7) between the noise contaminated windows as follows:

[TABLE]

For simplicity, we remove the absolute value of the input to the inner 2D DFT as we noticed through experimentation that it does not make any significant difference in the performance of the proposed method.

Using the linearity property, we can express the inner 2D DFT as follows

[TABLE]

Thus, the 2D DFT coefficients are the same as those of the noise-free case $\mathscr{F}\left\{\Delta_{W}\right\}$, except they are contaminated with the noise signal $\mathscr{F}\left\{\Delta_{n}\right\}$. According to Schoukens and Renneboog, (1986), we can express this noise signal as

[TABLE]

where each summation is performed over the square neighborhood window size of $(2n+1)\times(2n+1)$, and $a$ and $b$ are random variables. The noise signals affecting the real and imaginary parts of Fourier coefficients are derived from the sum of a number of random variables. Therefore, the original Fourier coefficients are subjected to additive random noise. According to Schoukens and Renneboog, (1986), the covariance matrix of the real and imaginary parts is nearly diagonal, which means that the real and imaginary coefficient noise components are independent.

Both $\mathscr{F}\left\{\Delta_{W}\right\}$ and $\mathscr{F}\left\{\Delta_{n}\right\}$ are complex. As a result, the absolute value of $\mathscr{F}\left\{\Delta_{W}\right\}+\mathscr{F}\left\{\Delta_{n}\right\}$ can be determined by vector algebra as shown in Figure LABEL:fig:sn_vectors, where $\vec{s}$ refers to $\mathscr{F}\left\{\Delta_{W}\right\}$ and $\vec{n}$ refers to $\mathscr{F}\left\{\Delta_{n}\right\}$. If the signal to noise ratio is much larger than 1, then the angle $\theta$ in Figure LABEL:fig:sn_vectors is very small. As a result, noise impact on the absolute value of Fourier coefficients can be approximated by simply projecting the noise component $\vec{n}$ on the original signal $\vec{s}$ vector direction. Consequently, such noise can be approximated as additive random noise $n_{m}$ as follows.

[TABLE]

where $\phi$, which is a random variable, denotes the angle between the signal and noise vectors in Figure LABEL:fig:sn_vectors.

Consequently, $|\mathscr{F}\left\{\Delta_{W}\right\}+\mathscr{F}\left\{\Delta_{n}\right\}|$, which is the input to the outer 2D DFT, suffers additive random noise $n_{m}$. Similarly, applying the same derivation for the magnitude of the outer 2D DFT, we conclude that the effect of noise on the magnitude of the outer 2D DFT coefficients suffer from additive random noise as well, say $n_{d}$. Since the dissimilarity measure $d()$ is the mean value of the the outer 2D DFT magnitude, then $d()$ is shifted by the mean of $n_{d}$, that is $\mu_{d}$.

Theoretically, the additive random noise affecting the migrated images shifts the dissimilarity measure (and hence, GoT) by a constant value $\mu_{d}$, assuming the noise is stationary. Because we apply an adaptive threshold to GoT, we expect almost no effect on the detection results. Practically, however, the window sizes are finite, and therefore, the proposed method will still suffer under noise because the averaging would be performed for a few samples. We conclude that to make the method more robust under noise, it is important to drop the smaller window sizes and use larger window sizes.

2.6 Salt Dome Tracking

In the detection of salt dome boundaries, the accuracy of labeled boundaries depends greatly on the selection of certain parameters. Since salt dome structures commonly vary across seismic sections, interpreters have to adjust parameters so that detected boundaries in each section can accurately capture structural changes. However, extra labor and time involved in the tuning of parameters may lower interpretation efficiency. Therefore, to avoid the frequent adjustment of parameters, we group seismic sections into reference and predicted sections by borrowing the ideas of I- and P-frames from video-coding techniques and apply the proposed detection method only to reference sections that constitute the minority of seismic sections. By tracking detected boundaries in reference sections through a seismic volume, we can synthesize salt dome boundaries in predicted sections with limited human intervention. The lower part of Figure LABEL:fig:diagram illustrates the pipeline of the proposed tracking method, each step of which is going to be introduced in detail in the following subsections.

2.6.1 Tensors and Multi-linear Analysis

By observing salt domes in seismic sections as Figure LABEL:fig:SaltDome shows, we notice that local areas along salt dome boundaries commonly have similar textures. Therefore, by classifying boundary textures in the reference section, we can build texture tensors and employ the tensor-based subspace learning algorithm to acquire their texture features, which act as important constraints to the synthesis of salt dome boundaries in predicted sections. Before we introduce the proposed classification method, we briefly review the basic concepts of tensors and multi-linear analysis.

In multi-linear algebra, a tensor represents a multidimensional array. Conventionally, we denote tensors using calligraphic letters such as $\mathcal{A}$, and matrices using bold uppercase letters such as $\mathbf{U}$. An $N$th-order tensor $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ with entries denoted $\mathcal{A}(i_{1}i_{2}\cdots i_{N})$ has $N$ indices or modes. By unfolding tensor $\mathcal{A}$ along the $n$th mode, we obtain matrix $\mathbf{A}^{(n)}\in\mathbb{R}^{I_{n}\times(I_{1}\times\cdots\times I_{n-1}\times I_{n+1}\cdots\times I_{N})}$. Figure LABEL:fig:unfolding illustrates the unfolding of a third-order tensor along three modes, distinguished by specific colors. The inverse operation of n-mode unfolding is n-mode folding, which restores $\mathcal{A}$ from matrix $\mathbf{A}^{(n)}$. The n-mode product of tensor $\mathcal{A}$ by matrix $\mathbf{U}\in\mathbb{R}^{J_{n}\times I_{n}}$, denoted $\mathcal{A}\times_{n}\mathbf{U}$, defines new tensor $\mathcal{B}$ with entries calculated as

[TABLE]

In addition, the product above can also be implemented by folding $\mathbf{U}\cdot\mathbf{A}^{(n)}$ along the $n$th mode. Similar to the singular value decomposition (SVD) of matrices, the decomposition of tensors can be expressed as $\mathcal{A}=\mathcal{S}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}\cdots\times_{N}\mathbf{U}^{(N)}$, in which $\mathbf{U}^{(n)}$, $n=1,2,\cdots,N$, represents an $I_{n}\times I_{n}$ orthogonal matrix with column vectors spanning the column space of $n$-mode unfolding matrix $\mathbf{A}^{(n)}$ de Lathauwer et al., (2000). Since a tensor with high dimensions requires compact representation in real applications, we can apply multi-linear principal component analysis (MPCA) Lu et al., (2008) to extract a low-dimensional tensor subspace that captures the most variance of the original tensor. Such a tensor subspace, denoted $\mathcal{C}\in\mathbb{R}^{P_{1}\times P_{2}\times\cdots\times P_{N}}$, $P_{n}<I_{n}$, $n=1,2,\cdots,N$, can be obtained by projecting the original tensor $\mathcal{A}$ onto matrices $\widetilde{\mathbf{U}}^{(n)}$ as

[TABLE]

where $\widetilde{\mathbf{U}}^{(n)}\in\mathbb{R}^{I_{n}\times P_{n}}$ is composed of eigenvectors corresponding to the largest $P_{n}$ eigenvalues of $\mathbf{A}^{(n)}\cdot\mathbf{A}^{(n)^{T}}$.

2.6.2 Texture Tensor Classification

To extract texture features from salt dome boundaries detected in a reference section, we define pairs of patches centered at boundary points that are subimages containing the boundary textures of the reference section and its corresponding GoT map. Since matrices are particular third-order tensors with a third-mode dimension equal to one, we denote patch pairs in a form of tensor pairs as $\{\mathcal{S}_{p_{i}},\mathcal{G}_{p_{i}}\}\subset\mathbb{R}^{I_{1}\times I_{2}\times 1}$, $i=1,2,\cdots,N_{b}$, where $N_{b}$ represents the number of all boundary points and $I_{1}$ and $I_{2}$ define the dimensions of patches along the depth and crossline directions, respectively. Since textures have no bias on either the depth or crossline direction, we commonly assume that $I_{1}$ is equal to $I_{2}$. In addition, to ensure that boundary textures can be captured by square-shaped patches, we determine the edge length of patches as around one-tenth of the larger of the vertical and horizontal resolutions. Furthermore, by evaluating the similarity between patches, we can classify boundary textures and build a set of texture tensors, denoted $\{\mathcal{S}_{k},\mathcal{G}_{k}\}$, $k=1,2,\cdots,N_{t}$.

To clarify the iterative classification strategy, we depict its block diagram in Figure LABEL:fig:classification. Because of the roughly semi-circular shape of detected boundaries in the reference section, we traverse clockwise along the boundary points and define $p_{1}$ as the first point in the bottom-left corner satisfying the dimension constraint of patches. After initializing $\{\mathcal{S}_{1},\mathcal{G}_{1}\}$ with patch pair $\{\mathcal{S}_{p_{1}},\mathcal{G}_{p_{1}}\}$, we search for $p_{2}$ in the $3\times 3$ neighborhood of $p_{1}$ on the basis of the point priority shown in Figure 1(a), where a smaller number represents a higher priority. To facilitate the evaluation of the texture similarity between the patch pair of $p_{2}$, $\{\mathcal{S}_{p_{2}},\mathcal{G}_{p_{2}}\}$, and the initial tensor pair, $\{\mathcal{S}_{1},\mathcal{G}_{1}\}$, we append $\mathcal{S}_{p_{2}}$ and $\mathcal{G}_{p_{2}}$ to $\mathcal{S}_{1}$ and $\mathcal{G}_{1}$, respectively, along the third mode and obtain the updated tensor pair $\{\widetilde{\mathcal{S}}_{1},\widetilde{\mathcal{G}}_{1}\}$. In Figure LABEL:fig:classification, function $(\cdot|\cdot)$ describes an operation that extends tensors along the third direction. Furthermore, by applying the MPCA to $\{\widetilde{\mathcal{S}}_{1},\widetilde{\mathcal{G}}_{1}\}$, we obtain the projection matrix of each mode, denoted $\widetilde{\mathbf{U}}_{\widetilde{\mathcal{M}}_{1}}^{(n)}\in\mathbb{R}^{I_{n}\times P_{n}}$, $(n=1,2)$, $\widetilde{\mathbf{U}}_{\widetilde{\mathcal{M}}_{1}}^{(3)}\in\mathbb{R}^{(I_{1}\times I_{2})\times P_{n}}$, $\mathcal{M}=\{\mathcal{S},\mathcal{G}\}$, in which $[P_{1},P_{2},P_{3}]$ are determined empirically. Because of the high computational efficiency of the sequential Karhunen-Loeve (SKL) algorithm Levey and Lindenbaum, (2000); Hu et al., (2011) on updating projection matrices, we extract $\widetilde{\mathbf{U}}_{\widetilde{\mathcal{M}}_{1}}^{(3)}$, $\mathcal{M}=\{\mathcal{S},\mathcal{G}\}$ from the row spaces of 3-mode unfolding matrices instead of their column spaces. Since eigenvectors in the projection matrices reflect the texture features of tensors, to evaluate the similarity between patches and tensors, we define reconstruction error $e_{c}$ as follows:

[TABLE]

where $p_{i}$ and $k$ represent the indices of the current boundary point and tensor pair, respectively. If $e_{c}$ is less than threshold $T_{e}$, the strong similarity allows us to extend $\{\mathcal{S}_{1},\mathcal{G}_{1}\}$ by appending $\{\mathcal{S}_{p_{2}},\mathcal{G}_{p_{2}}\}$. Otherwise, the great deviation requires us to keep $\{\mathcal{S}_{1},\mathcal{G}_{1}\}$ unchanged and initialize another tensor pair $\{\mathcal{S}_{2},\mathcal{G}_{2}\}$ with $\{\mathcal{S}_{p_{2}},\mathcal{G}_{p_{2}}\}$. By repeating the step above, we obtain classified tensor pairs that contain the local texture information of the detected boundary. Figure 1(b) illustrates one local area of a seismic section and its corresponding GoT map, in which the red curve represents the labeled salt dome boundary. In addition, Figure 1(b) also shows the 3-mode unfolding matrices of a tensor pair extracted from this local area, in which texture patches have a strong correlation.

2.6.3 Tracked Boundary Synthesis

On the basis of texture features extracted from classified texture tensors, we can localize boundary points in predicted sections. We first project the detected boundary in the reference section onto the target section and keep the coordinates of all points unchanged. Then, to identify the tracked point of every projected point $p$, we search among candidate points located along the normal direction of the projected boundary within radius $R_{s}$, which is adaptively determined by the inline number difference between the reference and the current predicted sections. To define the textures of candidate points, we extract pairs of patches centered at these points from the predicted section and its corresponding GoT map, denoted $\{\mathcal{S}_{p,j},\mathcal{G}_{p,j}\}\subset\mathbb{R}^{I_{1}\times I_{2}\times 1}$, $j=1,2,\cdots,(2R_{s}+1)$, where $j$ represents the index of the candidate point. Since in the reference section the patch pair of projected point $p$ belongs to tensor pair $\{\mathcal{S},\mathcal{G}\}$, by comparing the similarity between $\{\mathcal{S}_{p,j},\mathcal{G}_{p,j}\}$ and $\{\mathcal{S},\mathcal{G}\}$, we can determine tracked point $p^{\ast}$. The block diagram of the proposed tracking process is shown in Figure LABEL:fig:localization. The reconstruction error $e_{t}$ is calculated as follows:

[TABLE]

where $\lambda_{\mathcal{S}}=1$ and $\lambda_{\mathcal{G}}=|\log\left(\mathbf{G}[i,j]\right)|$ represent the weights of texture features extracted from seismic sections and corresponding GoT maps, respectively. $[i,j]$ denote the coordinates of the candidate point. Since GoT maps have been normalized to a range of zero to one, a larger $\mathbf{G}[i,j]$ indicates a smaller weight, which can move tracked points towards the real salt-dome boundaries of predicted sections by lowering the reconstruction error. Finally, the candidate point with the smallest reconstruction error is selected as the tracked point. We apply the same process to all boundary points and obtain binary image $\mathbf{B}$, which contains initial tracked points.

Because of the complicated structures of salt domes, it is inevitable that tracked points involve some outliers that may harm the smoothness and continuity of tracked boundaries. To remove these outliers, we apply the median filter to $\mathbf{B}$ and obtain filtered image $\mathbf{\widehat{B}}$ as follows:

[TABLE]

where function $\underset{R_{m}\times R_{m}}{\mbox{median}(\cdot)}$ applies the median filter to the $R_{m}\times R_{m}$ analysis window of point $[i,j]$. Without involving new points generated by median filtering, Equation (26) removes outliers and ensures the accuracy of the remaining tracked points. Furthermore, we connect the remaining tracked points clockwise with straight lines and obtain the initial labeling of the salt dome boundary. Since the connection based on lines may result in a jagged shape of the labeled boundary, which is geologically unreasonable, to delineate the salt-dome boundary more accurately, we need to employ two post-processing steps, dilation and skeletonization. We dilate initially labeled result $\mathbf{L}$ using disk-shaped structural element $\mathbf{H}$ as follows:

[TABLE]

where $\mathbf{H}_{z}$ represents the structural element centered at 2D point $z$. The dilated region can be understood as the locus of the points covered by $\mathbf{H}$ when the center of $\mathbf{H}$ moves inside $\mathbf{L}$. Furthermore, by applying skeletonization Kong and Rosenfeld, (1996) to the dilated region, we can extract the tracked boundary.

2.6.4 Salt Dome Tracking based on Noise-adjusted MPCA

In our previous discussion, we commonly assume that noise does not exist in seismic datasets. However, in real cases, it is possible that migrated seismic sections contain various levels of noise. If we still depend on the MPCA-based tracking method, extracted features corresponding to the most variations of texture tensors may have low signal-to-noise ratios (SNRs). Therefore, the appearance of noise may reduce the accuracy of tracked boundaries. To deal with this problem, on the basis of Roger’s work Roger, (1994), we modify the proposed tracking strategy and introduce noise-adjusted MPCA (NA-MPCA) that ranks principal components (PCs) along different modes by the SNR rather than the variation. Since noise is commonly assumed to be additive, we define image patches extracted from seismic sections as $\mathbf{X}=\mathbf{S}+\mathbf{N}$, where $\mathbf{S}$ is the image patch without noise and $\mathbf{N}$ represents noise that are independent to $\mathbf{S}$. In this paper, we estimate the amplitude of noise using wavelet-based denoising method Chang et al., (2000). In addition, we denote the covariance matrices of $\mathbf{X}$, $\mathbf{S}$, and $\mathbf{N}$ as $\mathbf{\Sigma}_{\mathbf{X}}$, $\mathbf{\Sigma}_{\mathbf{S}}$, and $\mathbf{\Sigma}_{\mathbf{N}}$, respectively. To arrange PCs in the order of SNRs, we need to apply a noise-whitening process that transforms $\mathbf{N}$ to identity covariance matrix $\mathbf{\widehat{N}}$. Therefore, we define transform matrix $\mathbf{P}$ and analyze the noise-whitening process as follows:

[TABLE]

where orthonormal matrix $\mathbf{V}$ contains the eigenvectors of $\mathbf{\Sigma}_{\mathbf{N}}$. To satisfy the constraint in Equation (28), we yield transform matrix $\mathbf{P}=\mathbf{V}\mathbf{\Delta}_{\mathbf{N}}^{-1/2}$. Then, we implement the noise-whitening process in the proposed salt dome tracking method by replacing the image patch $\mathcal{S}_{p_{i}}$, $i=1,2,\cdots,N_{b}$, with $\mathcal{S}_{p_{i}}\times_{1}\mathbf{P}$ and keeping $\mathcal{G}_{p_{i}}$ unchanged because of the blurring effect of noise to GoT maps. The NA-MPCA removes the influence of different noise levels in the reference and predicted sections and ensures the robustness of the proposed tracking algorithm.

2.7 Similarity Index Measurement of Salt Dome Boundaries

By employing the proposed salt dome detection and tracking method, we can semi-automatically delineate salt-dome boundaries in seismic sections. To evaluate the performance of the proposed method, we introduce a salt dome similarity (SalSIM) index that measures the distances between semi-automatically labeled results and the ground truth provided by experienced interpreters. To explain the definition of the SalSIM index, in Figure LABEL:fig:localDist we draw two types of curves for ilustration. The proposed SalSIM index is developed based on the Fréchet distance Alt and Godau, (1995), which can more accurately measure the deviation between two curves than the Haudorff distance Huttenlocher et al., (1993) by taking the continuity of curves into account.

On the basis of the Fréchet distance, we attempt to numerically describe local and global deviations between labeled results and the ground truth using local and global items in the proposed SalSIM index, respectively. To compare the local details of labeled boundaries and the ground truth, we define an analysis window that identifies a pair of curve segments, as the red rectangular shown in Figure LABEL:fig:localDist. Therefore, by moving the analysis window along the ground truth and calculating the Fréchet distance of every pair of local segments, we can obtain a sequence of distances, $\mathbf{d}=\left[d_{i}\right]$, $i=1,2,\cdots,N_{d}$, where $N_{d}$ is the total number of pairs. The mean and standard deviation of $\mathbf{d}$, denoted $\mu_{d}$ and $\sigma_{d}$, respectively, constitute the local item of the proposed SalSIM index. In addition, the global item involves the Fréchet distance of two entire curves, denoted $d_{max}$, which equals the largest distance in $\mathbf{d}$. Therefore, on the basis of obtained statistics, we introduce the SalSIM index as

[TABLE]

where $\alpha$ and $\beta$ are normalization factors determined empirically. Since the exponential function defined on negative real numbers has a range of zero to one, we apply it in the SalSIM index for normalization. Therefore, according to the expression of SalSIM in Equation (29), a greater SalSIM value represents smaller deviation between the semi-automatically labeled salt dome boundary and the ground truth and vice versa.

3 Examples

In this paper, we evaluate the performance of our proposed detection and tracking method on a real seismic volume acquired from the Netherlands offshore F3 block with the size of $24\times 16\mbox{ km}^{2}$ in the North Sea Opendtect, (1987). We focus on a local volume containing a distinguishable salt dome that has the inline number ranging from #389 to #409, the crossline number ranging from #401 to #701, and the time interval between 1,300ms and 1,848ms sampled every 4ms. Figure LABEL:fig:Inline400 illustrates a seismic section, Inline #400, extracted from the local volume.

3.1 Salt Dome Detection

The proposed salt dome detection method delineates the boundaries of salt domes based on the GoT maps of seismic sections, in which the point value represents the dissimilarity of neighboring square windows. To capture the multi-scale texture contrast of neighborhoods, we calculate GoT values by varying the window size from $3\times 3$ to $11\times 11$ and averaging the corresponding dissimilarity measures, denoted $\mathbf{G}$ in Equation (3). Figure 2(a) shows the GoT map of Inline #400, in which likely boundary regions have higher GoT values. On the basis of GoT maps, we utilize the region-growing method to identify salt dome areas and their corresponding boundaries. In the growing process, we tweak threshold $T_{g}$ for each seismic section to yield the best growing result. In Figure 2(b), we compare the green detected boundary with the manually labeled red ground truth and notice that the former overlaps the most part of the latter except the bottom-left and -right corners.

3.1.1 Subjective Comparison of Detected Boundaries

Since GoT maps mainly characterize the texture variations of local neighborhoods in seismic sections, to demonstrate the benefit of the GoT attribute, we compare it with two conventional seismic attributes, the GLCM contrast attribute Haralick et al., (1973) and the gradient attribute Toriwaki and Yoshida, (2009), both of which have been introduced to detect salt dome boundaries Berthelot et al., (2013); Aqrawi et al., (2011). To ensure fair comparison, we apply the same region-growing method to extract salt dome boundaries from different attribute maps. The GLCM describes the distribution of co-occurring grayscale values at a given offset over an image. For each point $[i,j]$, within its $(2R_{d}+1)\times(2R_{d}+1)$ analysis window, by varying the directions and pixel distances of offset $[\Delta i,\Delta j]$, we can obtain a series of GLCMs, denoted $\mathbf{P_{[\Delta i,\Delta j]}}$. Therefore, on the basis of these GLCMs, we define the GLCM contrast attribute $\overline{\mathbf{C}}[i,j]$ as

[TABLE]

where $u$ and $v$ indicate the row and column indices of the GLCM, respectively, and $N_{g}$ represents the number of all possible offsets. In the tested local volume, we obtain the GLCM contrast maps of seismic sections by setting $R_{d}=4$. Figure 2(c) illustrates the GLCM contrast attribute of Inline #400, in which higher contrast values indicate likely boundary areas. By applying the region-growing method to Figure 2(c), we can label the salt dome boundary as Figure 2(d) shows. Similar to the GLCM contrast attribute, the gradient attribute estimated based on 3D Sobel filter can also identify the boundaries of salt domes. To approximate the partial derivatives of different directions, we convolute three $3\times 3\times 3$ kernels of the Sobel filter with the local volume. For each point, the magnitude of the gradient attribute is calculated as follows:

[TABLE]

where $F_{x}$, $F_{y}$, and $F_{z}$ represent the partial derivatives of the crossline, depth, and inline directions, respectively. Figure 2(e) illustrates the gradient map of Inline #400, from which we extract the green salt dome boundary as Figure 2(f) depicts.

By comparing the GoT and GLCM contrast of Inline #400 in Figure 2, we notice that it is not easy to determine which one leads to the more accurate detection of salt dome boundaries because of their comparable performance on highlighting boundary areas. However, without multiscale-based description, the gradient map depending only on $3\times 3$ local regions shows noisy and discontinuous stripes around boundary areas, which limits the accuracy of detected boundaries. Figure 3 illustrates the local regions of Figures 2(b), (d), and (f), in which every column contains the same local regions of salt dome boundaries detected from various attribute maps and every row contain the different local regions of a detected salt dome boundary. The first column illustrates the bottom-left corner of detected boundaries, in which the second one is slightly more similar to the ground truth than the other two local boundaries. However, in second to fourth columns, boundaries detected by the proposed method show the highest accuracy.

3.1.2 Objective Comparison of Detected Boundaries

To quantize subjective perception, we employ the SalSIM index to evaluate the similarity between the ground truth and the detected boundary. Figure LABEL:fig:detection_cmpr shows the SalSIM indices of the salt dome boundaries of Inline #389 to #409 detected by various methods. Boundaries extracted from gradient maps have the lowest accuracy because of noise around boundary areas. In most of the tested seismic sections, boundaries detected by the proposed method have SalSIM indices higher than or equal to those obtained from GLCM contrast maps. Since the computation of GLCMs is time consuming, the proposed method outperforms the GLCM-based detection method on both efficiency and accuracy. Table 1 lists several statistical measures of SalSIM indices, in which the averaged maximum distance (AMD) in pixels represents the mean of the $d_{max}$ of all detected boundaries. The SalSIM indices of boundaries detected by the proposed method has the largest mean, the smallest standard deviation, and the shortest AMD, which proves the superiority of the proposed method.

To further verify the robustness of the proposed method, we apply the GoT-based detection method on the ten sections of the F3 block with crossline number ranging from #834 to #842, inline number ranging from #279 to #600, and the time interval between 1,300ms and 1,848ms sampled every 4ms. Figure LABEL:fig:crossline illustrates one section of this local dataset, and Table 2 shows the performance of various detection methods, in which the proposed GoT-based method still achieves the highest accuracy.

3.1.3 Robustness to Noise

To verify the robustness of the proposed method, we add Gaussian noise to seismic sections and compare the accuracy of salt dome boundaries detected from the corresponding GoT and GLCM contrast maps. The added zero-mean Gaussian noise has standard deviations ranging from $0.01$ to $0.05$, and Figure 4 illustrates the change of SalSIM indices according to the increasing of standard deviations. For one noise level, we extract the salt dome boundaries from noisy seismic sections using either the proposed method or the GLCM-based detection method, and the mean of corresponding SalSIM indices is shown in Figure 4. The error bars of means correspond to an uncertainty equal to one standard deviation, and the dashed curve represents the trend of means. Although the averaged accuracy of salt dome boundaries detected by the proposed method decreases slowly with the increasing of standard deviations, it is still greater than those of boundaries detected by the GLCM-based method at the most of noisy levels. We notice that with the increasing of noise the latter fluctuates around a certain SalSIM index rather than keep decreasing, the reason for which is that the quantization step in the calculation of the GLCM, functioning as a built-in noise filter, can weaken the influence of noise on the GLCM contrast maps. To fairly compare the robustness of these two methods, we apply an edge-preserving smoothing filter, the bilateral filter Tomasi and Manduchi, (1998), to noisy seismic sections and obtain the statistical measures of corresponding SalSIM indices in Figure 4(a). We notice that with the extra denoising operation, the GoT-based method shows higher accuracy than the GLCM-based method. If the noise level keeps increasing, depending only on the denoising operation we may not be able to label satisfied boundaries from migrated seismic sections. Therefore, geophysicists commonly prefer to regenerate migrated seismic sections with higher quality from raw seismic data by reducing noise that appears in every preprocessing step Chopra and Marfurt, (2014).

3.2 Salt Dome Tracking

In the local seismic volume we define Inline #400 as a reference section, the salt dome boundary of which has been labeled by experienced interpreters or computer-aided detection methods. We can track the reference boundaries through the seismic volume and synthesize salt dome boundaries in the neighboring twenty predicted sections ranging from Inline #389 to #409. Each point at the reference boundary corresponds to a pair of $31\times 31$ image patches extracted from the reference section and its corresponding GoT map. To acquire texture features from the local regions of the reference boundary, we group all these patch pairs into tensor pairs on the basis of the block diagram shown in Figure LABEL:fig:classification. The dimensions of texture features in three modes are $[15,15,5]$, and threshold $T_{e}$ on the reconstruction error is $2.3$. Furthermore, in predicted sections the proposed tracking method searches along the normal direction of the projected point and identifies the position of the tracked point by comparing the similarity between the patch pairs of candidate points and tensor pairs built from the reference section. The localization of tracked points is implemented automatically, which improves interpretation efficiency. Finally, we remove noisy points in predicted sections with $2\times 2$ median filters and connect remaining points to label the salt dome boundary under the shape constraint of salt domes. Figure 5(a) compares the green tracked salt dome boundary in Inline #391, synthesized based on the manually labeled reference boundary in Inline #400, with the red ground truth. We notice that these two curves almost overlap except for several local regions.

3.2.1 Subjective Comparison of Tracking Methods

To prove that tracking accuracy can be increased by involving the GoT attribute and the texture features of all modes, we need to compare the proposed method with three other tracking strategies. The first one is to implement the tracking process depending only on texture features extracted from vectorized patches rather than those from third-order tensors. The second one is to synthesize tracked boundaries within the framework of the proposed method, but without involving the GoT attribute. The third strategy is the same as the proposed method except for using GLCM contrast maps instead of GoT maps. We rename these three tracking strategies as the tracking method based on vectorization, the tensor-based tracking method without GoT maps, and the tensor-based tracking method with GLCM contrast maps. Figures 5(b) to (d) compare the red ground truth labeled manually with green salt dome boundaries labeled by the three tracking strategies mentioned above, and the tracked boundary in Figure 5(c) shows the greatest deviation from the ground truth. Figure 6 illustrates the local regions of Figures 5(a) to (d), in which every column contains the same local regions of salt dome boundaries synthesized based on various tracking methods and every row contains the different local regions of a tracked salt dome boundary. By comparing the local regions of tracked boundaries with the ground truth, we notice the boundary in the first row synthesized by the proposed tracking method is the most similar to the ground truth.

3.2.2 Objective Comparison of Tracking Methods

To objectively verify our conclusion, we synthesize tracked boundaries in seismic sections ranging from Inline #389 to #409 using various tracking strategies and plot their SalSIM indices in Figure LABEL:fig:tracking_cmpr. The horizontal and vertical axes represent the inline number and the SalSIM index, respectively. Although tracking method based on vectorization can achieve comparable accuracy with the proposed method in Inline #401 to #409, the SalSIM indices of salt dome boundaries in Inline #395 to #389 synthesized by the former one decrease more quickly than those obtained by the latter one. In addition, dashed curves in Figure LABEL:fig:tracking_cmpr shows the distribution trends of SalSIM indices. Although the SalSIM indices of four tracking strategies have similar trends that drop with the increasing offsets between predicted sections and the reference section, the proposed tracking method outperforms other tracking strategies particularly in Inline #389 to #399. Table 3 contains several statistical measures of SalSIM indices in Figure LABEL:fig:tracking_cmpr. The mean of SalSIM indices and the AMD evaluate the accuracy of the tracking method, and the standard deviation determines robustness. The proposed method has the greatest mean, the smallest standard deviation, and the shortest AMD, which proves the superiority of the proposed method.

3.2.3 Combination of Proposed Detection and Tracking Methods

In the previous sections, tracked boundaries are synthesized based on the boundary labeled by interpreters in the reference section. To further reduce human intervention in the tracking process, in the local seismic volume, we can synthesize the tracked boundaries of Inline #389 to #409 on the basis of the reference boundary in Inline #400 labeled by the proposed detection method. Figure 7 compares the tracked and detected green boundaries with the red ground truth. In Figure 7(a) and (b), the tracked boundary of Inline #392 is more similar to the ground truth than the detected result. In contrast, in Figure 7(c) and (d), since the bottom-left corner of the tracked boundary has great deviation to the ground truth, the detected boundary in Inline #408 seems to be more accurate. However, the appearance of the great deviation is caused by the deviation in the reference boundary shown in Figure 2(b) rather than our tracking method. To objectively compare the tracked and detected boundaries, we illustrate the corresponding SalSIM indices in Figure LABEL:fig:tracking_onDetection_cmpr. Since we employ $\lambda_{\mathcal{C}}$ as an important weight that helps move tracked points to the boundary area, the tracked boundaries may have higher SalSIM indices than the detected ones, particularly in Inline #395 to #399. The mean of the SalSIM indices of the tracked boundaries is $0.9362$, which is almost the same as those of the detected boundaries, $0.9348$. It shows that the proposed tracking method is robust and efficient enough to delineate salt dome boundaries.

3.2.4 Influence of Noise on Proposed Tracking Method

Since in some cases the noise level in the reference section may be different from those of predicted sections, to eliminate the influence of noise on the tracked process, we propose the noise-adjusted tracking method, which selects PCs corresponding to the greatest SNR other than the greatest variation. In the local volume, we manually add gaussian noise with different variances to seismic sections. In the reference section Inline #400, the added zero-mean gaussian noise has a variance of $0.01$. In contrast, the noise level increases to $0.05$ in predicted sections. We utilize the manually labeled reference boundary to synthesize tracked boundaries in Inline #389 to #409. In Figure LABEL:fig:noiseAdjusted, We compare boundaries synthesized by the tensor-based tracking method with or without noise adjusted. With the increasing offset between the reference and predicted sections, the appearance of noise magnifies the error in the localization of tracked points. Since the noise-adjusted tracking method is robust to noise, the difference of two curves in Figure LABEL:fig:noiseAdjusted grows with the increasing offset. The tracked boundaries obtained by the former one have a mean of SalSIM indices, $0.9468$, which is greater than that of the boundaries synthesized by the latter one, $0.9396$. It shows that the proposed noise-adjusted tracking method can alleviate the influence of noise and enhance the robustness of the tracking process.

4 Conclusion

In this paper, we propose noise-robust algorithms for detection and tracking of salt domes in post-migrated seismic volumes. By utilizing a texture-based gradient and a 2D DFT-based dissimilarity measure, we are able to accomplish noise-robust salt dome detection with both high accuracy and excellent efficiency. Then, by tracking a detected boundary from an initial seismic section through a 3D volume, we further improve the efficiency of the labeling process. The tensor-based noise-adjusted subspace learning approach proved very useful for the tracking of salt dome boundaries, especially when noise of varying levels is present. Combining together, the detection and tracking algorithms yield a practical salt dome labeling system that is effective, efficient, and robust.

5 ACKNOWLEDGMENTS

This work is supported by the Center for Energy and Geo Processing (CeGP) at Georgia Tech and King Fahd University of Petroleum and Minerals (KFUPM).