RESEARCH ARTICLE Deep Learning for Characterizing Paleokarst Collapse

10.1029/2020JB019685
Features in 3-D Seismic Images
Special Section:
Machine Learning for Solid Earth Xinming Wu1 , Shangsheng Yan1 , Jie Qi2 , and Hongliu Zeng3
Observation, Modeling and
1 School
of Earth and Space Sciences, University of Science and Technology of China, Hefei, China, 2 The University of
Understanding
Oklahoma, Norman, OK, USA, 3 Bureau of Economic Geology, University of Texas at Austin, Austin, TX, USA
Key Points:
• We propose a convolutional neural
network for collapsed paleokarst Abstract Paleokarst systems are found extensively in carbonate-prone basins worldwide. They can
detection in seismic images, which form large reservoirs and provide efficient pathways for hydrocarbon migration, but they can also create
is superior to conventional methods
• We propose a general workflow to
serious engineering geohazards. The full delineation of potentially buried paleokarst systems plays an
generate synthetic training data important role for reservoir characterization, oil and gas production, and other engineering tasks. We
sets with realistic structures and propose a supervised convolutional neural network (CNN) to automatically and accurately characterize
collapsed paleokarst features
• The results of our CNN-based
paleokarst and associated collapse features from 3-D seismic images. To avoid time-consuming manual
method provide an efficient way labeling for training the CNN, we propose an efficient workflow to automatically generate numerous 3-D
to fully and quantitatively analyze training image pairs including synthetic seismic images and the corresponding label images of the
the geometric details of paleokarst
systems
collapsed paleokarst features simulated in the seismic images. With this workflow, we are able to simulate
realistic and diverse geologic structures and collapsed paleokarst features in the training images from
which the CNN can effectively learn to recognize the collapsed paleokarst features in real field seismic
Correspondence to:
X. Wu,
images. Two field examples from the Fort Worth Basin demonstrate that our CNN-based method is
[email protected] superior to conventional automatic methods in delineating paleokarst collapse features from seismic
images. From the CNN-based paleokarst characterization, we can further automatically extract 3-D
Citation:
collapsed paleokarst systems and quantitatively measure their geometric parameters. Our CNN-based
Wu, X., Yan, S., Qi, J., & Zeng, H. method is highly efficient and takes only seconds to classify collapsed paleokarst features a 3-D seismic
(2020). Deep learning for image with 320 × 1, 024 × 1, 024 samples (approximately 268 km2 ) by using one graphics processing unit.
characterizing paleokarst collapse
features in 3-D seismic images.
Journal of Geophysical Research: Solid
Earth, 125, e2020JB019685. https:// 1. Introduction
doi.org/10.1029/2020JB019685
Karst is a typical type of carbonate terrain that has gone through significant diagenetic processes, which
often generate associated structural features of joints, caves, faults, and collapses (Qi et al., 2014). Buried,
Received 28 FEB 2020 collapsed paleokarst systems can form carbonate petroleum reservoirs including those explored in the Fort
Accepted 25 AUG 2020
Worth Basin (FWB) (e.g., Dou et al., 2011; Kerans, 1988; Loucks & Anderson, 1985), the Lower Cretaceous
Accepted article online 9 SEP 2020
Golden Lane field in eastern Mexico (Coogan et al., 1972; Viniegra & Castillo-Tejero, 1970), the Tarim Basin
in China (e.g., Desheng et al., 1996; Maoshan et al., 2011; Qi et al., 2014; Zeng et al., 2011), and other
paleokarst-related fields reported in previous works (e.g., Andre & Doulcet, 1970; Loucks, 1999; Mazzullo &
Mazzullo, 1970; Zhai & Zha, 1982). In addition, collapsed paleokarst systems pose potential drilling hazards
due to the weakness of the paleokarst zones (Qi et al., 2014; Soriano et al., 2019; Zhao et al., 2018). Therefore,
delineating subsurface collapsed paleokarst systems is important for petroleum reservoir characterization
and production.
Three-dimensional seismic images are widely and commonly used to map paleokarst systems and to study
paleokarst related features. To facilitate the paleokarst interpretation in 3-D seismic images, some seismic
attributes including coherence (Bahorich & Farmer, 1995; Li & Lu, 2014; Marfurt et al., 1999), structural
curvature (Al-Dossary & Marfurt, 2006; Di & Gao, 2016; Roberts, 2001), reflector rotation (Marfurt &
Rich, 2010), and spectral decomposition (Chen, 2016; Qi & Castagna, 2013) are calculated to highlight the
paleokarst features in 3-D seismic images. Automatic paleokarst interpretation based on such attributes
alone, however, remains challenging because these attributes are typically sensitive to noise and other types
of structural discontinuities or variations that are unrelated to the collapsed paleokarst features. Significant
human interactions are still required to delineate the collapsed paleokarst systems by integrating multiple
types of seismic attributes (Khatiwada et al., 2013; Qi et al., 2014; Sullivan et al., 2006; Zhao et al., 2018)

©2020. American Geophysical Union. We propose a convolutional neural network (CNN) to fully detect paleokarst collapse features in 3-D seismic
All Rights Reserved. images and extract 3-D paleokarst chimney tubes all at once, which would provide a clear 3-D perspective of

WU ET AL. 1 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

the collapsed paleokarst systems in the subsurface. CNN-based methods are powerful for multidimensional
image processing tasks including image classification (e.g., He et al., 2016; Krizhevsky et al., 2012; Simonyan
& Zisserman, 2014; Szegedy et al., 2015), object detection (e.g., Lin et al., 2017; Liu et al., 2016; Redmon
et al., 2016; Ren et al., 2015; Sermanet et al., 2013), and image segmentation (e.g., Badrinarayanan et al., 2017;
Chen et al., 2017; Long et al., 2015; Ronneberger et al., 2015). Recently, CNNs have been increasingly applied
in geoscience problems (Bergen et al., 2019; Bianco et al., 2019) including the interpretation of geologic faults
(e.g., Di et al., 2019; Lu et al., 2018; Wu, Liang, et al., 2019; Wu, Shi, et al., 2019; Zhao & Mukhopadhyay, 2018),
horizons (e.g., Geng et al., 2019; Wu, Zhang, et al., 2019; Wu et al., 2020), salt bodies (e.g., Di & AlRegib, 2020;
Di et al., 2018; Guillen et al., 2015; Shi et al., 2019), and channels (Pham et al., 2019) in seismic images.
In this paper, we consider the characterization of paleokarst collapse features in a 3-D seismic image as
an image segmentation problem and solve this problem by using a supervised CNN. The architecture of
our CNN is modified from the U-net (Ronneberger et al., 2015) by upgrading the convolutional layers from
2-D to 3-D and reducing the number of layers and features at each layer, which significantly saves GPU
memory and makes it applicable to our task of processing 3-D large seismic images. Training a CNN typically
requires a large amount of training data sets including input images and the corresponding label images. In
our problem of collapsed paleokarst interpretation in seismic images, however, we lack the input seismic
images, especially labeled images with interpreted paleokarst features. Manually interpreting or labeling the
paleokarst in seismic images is time-consuming, which makes it hard to prepare a large amount of reliable
or diverse label images for training a CNN.
To solve this problem of lacking training data sets, we propose a workflow to automatically generate
numerous synthetic seismic images, where a variety of realistic structural patterns and collapsed paleokarst
features are simulated with well-defined functions. The collapsed paleokarst features in these generated
synthetic seismic images are well defined, and we are able to automatically obtain the corresponding label
images with the ground-truth of collapsed paleokarst features. With this workflow, we automatically gen-
erate 100 pairs of 3-D synthetic seismic images and the corresponding labeled paleokarst images. These
100 pairs, with further data augmentation, are proven sufficient to train our CNN for collapsed paleokarst
characterization in 3-D seismic images. Although trained only by synthetic seismic images, our CNN shows
powerful performance in delineating collapsed paleokarst features in field seismic images. The detected
paleokarst results are consistent with previous manual interpretations and are superior to the collapsed
paleokarst features characterized by conventional seismic attributes of curvature and coherence. From the
images of collapsed paleokarst features computed by the trained CNN, we are able to automatically extract
the boundaries of 3-D collapsed paleokarst systems, which may be extensively and complicatedly devel-
oped in a 3-D seismic image. Based on the extracted 3-D collapsed paleokarst systems, we can further
automatically measure their geometric parameters for quantitative analysis of paleokarst development.
The structure of the paper is organized as follows: We start with discussing two 3-D field seismic surveys
in the FWB, the related geological background, and observations of the collapsed paleokarst features in
the seismic images. We then present a novel workflow to numerically simulate realistic folding structures
and collapsed paleokarst systems in 3-D synthetic seismic images. We further discuss the architecture of
the proposed CNN and train it with various synthetic data sets that are automatically generated by using
the proposed numerical simulation workflow. We finally demonstrate the applicability of the trained CNN
with two 3-D field seismic images and compare the results with those from both manual interpretation by
previous researchers and conventional automatic methods of curvature and coherence.

2. Field Seismic Data and Geologic Background

As denoted by the blue and red blocks in Figure 1, the two time-migrated 3-D seismic data sets (Figures 2a
and 3a) used in this paper were both acquired at the FWB but are located at different structural settings
within the basin. The first seismic data (Figure 2a) was acquired by the Bureau of Economic Geology at the
University of Texas at Austin and its industry partners including Arch Petroleum, Enserch, and Oxy USA
(Hardage, 1996). As denoted by the blue block in Figure 1a, the corresponding seismic survey covers an area
of 67 km2 at the Boonsville field in the northern FWB. The Boonsville field (colored cyan in Figure 1a) is
one of the largest gas fields in the United States. The second seismic image (Figure 3a) is from the central
FWB by the Marathon Oil Company with 16 live receiver lines, forming a 324 km2 wide-azimuth survey (red
block in Figure 1a) with a fine 16 × 16 m2 Common Depth Point (CDP) bin size (Qi et al., 2014).

WU ET AL. 2 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a)

(b)

Figure 1. In the general structure map of the Fort Worth Basin (modified after Wen et al., 2016), the red and cyan
blocks denote the survey locations of the two 3-D seismic data sets used in this paper. The diagram in (b) shows a
stratigraphic section (modified after Alhakeem, 2013) crossing the Fort Worth Basin from north to south as denoted by
the magenta line in (a).

WU ET AL. 3 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a) (b)

(c)

(d) (e)

Figure 2. The 3-D seismic image (a) is from the northern Fort Worth Basin as denoted by the blue block in Figure 1a.
In the paleokarst probability map (b), computed by our CNN-based method, the vertically extended paleokarst
chimney tubes are highlighted by relatively high probability values (colored by red). The vertical seismic section of
Line A (denoted in blue in a and b) is displayed with manual interpretation (c) by McDonnell et al. (2007), the
CNN-based paleokarst probability image (d), and the boundaries (cyan curves in e) that are automatically extracted
from the probability image.

The FWB is a major petroleum producing geological system which has yielded over 120 billion of gross
production since 2001 (George, 2016). The FWB is formed along the Ouachita thrust front as shown in
Figure 1a. It is a wedge-shaped (Figure 1a), asymmetrical, and northward deepening (Figure 1b) depression
with about 12,000 feet (3,657.6 m) of strata in its deepest northeast area near the Ouachita thrust belt and
Muenster Arch (Figure 1a). Previous studies suggested that the FWB was formed during the late Paleozoic
Ouachita Orogeny, which is a major tectonic event of thrust-fold deformation associated with the oblique
lithospheric convergence of the North American and South American plates (Walper, 1982).
Figure 1b shows a stratigraphic cross section of the FWB from north to south as denoted by the magenta
line in Figure 1a. The lower Ellenburger group of carbonate strata is unconformably overlain by the Ordovi-
cian Viola and Simpson formations, which pinch out just east of the first seismic survey (Figure 2 and
denoted by the blue block in Figure 1a) (McDonnell et al., 2007) and are also missing in the second seis-
mic survey (Figure 3 and denoted by the blue block in Figure 1a). Therefore, in the two study areas, the
lower Ellenburger group carbonate strata are unconformably overlain by the organic-rich Barnett Forma-
tion (as shown in Figure 1b), which recorded deposition of the deep-water foreland basin during the late
Mississippian (Sullivan et al., 2006) and plays a critical role in forming multiple gas fields in northern
Texas (Pollastro et al., 2007). The Barnett group is overlain by the Pennsylvanian formation, which includes
the Bend (lower Pennsylvanian), Strawn, and Canyon groups. As marked in Figure 2c, the Bend group

WU ET AL. 4 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a) (b)

(c) (d)

Figure 3. The 3-D seismic image (a) is from the central Fort Worth Basin as denoted by the red block in Figure 1a. The
′ ′
vertical seismic sections in (b) are extracted at Lines AA and BB , which cross a large paleokarst chimney as denoted
by the yellow lines in (a). From the 3-D seismic image, a paleokarst probability image (c) is computed by using our
CNN-based method to highlight the paleokarst chimneys, which are mostly consistent with the manually interpreted
′ ′
paleokarsts in the vertical sections in (a). The 2-D probability sections (d), extracted at Lines AA and BB , provide a
close-up view of our CNN-based detection of the paleokarst chimney, which is consistent with the manual
interpretation denoted by the dashed cyan lines.

(Morrowan and Atokan stages of the lower Pennsylvanian) consists of Marble Falls, Vineyard, Runawary,
Davis, and Caddo stratigraphic units (McDonnell et al., 2007).
Paleokarst caves and solution collapses are extensively generated in the Ellenburger group because it went
through a series of karst events ranging between the post Ellenburger and Early Pennsylvanian (Canter
et al., 1993). In addition, the strata of the Silurian and Devonian strata are missing in the two study areas,
which provides additional time for generating paleokarst in the Ellenburger group (Sullivan et al., 2006;
Zhao et al., 2018). Figure 4a shows a model (Kerans, 1988) of the paleokarst system and the associated
developments in the Ellenburger group. As shown in this model, the collapsed paleokarst system typically
contains the three basic components of a paleocave floor, fill, and roof. During burial, paleocaves often
collapse and sequentially produce a collapsed zone with some unique associated features as shown in the
diagram (Figure 4b) of the coalesced collapsed-paleocave system proposed by Loucks (1999). According to
the diagram, a coalesced collapsed paleokarst system typically contains fill deposits (breccias and stratified

WU ET AL. 5 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a)

(b)

Figure 4. The 3-D paleokarst model (modified after Kerans, 1988) in (a) shows the general features of a paleokarst
system where the sag structures or paleokarst chimneys may vertically extend hundreds of meters. The 2-D diagram
(b) of a single collapsed paleokarst system (modified after Loucks, 1999) shows features of the associated collapse and
extensive brecciation, where deformations (faults, fractures, and sag structures) are developed in postkarst-deposited
strata.

WU ET AL. 6 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a) (b) (c)

(d) (e) (f)

Figure 5. The proposed workflow for generating synthetic training data sets. In this workflow, we start with a flat reflectivity model (a) in which we then
simulate realistic folding structures (b) and collapsed paleokarst features (c). We further convolve the reflectivity model with a Ricker wavelet to compute a 3-D
synthetic seismic image (d). As the geometry of the collapsed paleokarst systems is well defined in the simulation, we are able to automatically obtain a binary
label image (e) with the ground-truth of the collapsed paleokarst areas labeled with ones. From the labeled binary paleokarst image (e), we can automatically
extract the 3-D collapsed paleokarst systems shown in (f).

deposits), a broader damage zone of crackle fractures, circular faults, and chimney or sag structures
(Figure 4a) in younger formations. Such collapse chimneys or sags in our seismic images vertically extend
upward from the Ellenburger group through the Mississippian (e.g., Barnett Shale) and Pennsylvanian strata
(as interpreted by McDonnell et al., 2007 in Figure 2c) over a distance of almost 800 m (Hardage et al., 1996).
Therefore, interpreting the collapsed paleokarst-related features, especially the sags with long vertically
upward extents, is important for reducing the risks of drilling and hydraulic fracturing in exploring the
main oil and gas reservoirs in the Barnett Shale (deposited during the Mississippian period) and Bend group
(deposited during the Middle Pennsylvanian period).
Three-dimensional seismic imaging provides an effective way to interpret the paleokarst system in the 3-D
subsurface. As shown in Figures 2a and 3a, collapsed paleokarst features, especially the large-scale fault
and sag structures, can be extensively observed. Manually interpreting such features, however, remains
a time-consuming task, while accurately and fully delineating the paleokarst-related features in the 3-D
subsurface space is challenging. We propose a CNN to efficiently and accurately segment or detect the
paleokarst-collapse features from the 3-D seismic images by computing a paleokarst collapse probability
volume as shown in Figures 2b, 2d, 3c, and 3d. From the probability volumes, we are able to further automat-
ically obtain the 3-D bodies of the paleokarst chimneys or sags in the 3-D subsurface space, which provides
a convenient way to quantitatively analyze collapsed paleokarst systems.

3. Seismic Simulation of Paleokarst Collapses

To train a CNN for segmenting the paleokarst-related features from a 3-D seismic image, we propose a work-
flow (as shown in Figure 5) to automatically generate many 3-D synthetic seismic images where the folding
structures and collapsed paleokarst features are diverse and realistic. The well-defined collapsed paleokarst
features can be automatically and perfectly labeled to obtain the target images for training the CNN.
3.1. Simulating Folding Structures
In this workflow, we begin with an initial 3-D reflectivity model r 0 (x, y, z) with flat layers (as shown in
Figure 5a), where the reflectivity values at each layer are randomly generated and smoothly vary in space.

WU ET AL. 7 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

We then randomly generate folding structures by vertically shearing the flat model, where the shearing field
s(x, y, z) is defined as a combination of two functions:

s(x, 𝑦, z) = s1 (x, 𝑦, z) + s2 (x, 𝑦, z), (1)

where the first function s1 is a combination of N Gaussian functions scaled by a vertically damping function
as suggested by Wu, Liang, et al. (2019) and Wu et al. (2020):
2 2
1.5 ∑
k=N (x−bk ) +(𝑦−ck )
2𝜎 2
s1 (x, 𝑦, z) = ak e k . (2)
zmax k=1

In defining this function, the parameters ak , bk , ck , and 𝜎 k are all randomly chosen from some predefined
ranges as discussed by Wu et al. (2020). By using the damping scalar function z1.5 , we gradually decrease
max
the curvature (or bending extent) vertically upward from bottom to top in the model. The second function
is simply a linear function as follows:
s2 (x, 𝑦, z) = px + q𝑦 + c0 , (3)

where the parameters p, q, and c0 are constant values that are randomly chosen for each reflectivity model.
This linear function is used to generate the purely dipping structures in the reflectivity model. The slopes
of the dipping structures in the x− and y− directions, respectively, are defined by the parameters p and q,
which are randomly chosen from a limited range of [− 0.25, 0.25] to avoid generating layers with extremely
high slopes after combining with s1 (Wu et al., 2020). We set c0 = −pxc − q𝑦c to ensure that the center trace
r(xc , yc , z) of the reflectivity model is not shifted.
With the defined vertical shearing field s(x, y, z) = s(x, y, z) + s2 (x, y, z), we are able to compute a folded
reflectivity model r(x, y, z) from the initial model r 0 (x, y, z) as follows:

r(x, 𝑦, z) = r0 (x, 𝑦, z + s(x, 𝑦, z)), (4)

where a classic sinc interpolation is applied to map the reflectivity model from the initially flat space to
the folded or sheared space. As the parameters for defining the shearing functions are randomly chosen
in generating the folded reflectivity model, the folding structures in each model are unique. In addition,
each parameter has many options; therefore, we have numerous possible combinations of parameters for
generating numerous folded reflectivity models with various structures.
3.2. Simulating Collapsed Paleokarst Structures
After generating the folded reflectivity model, the next step is to further simulate the collapsed paleokarst
structure features in the model. Based on the observations from the field seismic images in Figures 2 and 3,
the most dominant features of the collapsed paleokarst structures apparent in the seismic images are the
chimneys or sags, the fractures or faults within the chimney zones, and the steep circular or cylindrical faults
bounding the sags. As shown in the vertical seismic sections in Figures 2 and 3, each collapse chimney or sag
is apparent as a vertically elongated tube with an approximately ellipsoidal shape. The vertically extended
chimney or sag structures are cut by horizontal seismic slices, where we observe obvious circular features
or onion rings (as denoted by the red and magenta arrows in Figures 2a and 3a). Such circular features
indicate lateral changes in reflection time or depth due to the sag structures. Within some of the sag zones,
the seismic reflections are highly chaotic or discontinuous (as denoted by the magenta arrows in Figure 3
and magenta ellipses in Figure 3b), which indicates that the deposits within these sag zones may be highly
fractured or faulted.
To simulate a collapse chimney tube in the 3-D reflectivity model, we construct a 3-D vertically elongated
ellipsoid
𝑓 (p) = (p − c)⊤ R⊤ AR(p − c), (5)

with which we define the 3-D tube area of a chimney as follows:

{
if 𝑓 (p) ≤ 1 ∶ inside chimney cube
. (6)
if 𝑓 (p) > 1 ∶ outside chimney cube

WU ET AL. 8 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

In the ellipsoidal function (Equation 5), p = (x, 𝑦, z) represents the coordinates of a point in the 3-D reflec-
tivity model. c = (cx , c𝑦 , cz ) represents the center of the ellipsoid, which is randomly chosen within the 3-D
space of the model. A is a diagonal matrix defined by the three radii of the ellipsoid as follows:
1
⎡ r2 0 0 ⎤
⎢ 0x 1 0 ⎥
A=⎢ r𝑦2 ⎥, (7)
⎢0 0 1 ⎥
⎣ rz2 ⎦

where the values of r x , r y , and r z , respectively, are randomly chosen from the predefined ranges [1, 12], [1, 12],
and [4, 80] to construct ellipsoids of different sizes. Here, we set a wider range for r z to obtain primarily
vertically elongated ellipsoids. In addition, we set r x > 0.1r z and r y > 0.1r z to avoid generating extremely long
but narrow ellipsoids that are unrealistic in practice. R is a rotation matrix that rotates the ellipsoid around
the x− and y− axes as follows:

⎡1 0 0 ⎤ ⎡ cos 𝛽 0 sin 𝛽 ⎤
R = ⎢ 0 cos 𝛼 − sin 𝛼 ⎥ ⎢ 0 1 0 ⎥, (8)
⎢ ⎥⎢ ⎥
⎣ 0 sin 𝛼 cos 𝛼 ⎦ ⎣ − sin 𝛽 0 cos 𝛽 ⎦

where the rotation angles 𝛼 and 𝛽 are randomly chosen from a narrow range [− 10◦ , 10◦ ] to construct a
slightly dipping or vertically aligned ellipsoid. By randomly choosing all the parameters in Equation 5, we
are able to construct numerous possible ellipsoids with various shapes, sizes, orientations, and locations.
However, in reality, a paleokarst chimney tube is not a perfect ellipsoid; we therefore randomly add some
smooth perturbations to obtain irregular ellipsoids as shown in Figure 5f.
After defining the 3-D tube areas of the paleokarst chimneys, we then define the sag structures within
the chimney cubes. The sag structures are typically downward bending. We therefore vertically shift the
reflectors (inside a chimney cube) downward as follows to obtain a reflectivity model r k (x, y, z):

rk (x, 𝑦, z) = r(x, 𝑦, z + sk (x, 𝑦, z)), (9)

where the vertical shifts sk (x, y, z) are defined as

{
0, if 𝑓 (x, 𝑦, z) > 1
sk (x, 𝑦, z) = . (10)
𝛾( 𝑓 (x, 𝑦, z) − 1) + 𝜖(x, 𝑦, z), if 𝑓 (x, 𝑦, z) ≤ 1

In this equation, f (x, y, z) is the ellipsoidal function defined in Equation 5. 𝛾 is a positive scalar that is ran-
domly chosen from the range [10, 20]. 𝜖 (x, y, z) is a random perturbation field implemented to simulate
potential fractures or faults that dislocate the reflectors within the chimney tube. The perturbation field is
relatively small compared to the first term of 𝛾 (f (x, y, z) − 1). When the perturbation field is close to 0, the
shifts sk (x, y, z) ≈ 𝛾 (f (x, y, z) − 1) will be nonpositive and smoothly decrease from 0 at the chimney bound to
the most negative at the center of the chimney tube. In this case, the shifts will vertically shear the reflectors
so that they are smoothly curved downward, as shown between the cyan dashed lines (cylindrical faults) in
Figure 5c. Such generated sag structures will produce clear circular features or onion rings on the horizontal
slice as denoted by the red arrows in Figure 5c. If the perturbation field is significant, the shifts sk (x, y, z) will
no longer be smooth but noisy, resulting in the generation of curved and dislocated (fractured or faulted)
reflectors within the chimney tubes as denoted by the magenta arrows and dashed ellipses in Figure 5c. By
using this numerical method, we are able to automatically simulate various collapsed paleokarst features
(as in Figure 5c) that look realistic and are comparable to the real collapsed paleokarst features in the field
seismic images (Figures 2 and 3).
In summary, we have proposed a general workflow (Figure 5) to numerically simulate realistic folding and
collapsed paleokarst structure features in a reflectivity model. The parameters used to model the folding
and collapsed paleokarst features are summarized in Figure 6. A specific set of such parameters yields a
reflectivity model with unique folding and collapsed paleokarst features. These parameters are not fixed
but are randomly chosen to form numerous possible combinations, which allows us to generate numerous
reflectivity models with various structural and collapsed paleokarst features. This is critical for our next step
of training a CNN for collapsed paleokarst characterization which requires rich training data sets.

WU ET AL. 9 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 6. A summary of the parameters used to simulate folding and collapsed paleokarst structural features.

3.3. Training Data Sets

After constructing a reflectivity model (Figure 5c) with simulated folding structures and collapsed paleokarst
features, our next step is to simulate a synthetic seismic image (Figure 5d) by convolving a Ricker wavelet
with the reflectivity model. Instead of convolving in the vertical direction as in conventional methods, we
compute the convolution here in the direction perpendicular to the reflector structures as suggested by Wu,
Liang, et al. (2019) and Wu et al. (2020). The peak frequency of the Ricker wavelet is also randomly chosen
from a predefined range (10–35 Hz) to compute a synthetic seismic image with various frequency compo-
nents. As field seismic images are typically not as clean as the synthetic seismic image in Figure 5d, we
further add random noise or the noise from field images to make the synthetic one more realistic. Figure 5e
shows the corresponding label image that is automatically computed by labeling the paleokarst chimney
cubes with ones while the nonkarst areas with zeros. From the label image, we can further automatically
obtain the 3-D bodies of the paleokarst chimney cubes (Figure 5f) by simply calculating the isosurfaces (with
isovalue 0.5) of the label image using the method of marching cubes (Lorensen & Cline, 1987).
In training a CNN for detecting the collapsed paleokarst features, we need to prepare many training data
set pairs consisting of input seismic images (like the one in Figure 5d) and the corresponding label images
(like the one in Figure 5e). Fortunately, in our workflow, the parameters for simulating the synthetic seis-
mic images and the collapsed paleokarst features are not fixed. By randomly choosing the parameters, we
are able to obtain numerous combinations of the parameters to calculate numerous pairs of training image
with various and diverse structures and collapsed paleokarst features. In this work, we generate 120 pairs
of synthetic seismic and label images, 100 pairs for training and the rest for validation. Figure 7 shows eight
pairs of the automatically generated training data sets, where the first and third rows show the input seis-
mic images, while the second and fourth rows show the corresponding label images of collapsed paleokarst
features overlain with the seismic images. Notice that noise has been added to the synthetic seismic images.
The noise in the first row of images is extracted from various areas of the field seismic images. The noise
in the third row of images is randomly generated. In adding the noise, the noise-to-signal-ratio is randomly
selected from the range [0, 0.6]. The collapsed paleokarst zones are randomly scattered in the synthetic seis-
mic image, some are cut at the image boundaries. This is consistent with the prediction for a large field
seismic image, where we often need to cut a large volume into small subvolumes and make the prediction
for each subvolume.
To further increase the diversity of the training and validation data sets, we apply two types of simple data
augmentation to the original 100 pairs of training data sets. The first type of data augmentation involves
rotating the images around the vertical axis by 0◦ , 90◦ , 180◦ , and 270◦ , which increases the number of train-
ing data sets by a factor of 4. The second type of data augmentation involves randomly cutting smaller
subimages from the original images and use the subimages to train the CNN. The dimension size of each
of the originally generated images is 256 × 256 × 256 (samples). During the training process, we randomly
cut 128 × 128 × 128 subimages from the larger images, which ensures that the training data sets are mostly
different for different training epochs. By doing this, we are able to significantly increase the number and
diversity of the training data sets. In addition, training the CNN with smaller images can greatly reduce GPU
memory and computational costs during the training.

WU ET AL. 10 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 7. By using the proposed workflow (Figure 5), we are able to automatically generate numerous synthetic training data set pairs to train our CNN for
recognizing collapsed paleokarst features in seismic images. Shown here are eight pairs of synthetic seismic images (first and third rows) and the corresponding
binary label images (second and fourth rows) with the collapsed paleokarst features denoted by ones.

4. Deep learning for Detecting Collapsed Paleokarst Features

We design a CNN, as shown in Figure 8, to detect the collapsed paleokarst features in a 3-D seismic
image. The architecture of the designed CNN consists of a U-shaped network followed by a residual block
(Figure 9). The architecture of the 3-D U-shaped network is modified from the original 2-D U-net proposed
by Ronneberger et al. (2015) for 2-D medical image segmentation, making it applicable to our 3-D problem.
Compared to the original U-net, we reduce the number of layers and features at each layer to significantly
save memory and computational costs which is especially important for our 3-D seismic image segmenta-
tion problem. In addition, we reduce the number of downsamplings or poolings from four to three because
the dimension size (128 × 128 × 128) of our 3-D training images is relatively small and the downsampled
images after four poolings would be too small (only 8 × 8 × 8) to effectively preserve the spatial or structural
features in the original images.
Similar to the original U-net (Ronneberger et al., 2015), the modified U-shaped network still preserves the
general architecture of a bottleneck at the middle and symmetric left contracting and right expansive paths,
each of which contains three steps (reduced from four steps in the original U-net). In the left contracting

WU ET AL. 11 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 8. To detect the collapsed paleokarst features from an input seismic image, we design a CNN with the
architecture of a U-net followed by two residual blocks (Figure 9).

path, each step consists of two 3 × 3 × 3 convolutional layers (each followed by a ReLU activation func-
tion), and a 2 × 2 × 2 max-pooling layer with stride 2. In each convolutional layer, multiple 3 × 3 × 3 filters
are applied to convolve with the input image to obtain multiple output feature maps. The ReLU activation
function is a nonlinear function (𝑓 (x) = max(0, x)) applied to the feature maps to increase the nonlinear-
ity of the network (Krizhevsky et al., 2012). The max-pooling layer reduces the size of the feature map by
a factor of 2 and keeps only the maximum value with a 2 × 2 × 2 sliding window. At the first, second, and
third steps, the numbers of feature maps at the conventional layers are 16, 32, and 64, respectively, which
are significantly reduced with respect to those of the original U-net. The bottleneck at the middle contains
two 3 × 3 × 3 convolutional layers, each of which contains 512 features and is, again, followed by a ReLU
activation function.
The three steps in the right expanding path gradually upsample the downsampled features back to the orig-
inal size. Each step consists of a 2 × 2 × 2 upsampling operation with stride 2, a concatenation to combine
the features from the left path, and two 3 × 3 × 3 convolutional layers. Each convolutional layer is followed
by a ReLU activation function. The upsampling operation is implemented by using the UpSampling3D layer
defined in Keras (Chollet, 2015), which increases the size of an input image by a factor of 2. After the
U-shaped network, we further add another two residual blocks whose architecture is shown in Figure 9.
Each residual block contains two 3 × 3 × 3 convolutional layers, each of which contains 16 features and is
followed by a ReLU activation function. The final output layer is a 1 × 1 × 1 convolutional layer followed by
a sigmoid activation function that generates a 3-D probability image with values in the range [0, 1] as shown
on the right of Figure 8.
In constructing the network (Figure 8), the parameters of each layer are randomly initialized, which are
required to be further updated to create a good mapping of an input seismic image (left of the network) to
an output paleokarst image (right of the network). Updating the network parameters involves a training
process that uses an optimization algorithm to iteratively update the parameters until the error between the
outputs and the label images converges on the training data set. As we consider the detection of collapsed
paleokarst features a binary image segmentation problem, we use the following loss of binary cross-entropy
to measure the error between the outputs and the label images:

1∑
i=N
=− 𝑦 log(pi ) + (1 − 𝑦i ) log(1 − pi ), (11)
N i=1 i

Figure 9. The structure of a residual block.

WU ET AL. 12 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 10. The training (orange) and validation (blue) curves over 25 epochs. The loss (a) for both the training and
validation data sets quickly converges to a small value (nearly 0.03), while the prediction accuracy (b) significantly
increases to nearly 0.99.

where N is the number of samples in a 3-D output or label image and yi and pi , respectively, represent the
binary label and predicted values at the ith image sample. In the training process of optimizing the network
parameters, we use the Adam method (Kingma & Ba, 2014) with a learning rate of 0.0001. We iteratively train
the network with 25 epochs. One epoch involves passing the entire training data set (400 pairs, including the
original and rotated images) forward and backward through the neural network once. At each epoch, the
training images with a smaller size 128 × 128 × 128 are cut from larger generated images (256 × 256 × 256 as
in Figure 7) at random locations. Therefore, the training images at different epochs are mostly different.
In addition, each seismic image is subtracted by its mean and then divided by its standard deviation to obtain
a normalized image. This normalization process is applied to ensure the amplitudes of all the seismic images
are consistently distributed. Such normalizations modify the seismic amplitudes but do not change the geo-
metric features or structural patterns (relative amplitude variations), from which the collapsed paleokarst
features are distinguished from the nonkarst features. At each epoch of the training process, we evaluate the
updated network by computing the losses and the accuracies over both the synthetic training and validation
data sets. The accuracy a is defined as follows:
{
1∑
i=N
1, if 𝑦i = ⌊pi ⌉
a= e(𝑦 , p ), e(𝑦i , pi ) = , (12)
N i=1 i i 0, if 𝑦i ≠ ⌊pi ⌉

where, again, N is the number of samples in a 3-D output or label image and yi represents the binary label
value at the ith image sample. ⌊pi ⌉ represents the nearest integer to the predicted value pi (0 ≤ pi ≤ 1). We
stop the training at the 25 epoch as both the training and validation losses converge to nearly 0.33, while the
prediction accuracies over both the training and validation data set increase to nearly 0.99 (Figure 10).

5. Applications
To verify the effectiveness and generalization of the CNN trained with 25 epochs over the synthetic data
sets, we apply the trained CNN to both synthetic and field 3-D seismic images. Although the CNN is trained
by images of size 128 × 128 × 128, the size of the image input into the CNN for prediction is not fixed. This
means that the trained CNN can be flexibly applied to images of arbitrary dimension sizes, so long as each
of the dimensions can be divided by 8 due to the three downsamplings included in the CNN architecture
(Figure 8).
5.1. Synthetic Examples
Figure 11 shows the application of the coherence-based, curvature-based, and our CNN-based methods to
four synthetic seismic images (first column in Figure 11) for collapsed paleokarst characterization. These
four images (each with 256 × 256 × 256 samples) were generated by using the proposed numerical workflow
(Figure 5) and are different from the images in the training data set. The first two images contain noise
that is extracted from the field seismic images (Figures 2a and 3a). The second two images contain heavier
random noise. The second row of Figure 11 shows the paleokarst probability images pc (x) computed from
seismic coherence (Marfurt et al., 1999) as follows:

pc (x) = 1 − c(x)6 , (13)

WU ET AL. 13 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 11. Characterization of paleokarst collapse features in four synthetic validation seismic images (first row) by using coherence-based (second row),
curvature-based (third row), and our CNN-based (fourth row) methods. Our CNN-based method shows the best performance in obtaining accurate paleokarst
characterizations that are highly consistent with the ground-truth paleokarst features (last row).

where c(x) represents a seismic coherence map that typically highlights seismic discontinuities with rela-
tively lower values. The value of the exponent, 6, is somewhat arbitrarily chosen and is used to increase the
contrast between image samples with low and high coherences. The coherence-based paleokarst probability
images (second row in Figure 11) highlight most seismic discontinuities, including those in the collapsed
paleokarst zones and noise.
The third row of Figure 11 shows the paleokarst probability images computed from the most positive seismic
curvature Al-Dossary and Marfurt (2006) as follows:
( ) {
𝜅c (x) − 𝜅c (x)min 6 𝜅(x), if 𝜅(x) ≤ 0
p𝜅 (x) = 1 − , 𝜅c (x) = . (14)
−𝜅c (x)min 0, if 𝜅(x) > 0

In this equation, 𝜅 (x) represents the most positive curvature map computed from an input seismic image.
The structures within paleokarst collapse zones are typically curved downward and accordingly show

WU ET AL. 14 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a) (b)

Figure 12. (a) Precision-recall and (b) ROC curves are calculated to quantitatively evaluate paleokarst detection on the
synthetic validation volumes shown in Figure 11. Our CNN-based method (red curves) provides much more accurate
results than the coherence-based (black curve) and curvature-based (blue curve) methods.

negative values in the curvature map 𝜅 (x); we therefore clip 𝜅 (x) to obtain 𝜅 c (x) with only negative values.
We further define the curvature-based paleokarst probability by using the normalized and clipped curva-
ture map 𝜅 c (x) as shown in the above equation, where the power 6, again, increases the contrast between
image samples with low and high curvatures. The curvature-based paleokarst probability images (third row
of Figure 11) highlight some of the collapsed paleokarst zones as well as some downward folding structures
that are unrelated to the paleokarst collapses.

(a) (b)

(c) (d)

Figure 13. From a 3-D field seismic image (a), we calculate the paleokarst probability images by using
coherence-based (b), curvature-based (c), and our CNN-based (d) methods. The paleokarst probability image computed
by using our CNN-based method displays much cleaner and more accurate collapsed paleokarst features than those
computed by using the other two methods.

WU ET AL. 15 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a) (b)

Figure 14. (a) Precision-recall and (b) ROC curves are calculated to quantitatively evaluate collapsed paleokarst
detection on the 2-D vertical sections of the field seismic image in Figure 3. Our CNN-based method (red curves)
provides much more accurate results than the coherence-based (black curve) and curvature-based (blue curve)
methods. In calculating the evaluation curves, all the collapsed paleokarst detection results are compared to the
ground-truth collapsed paleokarsts that are manually interpreted on the two vertical sections in Figure 3a.

The fourth row of Figure 11 shows that our CNN-based paleokarst probability images visually display much
more accurate collapsed paleokarst characterizations than the coherence-based and curvature-based prob-
ability images. The collapsed paleokarst features highlighted (colored by red) in our CNN-based probability
images are consistent with collapsed paleokarst ground truth in the label images (last row of Figure 11). In
addition, our CNN-based method is highly efficient. It took only milliseconds to compute such a probability
image from an input seismic image with 256 × 256 × 256 samples by using a Quadro P6000 GPU.
In order to quantitatively evaluate the collapsed paleokarst characterization results, we further compute
precision-recall (Martin et al., 2004) and receiver operating characteristic (ROC) (Provost et al., 1998) curves
as shown in Figure 12. From the precision-recall curves (Figure 12a), we observe that our CNN-based
method shows the highest precision for all choices of recall. From the ROC curve, our CNN-based method
shows the highest true positive rates but the lowest false positive rates.
5.2. Field Examples in the FWB
In addition to the synthetic examples in Figure 11, we further use two field examples acquired from the FWB
(Figure 1a), to verify the generalization of the trained CNN to field seismic images. Compared to those in
the synthetic seismic images (Figure 11), the structures and features in the field images may be significantly
different.
Figure 13a shows the first 3-D seismic image (the same one in Figure 3a), which was acquired at the sur-
vey denoted by the red block in the FWB (Figure 1a). In this example, the paleokarst probability image
(Figure 13b) computed by using the coherence-based method (Equation 13) highlights some collapsed pale-
okarst zones but is highly sensitive to noise. The curvature-based method (Equation 14) yields a better
paleokarst probability image that detects more collapsed paleokarst zones but still contains some noisy
features (e.g., the area denoted by the cyan circle) and fails to completely delineate some true paleokarst
chimneys (as denoted by the magenta arrows in Figure 13c). Compared with the coherence-based and
curvature-based methods, our CNN-based method yields the best probability image (Figures 3c and 13d),
where the highlighted collapsed paleokarst zones are most consistent with the manual interpretation in
Figure 3a.
As shown in the vertical sections in Figure 3c, the boundaries of the highlighted paleokarst chimneys are
consistent with the manually interpreted cylindrical faults (cyan dashed lines) bounding the chimney tubes.
In Figure 3d, we display two seismic Sections AA′ and BB′ that vertically slice a large paleokarst chimney
tube, where we are able to more clearly view the details of the result. In these two sections, we observe that
the highlighted chimney tube starts appearing in the Ellenburger group (below the blue horizon) and con-
tinuously extends upward into the Pennsylvanian strata at the top, which is consistent with the conclusions
of previous studies (e.g., Hardage et al., 1996; McDonnell et al., 2007) on the collapsed paleokarst systems
in the FWB. In addition, the boundaries of the detected paleokarst chimney are highly consistent with the
manually interpreted bounding cylindrical faults as denoted by the cyan dashed lines in Figures 3b and 3d.
The precision-recall and ROC curves in Figure 14 show quantitative evaluations of the CNN-based (red),

WU ET AL. 16 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 15. The 3-D paleokarst chimney tubes in (a)–(c) are automatically extracted from the paleokarst probability
image (Figure 13b), the most negative curvature image (Figure 13c), and the most positive curvature image
(Figure 13d), respectively.

WU ET AL. 17 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 16. The image in (a) shows a different view of the 3-D seismic image in Figure 2a. From this 3-D seismic image,
we compute a 3-D paleokarst probability image overlain with the seismic image in (b) as well as in Figure 2b. From the
paleokarst probability image, we are able to automatically reconstruct the 3-D paleokarst chimney tubes (c) by simply
extracting the isosurfaces (with an isovalue of 0.5) of the probability image.

curvature-based (black), and coherence-based (blue) paleokarst detections compared to the manual inter-
pretation in the two vertical sections in Figure 3a. Our CNN-based method achieves the best performance
while the coherence-based method shows the worst performance on this field example.
From the paleokarst probability image (Figures 3c and 13b), we can further extract the 3-D paleokarst chim-
ney tubes (Figure 15a) by simply extracting the isosurfaces (with an isovalue of 0.5) using the method of
marching cubes (Lorensen & Cline, 1987). As comparisons, Figures 15a and 15b, respectively, display the
chimney tubes extracted from the coherence-based (Figure 13b) and curvature-based (Figure 13c) proba-
bility images. We observe that the tubes extracted from the coherence-based and curvature-based images
contains many noisy and irregular blocks, most of which are not true paleokarst chimney tubes. In addi-
tion, the shapes of the extracted tubes are highly noisy as well, which makes it difficult to quantitatively
analyze their geometries. In contrast, the 3-D chimney tubes extracted from our paleokarst probability
images are geologically more reasonable and the boundaries of the tubes are more clearly presented. Most
of these extracted chimney tubes narrow upward, which is consistent with the observations by McDonnell
et al. (2007). From each of these extracted chimney tubes, we can automatically and accurately measure the
3-D geometric parameters (e.g., diameter and vertical length) of each tube and analyze the 3-D distribution
of the collapsed paleokarst systems.
The second field seismic image shown in Figures 16a and 2a is from the survey denoted by the blue block
in the FWB (Figure 1)a. The collapsed paleokarst systems in this seismic image have been extensively
analyzed and well interpreted by previous studies (e.g., Alhakeem, 2013; Hardage et al., 1996; McDonnell
et al., 2007). In this work, we compare our automatically interpreted results using the deep learning method
with collapsed paleokarst systems that are manually interpreted by McDonnell et al. (2007). Figures 2b and
16b show the automatically estimated paleokarst probability images obtained by using the our CNN. We
observe that the major paleokarst chimneys are clearly highlighted in red, which correspond to relatively
high probabilities.
To closely visualize the details of the detected collapsed paleokarst features and compare them with the
manual interpretation results, we extract two seismic sections marked by Line A and Line B in Figures 2a
and 2b. As shown in Figure 2b, Line A passes through five circular anomalies with high paleokarst proba-
bility. The extracted seismic section at Line A is shown in Figure 2c, where the paleokarst chimneys (sags)
and the steep concentric or cylindrical faults were manually interpreted by McDonnell et al. (2007). Based
on the manual interpretation, the paleokarst chimneys originate within the Ellenburger Group (the top of
Ellenburger is marked by the blue line in Figure 2c) and extends vertically upward through the lower Atokan
to the Caddo Pool Formation (McDonnell et al., 2007). Figure 2d shows the corresponding paleokarst prob-
ability section (overlain on the seismic section of Line A), where the areas of the five paleokarst chimneys
are clearly highlighted. The contours (cyan curves), extracted from the probability image at the probability
value of 0.5, delineate the boundaries of the paleokarst chimneys, which are consistent with the manually
interpreted concentric faults (dashed black lines in Figure 2c) that bound the chimneys. Line B passes one
circular anomaly with high paleokarst probabilities as shown in Figure 2b. Figure 17 shows the seismic

WU ET AL. 18 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Figure 17. The vertical seismic section at Line B (denoted by the magenta line in Figures 2a and 2b) is displayed with the manual interpretation (a) by
McDonnell et al. (2007), our CNN-based paleokarst probability image (b), and the paleokarst boundaries (cyan curves in c) that were automatically extracted
from the probability image.

section of Line B with manual and our automatic interpretation of the paleokarst chimney in the middle.
Again, our automatic interpretation of the paleokarst chimney (Figures 17b and 17c) in this seismic section
is consistent with the manual interpretation (Figure 17a) by McDonnell et al. (2007).
Figure 18a shows a horizontal coherency slice of the 3-D seismic image, where the chimney structures are
indicated as relatively small coherence anomalies (colored in dark gray). On this coherence slice, 10 circular
chimney structures are manually marked as white circles according to McDonnell et al. (2007). Figure 18b

(a) (b)

Figure 18. A coherency slice (a) of the 3-D seismic image is extracted at 1,140 ms TWT, where the circular paleokarst
features are characterized as dark gray colors and the boundaries of ten circular paleokarsts (white curves) are
manually interpreted by McDonnell et al. (2007). The CNN-based paleokarst probability slice (b) at the same location
shows a much clearer characterization of the circular paleokarst features (highlighted by red colors), from which we
can automatically extract the boundaries of all the circular paleokarst features as denoted by the cyan curves. The
automatically extracted boundaries of the 10 paleokarsts (denoted by yellow numbers in b) are quite consistent with
the manual interpretation (denoted by black numbers in a).

WU ET AL. 19 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

(a) (b) (c)

(d) (e) (f)

Figure 19. To closely view our CNN-based paleokarst detection, we extract a small subseismic image (a) containing the
third circular paleokarst denoted by the cyan arrow in Figures 18b and 2b. The corresponding subvolume of the
paleokarst probability image is overlain with the seismic image in (b), from which we automatically extract the 3-D
paleokarst chimney tube (c) and automatically calculate the geometric parameters of the paleokarst. Horizontal slices
of the subseismic image (a), subpaleokarst image (b), and the extracted paleokarst tube (c) are shown in (d)–(f),
respectively.

shows the corresponding horizontal slice of the paleokarst probability map overlain with the coherence slice.
The circular red areas in Figure 18b denote all the automatically detected chimney structures including the
10 manually interpreted chimneys in Figure 18a. Of course, some of these detected paleokarst chimneys still
need to be verified although most of them look reasonable as shown in the vertical sections in Figures 16b
and 16c. By extracting the contours (with a probability value 0.5) from the paleokarst probability slice, we are
able to automatically obtain the boundaries (cyan curves in Figure 18b) of all the chimneys, which provide
a convenient way to quantitatively analyze the geometry of each chimney, including its shape and size.
Figure 19 shows a subvolume containing the third paleokarst chimney denoted by the cyan arrows in
Figures 18b and 2b. From the subvolume of the paleokarst probability map (Figure 19b), we can automat-
ically extract the paleokarst chimney structure as a vertically extended 3-D tube, as shown in Figure 19c.
Figures 19d–19f, respectively, show horizontal slices of the seismic amplitude, paleokarst probability, and
extracted chimney tube. From the extracted 3-D chimney tube, we can automatically estimate the geometric
parameters of the paleokarst system such as the length and width as marked in Figures 19c and 19f, which
are consistent with the manual measurements provided by McDonnell et al. (2007).

6. Conclusions
We have proposed a CNN-based method to automatically delineate collapsed paleokarst features in 3-D
seismic images, where the paleokarst detection is considered a problem of 3-D binary image segmentation.

WU ET AL. 20 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

The architecture of the CNN is modified from the U-net, which was originally proposed for 2-D medical
image segmentation and is today the most widely used network for common image segmentation problems.
Specifically, we upgraded the original 2-D convolutional layers into 3-D layers and reduced the number of
layers and features at each layer to make it applicable to our problem of collapsed paleokarst interpretation
in 3-D seismic images.
The main challenge or limitation of applying supervised CNNs to geoscience problems, including our col-
lapsed paleokarst interpretation, is the lack of training data sets, especially the label images or interpretation
results. To solve this problem, we proposed an efficient and effective workflow to automatically generate
reflectivity models with realistic folding structures and collapsed paleokarst features and then simulated
the corresponding synthetic seismic images by convolving the reflectivity models with Ricker wavelets and
adding some extent of random noise or real noise from field seismic images. As the folding structures and
collapsed paleokarst features are well defined by functions in generating the reflectivity models, we were
able to automatically obtain the corresponding label images with ground-truth collapsed paleokarst fea-
tures. In addition, the parameters of functions are not fixed but randomly chosen, which makes it possible
to generate numerous synthetic training data sets with various folding and collapsed paleokarst features.
Although trained with only synthetic data sets, the CNN can be applied to characterize the collapsed pale-
okarst features in field seismic images that are not included in the training data sets. The application of
two field examples demonstrated that the collapsed paleokarst interpretation of our CNN-based methods
is consistent with manual interpretations and is significantly superior to conventional automatic methods.
From the CNN-based collapsed paleokarst interpretation results, we were able to automatically extract the
3-D boundary of each paleokarst chimney tube in the seismic images, which enabled us to automatically
measure the geometric parameters of the paleokarst tubes for quantitative analysis of the paleokarst system
development.
The features of collapsed paleokarsts formed in different areas can be highly variable. However, the current
trained CNN model may be limited to characterizing collapsed paleokarst systems in 3-D seismic images as
we mainly simulated the collapsed paleokarst features in our synthetic training data sets. Fortunately, the
current CNN model can be easily extended to delineate any type of paleokarst system by simply including
more training data sets that contain the corresponding paleokarst features.

Data Availability Statement

These two seismic images are available through Hardage et al. (1996) and Qi et al. (2014). The 120 pairs
of synthetic data sets, used for training and validating our CNN, are uploaded to Zenodo and are freely
available through the DOI link (https://doi.org/10.5281/zenodo.3690252). The codes related to this work is
be available through GitHub (https://github.com/xinwucwp/KarstSeg3D) or Zenodo (https://zenodo.org/
record/4013145).

Acknowledgments References
The authors would like to thank the
Al-Dossary, S., & Marfurt, K. J. (2006). 3D volumetric multispectral estimates of reflector curvature and rotation. Geophysics, 71, P41–P51.
Bureau of Economic Geology and
Alhakeem, A. A. (2013). 3D seismic data interpretation of Boonsville Field, Texas (PhD thesis), Missouri University of Science and
Marathon Oil, respectively, for kindly
Technology.
providing the first 3-D field seismic
Andre, P., & Doulcet, A. (1970). Rospo Mare Field, Italy, Apulian Platform, Adriatic Sea, Stratigraphic traps II: AAPG Treatise of
image shown in Figures 2 and 16 and
petroleum geology (pp. 29–54). Rospo Mare, Italy: Atlas of Oil and Gas Fields.
the second one shown in Figures 3a
Badrinarayanan, V., Kendall, A., & Cipolla, R. (2017). Segnet: A deep convolutional encoder-decoder architecture for image
and 13a.
segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 39, 2481–2495.
Bahorich, M., & Farmer, S. (1995). 3-D seismic discontinuity for faults and stratigraphic features: The coherence cube. The Leading Edge,
14, 1053–1058.
Bergen, K. J., Johnson, P. A., de Hoop, M. V., & Beroza, G. C. (2019). Machine learning for data-driven discovery in solid Earth
geoscience. Science, 363, eaau0323.
Bianco, M. J., Gerstoft, P., Traer, J., Ozanich, E., Roch, M. A., Gannot, S., & Deledalle, C.-A. (2019). Machine learning in acoustics: Theory
and applications. The Journal of the Acoustical Society of America, 146, 3590–3628.
Canter, K. L., Stearns, D. B., Geesaman, R. C., & Wilson, J. L. (1993). Paleostructural and related paleokarst controls on reservoir
development in the Lower Ordovician Ellenburger Group, Val Verde Basin, Paleokarst related hydrocarbon reservoirs: SEPM Core
Workshop (Vol. 18, pp. 61–99). Tulsa, OK: SEPM Society for Sedimentary Geology.
Chen, Y. (2016). Probing the subsurface karst features using time-frequency decomposition. Interpretation, 4, T533–T542.
Chen, L.-C., Papandreou, G., Kokkinos, I., Murphy, K., & Yuille, A. L. (2017). Deeplab: Semantic image segmentation with deep
convolutional nets, atrous convolution, and fully connected CRFS. IEEE Transactions on Pattern Analysis and Machine Intelligence, 40,
834–848.

WU ET AL. 21 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Chollet, F. (2015). Keras. https://github.com/fchollet/keras

Coogan, A., Bebout, D., & Maggio, C. (1972). Depositional environments and geologic history of Golden Lane and Poza Rica trend,
Mexico, an alternative view. AAPG Bulletin, 56, 1419–1447.
Desheng, L., Digang, L., Chengzao, J., Gang, W., Qizhi, W., & Dengfa, H. (1996). Hydrocarbon accumulations in the Tarim Basin, China.
AAPG bulletin, 80, 1587–1603.
Di, H., & AlRegib, G. (2020). A comparison of seismic saltbody interpretation via neural networks at sample and pattern levels.
Geophysical Prospecting, 68, 521–535.
Di, H., & Gao, D. (2016). Efficient volumetric extraction of most positive/negative curvature and flexure for fracture characterization from
3D, seismic data. Geophysical Prospecting, 64, 1454–1468.
Di, H., Shafiq, M. A., Wang, Z., & AlRegib, G. (2019). Improving seismic fault detection by super-attribute-based classification.
Interpretation, 7, SE251–SE267.
Di, H., Wang, Z., & AlRegib, G. (2018). Deep convolutional neural networks for seismic salt-body delineation. Presented at the AAPG
Annual Convention and Exhibition.
Dou, Q., Sun, Y., Sullivan, C., & Guo, H. (2011). Paleokarst system development in the San Andres Formation, Permian Basin, revealed by
seismic characterization. Journal of Applied Geophysics, 75, 379–389.
Geng, Z., Wu, X., Shi, Y., & Fomel, S. (2019). Relative geologic time estimation using a deep convolutional neural network. 88th Annual
International Meeting, SEG, Expanded Abstracts, 2238–2242.
George, M. C. (2016). The Muenster Uplift of North Texas: The Easternmost Expression of the Pennsylvanian Ancestral Rockies (PhD
thesis), The University of Texas at Dallas.
Guillen, P., Larrazabal, G., González, G., Boumber, D., & Vilalta, R. (2015). Supervised learning to detect salt body, SEG Technical
Program Expanded Abstracts 2015 (pp. 1826–1829). Tulsa, OK: Society of Exploration Geophysicists.
Hardage, B. A. (1996). Boonsville 3-D, data set. The Leading Edge, 15, 835–837.
Hardage, B., Carr, D., Lancaster, D., Simmons, J., Elphick, R., Pendleton, V., & Johns, R. (1996). 3-d seismic evidence of the effects of
carbonate karst collapse on overlying clastic stratigraphy and reservoir compartmentalization. Geophysics, 61, 1336–1350.
He, K., Zhang, X., Ren, S., & Sun, J. (2016). Deep residual learning for image recognition. In 2016 IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), Las Vegas, NV (pp. 770–778). https://doi.org/10.1109/CVPR.2016.90
Kerans, C. (1988). Karst-controlled reservoir heterogeneity in Ellenburger group carbonates of West Texas. AAPG bulletin, 72, 1160–1183.
Khatiwada, M., Keller, G. R., & Marfurt, K. J. (2013). A window into the Proterozoic: Integrating 3D seismic, gravity, and magnetic data to
image subbasement structures in the southeast Fort Worth Basin. Interpretation, 1, T125–T141.
Kingma, D. P., & Ba, J. (2014). Adam: A method for stochastic optimization. CoRR, https://arxiv.org/abs/1412.6980
Krizhevsky, A., Sutskever, I., & Hinton, G. E. (2012). Imagenet classification with deep convolutional neural networks. In F. Pereira,
C. J. C. Burges, L. Bottou, & K. Q. Weinberger (Eds.), Advances in Neural Information Processing Systems (pp. 1097–1105). New York:
Curran Associates, Inc. Retrieved from http://papers.nips.cc/paper/4824-imagenet-classification-with-deep-convolutional-neural
-networks.pdf
Li, F., & Lu, W. (2014). Coherence attribute at different spectral scales. Interpretation, 2, SA99–SA106.
Lin, T., Dollár, P., Girshick, R., He, K., Hariharan, B., & Belongie, S. (2017). Feature pyramid networks for object detection. In 2017 IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), Honolulu, HI (pp. 936–944). https://doi.org/10.1109/CVPR.2017.106
Liu, W., Anguelov, D., Erhan, D., Szegedy, C., Reed, S., Fu, C.-Y., & Berg, A. C. (2016). SSD: Single shot multibox detector, European
conference on computer vision (pp. 21–37). Springer.
Long, J., Shelhamer, E., & Darrell, T. (2015). Fully convolutional networks for semantic segmentation. In 2015 IEEE Conference on
Computer Vision and Pattern Recognition (CVPR), Boston, MA (pp. 3431–3440). https://doi.org/10.1109/CVPR.2015.7298965
Lorensen, W. E., & Cline, H. E. (1987). Marching cubes: A high resolution 3D surface construction algorithm. ACM Siggraph Computer
Graphics, 21, 163–169.
Loucks, R. G. (1999). Paleocave carbonate reservoirs: Origins, burial-depth modifications, spatial complexity, and reservoir implications.
AAPG bulletin, 83, 1795–1834.
Loucks, R. G., & Anderson, J. H. (1985). Depositional facies, diagenetic terranes, and porosity development in Lower Ordovician
Ellenburger Dolomite, Puckett Field, West Texas, Carbonate petroleum reservoirs (Vol. 18, pp. 19–38). Berlin: Springer-Verlag.
Lu, P., Morris, M., Brazell, S., Comiskey, C., & Xiao, Y. (2018). Using generative adversarial networks to improve deep-learning fault
interpretation networks. The Leading Edge, 37, 578–583.
Maoshan, C., Shifan, Z., Zhonghong, W., Hongying, Z., & Lei, L. (2011). Detecting carbonate-karst reservoirs using the directional
amplitude gradient difference technique. In SEG Technical Program Expanded Abstracts 2011 (pp. 1845–1849). Tulsa, OK: Society of
Exploration Geophysicists.
Marfurt, K. J., & Rich, J. (2010). Beyond curvaturevolumetric estimates of reflector rotation and convergence. In SEG Technical Program
Expanded Abstracts 2010 (pp. 1467–1472). Tulsa, OK: Society of Exploration Geophysicists.
Marfurt, K. J., Sudhaker, V., Gersztenkorn, A., Crawford, K. D., & Nissen, S. E. (1999). Coherency calculations in the presence of
structural dip. Geophysics, 64(1), 104–111.
Martin, D. R., Fowlkes, C. C., & Malik, J. (2004). Learning to detect natural image boundaries using local brightness, color, and texture
cues. IEEE Transactions on Pattern Analysis and Machine Intelligence, 26, 530–549.
Mazzullo, S. J., & Mazzullo, L. J. (1970). Paleokarst and karst-associated hydrocarbon reservoirs in the Fusselman Formation, West Texas,
Permian Basin. In Paleokarst, karst related diagenesis and reservoir development: Examples from ordoviciandevonian age strata of west
Texas and the mid-continent (Vol. 92-33, pp. 110–120). Tulsa, OK: Permian Basin Section, SEPM Publication.
McDonnell, A., Loucks, R. G., & Dooley, T. (2007). Quantifying the origin and geometry of circular sag structures in northern Fort Worth
Basin, Texas: Paleocave collapse, pull-apart fault systems, or hydrothermal alteration? AAPG Bulletin, 91, 1295–1318.
Pham, N., Fomel, S., & Dunlap, D. (2019). Automatic channel detection using deep learning. Interpretation, 7(3), SE43–SE50.
Pollastro, R. M., Jarvie, D. M., Hill, R. J., & Adams, C. W. (2007). Geologic framework of the Mississippian Barnett Shale,
Barnett-Paleozoic total petroleum system, Bend arch-Fort Worth Basin, Texas. AAPG Bulletin, 91, 405–436.
Provost, F. J., Fawcett, T., & Kohavi, R. (1998). The case against accuracy estimation for comparing induction algorithms. ICML, 445–453.
Qi, J., & Castagna, J. (2013). Application of a PCA fault-attribute and spectral decomposition in Barnett Shale fault detection. In SEG
Technical Program Expanded Abstracts 2013 (pp. 1421–1425). Tulsa, OK: Society of Exploration Geophysicists.
Qi, J., Zhang, B., Zhou, H., & Marfurt, K. (2014). Attribute expression of fault-controlled karst—Fort Worth Basin, Texas: A tutorial.
Interpretation, 2, SF91–SF110.
Redmon, J., Divvala, S., Girshick, R., & Farhadi, A. (2016). You only look once: Unified, real-time object detection. In 2016 IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), Las Vegas, NV (pp. 779–788). https://doi.org/10.1109/CVPR.2016.91

WU ET AL. 22 of 23
 Journal of Geophysical Research: Solid Earth 10.1029/2020JB019685

Ren, S., He, K., Girshick, R., & Sun, J. (2015). Faster R-CNN: Towards real-time object detection with region proposal networks. Advances
in Neural Information Processing Systems, 39, 91–99.
Roberts, A. (2001). Curvature attributes and their application to 3D interpreted horizons. First Break, 19, 85–100.
Ronneberger, O., Fischer, P., & Brox, T. (2015). U-net: Convolutional networks for biomedical image segmentation. In International
Conference on Medical image computing and computer-assisted intervention (pp. 234–241). New York: Springer.
Sermanet, P., Eigen, D., Zhang, X., Mathieu, M., Fergus, R., & LeCun, Y. (2013). Overfeat: Integrated recognition, localization and
detection using convolutional networks. arXiv preprint arXiv, 1312, 6229.
Shi, Y., Wu, X., & Fomel, S. (2019). Saltseg: Automatic 3D salt segmentation using a deep convolutional neural network. Interpretation,
7(3), SE113–SE122.
Simonyan, K., & Zisserman, A. (2014). Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv, 1409,
1556.
Soriano, M. A., Pocoví, A., Gil, H., Perez, A., Luzón, A., & Marazuela, M. Á. (2019). Some evolutionary patterns of palaeokarst developed
in Pleistocene deposits (Ebro Basin, NE Spain): Improving geohazard awareness in present-day karst. Geological Journal, 54, 333–350.
Sullivan, E. C., Marfurt, K. J., Lacazette, A., & Ammerman, M. (2006). Application of new seismic attributes to collapse chimneys in the
Fort Worth Basin. Geophysics, 71, B111–B119.
Szegedy, C., Liu, W., Jia, Y., Sermanet, P., Reed, S., Anguelov, D., et al. (2015). Going deeper with convolutions. In 2015 IEEE Conference
on Computer Vision and Pattern Recognition (CVPR), Boston, MA (pp. 1–9). https://doi.org/10.1109/CVPR.2015.7298594
Viniegra, F., & Castillo-Tejero, C. (1970). Golden Lane Fields, Veracruz, Mexico. In Geology of giant petroleum fields (Vol. 14, pp. 309–325).
Tulsa, OK: AAPG Memoir.
Walper, J. L. (1982). Plate tectonic evolution of the Fort Worth Basin. In Petroleum geology of the Fort Worth Basin and Bend arch area
(pp. 237–251). Dallas, TX: Dallas Geological Society.
Wen, T., Castro, M. C., Nicot, J.-P., Hall, C. M., Larson, T., Mickler, P., & Darvari, R. (2016). Methane sources and migration mechanisms
in shallow groundwaters in Parker and Hood Counties, Texas—A heavy noble gas analysis. Environmental science & technology, 50,
12,012–12,021.
Wu, X., Geng, Z., Shi, Y., Pham, N., Fomel, S., & Caumon, G. (2020). Building realistic structure models to train convolutional neural
networks for seismic structural interpretation. Geophysics, 85, WA27–WA39.
Wu, X., Liang, L., Shi, Y., & Fomel, S. (2019). FaultSeg3D: Using synthetic datasets to train an end-to-end convolutional neural network
for 3D seismic fault segmentation. Geophysics, 84(3), IM35–IM45.
Wu, X., Shi, Y., Fomel, S., Liang, L., Zhang, Q., & Yusifov, A. Z. (2019). Faultnet3D,: Predicting fault probabilities, strikes, and dips with a
single convolutional neural network. IEEE Transactions on Geoscience and Remote Sensing, 57, 9138–9155.
Wu, H., Zhang, B., Lin, T., Cao, D., & Lou, Y. (2019). Semiautomated seismic horizon interpretation using the encoder-decoder
convolutional neural network. Geophysics, 84, B403–B417.
Zeng, H., Wang, G., Janson, X., Loucks, R., Xia, Y., Xu, L., & Yuan, B. (2011). Characterizing seismic bright spots in deeply buried,
Ordovician Paleokarst strata, Central Tabei uplift, Tarim Basin, Western China. Geophysics, 76, B127–B137.
Zhai, G., & Zha, Q. (1982). Buried-hill oil and gas pools in the North China Basin, in The deliberate search for the subtle trap. AAPG
Memoir, 32, 317–335.
Zhao, T., Li, F., & Marfurt, K. J. (2018). Seismic attribute selection for unsupervised seismic facies analysis using user-guided
data-adaptive weights. Geophysics, 83, O31–O44.
Zhao, T., & Mukhopadhyay, P. (2018). A fault detection workflow using deep learning and image processing. In SEG Technical Program
Expanded Abstracts 2018 (pp. 1966–1970). Tulsa, OK: Society of Exploration Geophysicists.

WU ET AL. 23 of 23