Characterization of migrated seismic volumes using texture attributes: a comparative study

TL;DR

This study compares various texture attributes derived from image processing techniques to assess their effectiveness in characterizing migrated seismic volumes within an image retrieval framework.

Contribution

It provides a comprehensive comparison of frequency and space domain texture attributes for seismic volume characterization, aiding seismic interpretation.

Findings

Steerable pyramid and curvelet attributes show high retrieval accuracy.

Local binary pattern offers computational efficiency.

The study enhances understanding of texture attribute effectiveness in seismic analysis.

Abstract

In this paper, we examine several typical texture attributes developed in the image processing community in recent years with respect to their capability of characterizing a migrated seismic volume. These attributes are generated in either frequency or space domain, including steerable pyramid, curvelet, local binary pattern, and local radius index. The comparative study is performed within an image retrieval framework. We evaluate these attributes in terms of retrieval accuracy. It is our hope that this comparative study will help acquaint the seismic interpretation community with the many available powerful image texture analysis techniques, providing more alternative attributes for their seismic exploration.

Figures

Tables

Table 1. Table 1: Retrieval performance for different attributes.

	P@20	P@50	MAP	RA	AUC
SP	1.000	0.965	0.928	0.866	0.965
CT	1.000	0.968	0.954	0.914	0.988
LBP	0.999	0.953	0.932	0.871	0.967
LRI	0.997	0.977	0.953	0.896	0.968
SeiSIM	1.000	0.992	0.974	0.943	0.990

Table 2. Table 2: Comparison of computation time in seconds.

SP	LBP	LRI	CT	SeiSIM
0.2137	0.1638	2.3946	0.1201	1.1872

Equations

$$J=\lceil {\log }_2\min (N_1,N_2)-3\rceil ,$$

$$K(j)=16\times 2^{\lceil (j-1)/2\rceil }.$$

$$LB{P}_{R,P}\left[i_c,j_c\right]=\sum _{p=0}^{P-1}s\left(I_c-I_p\right)\cdot 2^p,$$

$$d_{SC}(H^f,H^g)=\sum _{i=1}^M{\left(\sqrt{H^f(i)}-\sqrt{H^g(i)}\right)}^2,$$

$$D(f,g)=\sum _{j=1}^{J-1}{d}_{SC}(H_j^f,H_j^g)$$

$$Q_{STSIM-1}(\text{x},\text{y})=[l(\text{x},\text{y}){]}^{\frac{1}{4}}[c(\text{x},\text{y}){]}^{\frac{1}{4}}[a_h(\text{x},\text{y}){]}^{\frac{1}{4}}[a_v(\text{x},\text{y}){]}^{\frac{1}{4}}$$

$$Q_{DM}(D{M}_1,D{M}_2)=[a_{DM}^h(D{M}_1,D{M}_2){]}^{\frac{1}{2}}[a_{DM}^v(D{M}_1,D{M}_2){]}^{\frac{1}{2}}$$

$$SeiSIM=[Q_{STSIM-1}{]}^{\frac{1}{2}}[Q_{DM}{]}^{\frac{1}{2}}.$$

Full text

1 Introduction

In image processing, texture attributes are quantities generated from a texture pattern that capture the unique spatial distribution of the pixel intensities Gonzalez and Woods, (2006). Such texture attributes have been employed in some seismic interpretation applications. For instance, the gray level co-occurrence matrix (GLCM) was applied to salt dome detection Berthelot et al., (2013) and deep-marine facies discrimination Gao, (2007); Hilbert transform features were utilized for seismic image segmentation Pitas and Kotropoulos, (1992); Gabor filters were adopted for seismic image segmentation as well Röster and Spann, (1998).

Although these applications were successful, we believe texture attributes can be even more useful for seismic exploration. Given that migrated seismic volumes are textural in nature, texture attributes have the potential of serving as local descriptors that characterize the migrated data. Such descriptors are essential for a computer-assisted understanding of the ”subsurface scene,” which can help pinpoint spots that are of more interest to an interpretor. In recent years, many powerful texture analysis techniques have been developed in the image processing community. Most of them have not been exposed adequately to exploration geophysicists yet. Therefore, in this paper, we conduct a comparative study examining several typical texture attributes with respect to their capability of characterizing a seismic volume. Please note that, although some of the techniques discussed in this paper may have been introduced to the community before, they will be explored in a different context, i.e., as texture attributes for general description of migrated seismic volumes.

The attributes we examine here can be categorized into frequency or space domain techniques. For the frequency domain, we explore the steerable pyramid (SP) Simoncelli et al., (1992), and the curvelet transform (CT) Candés et al., (2005). These are two typical extensions of the standard wavelet transform (WT), which is probably the most popular technique for spectral content-based image and texture analysis. SP improves over WT Chui, (2014) in that it achieves translation-invariance and orientation-invariance by dropping the orthogonality constraint. In addition to SP, several other multi-scale techniques have been developed more recently to overcome an inherent shortcoming of WT, which is the lack of directionality. Among such techniques, CT gained popularity particularly for seismic data processing, mainly because it describes with high efficiency signals of smooth curves.

For the space domain, we study the local binary pattern (LBP) Ojala et al., (2002) and the local radius index (LRI) Zhai et al., (2013). LBP generates binary numbers in local areas according to gray-scale intensity differences among image pixels. The pattern of the binary numbers is captured for a local area as an occurrence histogram of the numbers. Since its introduction, LBP (and its variants) has been widely used for very successful texture analysis due to its robustness and computational efficiency. LRI is one of the most recently developed techniques following the concept of local patterns. It measures distances between local edges along all directions. The resulting histograms serve very well as a local texture discriptor. We note that although some of the older techniques such as GLCM have been used for seismic applications for a long time, they are not included in this paper. The major reason is that, first, they are already familiar to the community; and second, the newer techniques are more powerful.

To determine the capability of the attributes for characterizing migrated seismic data, we perform retrieval experiments, where a dataset is searched to identify images that are similar to a given query image in content, or in seismic structure when our dataset is concerned. The structural similarity between seismic images is measured on the basis of the attributes. Thus, the attributes’ capability of characterizing the data (or the seismic structures contained therein) can be evaluated in terms of retrieval accuracy. Higher accuracies indicate that the associated attributes are better able to distinguish between different seismic structures. To achieve a more informed comparison, we choose SeiSIM Long et al., (2015) as our benchmark. This method considers both the frequency domain statistics obtained using SP, and the space domain statistics calculated from discontinuity maps instead of original seismic images. One of very few techniques designed for evaluating similarity between seismic images, SeiSIM was demonstrated to be capable of capturing differences in geological structures with reliability. We hope this comparative study will help acquaint the seismic interpretation community with the many available powerful image texture analysis techniques, providing more alternative attributes for their seismic exploration.

2 Methods

2.1 Frequency domain techniques

2.1.1 Steerable pyramid (SP)

SP is a multiresolution image representation developed by Simoncelli et al., (1992). As illustrated in Figure LABEL:fig:pyramid, the technique first decomposes a given image into a highpass subband and a lowpass subband. Then it processes the lowpass subband, obtaining a series of bandpass subbands and another lowpass subband. The bandpass subbands reveal image details along various orientations. The newly obtained lowpass subband is subsampled and then further processed in a similar manner to yield orientational details at a coarser spatial scale. Such recursive decomposition eventually yields a pyramid of subbands, representing the original image along different orientations at different scales. Histograms of the coefficients from the decomposition can be established for each subband, which capture the statistical characteristics of the coefficients. The histograms are further examined for retrieval purpose, details of which will be discussed later in the experiments.

\plot

pyramidwidth=0.48 Illustration of a SP decomposition with $K$ scales and $N$ orientations at each scale. In this paper, we set $K=4$ and $N=8$.

2.1.2 Curvelet transform (CT)

CT is also a directional multiscale decomposition, first introduced by Candés et al., (2005). It works by first applying the two-dimensional fast Fourier transform (2-D FFT) to an image, and then dividing the frequency plane into small sections (or wedges) corresponding to multiple scales and orientations. The total number of scales in the curvelet tiling, $J$, is dependent on the size of the image as

[TABLE]

where $N_{1}$ and $N_{2}$ are the image height and width in pixels, respectively; and $\lceil\cdot\rceil$ is the ceiling function. The number of orientations at scale $j\geq 1$, $K(j)$, is given by:

[TABLE]

Once the frequency plane is partitioned (see Figure 1 for an example), curvelet coefficients are generated by applying the 2-D IFFT to each wedge (after smoothing). Since the FFT of real images is symmetric around the origin, only two quadrants of the Fourier spectrum are necessary for obtaining the coefficients. Again, histograms are formed for coefficients in each subband and used in the retrieval experiments.

2.2 Space domain techniques

2.2.1 Local binary pattern (LBP)

LBP Ojala et al., (2002) describes the local spatial structure of textures by thresholding the neighborhood of each pixel and defining the result as a binary number. Mathematically, the LBP operator is expressed as

[TABLE]

where $P$ represents the number of points in the neighborhood with radius $R$, $\left[i_{c},j_{c}\right]$ indicates the coordinates of the center point, and $I_{c}$ and $I_{p}$ denote the intensity of the center and neighboring points, respectively. Function $s(\cdot)$ has a value of $1$ if $I_{c}\geq I_{p}$. Otherwise, the value of $s$ is [math]. Since the LBP operator encodes only the signs of the difference between the center and neighboring points, however, the information of difference magnitude has been discarded. To overcome this problem, Guo and Zhang, (2010) proposed completed LBP (CLBP), where three components are considered as follows. First, CLBP_C encodes the center pixel intensity into a binary number. Then, CLBP_S and CLBP_M are generated using the difference between the center and its neighbors, with the former encodes the sign of the difference and the latter the magnitude. Histograms of the three components are concatenated into one feature vector to describe the local texture pattern. In fact, CLBP_S is exactly the same as LBP. In this paper, we use CLBP instead of the original LBP. We set $P=20$ and $R=3$.

2.2.2 Local radius index (LRI)

LRI characterizes a texture pattern by the distribution of distances between adjacent edges along a certain orientation Zhai et al., (2013). A local index can be computed for each image pixel in two different ways, resulting in two variations of LRI. For LRI-A, inter-edge distance (i.e., width of adjacent smooth regions) in each given direction is calculated; while for LRI-D, the distance from pixels to the nearest edge (i.e., boundary of next smooth region) is adopted. In Figure 2, an example is given illustrating how to compute LRI-A and LRI-D. Histograms calculated from LRI-A and LRI-D are concatenated to form one vector for the retrieval experiments.

There are two key parameters for LRI calculation. One is a threshold $T$ to determine edge. The other is a range $K$ to control how far to search for edges. In our experiments, we use $T=\sigma/2$ (where $\sigma$ is the local standard deviation) and $K=3$. We obtain a $(2K+1)$-bin histogram in each of eight directions.

3 Experiments

3.1 Data

We created a dataset consisting of 400 images extracted from the public dataset of Netherlands offshore F3 block with the size of $24\times 16km^{2}$ in the North Sea. The images are grouped into 4 classes according to their geological structures, namely, clear horizon, chaotic horizon, faults, and salt dome. Each group includes 100 images. All images are $150\times 300$ in pixels. Example images are given in Figure 3.

3.2 Experimental procedure

We conduct retrieval experiments on the dataset where the similarity is measured based on the various texture attributes discussed above. For a complete evaluation, every image in the dataset is selected once as the query image. It is then compared with the rest of images, ranking them according to the similarity score in a descending order. Such ranking is obtained for all images. Then the performance can be assessed using all rankings.

The images are first normalized. Then for each image, a specific texture attribute is computed and histograms are formed for the attributes. The images are compared in terms of a certain distance measure calculated between corresponding histograms. We used squared chord distance (SCD) Cha, (2007) to generate all the results presented in this paper. We also tested other distance measures such as the Kullback-Leibler divergence (KLD) Kullback and Leibler, (1951), and noticed no significant difference in the results. SCD is defined as:

[TABLE]

where $H^{f}$ and $H^{g}$ are the two histograms to be compared, from images or subbands $f$ and $g$, respectively. $H^{f}(i)$ (or $H^{g}(i)$) is the $i^{th}$ bin in the histogram $H^{f}$ (or $H^{g}$). $M$ is the number of bins in the histogram.

The overall distance between the two images is given by combining SCD calculated over all $J$ pairs of histograms:

[TABLE]

The overall distance can then be easily converted to similarity values bounded between [math] and $1$.

3.3 Evaluation metrics

We adopt the following metrics that are typically used to evaluate retrieval performance (Zujovic et al., (2013); Zhai et al., (2013)). They are all in the range of $0-1$, with the higher value indicating better performance.

• •

Precision at $n$ (P@$n$)

This index presents the precision up to the $n^{th}$ retrieved image.

• •

Mean average precision (MAP)

This index accounts for the case when multiple matching images are in existence in the dataset. For each query image, the database is retrieved till the last matching image is identified. For a certain matching image of which the rank is $n$, the fraction of total matching images among the first $n$ retrieved images is used as the associated precision. The average retrieval precision for a query is calculated using the precision associated with each of its matching images. Then, the MAP is obtained by averaging across all query images.

• •

Retrieval accuracy (RA)

This is the overall retrieval accuracy.

• •

Area under curve (AUC)

This index calculates the area under the curve of the receiver operating characteristic (ROC), which is typical for binary detection problem considering both the rate of true detection and the rate of false positive at the same time.

3.4 The benchmark: SeiSIM

SeiSIM is a metric that was proposed to evaluate similarity between two migrated seismic images Long et al., (2015). It combines frequency domain texture attributes with space domain geological attributes. In the frequency domain, texture similarity is evaluated using statistics calculated on subbands decomposed from SP, an approach denoted as STSIM-1 Zujovic et al., (2013).

[TABLE]

where $l(\textbf{x},\textbf{y})$ is the term accounting for luminance similarity, $c(\textbf{x},\textbf{y})$ represents contrast similarity, $a_{h}(\textbf{x},\textbf{y})$ is the structural similarity based on horizontal autocorrelation, and $a_{v}(\textbf{x},\textbf{y})$ gives the structural similarity according to the vertical autocorrelation, all calculated between subband images x and y.

In the space domain, geological similarity, denoted as DM Long et al., (2015), is assessed using statistics obtained from discontinuity maps associated with the original images. The similarity between two discontinuity maps $DM_{1}$ and $DM_{2}$ is then defined as

[TABLE]

where $a_{DM}^{h}(DM_{1},DM_{2})$ finds the similarity based on horizontal autocorrelation calculated in each map, and $a_{DM}^{v}(DM_{1},DM_{2})$ determines the similarity using vertical autocorrelation from each map.

Consequently, the seismic similarity is expressed as

[TABLE]

When applied to the retrieval problem in this paper, SeiSIM does not follow the framework used for the other texture attributes discussed earlier. It extracts simple statistics such as means, variances, and autocorrelations from the attributes and directly compares them between images. In contrast, the retrieval framework in this paper for all other techniques forms histograms of the attributes and compares them instead. Obviously, SeiSIM is not exactly comparable to the others. Therefore, we use it in this study as a benchmark, since it has been demonstrated to be a good measure of the seismic similarity Long et al., (2015).

3.5 Results

The retrieval results are shown in Table 1. The overall best performance is observed with SeiSIM, for which we believe the incorporation of the spatial geological attributes (discontinuity) plays an important role. As we pointed out earlier, SeiSIM is not exactly a comparable method within the framework adopted in this paper, but used as a benchmark. Comparing with SeiSIM, generally the retrieval performance was excellent for the four attributes being examined. P@20 values are all perfect. P@50 values are close to each other among the four, and still close to SeiSIM. The same observation is made with AUC. For MAP, CT and LRI yield higher values than SP and LBP. While for RA, CT is the best, followed by LRI. RA results for SP and LBP are very close and much lower than CT. Among the four attributes, CT yields the best overall performance, closest to SeiSIM. LRI is the second best when all metrics are considered. We conclude that the attributes all demonstrate great potential of effectively capturing the characteristics of the different geological structures, with CT and LRI being the most promising ones.

We also examined the computational time required for comparing a pair of images using the techniques. The comparison was performed on a computer with the following configuration: Intel Core i7, 3.4 GHz, with 32GB of RAM, and running on 64-bit windows 7. The results are shown in Table 2. CT is the fastest among all techniques.

4 Conclusions

In this paper, we conducted a comparative study of several typical examples of texture attributes developed in the image processing community in recent years, covering different techniques in frequency and space domain. Within a framework for image retrieval, the attributes were examined in terms of retrieval accuracy. The study demonstrated that texture attributes are generally capable of characterizing a migrated seismic volume according to its geological structure, thus can be effective for computer-assisted understanding of subsurface structures. We hope this study will also introduce to the exploration geophysicists the existing powerful image texture analysis methods, which have the potential to provide useful attributes for various seismic exploration applications.

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).