Tensor-based subspace learning for tracking salt-dome boundaries

TL;DR

This paper introduces a tensor-based subspace learning method for tracking salt-dome boundaries in seismic volumes, improving detection accuracy by leveraging texture information and tensor analysis.

Contribution

It presents a novel tensor-based approach that enhances salt-dome boundary tracking in seismic data, outperforming existing methods.

Findings

Outperforms state-of-the-art salt-dome detection methods

Utilizes texture tensors for improved boundary identification

Achieves higher accuracy in seismic volume analysis

Abstract

The exploration of petroleum reservoirs has a close relationship with the identification of salt domes. To efficiently interpret salt-dome structures, in this paper, we propose a method that tracks salt-dome boundaries through seismic volumes using a tensor-based subspace learning algorithm. We build texture tensors by classifying image patches acquired along the boundary regions of seismic sections and contrast maps. With features extracted from the subspaces of texture tensors, we can identify tracked points in neighboring sections and label salt-dome boundaries by optimally connecting these points. Experimental results show that the proposed method outperforms the state-of-the-art salt-dome detection method by employing texture information and tensor-based analysis.

Equations

$$\overline{C}[x,y]=\sum _{\left[\Delta x,\Delta y\right]}\left(\frac{1}{N_g}\sum _i\sum _j(i-j{)}^2{G}_{\left[\Delta x,\Delta y\right]}[i,j]\right).$$

Full text

**Index Terms— ** salt-dome tracking, contrast attribute, texture tensors, and subspace learning

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 estimate the possible positions of reservoir regions, experienced interpreters need to accurately delineate the boundaries of salt domes in collected 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 computer-aided methods that detect salt-dome boundaries using graph theory and image-processing techniques. Lomask et al. [1] defined seismic sections as weighted undirected graphs and applied the normalized cut image segmentation (NCIS) to globally optimize the delineation of salt-dome boundaries. The method in [2], an extension of [1], employed local dips in the weight matrix and utilized bound constraints to remove boundary artifacts. Similarly, Harpert et al. [3] introduced the modified NCIS by combining multiple seismic attributes with adaptive weights. Although these NCIS-based methods can be implemented in parallel, their high computational cost limits their future application on high-resolution or three-dimensional (3D) seismic data. To improve the efficiency of global segmentation, Harpert et al. [4] proposed to detect the boundaries of salt domes using the pairwise region comparison based on the minimum spanning tree, which reduces the algorithm complexity from $\mathcal{O}\left(n^{2}\right)$ to $\mathcal{O}\left(n\log n\right)$. In [5], Aqrawi et al. applied a dip-guided 3D Sobel filter to detect salt-dome boundaries. Recently, Berthelot et al. [6] have introduced to segment salt domes by testing the combination of multiple seismic attributes–texture attributes, frequency-based attributes, and dip attributes–in the supervised Bayesian classification model. In [7], to characterize texture differences between salt domes and surrounding geological structures, Hegazy and AlRegib proposed the directionality attribute using the moment of inertia tensors of gradient components

In this paper, we propose a new algorithm for tracking the boundaries of salt domes through 3D seismic volumes. We first derive the contrast map from an initial image, the salt-dome boundary of which has been accurately labeled by experienced interpreters or practical detection methods. On the basis of the initial image and its corresponding contrast map, we build third-order texture tensors along the labeled boundary. Then, by employing features extracted from the subspace of texture tensors, we identify the positions of tracked boundary points in neighboring images. Finally, the optimal connection of tracked points synthesizes tracked boundaries. To maintain consistency with the terminology in video coding, we rename initial and neighboring images as reference and predicted sections, respectively, in the rest of this paper.

2 The Proposed Method

The block diagram of the proposed method is shown in Fig. 1. In the following subsections, we explain each step of the proposed pipeline in detail.

2.1 Texture Attribute Extraction

Since rock strata surrounding salt domes are commonly formed by sedimentary rocks such as shale and limestone, the various properties of salt and sedimentary rocks result in their distinctive appearances in seismic sections. Fig. 2(a) illustrates a reference section obtained from the North Sea dataset [8] with inline number 399, in which a blue curve indicates the boundary of a salt dome. By observing textures along the labeled boundary, we notice that in contrast to the homogeneous textures of the salt body, neighboring rock strata have varying, more complicated textures. Such a prominent texture contrast represents an important feature of salt-dome boundaries.

To accurately characterize boundary textures, we employ a gray level co-occurrence matrix (GLCM) [9] that describes the distribution of co-occurring grayscale values at a given offset over an image. For each point in the reference section, we calculate its corresponding GLCM within an $(2R_{d}+1)\times(2R_{d}+1)$ analysis window as

$G_{\left[\Delta x,\Delta y\right]}[i,j]=\mbox{Pr}\left(S[x,y]=i,S\left[x+\Delta x,y+\Delta\right]=j\right)$

, where $x$ and $y$ correspond to the indices along the crossline and depth directions, respectively, and $S[x,y]$ represents the intensity of seismic signals. $\mbox{Pr}(\cdot)$ denotes the probability of an event. In addition, $i$ and $j$ correspond to grayscales under specified quantization, and $\left[\Delta x,\Delta y\right]$ defines the predetermined offset, which is limited to four directions, $\left\{0^{\circ},45^{\circ},90^{\circ},135^{\circ}\right\}$ and has a pixel distance within the range of $1$ to $R_{d}$. By selecting various directions and pixel distances, we can derive $N_{g}$ GLCMs for each point that contain multi-directional and multi-scale texture information. From obtained GLCMs, we extract the contrast attribute that describes the variation of textures as follows:

[TABLE]

Fig. 2(b) illustrates the contrast map of Inline #399. Blue areas with contrast values close to zero represent homogeneous textures. In contrast, green and red stripes with contrast values close to one indicate great texture variations around the salt-dome boundary.

2.2 Texture Tensor Classification

2.2.1 Tensors and Multi-linear Analysis

In the field of multi-linear algebra, tensors are commonly used to describe high-dimensional ($N\geq 3$) data. For an $N$-th order tensor $\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$, each order represents a mode of $\mathcal{A}$. We can unfold $\mathcal{A}$ along the $n$-mode and 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})}$ [10]. 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}$, is new tensor $\mathcal{B}\in\mathbb{R}^{I_{1}\times\cdots\times I_{n-1}\times J_{n}\times I_{n+1}\cdots\times I_{N}}$, the entries of which are calculated as $\mathcal{B}(i_{1}\cdots i_{n-1}j_{n}i_{n+1}\cdots i_{N})=\sum_{i_{n}}\mathcal{A}(i_{1}\cdots i_{N})\cdot\mathbf{U}(j_{n},i_{n})$. Tensors can be decomposed 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 and its column vectors span the column space of $n$-mode unfolding matrix $\mathbf{A}^{(n)}$. Since image and video tensors commonly have high dimensions along each mode, we apply multi-linear principal component analysis (MPCA) [11] to extract a tensor subspace with lower dimensions that captures the most variation of input data. The projection operation can be expressed as $\mathcal{C}=\mathcal{A}\times_{1}\widetilde{\mathbf{U}}^{(1)^{T}}\times_{2}\widetilde{\mathbf{U}}^{(2)^{T}}\cdots\times_{N}\widetilde{\mathbf{U}}^{(N)^{T}}\in\mathbb{R}^{P_{1}\times P_{2}\times\cdots\times P_{n}}$. $\widetilde{\mathbf{U}}^{(n)}$, $n=1,2,\cdots,N$, represent $I_{n}\times P_{n}$ matrices, the column vectors of which are orthonormal and composed of eigenvectors corresponding to the largest $P_{n}$ eigenvalues of $\mathbf{A}^{(n)}\cdot\mathbf{A}^{(n)^{T}}$.

2.2.2 Adaptive Classification of Texture Tensors

To fully utilize texture information along salt-dome boundaries, we extract pairs of patches centered at boundary points, $p_{i}$, $i=1,2,\cdots,N_{b}$, from the reference section and its corresponding contrast map. $N_{b}$ defines the number of all boundary points. Since matrices are special third-order tensors with $I_{3}=1$, patch pairs, denoted $\left\{\mathcal{S}_{p_{i}},\mathcal{C}_{p_{i}}\right\}\subset\mathbb{R}^{I_{1}\times I_{2}\times 1}$, contain two types of third-order texture tensors. $I_{1}$ and $I_{2}$ represent the dimensions of patches along the crossline and depth directions, respectively, which are empirically determined based on the complexity of boundary textures. Since the formation of salt domes lasts hundreds of years, local areas commonly have strong uniformity, which is an important constraint to classifying textures along the boundaries of salt domes.

In this section, we attempt to build pairs of tensors, $\{\mathcal{S}_{k},\mathcal{C}_{k}\}$, $k=1,2,\cdots,N_{t}$, constructed by the grouped patches of boundary points containing similar textures. The pseudocode of the proposed classification method is listed in Algorithm 1.

We first initialize $\{\mathcal{S}_{1},\mathcal{C}_{1}\}$ using patches $\{\mathcal{S}_{p_{1}},\mathcal{C}_{p_{1}}\}$. Since in seismic sections the boundaries of salt domes commonly appear in roughly semi-circular shapes, to transverse all boundary points clockwise, we define $p_{1}$ as the first point in the bottom-left corner satisfying the dimension constraint of patches. To identify the next boundary point $p_{2}$, we search clockwise in the $3\times 3$ neighborhood of $p_{1}$. Fig. 3(a) illustrates the priority of neighboring points, in which a smaller number represents higher priority. After determining the position of $p_{2}$, we append the corresponding patches $\{\mathcal{S}_{p_{2}},\mathcal{C}_{p_{2}}\}$ to the current tensor pair $\{\mathcal{S}_{1},\mathcal{C}_{1}\}$ along the $z$ direction and denote updated tensors as $\widetilde{\mathcal{S}}_{1}$ and $\widetilde{\mathcal{C}}_{1}$, respectively. By applying the MPCA on tensors $\{\widetilde{\mathcal{S}}_{1},\widetilde{\mathcal{C}}_{1}\}$, we obtain the projection matrices of all modes, 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{C}\}$, in which $[P_{1},P_{2},P_{3}]$ are determined empirically. For 3-mode unfolding matrices, to utilize the high computational efficiency of the sequential Karhunen-Loeve (SKL) algorithm [12, 13], we extract the projection matrices of row spaces. On the basis of reconstruction error $e$, we can evaluate the similarity between $\{\mathcal{S}_{p_{2}},\mathcal{C}_{p_{2}}\}$ and $\{\mathcal{S}_{1},\mathcal{C}_{1}\}$. If $e$ is less than threshold $T_{e}$, high similarity leads to the extension of current tensors with $\{\mathcal{S}_{p_{2}},\mathcal{C}_{p_{2}}\}$ appended. Otherwise, $\{\mathcal{S}_{p_{2}},\mathcal{C}_{p_{2}}\}$ initializes the next new tensor pair $\{\mathcal{S}_{2},\mathcal{C}_{2}\}$. By repeating the steps above, we obtain classified texture tensors. Fig. 3(b) illustrates the patches of two tensor pairs, the strong correlation among which can be better captured by tensor-based analysis.

2.3 Tracked Boundary Synthesis

On the basis of projection matrices extracted from classified tensors, we can localize tracked boundary points in predicted sections. We first project the labeled boundary onto the target predicted section and keep the coordinates of the boundary points unchanged. Then, at each projected boundary point $p$, we search for the tracked point along its normal direction within a $(2R_{s}+1)$ window, in which $R_{s}$ represents the inline number difference between the reference section and the predicted section. Each candidate point corresponds to a pair of texture patches, denoted $\{\mathcal{S}_{p,j},\mathcal{C}_{p,j}\}\subset\mathbb{R}^{I_{1}\times I_{2}\times 1}$, $j=1,2,\cdots,(2R_{s}+1)$, extracted from the predicted section and the corresponding contrast map. We assume that the texture patches of current projected point $p$ belong to tensor pair $\{\mathcal{S},\mathcal{C}\}$. The localization of tracked point $p^{\ast}$ is explained in the following pseudocode:

By evaluating the similarity between current tensors $\{\mathcal{S},\mathcal{C}\}$ and the patch pairs of the candidate points, we can localize tracked point $p^{\ast}$. $\lambda_{\mathcal{S}}$ and $\lambda_{\mathcal{C}}$ represent the weights of $\mathcal{S}$ and $\mathcal{C}$ in reconstruction error $e$. Since the ranges of seismic sections and contrast maps have been normalized between [math] and $1$, we define $\lambda_{\mathcal{S}}=1$ and $\lambda_{\mathcal{C}}=|\log(\overline{C})|$. Therefore, candidate points with higher contrast values have lower weights in $e$. Fig. 4(a) illustrates tracked points in Inline #409. After removing noisy points using the $2\times 2$ median filter, we connect the remaining tracked points to synthesize the tracked boundary. In Fig. 4(b), green and blue curves represent the tracked and projected boundary in Inline #409, respectively.

3 Experimental Results

In this paper, we apply the proposed salt-dome tracking method on a 3D real seismic dataset acquired from the Netherlands offshore F3 block with the size of $24\times 16$ km2 in the North Sea [8]. To illustrate the performance of the proposed method, we focus on a local volume of the dataset containing discernible salt-dome structures. The tested volume has an inline number ranging from #151 to #501, a crossline number ranging from #401 to #701, and a time direction ranging from 1,300ms to 1,848ms with a step of 4ms. Both Figs. 2(a) and 4(b) illustrate seismic sections extracted from the local volume.

As we mentioned in previous sections, on the basis of the labeled salt-dome boundary in the reference section Inline #399, we attempt to synthesize tracked salt-dome boundaries in neighboring predicted sections ranging from Inline #389 to #409. To acquire the contrast map of each seismic section, we derive the GLCMs of points from $9\times 9$ neighborhoods by selecting various directions and pixel distances. Then, we extract $31\times 31$ patches from the boundary areas of the reference section and its corresponding contrast map and group similar patches into texture tensors using the proposed classification method in Algorithm 1. The dimension of the subspace is $[15,15,5]$, and threshold $T_{e}$ for the reconstruction error in the classification of texture tensors is $3$. Furthermore, in the predicted sections, we search along the normal directions of projected points and localize tracked points by comparing the similarity between the patches of candidate points and texture tensors built from the reference section. The searching process and the connection of tracked points are implemented automatically, which improves interpretation efficiency.

To evaluate the necessity of the contrast attribute and the tensor-based learning process in salt-dome tracking, we track the manually labeled boundary through predicted sections using three methods, the proposed method, the proposed method without the contrast attribute, and the vector-based subspace learning method. The last one refers to a method that extracts features only from the subspace of vectorized texture patches, which constitute the 3-mode unfolding matrices of tensors in the proposed method. To objectively evaluate the similarity between tracked boundaries and the ground truth, we calculate Fréchet distances [14] between the segments of curves and normalize the averaged distance as a similarity index. Fig. 5 illustrates the similarity indices of salt-dome boundaries labeled by the three methods. The proposed method synthesizes tracked boundaries with the highest accuracy, particularly in predicted sections with larger inline offsets. For example, in Fig. 7(a), the green tracked boundary in Inline #409 almost coincides with the ground truth, labeled by the red curve. In addition, to verify the robustness of the proposed method in practical cases, we detect the salt-dome boundary of the reference section by Aqrawi’s method [5] and synthesize tracked boundaries in predicted sections. In Fig. 6, tracked boundaries in most of the predicted sections have higher similarity indices than those detected by Aqrawi’s method, which indicates the higher reliability of the proposed method. In Figs. 7(b) and (c), we compare tracked and detected boundaries (green) with the ground truth (red), respectively, and find that the former has more deviations.

4 Conclusion

In this paper, we developed the tracking of salt-dome boundaries using a tensor-based subspace learning algorithm. We built texture tensors by classifying image patches along the boundary areas of the reference section and its corresponding contrast map. With features extracted from the subspace of texture tensors, we identified the positions of tracked points by evaluating their similarity to texture tensors. The connection of tracked points synthesizes tracked boundaries. Experimental results showed that to obtain more accurate labeling, the tracking method needs to employ contrast maps and tensor-based analysis. In practical cases, the proposed method also showed higher reliability than the state-of-the-art detection method.

5 Acknowledgements

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