Structural smoothing on mixed instantaneous phase energy for automatic fault and horizon picking: case study on F3 North Sea

Abstract

Today, with the development of 3-D studies and the increase in seismic data volume, there is a growing need to expand interpretation techniques for achieving higher speed and accuracy of interpretation tasks. Determining seismic faults and horizons is vital to accomplish the process as one of the essential stages of data interpretation. With the recent development of computational methods in seismic interpretation and their benefits, different approaches have been promoted. The specialist can make the understanding much faster with higher accuracy. In this research, a fully automated dual horizon and fault selection approach in the presence of semi-vertical faults is presented using a structural smoothing condition. Geological faults make it challenging to map sedimentary layers appropriately which is targeted for review in this work. Unlike Image processing techniques that determine the location of faults only, the proposed approach gives the benefit of the estimated fault displacement. In this method, faults are modeled as a displacement vector field. Despite traditional methods (such as similarity and coherence), in this method, the vector field of the estimated fault displacement determines the displacement and its location. This vector field can be used for auto-determination of fault-related layers displacement. As a result, automatic horizon picking in the presence of such faults is possible, thereby simplifying the mapping of sedimentary layers.

Introduction

Geological modeling, seismic properties, and petrophysical modeling using logs, interval velocities, cores, and 3D seismic data are utilized to create structural elements on seismic data for reservoir modeling (Anees et al. 2022a, b). In the case of automatic horizon and fault picking, knowing about seismic objects is essential. For example, fluvial paleo-channels can be detected easily by logs data, core, and outcrop. Still, due to the outreach of several faults and limited seismic separability in the study region, it is difficult to recognize them by conventional seismic methods(Anees et al. 2019).

Due to the character of seismic data, fault identification is a complicated process, and using seismic data to determine the exact location of the fault is rarely possible (Taheri and Hadadi 2020). Seismic image processing usually involves locating faults (Hale 2009). Faults play an essential role in the zoning of the production reservoir so that the results are visible in three-dimensional modeling (Taheri and Morshedy 2017). To determine the horizons in the process of further interpretation and processing, it is necessary to capture the travel time of this five-dimensional data on an n-dimensional subsurface (n <  = 5) (Bienati et al. 1999). Neural network cells have the property of local communication. This property is used to connect the horizons. From the experimental results of the seismic bright spot pattern, the compatibility of the determined horizons with the reality in the image has been proven (Huang et al. 2005). Another method to determine seismic horizons automatically is using pattern recognition, which eliminates the interference of peak points on the same horizon (Tsiotsios and Petrou 2013). The minimum spanning tree (MST) is an algorithm used to guide and optimize the selection of trays to minimize interference points and completely outcoming the horizon. Another method of automatically determining seismic horizons is a minimum spanning tree (MST) algorithm, which calculates the sum of edge values by searching between image edges and finding the spanning tree with a low value. The automatic determinant of seismic horizons in this method, by using MST algorithm, ensures the optimal horizon with the minimum discount and consequently the maximum score. After creating directional vectors, the MST algorithm is applied to search horizontal peak points in three dimensions (Yu et al. 2011). In a study by Faraklioti and Petro, a one-dimensional linear filter in the vertical direction collected candidate horizon points (Faraklioti and Petrou 2004). At the same time, its previous methods have used combinations of filtering and two-dimensional horizon tracking. A fault displacement area is estimated if the pattern analysis is omitted in the image cropping step to virtualize the faults (Kadlec et al. 2008). Most proposed seismic image processing techniques trace the fault location but cannot estimate the resulting displacement (Liang et al. 2010). In an ideal two-dimensional seismic image, faults appear as discontinuities without any random or coherent noise. Commentators can quickly identify these discontinuities (Taheri et al. 2018). The fault zone parameter is usually stands for faults that include fractures (Taheri and Mohammad Torab 2017). This study's fault is a single discontinuity.

Wu and Fomel (2018) performed automatic fault interpretation using the optimal level. By calculating the weighted average of the optimal surface directions, they could accurately estimate the fault orientations (slopes and impacts). In 2019, a study was conducted on constructing a 3-D geological model of a fractured reservoir. Geological and engineering data from fractured basement formation for the 3D geological model were used. The results showed that the calculated storage capacity for the fractured reservoir is limited from 7 to 99 million metric tons. These estimated results indicate that fractured reservoir has a combined potential for storage (Vo Thanh et al. 2019).

Wu et al. (2019) Conducted a study on seismic image processing to detect faults and smooth their structures. This study was performed using a neural network. Its field studies showed that this processing method shows all three stages of accurate fault detection, smoothed seismic volume with structure, and reflection slopes more accurately than conventional methods.

A study was conducted in the Hangjinqi region in China. The base target of this research is to identify the horizontal geometry of the numerous channels. Based on the study of unified analyses of geological cores, geophysical well-logs, and seismic properties, it is confirmed that multiple channels are available in the study area. Furthermore, channel dispensation based on seismic properties and its overlapping on the paleogeographic map showed that they are coming from the North and are merging in the South direction (Anees et al. 2019).

Ashraf et al. (2020) Conducted a seismic study to identify the fault network and fractures. This study was performed using machine learning techniques. The density, slope, and curvature of fractures were calculated through the neural network. This study showed that unsupervised vector quantization (UVQ) shows the field's direction and spread of fractures. In addition, ant colony optimization (ACO) shows the faults' fracture length, slope angle, azimuth, and surface area.

An et al. (2021) used neural networks to detect three-dimensional seismic data faults automatically. This study showed that this method is comparable to human performance in analyzing image data and has even identified subtle faults. It also reduces fault interpretation time and ensures results. Another study purposed to delineate subsurface geobodies by extending a facies model to study the sedimentation processes and facies dispensation that have been neglected formerly. This study aimed to delineate subsurface geobodies by developing a facies model to study the depositional processes and facies distributions that have been ignored previously. The result of this study showed that the used methodology for the produced facies model could be extended to various basins with equal geological settings, and it can be used for prospect evaluation, future drilling, and development plans within the Sawan gas field in Pakistan (Ashraf et al. 2021).

In 2022 a study was conducted on sweet spots prediction through fracture genesis. This research engaged structural smoothing and data analysis using directional filtering to sharpen the structural discontinuities as well as 3D visualization and automatic tracking were conducted to expound on the faults and horizons (Jiang et al. 2022). In 2022, a study was conducted in the Northern Basin of China to identify gas accumulation in appropriate areas using faults and sedimentary facies. This study is conducted by incorporating the 3D seismic grid, well logs, and several cores using seismic stratigraphy, geological modeling, seismic attribute analysis, and well logging to identify gas accumulation zones. The unified results showed that the North-Western sector was uplifted compared to the southern sector. The natural gas accumulated in the southern region was migrated through fault into the northern zone and showed that the favorable zones of gas accumulation lie toward the northern region (Anees et al. 2022a, b).

This study sought to replace the automated process using the structural smoothing method with the manual method to show the displacement faults to determine the seismic horizons in the presence of a fault. In this process, fault curves are extracted from two-dimensional images of a mixed seismic attribute (mapped instantaneous phase energy abbreviated as MIPE). The displacements resulting from faults around these curves are obtained. Theoretically, the automatic determination of seismic horizons was thoroughly investigated using the structural smoothing method. It was tested in computer programs using various artificial and accurate data. This method fully includes ideal conditions, assuming no random and non-random noise, turbulent zone conditions, convergent horizons in the presence of a strong reflector, and the use of other seismic indicators. Matlab software is used to develop computer codes.

In continuation, the methods and applications are presented. Then, the mathematics of the method, which includes fault modeling, filtering, shearing, and determination of structural tensors, was discussed. In the following, the automatic determination of seismic horizons in the presence of semi-vertical faults was discussed using the structural smoothing method in the MIPE section. The necessary processes for determining the position of the fault and then selecting the seismic horizon before and after its location were discussed in detail. Finally, the horizon and faults are imported to Opendtect.

Geological setting

The sedimentary basin of the North Sea is under the mechanism of rifting in most of the Mesozoic periods with a subsidence phase after the Cenozoic rift sag. As well as rifting began in the Triassic period and reached its maximum in the Jurassic and Early Cretaceous. In this sedimentary basin, the main source rocks for hydrocarbons are Westphalian coal layers and Lower Jurassic Posidonia shales, which contain gas and oil, respectively. The deposits related to the gas, pertain to the clastic sedimentary successions deposited after the Mid-Miocene. The ice build-up had a significant impact on the faults and tectonics movement of salts in this sedimentary basin and caused it to disrupt the sedimentation channels and create routes for hydrocarbon. Obstruction by the ice build-up in the North Sea region caused the diversion of rivers that used to flow from the north to the west.

Block F3 is in the Dutch part of the North Sea and is covered by 3-D seismic operations to explore gas and oil in the layers above the Jurassic and below the Cretaceous. In this block, large-scale sedimentary layers resulting from the burial of deltaic systems in the Baltic Sea region are visible. This delta package contains sand and shale with porosity above 20–33%. Also, cemented carbonate veins are available. This package has exciting features, including large-scale book layers. Also, in this block, Bright Spots are the result of biogenic gas. In the seismic section, some transparent, jagged, linear, and, hollow outcrops are recognizable. Well logs show that evident outcrops have a single lithology and can be sand or shale. Irregular and cluttered outcrops represent faults and shrouds containing sandy turgidities. The data is accessible in the public domain, and many publications use this dataset.

Methods and applications

Structural smoothing method

Structural smoothing coefficients are extracted from the seismic/attribute input image due to the related structures. Smoothing along structures apparent in seismic images can enhance these structural features while preserving significant discontinuities such as faults or channels. (Ashraf et al. 2019). The selected directions are obtained using the scan estimates (over a group of input directions) and choosing the most similar ones or among those whose comparison markers are maximal. Directional scans can provide the similarities and hints needed for structural smoothing. Searching for the optimal direction before scanning is costly if the information is unavailable on the linearity or plane of the imaged structures. Buried channels are the properties of curved lines in seismic images. In comparison, geological horizons appear as curved planes. Both structural smoothing and similarity methods are expected to consider the different dimensions of apparent structures in seismic images (Jiang et al. 2022). Therefore, directional and dimensional scanning is computationally expensive for 3D images.

Coherency is a guide for geological layers in seismic images of the Earth's subsurface (Fig. 1a). Sparse discontinuities are usually seen, that is the idea of picking faults (Fig. 1b) shows two apparent discontinuities in the center image. This image shows the same section showing the vertical displacements caused by the fault. Displacement signs indicate that the geological layers have shifted from left to right or down (positive) or up (negative). These types of discontinuities are consistent with geological faults. The fault is more perpendicular to the layers, so the rocks on one side of the fault tend to move down or up compared to the rocks on the other side.

Fig. 1

a A two-dimensional seismic section, b extracted faults on section(Liang et al.)

In this paper, the fault is considered as a single image roughness. Seismic image processing usually involves locating faults (Salom et al. 2009). However, this processing cannot yet be a reliable method for estimating displacement along faults. Today, faults' location is done manually, using the complex process of viewing seismic images, and determining the corresponding points on both sides of the fault. This study aims to replace this manual method with an automated process that shows fault displacements. This process extracts the fault curves from the two-dimensional images and obtains the fault's displacement around these curves.

Fault modeling

A common feature that most fault detection techniques have is calculating markers that identify seismic image discontinuities. The structural smoothing method, calculated by the symbol, shows the fault zone. This method is quite relatively in identifying fault curves and fault positions. Geologists and geophysicists performed mathematical modeling of faults based on seismic images and well survey data (Barnett et al. 1987). In this study, the focus is on two-dimensional sections of seismic images. Therefore, using the two-dimensional scale matrix $f({x}_{1},{x}_{2})$ where ${x}_{1}$ and ${x}_{2}$ are the vertical and horizontal coordinates of the sampled points, and the faults are as discontinuities in the seismic image. The fault displacement estimation is a two-dimensional vector which is: $u=u({u}_{1},{u}_{2})$ (see Fig. 2. a) where ${u}_{1}={u}_{1}({x}_{1},{x}_{2})$ and$,{u}_{2}={u}_{2}({x}_{1},{x}_{2})$, respectively. Vertical and horizontal displacement components are the result of faults. The vertical displacement of the sequence to the sequence usually corresponds to the direction of the geological layers in the seismic images. Still, it does not reasonably describe fault displacement problems such as insufficient use of two sequences, non-vertical faults, and lack of faults in sloping layers. These are the three weak points factors that will improve in the following sections.

Fig. 2

a Fault displacement vector u(x) (Liang et al.), b Cut a sloping layer(Liang et al.)

The Fault position assessment is by two-stage pattern analysis. The first step in the study is to estimate the fault position. The next step is to remove the fault from the fault zone with one-pixel precision. In seismic images, the vertical displacements between sequences, before and after cutting, are zero in horizontal geological layers. Therefore, the cut can change the removal between the sequences when the layers are sloping. As shown in Fig. 2b, $f$ and $f^{\prime}$ are the pre-and post-cut images.

A sloping geological layer is the modeling of a linear cross-section at a specific point. Thus,$f$ becomes a function of a variable $g$ as follows:

$$ f\left( {x_{1} ,x_{2} } \right) = g\left( {x_{1} - \rho x_{2} } \right) $$

(1)

where $\rho$ Is the slope value.

Therefore:

$$ \begin{aligned} f^{\prime}\left( {x_{1} ,x_{2} } \right) & = g\left[ {x_{1} - \rho \left( {x_{2} - sx_{1} } \right)} \right] \\ & = g\left[ {\left( {1 + s\rho } \right)x_{1} - \rho x_{2} } \right] \\ \end{aligned} $$

(2)

Similarly, we will have:

$$ \rho^{\prime} = - \frac{{\frac{\partial f^{\prime}}{{\partial x_{2} }}}}{{\frac{\partial f^{\prime}}{{\partial x_{1} }}}} $$

(3)

where $\rho^{\prime}$ Is the slope of the layer after cutting and substituting $f^{\prime}$ in Eqs. (2) using (3) we have:

$$ \rho^{^{\prime}} = - \frac{{\frac{{\partial g\left[ {\left( {1 + s\rho } \right)x_{1} - \rho x_{2} } \right]}}{{\partial x_{2} }}}}{{\frac{{\partial g\left[ {\left( {1 + s\rho } \right)x_{1} - \rho x_{2} } \right]}}{{\partial x_{1} }}}} = - \frac{{ - \rho \frac{\partial g}{{\partial x_{2} }}}}{{\left( {1 + s\rho } \right)\frac{\partial g}{{\partial x_{1} }}}} = \frac{\rho }{1 + s\rho } $$

(4)

where ρ' is a function of the shear value of s, the slopes $\rho$ and $\rho^{\prime}$ are vertical displacements, precisely at the selected point, before and after cutting. When $\rho^{\prime}\left( s \right) = \rho = 0$ and $\rho = 0$, which means that the cut will not affect the horizontal layer The innovation in our methodology is using a seismic attribute that is the map of the instantaneous phase in the energy section. Mixed instantaneous phase energy (MIPE) is simply a multiplication of normalized instantaneous phase and energy attributes that highlights both the significance of amplitudes and phase characteristics of seismic reflected waveshape simultaneously.

Figure 3 shows the different functional relationship patterns where the apex of the black curve in Fig. 3b corresponds to the black pixels in Fig. 3a, and this pixel is located precisely on the fault. Figure 3b shows the horizontal coordinates of the peak point $l_{2}$, which cuts the image vertically (with a ratio of l $\frac{{l_{2} }}{{l_{1} }}$). The curves in Fig. 3c correspond to the relationships shown in the equation.

Fig. 3

Investigation of the relationship between vertical displacement and ${l}_{2}$: a seismic artificial image, b five corresponding curves with five selected points in the seismic artificial image, c shape diagrams for different sections (with lower values) (Liang et al.)

These vectors should be zero or close to zero in areas with no faults. Near-horizontal faults are not easily discernible even by commentators. Therefore, the faults considered in this study are more vertical, and it is assumed that the angle between the fault curve and the vertical line is less than 45°

Filtering

Many different studies describe filtering techniques for seismic images to help interpreters. Many of these techniques add features such as structural layers when reducing noise. One of the techniques for filtering seismic images is the use of the Van Gogh filter, which increases the dispersion coherence in seismic images. This method solves partial differential equations with the help of a field scattering sensor. Another seismic image filter is the edge protection filter, which extends to three dimensions (Liu et al. 2022). The most important advantage of this filter is its efficiency and simplicity. However, none of these filters provides the required correlation (Verschuur 2013). Here, designing a filter is necessary to create a cross-correlation between both sides of a fault. This filter collects information about the input of adjacent effects but cannot collect them along the fault. The designed filter has two input effects that produce a weighted average impact on each side of the fault, and its preparation is as follows.

It is smooth and non-edged from left to right and from right to left, and these filters are created by inverting the direction. This method is a powerful and flexible one-way filter. The input of this filter is a two-dimensional image that acts as a one-dimensional set of vertical effects. This image (sequentially) is processed from left to right and from right to left. Changes in the one-way power filter resulted in more fault retention during smoothing. This filter is suitable due to its low cost and little smoothing in the presence of discontinuities. In the tests performed, the scan of the displacement fields is the selection of a sequence of the smoothing image (from left to right), and the connection of the created differences is the unevenness of the fault.

Shearing

The shear wave reflectivity of the fault zone is taken from one-dimensional normal occurrence modeling of accurate geologic fault zone models. The shear wave reflectivity of the fault zone, which originates from fine compositional layering, is moderate and can be directly associated with singular geologic formations. Within the fault zone, variable shear wave reflectivity, both horizontally and vertically, originates from geologic dissimilarity(McCaffree and Christensen 1993). Seismic anisotropy, especially shear wave splitting, provides strong evidence for coherent deformation over several tens of km wide in the lithospheric mantle under major transcurrent faults. Yet it cannot find narrow strain localization zones or shallowly dipping faults (Vauchez et al. 2012).

When the faults are not vertical, non-vertical sections connect the pixels. Here $l_{1 }$ is defined as the maximum vertical displacement of the fault, and the value of the variable $s$ is defined as the cut-off image value that controls the number of displaced samples. Using a mathematical approach, cutting a two-dimensional vector field $f\left( {x_{1} ,x_{2} } \right)$ means creating a new two-dimensional vector field $f^{\prime}\left( {x_{1} ,x_{2} } \right)$ as follows:

$$ f^{\prime}\left( {x_{1} ,x_{2} } \right) = f\left( {x_{1} ,x_{2} - sx_{1} } \right) $$

(5)

where $s > 0$ indicates the rows' movement to the right and $s < 0$ shows the movement of the rows to the left. Assume the value of $s \in \left( { - 1, + 1} \right)$ because, in this model, the slope of the fault must be less than 45 degrees relative to the vertical line. To sample the cut-off value range $\left( { - 1, + 1} \right)$, we define another value of $l_{2}$ with the assumption $l_{2} \in \left( { - l_{1} ,l_{2} } \right)$ and $s = \frac{{l_{2} }}{{l_{1} }}$. Since s is not an integer, the sequence is a floating value using interpolation. Shear displacement also involves the protrusion and cutting of a portion of the seismic image, which must eject pre-cut information for the lost image so that post-cutting information is not lost.

Fault coefficient (c)

According to Fig. 3, the difference between the fault location pattern and the fault-free location on the MIPE section is easily visible. The result is the selection of the value of C over an extensive range in this image with artificial data. In actual data, however, the fault pattern is less detectable. In this case, the scope for determining the value of c is minimal, resulting in incorrect fault detections, which is due to the selection of small values of this coefficient. When a value increases, these incorrect choices disappear. Therefore, since seismic images usually have random noise and slopes change in the direction of the geological layers, it is necessary to choose a more appropriate value or a better analytical model for processing seismic images with lower quality.

In the first step, samples around the fault are detected from distant samples to the fault according to the following criteria: If the displacements change slowly according to different shear values, there is no fault, in other words, if a significant point in the displacement-shear curve, as in the black curve shown in Fig. 3, there is a fault. In practice, η is the fault threshold. The calculation of the mean value together with the maximum number of vertical displacements with different sections at a given point is as follows:

$$ {\text{Max}}\left( {\rho^{\prime}\left( s \right)} \right) > \eta $$

(6)

and

$$ {\text{Max}}\left( {\rho^{\prime}\left( s \right)} \right) > c\left| {{\text{ave}}\left( {\rho^{\prime}\left( s \right)} \right)} \right| $$

(7)

The use of genetic algorithms and multi-scale models of Bissin are new methods to solve this problem (Admasu and Tönnies 2006). After obtaining the vertical displacement V₁ and the sheer value S, the calculation of the displacement vector is as follows:

$$ u_{1} \left( {x_{1} ,x_{2} } \right) = v_{1} \left( {x_{1} ,x_{2} - sx_{1} } \right) $$

(8)

$$ u_{2} \left( {x_{1} ,x_{2} } \right) = - su_{1} \left( {x_{1} ,x_{2} } \right) $$

(9)

The continuity of the fault curve depends on the maximum amount of vertical displacement $l_{2}$. $l_{2}$ Controls the vertical size of the window and, therefore, cannot be too small. Larger values of $l_{2}$ cause the fault curves to be continuous, while the use of smaller values of $l_{2}$ with a sudden change of direction shows more minor faults. This is because larger values of $l_{2}$ help to detect more significant faults. Figure 4 shows these relationships in actual data where the blue curve is still distinguishable from other angles.

Fig. 4

Investigation of the relationship function between vertical displacement and ${l}_{2}$: a five selected points of the image with artificial data, b the five curves equivalent to the five selected points with the same colors (Liang et al.)

Selecting different values of the η fault threshold creates different levels of detail for the displacement field. The logical value of the fault threshold corresponds to the sampling intervals of the image. The interpreter can get the best result by constantly changing the threshold number value. The workflow of the methodology is shown in Fig. 5.

Fig. 5

Workflow of the methodology of structural smoothing on MIPE attribute yields two horizon and fault cubes

Results and discussion

Combining the structural smoothing method on the MIPE section with seismic markers to determine the horizon and other geological and stratigraphic features is done. Determination of seismic horizon and geological and stratigraphic characteristics of the F3 block located in the North Sea, using Opendtect software and structural smoothing program, is targeted. Data were analyzed using seismic software indicators. To Identify the functionality of the proposed workflow inline No. 250 is selected (Fig. 6). The abnormal high values are even evident clearly on mixed instantenous phase energy (MIPE) section (Fig. 7).

Fig. 6

Seismic section inline direction (No 250) of F3 Dutch Sector North Sea data

Fig. 7

Mixed instantaneous phase energy (MIPE) attribute for seismic inline number 250. Bright spots and abrupt changes are modified due to the mapping of normalized energy. The coherency is enhanced as a result of instantaneous phase attribute functionality

An energy array identified bright spots that are suitable for finding seismic targets. The energy indicator detects bright spots in the upper left of the section.

MIPE improved the balance of amplitude spectrum (the observable weighting in Fig. 8). The contribution of both the high and low frequencies of the seismic frequency band by preparing and performing this indicator is visible. The amplitude of the MIPE section showed the first peak as the minimum frequency and the last peak as the maximum frequency.

Fig. 8

Amplitude balanced spectrum of MIPE attribute and its frequency distribution diagram in terms of frequency

As a result, the blow frequencies are reported for the MIPE section,

Min freq. = 8 HZ,

Max freq. = 75 HZ,

Mid freq. = 41 HZ.

The spectral analysis index defined these three frequencies. By applying this indicator and assigning blue to the high-frequency indicator, green to the medium frequency indicator, and red to the low-frequency indicators, and also flipping the colors, the following figure was obtained in which you can easily preserve channels in a sedimentary environment.

The image with this condition (Fig. 9) contains an additional similarity marker in which, in addition to the channel, there is also a fracture. RGB refers to a system for representing the colors to be used on a computer display. Red, green, and blue can be combined in various proportions to obtain any color in the visible spectrum. Levels of R, G, and B can each range from 0 to 100% of full intensity. Based on the MIPE section, RGB frequency mixing methods are applied to break down the frequency spectrum and map the dispensation of the Channel system and micro-fractures in the rock. Compared with the conventional RGB color mixing manner, the RGB color mixing technique based on the relict impedance inversion not only can better delineate the relationship between the large-scale fractures but also can effectively predict the distribution of those small-scaled fractures in the region (Ruizhao et al., 2017). This information will later be considered as directions for structural smoothing optimization in our algorithm.

Fig. 9

Channel system and micro-fractures using RGBA colors and similarity markers

The fault model for seismic inline number 250 is shown in Fig. 10. These pink lines are automatically extracted faults. Seismic horizons are also shown for a part of this inline in Fig. 11 in a seismic line. The standard tracking horizon in Opendtect(pink lines in Fig. 11) are located in a different view from our method (yellow lines in Fig. 11). Among three options (positive, negative, or zero-crossing), the MAX option is chosen. Hence the peak of positive reflectors is as reported in this figure for both the standard Opendtect tracking algorithm and our methodology, i.e., structural smoothing on the MIPE attribute. Although the yellow and pink horizons in Fig. 11 approve each other in most places, in the vicinity of the faulted area, the yellow one has a higher correlation with the reflector character.

Fig. 10

Determining the position of major and minor faults. The faults are tracked in both vertical and semi-vertical directions

Fig. 11

Comparing the results of running two methodologies, Opendtect standard horizon tracking (pink) and our proposed strategy (yellow), in the presence of a semi-vertical auto-picked fault

Summary and conclusions

• 1. The results showed that structural smoothing filters on MIPE (modified energy and phase in one section) are easy and efficient ways to simultaneously delineate semi-vertical and seismic horizons. The parameters are derived from structural tensor fields from a modified seismic attribute. Unlike conventional smoothing filters, whose parameters are determined based on the slope of the image features, this filter allows the interpreter to smooth 2D and 3D images in the desired directions and dimensions.

• 2. Also, the computed structural-directional similarity is defined as the ratio of squared image squares to image squares smoothing, which is a particular example of weighted similarity (in the presence of seismic objects like fractures and channels in this study). The direction of these features is arbitrary and has fewer vertical restrictions, and unlike the weights used in conventional similarity markers, the weight used in the structural-directional similarity markers has a slow and continuous attenuation to zero.

• 3. The estimated displacement field provides pixel accuracy in horizontal-vertical characteristics, but this accuracy at the fault location is usually computed between two nearby pixels. This fault displacement vector field could improve the estimated tensors in the areas around the fault and thus smooth the seismic images. Besides, fault displacement vectors could use to smooth seismic images later, like a median fault enhancement filter.

• 4. This algorithm assumes the relative continuity and detectability of the horizons, which in many cases is not the case, and therefore makes this method unreliable when considering all seismic horizons. This method effectively measures and determines the trend of horizons and faults simultaneously, and its application is in the analysis of seismic images for better structural interpretation models.

Appendix A: a computer program

This study used Matlab software for the program file. The code is not planned to be public as it is part of a bigger toolbox at the University of Tehran.