The application of volume texture extraction to three-dimensional seismic data – lithofacies structures exploration within the Miocene deposits of the Carpathian Foredeep

There are numerous conventional fields of natural gas in the Carpathian Foredeep, and there is also evidence to suggest that unconventional gas accumulations may occur in this region. The different seismic signatures of these geological forms, the small scale of amplitude variation, and the large amount of data make the process of geological interpretation extremely time consuming. Moreover, the dispersed nature of information in a large block of seismic data increasingly requires automatic, self-learning cognitive processes. Recent developments with Machine Learning have added new capabilities to seismic interpretation, especially to multi-attribute seismic analysis. Each case requires a proper selection of attributes. In this paper, the Grey Level Co-occurrence Matrix method is presented and its two texture attributes: Energy and Entropy. Haralick’s two texture parameters were applied to an advanced interpretation of the interval of Miocene deposits in order to discover the subtle geological features hidden between the seismic traces. As a result, a submarine-slope channel system was delineated leading to the discovery of unknown earlier relationships between gas boreholes and the geological environment. The Miocene deposits filling the Carpathian Foredeep, due to their lithological and facies diversity, provide excellent conditions for testing and implementing Machine Learning techniques. The presented texture attributes are the desired input components for self-learning systems for seismic facies classification.


INTRODUCTION
Modern, high-resolution 3D seismic volumes make a significant contribution to the quality of subsurface geological imaging and are widely used in the oil industry for hydrocarbon exploration. However, due to the huge amount of data obtained by the method of three-dimensional seismic surveying, the effective interpretation of seismic volumes increasingly requires the automation of cognitive processes through the use of the technique of self-learning systems. Recent developments with Machine Learning (ML) have added new capabilities to seismic interpretation, especially to multi-attribute seismic analysis. Moreover, at present ML is perceived as a method that will undoubtedly affect the interpretation of seismic data in the future (Wrona et al. 2018). These processes can proceed in a supervised and unsupervised form, and the key to both cases is the correct selection of seismic attributes (Marfurt 2018).
The results of the latest research prove that instead of the full range of available seismic attributes, only the selection of those attributes which most distinguish elements from the seismic data should be used in the process of correlation with identified seismic facies (Roden et al. 2015, Infante-Paez & Marfurt 2019. The application of texture attributes for seismic facies classification using ML techniques can provide seismic facies information that is not provided by any conventional seismic approach (Marfurt & Chopra 2007). Hence, in this paper two texture attributes: Energy and Entropy will be presented. Both were derived by the use of the Grey Level Co-occurrence Matrix (GLCM) method, which emphasize (in a special way) the continuity of the identified seismic reflection patterns.
The seismic data used in the research comes from a modern 3D survey located in the NE part of the Carpathian Foredeep Basin, which belongs to the foreland basin system that surrounds the Carpathian orogenic belt (Fig. 1). The 3D seismic survey was acquired in order to develop the stratigraphic model that would match the production history as well as being consistent with the available well control.
The study area is located within the northern, outer part of the basin which is filled with thick pelagic and turbidite formations (Wysocka 2006). The main objective of the research was a fragment of sedimentation profile referred to as the Machów Formation (Myśliwiec 2004). The sedimentation process of this Late-Badenian and Sarmatian clastic succession was influenced by its close proximity to the Roztocze Hills ( Fig. 1), forming part of the Carpathian forebulge zone (Wysocka 2016). The results of the seismic data interpretation confirmed the significant tectonic activity in the Miocene basement, highlighted by the high variability of sedimentation environments during the temporal -spatial evolution of a basin fill. All of the above information makes the study area extremely interesting for testing modern tools and unconventional methods to optimally understand the variability of environments and to estimate the exploration potential of structures within the basin. There are numerous conventional fields of natural gas in the Carpathian Foredeep and there is also evidence to suggest that unconventional gas accumulations, related to non-structural types of traps, may occur in this region (Pietsch et al. 2010). Such exploration objects are, in particular, the external and internal elements of submarine channels and fans. The different seismic signatures of these geological forms and the small scale of amplitude variation and limited visibility, along with the large amount of data, make the process of geological interpretation extremely time consuming. The key stage of the interpretation workflow with the use of self-learning techniques is the selection of seismic attributes that are highly sensitive to the presence of the objects sought. They include, among others, amplitude attributes (the most common one), frequency, geometric or those calculated on the basis of seismic texture. The latter, describing the spatial relationships of seismic facies, are the subject of this research. GLCM attributes in the interpretation of seismic data, in addition to frequency analysis, have a significant impact on the understanding of the geology of the area under study, both at the early stage of geological reconnaissance across the seismic volumes and in the later, advanced detection of selected prospecting objects. These cognitive processes can take place in a traditional way, through the conventional work of the interpreter, or in an automated manner by the use of ML techniques which are supervised or not.

ATTRIBUTE ANALYSIS METHODOLOGY
It has been almost 50 years since Haralick and his co-authors published a series of articles covering ways to use the statistical Grey-Level Co-occurrence Matrix method (Haralick et al. 1973, Haralick 1979. His work was based on the results of previous mathematical research (authors quoted by Kupidura et al. 2015, Hall-Beyer 2017, however Haralick was the first to propose the use of texture attributes in the image classification process. The GLCM attributes can be divided into three groups (Hall-Beyer 2017): -contrast (contrast, homogeneity, dissimilarity), -orderliness (entropy, angular second moment, energy), -statistics (mean, variance, correlation).
The GLCM is a square matrix and its size depends on the number of grey levels Ng in the analyzed image. Neighborhood relationship analysis is performed for each pixel in the texel, a defined area (x, y) of the image (Fig. 2). These relations are usually determined for four directions (angles α: 0°, 45°, 90° and 135°) and for a specific distance d between the analyzed pixels (Haralick et al. 1973, Hall-Beyer 2017. Typically, the distance d between the reference pixel and the adjacent pixel is 1, but other values can also be used for calculations.

Fig. 2. Texel (x, y) and defined four directions of neighborhood analysis between pixels at a distance d
As already mentioned, the GLCM method was originally dedicated to texture analysis of digital two-dimensional images. Texture refers to the characteristic patterns in an image, defined by the value and variability of samples/pixels at a given location. In such an approach, the 2D seismic section can also be considered as a texture image (Fig. 3A). Then we are dealing with seismic textures formed by the set of seismic reflections, and the pixel values are determined by the size and variability of seismic samples located along the seismic trace (Gao 1999(Gao , 2002(Gao , 2003. Hence the concept of a seismic texel (Fig. 3B), within which the texture parameters (attributes) are calculated for each seismic sample by analyzing the arrangement of adjacent reflections ( Fig. 3C-F). The values of the GLCM matrix are defined by the frequency of coexistence of a reference pixel of a certain grey level i with an adjacent pixel of a different grey level j. The calculation process is schematically represented in Figure 3C, D by means of green and orange rectangles and squares. There is one coexistence of pixels with values 2 and 1 (green color) and two pixel neighborhoods with the values 3 and 4 (orange color). The GLCM matrix is the basis for further statistical analyzes, but before the texture parameters are determined, it should be normalized (Hall-Beyer 2017) (Fig. 3E). Normalization gives the matrix a probability distribution, and it is obtained by dividing the value of each of its elements (i, j) by the sum of all elements (Badurska, 2007).
At this point, it is worth emphasizing the increased informativeness of the presented vertical seismic section with the textural attribute ( Fig. 3F), where it is much easier to see channel-like objects than in the classic amplitude image (Fig. 3A).
In the literature of the subject, several previous articles can be found on the use of the GLCM method for the interpretation and visualization of 2D seismic data (Zhang & Simaan 1989, Vinther et al. 1996, Gao 1999, 2002. In later years, attempts were made to modify the GLCM method in order to determine the possibility of describing seismic structures within 3D seismic volumes. Two interesting papers presenting the results of research on the characteristics of seismic textures based on 3D volumes were published by Gao (2003) and Eichkitz et al. (2013). In this article, the assumptions of the Gao method will be presented, as his work, along with Haralick's, formed the basis for the development of an interpretative tool that was used by the author to study the GLCM attributes: Energy and Entropy.
The geometry of the 3D seismic volume is based on the perpendicular system of inlines (IL) and xlines (XL), in a defined vertical space (Z) (Fig. 4A). A similar geometry has a three-dimensional texel, a mini-cube composed of N IL × N XL × N Z volume picture elements (voxels) (Fig. 4B). The size of the mini-cube is usually adapted to the characteristics of the search objects, however, the optimal number of voxels in the N IL and N XL domains ranges from 3 to 9, while in the N Z domain ranges from 7 to 21 (Gao 2003). In some cases, the size of one of the 3D texel domains may be neglected comparing to the other two. We use this possibility, for example, to capture the horizontal variability of seismic amplitudes in the analyzed interval of the seismic volume, then N IL = N XL >> N Z .
To perform statistical analysis based on 3D texel, equivalent to the GLCM method used in the texture analysis of 2D images, Gao proposed the use of the Voxel Co-occurrence Matrix (VCM) (Gao 1999(Gao , 2003. The VCM is a statistical representation of the set of seismic reflections in a 3D texel, presented in tabular form. Following the Gao's concept, the VCM matrix is calculated per texel for each analyzed sample in the seismic volume. A square, symmetric matrix is constructed of N g × N g elements, where N g is the number of grey levels in the analyzed 3D seismic data. In the VCM matrix, each element E (i, j, α, β) determines the frequency with which, in a given 3D texel, the voxel of amplitude i (<N g ) is neighbored to the voxel of amplitude j (<N g ) in the direction specified by the angles α and β. The layered pattern of reflections in the seismic image makes the VCM different for each of the three perpendicular directions, determined by the domains: N IL (α = 0°), β = 0°), N XL (α = 90°, β = 0°) and N Z (β = 90°). For example, the VCM matrix for the domain N IL (α = 0°, β = 0°) takes the following form of a mathematical expression (Reed & Hussong 1989, modified by Gao 2003: where (x, y, z) represents the size of the 3D texel, S expresses the frequency of the co-occurrence relationship between voxels specified in curly braces, g(m, n, o) and g(p, q, r) correspond to the values of two voxels with the location (m, n, o) and (p, q, r), respectively.
In the presented research, the Petrel E&P software platform was used, which is applied to generate a number of seismic attributes (volume attributes) what also provides the GLCM method. The platform is dedicated here to the analysis of 3D seismic volumes and is used in the detection of such geological elements as mass transport deposits, channels, dewatering structures, salt forms and other structures separating from the seismic background. In the calculation of GLCM attributes, the application uses the VCM method described above, therefore the elements of the matrix are counts of the frequency which voxels with specific values (corresponding to seismic amplitudes) co-occur in the analyzed 3D texel. The variation of grey levels for each sample (voxel) in a 3D volume is calculated for three basic, perpendicular directions, corresponding to the geometry of inlines, xlines and vertically. In order to use the GLCM method, the interpreter defines four parameters necessary to build the matrix: the number of grey levels of the input volume (quantization level), the size of the 3D texel (here a moving window), the distance between the reference and the adjacent voxel and the type of texture attribute to perform statistics. There are two attributes to be chosen: Energy, which is used to determine texture homogeneity, and Entropy, which measures the amount of disorder between samples/voxels. As a result, for each of the calculated texture attributes, three seismic volumes are generated (one in each considered direction along the inline (IL), xline (XL) and vertically (Z); GLCM_IL, GLCM_XL and GLCM_Z) (Fig. 4C, D, E). The first two volumes are characterized by lateral variation in the amplitude values, from seismic trace to trace (intertrace), while the third volume relates to vertical changes within the seismic trace (intratrace).
In the last stage of the GLCM analysis, the obtained attribute volumes are subject to the interpretation process. It can run both with the use of classic vertical and horizontal sections, especially when the purpose of the texture attribute is to emphasize the continuity of seismic reflections or to emphasize structural objects from the seismic image. However, the most effective analysis of GLCM attribute volumes assumes the use of interactive interpretation tools, like RGB blending within a box probe (Fig. 5). Then, apart from interesting texture objects highlighted from the seismic background, the variability of textures within the analyzed element may reflect, for example, changes in lithology or thickness.

THE RESULTS OF ATTRIBUTE ANALYSIS
The analysis was conducted on one of the most interesting Miocene intervals, with a time thickness of 200 ms. Its upper limit was determined by the reference seismic horizon H, close to the gas-bearing level identified in several production wells. GLCM analyses were performed on the volume of 3D seismic data in the processing version with preserved real amplitude relations. As the GLCM method used for the analysis of seismic data uses the seismic amplitude as a sample, this volume choice should ensure appropriate geometric imaging in the results obtained for objects with similar acoustic parameters. During its geological history, the studied area underwent a continuous structural evolution. As a result of tectonic activity and compaction, the initial structural system of the Miocene deposits was distorted. Therefore, horizon flattening along the selected horizon H was applied to the volume of 3D seismic data, during the data preconditioning stage. Thanks to this transformation, in the course of further cognitive processes (using a box probe or RGB blending), genetically related structural elements from different sedimentary environments could be tracked.
Three seismic volumes were generated for each of the directions analyzed by the GLCM method, i.e. two horizontal components and one vertical component. During careful skimming of the composites of the three obtained volumes, for both the Energy and Entropy attributes, several interesting texture objects were noticed (Figs. 6, 7). The presented images come from the upper part of the analyzed interval, located 10 ms below the seismic horizon H (Fig. 8). The similarity of the results of using both attributes, Energy and Entropy, indicates a high degree of their correlation, as they belong to the same group of ordering attributes (Hall--Beyer 2017). However, in both cases, there is a perfect definition of a few anomalous textural features that stand out from the seismic background. Several exploration wells are located in their vicinity, four of which are important due to the existing production tests in the analyzed interval below the seismic horizon H. Well O-6 is negative, no gas flow was recorded as a result of the tests. On the other hand, O-7, O-12 and O-3 are gas production wells with very good reservoir parameters. The texture attributes indicate areas which calibrate well with above exploration wells. In the remaining boreholes, those located in the immediate vicinity of the obtained texture objects, no production tests were carried out for the analyzed interval. Fig. 7. A blended image using RGB color model. GLCM Entropy. Box probe (15 ms high) located 10 ms below the reference seismic horizon H. An arbitrary section I-I' location (see Fig. 8)   Fig. 6. A blended image using RGB color model. GLCM Energy. Box probe (15 ms high) located 10 ms below the reference seismic horizon H. An arbitrary section I-I' location (see Fig. 8)

DISCUSSION AND CONCLUSIONS
The legacy 2D seismic data available so far in the study area, apart from the lack of coherence of the seismic image, has at most provided indications of the presence of a potential gas-saturated zone in the form of an amplitude anomaly. Hypotheses of sedimentation environments based on borehole data, even with a large number of wells, are still characterized by the point character of the information and also do not allow the reconstruction of the trap geometry. This modern 3D seismic survey, supplemented by carefully chosen attributes, allowed us to create a stratigraphic model consistent with the borehole data. It replaced the ambiguous results obtained from several older 2D seismic surveys. The aim of the interpreter's work has now become the search for specific shapes and then their identification. As a result, it is possible to diagnose the geological form of a clear genesis and, consequently, to identify potential reservoirs. The texture attributes are a quantitative suite that aids the interpreter in recognizing and detecting those objects more effectively than from the seismic amplitude volume. The conducted research on the selected interval of the Sarmatian sediments brought a lot of relevant information related to this gas-bearing formation.
The results of the work can be summarized in the following conclusions: 1. Contemporary methods of 3D seismic data acquisition and processing provide a huge amount of information about subsurface geology in the form of multi-variant seismic data volumes. For their effective use, it is necessary to apply advanced interpretation tools, which, based on seismic attributes, provide the interpreter with a number of clues in the process of detecting potential prospecting objects. The use of the GLCM method texture attributes in the described case was one of the important elements of the preliminary exploration phase of the seismic data interpretation work. This mainly concerned the stage of the time-consuming analysis of seismic volumes, aimed at detecting potential unconventional reservoirs. 2. The resulting texture images were the starting point for further, now object-oriented, interpretation within a specific seismic layer, a target zone, covering only an interval of about 10-12 m thick. The details of the analyzed deposition environment are interestingly presented by the Entropy attribute of the GLCM_Z volume, calculated along the seismic trace. In Figure 9 geobodies are presented that have been interpreted within the chosen interval, using a box probe with manipulated transparency. For comparative purposes, the results of two other interpretation tools have also been posted for the same seismic layer. The first one (Fig. 10) is the popular attribute of the distribution of mean square RMS amplitude values. The second one (Fig. 11) is a composite image of iso-frequencies 38, 44 and 48 Hz using cyan, magenta and yellow blending (CMY). The texture attribute carries a lot of information about the complexity of the channel system present here. The observed texture changes (in the form of changing colors) most likely reflect the lithological variability of the studied complex. The application of volume texture extraction to three-dimensional seismic data -lithofacies structures exploration... 3. The indicated position within the seismic volume allowed for a quantitative interpretation of the discovered elements of submarine channels. Figure 12A presents the extraction of real seismic amplitudes performed along an offset horizon H', located 5 ms below the seismic horizon H. This map is presented using a dedicated color palette whose task is to show the geometry of the channels, probably crevasse splays and fans visible here. In the presented interpretation, the area of the color zone reflects potential gas-saturated zones. It is an interpreter's subjective assessment, based on information obtained from gas tests result of the analyzed level in boreholes. The information about the negative gas test result in the O-6 well was of key importance for limiting the potential area. The result of the interpretation provides an explanation of the field condition presented here -its genesis and geometry. It also indicates potential areas where the gas field may continue. This is an example of the presence of a lithofacial gas object in the Miocene formations, the geometry of which is not determined by the structural arrangement of the layers. Figure 12B shows an arbitrary seismic section I-I' (same seismic data volume) through the negative and positive wells associated with the channel system under test.

The Miocene deposits filling the Carpathian
Foredeep, due to their lithological and facies diversity, provide excellent conditions for testing and implementing self-learning systems, referred to in the literature as Machine Learning. The presented texture attributes are the desired input components for ML techniques for seismic facies classification. In the case under consideration, the Energy and Entropy attributes show a very high correlation with each other, at the level of 0.9. Therefore, to reduce the attribute space, only one of them will be used in further interpretation studies using ML. Undoubtedly, a helpful tool in the process of selecting attributes will be Principal Component Analysis.
The author wishes to thank PGNiG SA for its permission to publish this paper. He is also indebted to Paweł Pomianowski and Anna Świerczewska, as well as to the anonymous reviewers for their constructive comments which helped greatly to improve the manuscript.