Elastic and shear moduli of coal measure rocks derived from basic well logs

using fractal statistics and radial basis functions

* C. O ̈ zgen Karacan

* Tel.: +1412 386 4008; fax: +1412 386 6595.

E-mail address: cok6@cdc.gov

CDC/NIOSH, Pittsburgh Research Laboratory, Disaster Prevention and Response Branch, 626 Cochrans Mill Road, PO Box 18070, Pittsburgh, PA 15236, USA

abstract

Gamma ray, density, sonic and core logs obtained from two boreholes drilled over a longwall panel in

Southwestern (SW) Pennsylvania were analyzed for formation boundaries, log-derived porosities and

densities and for rock elastic properties from sonic transit times. Gamma ray (GR) and density logs (DL)

were analyzed using univariate statistical techniques and fractal statistics for similarity and ordering of

the log data in depth. A Fourier transformation with low-pass filter was used as a noise elimination

(filtering) technique from the original logs. Filtered data was tested using basic univariate and fractal

statistics, rescaled range (R/S) and power spectrum (PS) analysis to compare the information

characteristics of the filtered logs with the original data. The randomness of log data in depth was

analyzed for fractional Gaussian noise (fGn) or fractional Brownian motion (fBm) character.

A new prediction technique using radial basis function (RBF) networks was developed to calculate

shear and Young’s moduli of the formations when sonic logs are not available. For this approach, the

filtered logs were used as input to an RBF based upon a combination of supervised and unsupervised

learning. The network was trained and tested using rock elastic properties calculated from the sonic log

of one of the boreholes. The network was used to predict the elastic and shear moduli of the coalmeasure rocks over a longwall coal mine in SW Pennsylvania. This approach demonstrated that it could

be used for prediction of elastic and shear moduli of coal-measure rocks with reasonable accuracy.

1. Introduction

The properties of coal-measure rocks within the longwall

overburden are important because of their controlling effect on

fluid storage and flow before and after coal extraction. Methane

inflow into the mines from fractured strata during longwall

mining and production potential of the methane degasification

boreholes drilled from surface (gob gas ventholes or GGVs) to

capture these emissions are influenced by reservoir, elastic and

strength characteristics of the overlying strata. Thus, the ability to

accurately determine reservoir and elastic properties of the coal

measure rocks is extremely important for ground control and

methane control objectives in underground coal mining. The

elastic modulus, or Young’s modulus, and shear modulus are used

when deformations in underground mines need to be computed.

Experimentally, the elastic modulus can be determined from the

stress–strain response of a rock sample subjected to uniaxial

compression. However, this is a complicated and time consuming

experiment.

With the advent of new computing and experimental techniques, complicated and time consuming experiments and methods

to determine rock elastic and strength properties are being

replaced by quicker, and may be even more accurate, techniques.

For instance, the uniaxial unconfined compressive strength (UCS)

test has been largely replaced by simpler, faster and cheaper

‘‘indirect’’ tests such as point loading method [1] and Schmidt

hammer [2], although UCS was the main direct quantitative

method for rock strength determination for many years. Triaxial

testing of a core from a rock material to find rock elastic and

strength properties is another difficult, expensive and complicated

laboratory test, in which the stress state in the rock sample is

axisymmetric and defined by the axial stress, confining pressure

and pore pressure. The results of this test are usually studied in

terms of Terzaghi’s effective mean pressure and deviatoric stress

which correspond to volume strain and the deviatoric strain [3].

In order to overcome the difficulties related to determination

of elastic properties of rocks by laboratory tests, it has become

very popular to develop alternatives that would predict the elastic

moduli of rocks using theoretical and empirical approaches. In

these approaches, elastic and shear moduli are either related to

simple physical properties like total porosity or to other mechanical properties, indices, and even to mineral composition of rocks

that are determined in the laboratory or documented in the
literature. However, as shown by the equations and conversion

factors in the literature, it is clear that there are no unique factors

and relations which can be used, even for the same type of rock.

Also, experimental data available in the literature that sometimes

is used to generate predictive equations mainly concerns limestones and cemented sandstones. Poorly cemented sandstones are

difficult to characterize due to coring and heterogeneity problems.

Low permeability of shales and their sensitivity to water

complicate the tests by showing a strong interaction of the strain

rate and the saturating fluid. Thus, literature data on shales needs

to be handled carefully [3] for generating conversion factors and

relations.

Palchik [4] reported the relationships between Young’s

modulus (E) and sc are obtained from the stress–strain responses

of six carbonate rocks (dolomites, chalks and limestones)

collected from different locations in Israel. He reported that the

use of these empirical equations is not universal, but rather

limited, since empirical constants can only be used for the

carbonate rocks with which the tests were conducted. In addition,

he noted that although E generally increases with increasing sc,

these relations have different degrees of reliability with R2 ranging

from 0.57 to 0.98. Thus, he proposed a stress–strain model based

on Haldane’s distribution function to predict the elastic modulus

of intact carbonates. He concluded that this model is better

correlated with the experimental data for Israeli carbonate rocks

having a UCS o100 MPa.

Deriving the elastic properties of a rock from its mineralogical

content has been a tempting alternative approach [3]. The idea

behind this approach is that the mineralogical characteristics of

rocks influence the macromechanical properties and especially the

rock strength. Techniques such as variational methods [5] and the

theory of homogenization [6] are available for this objective.

However, these methods are far from practical and require powerful

numerical procedures. Thus, approaches based on a restrained

number of parameters accessible from drilling logs, which give

reasonable estimations of rock mechanical properties, were sought.

Mineral composition, grain size and microstructure are known

to be the most important parameters that influence rock strength

and elastic properties. Bemer et al. [3] also stated that it is difficult

to address sandstones since their composition is extremely

variable and their elastic properties depend on the stress state

due to their often nonlinear behavior. Sabatakikis et al. [7]

conducted a large number of laboratory tests on intact marlstones,

limestones and sandstones. Index properties such as porosity, unit

weight, mineralogical content, Schmidt hammer and point

loading indices, as well as the strength under uniaxial and triaxial

compression were determined. From the analyses of the data,

conversion factors relating rock strength to rock compositional

properties and index properties were calculated. They reported

that textural characteristics in limestones are more important

than mineral composition for rock strength. In sandstones, the

mineral composition was found to be more important.

Further attempts to determine UCS, shear and Young’s moduli

included data driven approaches. These used material properties

to generalize the relations between input and output space for

nonlinear approximations. Grima [8] used fuzzy modeling and

artificial neural networks (ANN) for the modeling of unconfined

compressive strength of rock samples. He compared results with

those of statistical models and the empirical relationship of

Verwaal and Mulder [9]. Density, porosity and equotip hardness

test (EHT) values were used as the input parameters for his

models. He reported that data driven models are more reliable

than statistical and empirical models. Sonmez et al. [10]

developed an artificial neural network-based chart which considers sc and unit weight as input to predict elastic modulus of

intact rock. Maji and Sitharam [11], on the other hand, applied

back propagation algorithm and radial basis function to predict

elastic modulus of jointed rocks from the elastic modulus of intact

rocks, confining pressure, joint frequency, joint inclination and

joint roughness. They concluded that these models could

accurately represent the effect of confining pressure and joint

properties.

One of the major disadvantages inherent in all laboratory tests,

and in the derived correlations is that the tested samples will not

be in the in situ conditions anymore when they are tested in the

laboratory. It is very difficult or nearly impossible to create exactly

the same underground stress and fluid saturation conditions in

the laboratory due to the many interacting parameters under in

situ conditions. In addition, the scaling effect in the laboratory

measurements due to sample size and the testing effects at the

‘‘end points’’ of the sample may be an issue to interpret the

results. Furthermore, the elastic moduli can be either dynamic or

static, and undrained or drained. Standard laboratory tests usually

give access to static and drained elastic moduli, which can then be

directly used to compute strains. However, dynamic undrained

tests can be conducted either ‘‘in situ’’ or with complex and time

consuming methods.

Geophysical well logs and seismic, such as time-lapse seismic,

techniques can measure the in situ properties of rocks and fluid

reservoirs. For instance, a well log interrogates the formations

around the wellbore with different techniques and obtains

response data as the tool travels downward or upward in a

borehole. Among the well log techniques that are most frequently

used by oil and gas companies are spontaneous potential (SP),

gamma ray (GR), density log (DL) and sonic log (SL) that calculate

basic formation properties such as porosity, density, shale content

and boundaries of the formations of interest. Sonic logs that

employ compression and shear waves to interrogate the formations can also be used as the source of data for an ‘‘in situ’’

geotechnical evaluation of the rock formations since acoustic

travel time is explicitly tied to the density and the elasticity of the

medium [12]. They also can give access to dynamic undrained

elastic moduli, which require data on the compressibility of the

saturating fluids to be turned into dynamic drained elastic moduli.

Gamma ray and density logs are two of the conventional

logging techniques frequently employed by mining companies to

determine formation boundaries, formation thicknesses and

formation types. On the other hand, full wave sonic logs are

seldom used since they are more expensive to run, data intensive

and more complicated in terms of operation procedures. However,

SL is capable of providing more detailed information about the

condition of the borehole, surrounding strata, their elastic

properties and how they may behave during underground mining.

This information is particularly important for ground control and

for effective control of methane from a fractured zone using gob

gas ventholes (GGV).

The objective of this study is to develop a new technique to

calculate shear and elastic (Young’s) moduli of the formation,

which can then be used to determine bulk modulus and Poisson’s

ratio, using radial basis function (RBF) networks. The method is

based on processing of GR and DL by Fourier transform, fractal

statistics and then modeling using RBF to predict elastic and shear

moduli, when full wave sonic logs are not available or too costly to

obtain. This approach evaluates the ‘‘in situ’’ rock elastic properties of coal-measure rocks and for prediction purposes with a

reasonable accuracy.

2. Borehole data evaluated in this study

Data from two exploration boreholes (EBH-1 and EBH-2) were

used in this study. These two boreholes were in the same mining
area but were 8 km apart from each other and were located

south and north of Waynesburg, Pennsylvania (Fig. 1).

Interval 

EBH-1 

EBH-2 

~8 km 

N 

S W 

E 

thicknesses 

~70.0m 

~25.0 m 

~12.0m 

Major layers with 

thicknesses 

Sandstone-~18 m 

Waynesburg-~1m 

Limestone-~30 m 

Sandstone-~30 m 

Shaly Sandstone-~6 m 

Sewickley-~1.5 m 

Shale-~1m 

Limestone-~8 m 

Shale-~10 m 

Shale+Sandstone-~3.0 m 

Sandstone Chan. ~0/12 m 

Shale-~0/12 m 

Pittsburgh-2 m 

Fig. 1. Exploration boreholes (EBH) studied in this paper, their locations in the field and a generalized stratigraphic column of this area.

EBH-1 and EBH-2 were drilled for exploration purposes to

characterize the mine roof and the mining coal thickness at the

two locations. They were drilled from surface to an approximate

depth of 12.2 m with rotary techniques and continued with a

5.8 cm core drill until the bottom of the Pittsburgh coal bed was

exceeded by about 6.1 m. Cores were fully recovered from EBH-1

and were marked for depths and lithology identification. The total

depths of EBH-1 and EBH-2 were 255.1 and 263 m in depth,

respectively. After completion of drilling, boreholes were logged

with gamma and density tools for the entire length of the

borehole. In addition, EBH-2 also was logged with a full wave

sonic tool.

In this area of the Northern Appalachian basin, the overburden

depths above Pittsburgh coal seam ranges between 152 and

274 m. A generalized stratigraphic section of the strata above the

Pittsburgh coal bed in the study area is shown in Fig. 1. At some

locations in this area, there is a sandstone paleo-channel complex

and associated shale unit overlying the Pittsburgh seam.

The Sewickley coal bed and the rider coals (not shown in this

figure) directly above the Pittsburgh coal bed are believed to be

the primary sources of gob gas during longwall mining. Gas

released from the Pittsburgh rider coals located in the caved zone

is expected to migrate to the mine ventilation system, while gas

from the Sewickley coalbed and from any other gas-bearing

horizons above the caved zone migrates to the pressure sink of the

operating gob gas ventholes in the fractured zone.

3. Characterization of coal-measure rocks using well logs

3.1. Characterization of rock formations using gamma ray and

density logs

In this section, gamma ray and density logs obtained from

EBH-1 and EBH-2 were evaluated with core logs to characterize

the coal-measure rocks above Pittsburgh coal bed. Fig. 2 shows

the lithological log, thicknesses of major layers defined from the

driller’s log, and raw gamma and density logs from EBH-2. At EBH2, the top of the Pittsburgh coal bed was at a depth of 252 m. The

general sequence between 204 and 241 m above the Pittsburgh

seam layer can be divided in four sections (dotted lines in Fig. 2).

These zones are mainly alternating sequences of limestone, shale,

and sandstone, as determined from gamma ray and density logs.

In each of these zones, thin shale layers (A, C–F) are sandwiched

between stronger limestone or sandstone formations. In this log,

the Sewickley coalbed is located at ‘‘B’’ (low gamma ray, low

density) and overlain by carbonaceous shale (A) which records a

high gamma ray reading accompanied by a decreased density.

Raw gamma ray and density logs, as well as lithological logs

and strata thicknesses, for EBH-1 are shown in Fig. 3. At this

location, the Pittsburgh coal bed is found at a depth of 254 m.

Similar sequence of layers, as in EBH-2, is observed at this location

too. Limestone layers are separated by shale layers as annotated

on gamma and density logs (A–G). Some of the low density

readings (A/ and B/

) are associated with carbonaceous and limy

shales and with highly porous shale units. At this borehole site,

the Sewickley coal seam is thicker than it is at EBH-2 location. This

might influence the amount of methane that will be experienced

once this location is undermined. Also, this location is missing the

thick Pittsburgh sandstone that overlays the Pittsburgh coal seam.

The absence of this layer may affect the fracturing and the amount

of methane associated with this layer.

3.2. Shale content and porosity of coal-measure formations at EBH-1

and EBH-2 using GR and DL

The GR log is a record of a formation’s radioactivity. During

deposition of coal measure shales, its clay constituents adsorb

ions of the heavy radioactive elements from mineralized waters

involved in the weathering process of igneous formations. Due to

adsorption of radioactive elements, fine grained sediments such

as clayey sands and shales are generally quite radioactive. On the

other hand, carbonate rocks have low radioactivity since they are

mainly derived from calcareous marine-life skeletal material.

Similarly, pure sandstones have low radioactivity due to the wellordered structure of quartz that exclude impurities during

crystallization. Since sandstones are usually porous, they may be

penetrated with radioactive ions and clays later [13].

The GR log is usually used to identify boundaries, primarily

shale units from other low radioactivity formations, and to

quantify shale volume of the rock formations. In fractured

formations, an increase in the gamma ray reading without

concurrently higher formation shaliness can be observed. This

increase has been explained by deposition of uranium salts along

the discontinuity surfaces of a fracture or within the crack itself.

Also, it is not uncommon to find as much as 0.01% uranium or

thorium in dark bituminous shale units, which increases the

gamma ray reading [13].

In this study, formation type and thicknesses were first

identified from the GR logs of the boreholes for calculating shale

volume. The clean-formation readings for each rock type were

determined by locating the clean sandstone and limestone units

in the logging data. Pittsburgh sandstone and limestone intervals

were suitable for this purpose since they were free from shales.

The GR reading for coal was determined based on the average GR

reading along Sewickley layer. For the shale GR reading, the

average reading of the shale intervals was determined and used in
the calculations. In shale volume calculations, each interval was

calculated separately using respective clean-formation readings.

The aim was to minimize errors from using a single formation GR

value for the entire borehole interval.

Strata Thickness (m) Gamma 

0 2 4 6 8 0 50 100 150 

200 

205 

210 

215 

220 

225 

230 

235 

240 

245 

Depth (m) 

200 

205 

210 

215 

220 

225 

230 

235 

240 

245 

Depth (m) 

Density 

Gamma Ray (cps) 

Density (g/cc) 

Pittsburgh 

sandstone 

A 

B 

C 

D 

E 

F 

0 1 2 3 

Fig. 2. Lithological log of EBH-2 with thicknesses of major layers defined from driller’s log and raw gamma and density logs.

Porosities were calculated based on the readings of the density

log which corresponds to bulk density. There are a variety of

factors affecting bulk density measurements. One of the most

important is the shaliness of the formation which affects the

measurement by the amount of its contribution to the total signal.

In this study, a correction was made by the volume of shale in

order to obtain a shale-corrected density measurement [13].

The following equation can be used to quantify the shale

volume by using GR log [14]:

GR GR Vshale

- 1⁄4 l clean rock (1) GRshale - GRclean rock

In this equation, GRl is the gamma ray reading (cps, count per

second) from the log. The other two terms are the GR readings of

the clean formations (clean sandstone and limestone) and the

pure shale [13,14]. These readings were obtained from the

formations and using the methodology mentioned in the previous

paragraph.

Bulk density (rb), on the other hand, is determined from

density log, from which porosity can be calculated (Eqs. (2) and

(3)). However, before porosity can be determined, the lithology of

the formation, the matrix density, and the fluid density filling the

pore space must be known. For a clean formation of known matrix

density, rma, with a fluid of average density, rf, the linear sum of

the contributions can be used to calculate porosity (f) as

described in Karacan [15] and by using below equations:

rb 1⁄4 frf þ ð1 - fÞrma (2)

rb 1⁄4 rbsh:free ð1 - VshaleÞ þ rshaleVshale (3)

Fig. 4 shows the shale content and porosity values for EBH-1.

The calculations show that the lowest shale percentage is found in

limestone bearing units at depths of 225.5 m, 230.1–234.7 m, with

an average shale amount of 5–10%. Within this interval, there are

layers of weak shale formations at 229.9 and 233.5m. These and

similar layers are potentially weak interfaces to be affected by

mining stresses and by the resultant bedding plane separations.

Outside of limestone intervals, the average shale content in the

strata varies from 40% to 50%.

Fig. 4 also shows the porosity values calculated in the EBH-1

borehole using the density log. This plot shows that porosity

values are generally low (0–0.1), except for some intervals where

values exceeding 0.5 are observed. These locations are at 210.3 m

(A), 217.9 m (B), 224.0–225.5 m (C–D), 233.1 m (E), 237.8 m (F) and

243.8 m (G), where 224 m corresponds to the Sewickley coal bed.

These intervals are generally associated with fractures, laminated

layers and shales inter-layered between limestone units. These

high porosity intervals can result in weak spots in the strata

leading to eventual fracturing and separation, which may also
promote methane flow towards any borehole or to the mine.

These areas may potentially be the reservoirs of free methane that

will flow into a pressure sink once the fracture connection is

established during mining.

200

205

210

215

220

225

230

235

240

245

250

255

0

Gamma 

Depth (m)

0

Density 

Gamma ray (cps)

Density (g/cc)

200

205

210

215

220

225

230

235

240

245

250

255

0

Strata Thickness (m)

Depth (m)

A 

B 

C 

D 

E 

F 

G 

Al

Bl

0.5 1 1.5 2 2.5 3

2 4 6 8 10 12 50 100 150

Fig. 3. Raw GR, DL, lithological log and strata thicknesses for EBH-1.

Similar observations can also be made for EBH-2 location

based on shale content and porosity calculations, which are

shown in Fig. 5. In this borehole, average shale content in the

sandy-shale layers is around 30–40%. Higher shale-content areas

are interbedded with more competent limestone sequences and

the Pittsburgh sandstone 217.9–228.6 m and below 236.2 m,

respectively. The layer that is almost exclusively shale at around

211.8 m is closely associated with thin carbonaceous shale and is

either a clay-filled fracture or a fracture surface that is carrying a

high amount of radioactive ions.

The calculated porosities for EBH-2, based on the density log,

are also shown in Fig. 5. The values indicate that most porosity

values are around 0.1. However, there are higher porosity sections

based on the strata intervals. The highest values are marked from

A to H on the figure and they are generally associated with shale

layers, especially where clay content is high and possibly where

some natural fractures within the strata or along the bedding

interfaces are present. In these sections, porosities were as high as

0.4–0.5 and in some instances even more. These high porosities

can be related to weak spots in the strata that may lead to

fracturing and bedding plane separation. They can also constitute

pockets for free methane and can serve as flow paths within the

formation.

3.3. Shear modulus and Young’s modulus calculation at EBH-2 using

sonic log

In sonic logging, ultrasonic frequencies are employed in the

form of compressional and shear waves. Compressional wave

measurement used in sonic logging exhibits longitudinal particle

motion and can be propagated in solids, liquids and gases. The
shear wave, on the other hand, is a transverse wave, in which the

direction of propagation is perpendicular to the direction of

particle displacement. In case of a compressional wave, the

acoustic wave alternately compresses the surrounding medium on

a forward movement and rarifies it on a backward movement [13].

220 

245 

Shale (%) Porosity 

0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 

200 200 

205 205

210 A 210 

215 215 Depth (m) 

220 

C 

225 D

B 

225

Depth (m)

230 230

E 

235 235

F 

240 240

G 

245 

250 250

255 255

Fig. 4. Shale content and porosity log calculated from raw gamma and density logs

obtained from EBH-1.

Shale (%) Porosity 

0 20 40 60 80 100 0.0 0.2 0.4 0.6 

200 200

205 205

A 

210 210

B

215 215

C 

D 

F

Depth (m)

Depth (m)

220 220

E 

225 225

230 230

G 

235 235

H 

240 240

245 245

Fig. 5. Shale content and porosity log(s) calculated from raw gamma and density

logs obtained from EBH-2.

Velocities obtained from sonic logs can be used as the

fundamental information for rock classification and geotechnical

evaluation. Acoustic travel time of any type of ultrasonic

disturbance is explicitly tied to the density and the elasticity of

the medium. Compressional waves are characterized by first

arrival times. Although slower than compressional waves (0.5–0.7

times), shear waves are usually stronger and can be identified in

the total wave train.

In this paper, shear modulus and Young’s modulus values were

calculated using the velocities of compressional and shear waves

determined from full wave sonic logs as described in [12].

Dynamic moduli were calculated using both compressional and

shear waves and a determination of dynamic Poisson’s ratio (n)

using Eqs. (4)–(6):

G 1⁄4 rV2

s (4)

n 1⁄4 ( V 2 2 p

),( V p 2 2 2 s)

(5) Vs V 

where r is the density, G is shear modulus, Vp and Vs are the

compressional and the shear wave velocities, respectively. The

Young’s modulus (E) and the bulk modulus (K) can then be

calculated using

E 1⁄4 2Gð1 þ nÞ (6)

E K 1⁄4 (7) 3ð1 2nÞ 

Young’s moduli and shear moduli for EBH-2 are given in Fig. 6.

This figure shows that highest Young’s and shear moduli are

associated with limestone and competent sandstone layers. The

values in these layers are as high as 4 and 9 GPa for shear modulus

and Young’s modulus, respectively. When the layer is shale, sandy

or limy shale, or coal, values decrease abruptly to values as low as

1–2 GPa, indicating weaker rock units that will be deformed easily

when subjected to high stress and strain conditions prevailing

during mining. These weak layers and their interfaces with

stronger rocks are candidates for fracturing and bedding plane

separations to form increased permeability pathways for methane

migration into longwall gobs.

Despite the crucial information obtained that is related to

possible rock behavior under certain stress and strain conditions,

sonic logging is more unconventional than both density and

gamma ray logging, especially in mining-related strata characterization. Part of the reason why it is not used as often as the other

two is the cost and time associated with recording and processing

sonic data. Processing sonic full wave profiles can also be a tedious

task since the arrival times of compressional and shear waves

should be determined from an entire wave train to calculate

Young’s modulus, shear modulus and sonic porosity.
200

205

210

215

220

225

230

235

240

245

0

Modulus (GPa)

Depth (m)

Shear M

Young M

2 4 6 8 10

Fig. 6. Shear and Young’s moduli of the formations calculated based on

compression and shear-wave sonic velocities at EBH-2 drill site.

4. Methodology development for log-data treatment and for

generating rock elastic properties from gamma and density logs

The aim of this section is to gain a better understanding of the

data ordering in available logs from EBH-1 and EBM-2 using

rescaled range (R/S) and power spectrum (PS) analyses to

determine their fractal (fGn/fBm) character. This information is

used to compare the original logs with the filtered ones generated

using Fourier transform (FT) to make sure that both sets of data

carry the same or similar information content and data ordering.

Next, a methodology to generate shear and Young’s moduli using

filtered gamma and density records is developed. The method is

based on using filtered GR and DL with radial basis function

networks to predict available Young’s and shear moduli in the

same borehole. An optimized network system is then used to

predict the elastic and shear moduli for EBH-1 using its gamma

and density readings, which unfortunately did not have a sonic log

record.

The flowchart in Fig. 7 graphically presents the details of

the methodology followed for this approach. The details of the

techniques will be given in the next sections and in the

appendices in this paper.

Testing/Verification

on EBH-2 

Prediction for EBH-1

Young’s Modulus 

Shear Modulus 

Fourier Transform 

(Low-Pass Filter) 

Inversion

 Radial Basis Function NN 

 Training (Optimization) 

 Cross Validation 

 Testing 

EBH-2 EBH-1

Gamma Ray 

(Original)

Density Log 

(Original)

Gamma Ray 

Density Log 

Young’s Mod 

Shear Mod. 

Fourier Transform 

(Low-Pass Filter) 

Inversion

Filtered

DL

Filtered

GR

Comparison with original logs 

Fractal statistics 

Univariate statistics 

Rescaled Range 

(R/S) Analysis

Power Spectrum 

(PS) Analysis

Comparison

satisfactory ? 

Yes

No

Fail Pass

Filtered

DL

Filtered

GR

Fig. 7. A flowchart representation of the details of the mathematical methodology

followed for the log-data analysis and generation of radial basis function network.

4.1. Stationarity of self-affine fractal character of log traces

In analyzing the self-affine fractal data sequences and the

ordering of data within, stationarity and ergodicity are important

characteristics to be determined. A stationary data ensemble is

one whose ensemble average properties are independent of depth

or time. Thus, if the logs show an increasing (or decreasing) trend

with depth, even if the ensemble averages are the same, the

ensemble is non-stationary and, therefore, it is also non-ergodic

[16,17]. These two properties usually determine which mathematical techniques, fractional Gaussian noise (fGn) or fractional

Brownian motion (fBm), are appropriate for an analysis.

Fractional Gaussian noise and fractional Brownian motion are

data sequences with self-affine fractal character that can be used

to describe time series and geological distributions as well as

heterogeneities. Self-affine fractal sequences are statistical in

character and mathematically described by using ‘‘Hurst coefficient’’, which is a dimensionless number between 0 and 1 and

represents the intermittency of the data. Unlike pure Gaussian

(white noise, or uncorrelated data), fGn exhibits spatial correlation between data points. fBm is the integral of fGn and contains

more low frequency components than high frequencies. The

characteristics of fGn and fBm are mathematically described in

Appendices A and B.

Histograms of fGn- and fBm-type data sequences present

differences, too. Histogram of fGn is Gaussian where the

histogram of fBm is box-shaped, or rectangular. In this study

Gaussians, which represented the data distributions with accuracies of 0.94 and 0.96, were fitted to those data in the depth

interval of interest (204.2–240.8 m) and a good representation of

histograms with these curves were observed. The Gaussian shape

of the histograms suggested that the data series could be

classified as fGn. Also, DL and GR of EBH-2 (Fig. 2) were locally

increasing or decreasing based on the formations within this

interval. However, there was not a consistent trend in GR and DL
within this interval, suggesting that both curves are stationary. A

similar stationarity property could be proposed for Young’s and

shear moduli by examining the sonic-log derived data shown in

Fig. 6. They both could be represented almost perfectly by

Gaussian distribution with an accuracy of 0.98 (Fig. 8).

0

0.81.62.43.2 

4

4.85.66.4

More 

0 

1 2 3 

4 5 

6 7 8 9

10 

0 

20 

40 

60 

80 

100 

120 

140 

160 

180 

200 

Frequency 

0 

20 

40 

60 

80 

100 

120 

Frequency 

Shear Modulus (GPa) Young's Modulus (GPa) 

Fig. 8. Histogram for shear and Young moduli data calculated from SL of EBH-2 and their Gaussian representation.

4.2. Noise reduction in log data (EBH-2) by Fourier transforms

method

Random noise (RN) can greatly affect the result of certain types

of analyses by shifting the location of the peak, which leads to

systematic overestimation of the height of the peak and to

amplification after differentiation. RN can also decrease the

prediction capability of any predictive methods by creating

additional nonlinearities and local minimums that may distort

the approach to convergence in iterative techniques. Since random

noise is distributed over all frequencies and the real signal is

typically limited to low frequencies, a reduction of high frequency

components improves the signal/noise (S/N) ratio. Before proceeding further in the elastic and shear moduli generation using

GR and DL in RBF, the original logs were filtered to remove the

high frequency signals, which correspond to possible noise. The

purpose of this step is to increase S/N and thus improve the

prediction capability of RBF on the true ‘‘signal’’.

There are different methods available to generate a smoother

data trace. For instance, a smoother log signal can be generated

from the original curve by averaging every ‘‘n’’ values of the curve

and plotting these averages at the middle of each of the averaged

intervals. However, this process significantly reduces the number

of data points, and thus the information content. Running average

is another option but it creates large lags at the ends of the data

sequence.

For filtering the data sequence, the aim is to generate a curve

with improved S/N that honors the properties of the original curve

and that is statistically similar. This can be done by generating a

curve with the same (or close) mean, standard deviation, the same

average values and self-affine fractal properties investigated by

different methods as the original curve. In this work, the objective

of ‘‘non-destructive’’ filtering on well logs was achieved by

employing a low-pass filter using Fourier transform methods.

Well logs are processed as discrete and finite signals and

transformed into frequencies for filtering of high-frequency

signals, while trying to preserve as much of the original

information content as possible.

During FT, selected parts of the frequency spectrum ‘‘G’’ can

easily be subjected to piecewise mathematical manipulations

based on the frequency of the signal. These manipulations result

in a modified or ‘‘filtered’’ spectrum. Therefore, signal smoothing

can be performed by completely removing higher frequency

components while information-bearing low frequency components are retained. By eliminating undesirable high-frequency

components, we can filter a data array to remove noise and other

artifacts. Taking an inverse Fourier transform converts the data

back into the depth domain [16]:

DoXN

2p g i j 1⁄4 Gme ðj-1Þðm-1Þ=N (8) 2p m1⁄41

where o is the angular frequency of a sine or cosine function, and

f 1⁄4 o/2p.

Fig. 9 compares the original GR and DL of EBH-2 with the curve

generated by averaging the data every 0.61 m and with the filtered

curve generated using the low-pass filter in Fourier transform

methods. This figure shows that the filtered curve better

represents the original curve in both cases without destroying

the overall information in high and low frequency regions.

The success of FT filtering techniques is further investigated by

comparing the descriptive statistics of the original log data with

the filtered one to see if similar mean, standard deviation and

median can be captured. The Hurst exponents and global scaling

parameters obtained from the two sets to characterize the fractal

nature of the original and filtered data using rescaled range and

power spectrum were also studied in the following sections.

4.3. Univariate measures to investigate the statistical similarity of

filtered and original logs of EBH-2

The filtered logs were investigated using univariate and fractal

statistics to see if they honored the properties of the original

curves. Comparing the histograms obtained from filtered GR and

DL logs to the original histograms showed that FT preserved the

Gaussian properties of the original logs with accuracies of 0.93

and 0.98 for GR and DL, respectively. A similar filtering and

comparative study for shear and Young’s moduli logs were also

performed and proved that these data could also be filtered

without affecting the character of the original moduli curves.

In addition to Gaussian checks of original and filtered log data,

statistical parameters of original and filtered sequences of GR, DL,

shear modulus and Young’s modulus from EBH-2 were compared.

This approach checked if the filtered data sequences honored the

basic statistical features of the original data sequence. Table 1

shows the results of these analyses. Values of mean, standard

error, median and standard deviation are very similar for the

original and filtered logs of EBH-2. These results encourage

the use of filtered data instead of noisier, original data in the

subsequent analysis and proxy model development for prediction

of elastic properties of rock strata for EBH-2 and EBH-1.

4.4. Self-affine fractal analysis on logs of EBH-2: power spectrum

and rescaled range analysis for determining Hurst (Hu) and global

scaling (H) coefficients

Univariate statistical analyses on the original and filtered

sequences of different logs obtained from EBH-2 show statistical
similarities between data sequences before and after FT. In order

to have a better understanding of the data ordering in specific

geophysical logs, multivariate measures, such as rescaled range

analysis and power spectrum analysis should be used (Appendices

C and D).

225 

245 

Gamma (cps) Density (g/cc) 

0 20 40 60 80 100 1.25 1.75 2.25 2.75 

200 200 

205 205 

210 210 

215 215 

220 220 

Gamma-Original 

Gamma-Averaged 

Gamma-Filtered 

Depth (m) 

230 230 

235 235 

240 240

Density-Original 

Density-Averaged 

Depth (m)

225 

245 

Fig. 9. Comparison the original GR and DL of EBH-2 with the curve generated by

Density-Filtered 

averaging the data in every 0.61 m and with the filtered curve generated using

Fourier transform methods.

Table 1

Basic univariate statistical

DL (B), shear modulus (C)

parameters of original and filtered

and Young modulus (D) logs from

sequences

EBH-2.

of GR (A),

Original Filtered

(A) GR (cps)

Mean

Standard error

Median

Standard deviation

(B) DL (g/cc)

Mean

Standard error

Median

Standard deviation

(C) Shear modulus (GPa)

Mean

Standard error

Median

Standard deviation

(D) Young modulus (GPa)

Mean

Standard error

Median

Standard deviation

36.242

0.534

36.690

18.744

2.473

0.004

2.500

0.147

2.213

0.018

2.260

0.631

5.318

0.043

5.330

1.529

36.262

0.495

35.935

17.234

2.473

0.004

2.499

0.138

2.213

0.016

2.279

0.558

5.319

0.040

5.284

1.411

Fractional Gaussian noise-type data can be analyzed by Hurst

coefficients (Hu), which can be evaluated numerically by rescaled

range analysis. R/S is considered a robust method for investigating

the presence of correlations in random events. This analysis is a

measure of how a sequence varies as the distance (lag) between

data points increases according to R=S 1⁄4 ðaNÞ

Hu, where N is the lag

and a is a constant. In general, the longer the interval studied (the

lag), the greater the rescaled range. In R/S analysis, the logarithm

of rescaled range is plotted against the logarithm of lag. If the

resulting plot gives a straight line for the region of long lags, the

slope of this line is denoted as ‘‘Hu’’. The values of Hu determined

from R/S are given in Table 2.

Table 2

Results of R/S analysis

and Young’s modulus

for original

(EBH-2).

and filtered data series of GR, DL, shear modulus

Hu (from original log) Hu (from filtered log)

R/S on EBH-2 logs

Gamma

Density

Shear moduli

Young’s moduli

0.420

0.299

0.401

0.248

0.405

0.285

0.447

0.285

Fig. 10 shows the R/S plots for original and filtered sequences

for GR and DL from EBH-2. These plots show that filtering changes

the behavior of R/S in the whole lag range. However, Hu at long

lags is not severely affected. The Hu values calculated for original

and filtered data using these graphs ranged between 0.2 and 0.4

for density and gamma logs (Table 2). These values are in

0oHuo0.5 range and represent anti-persistent nature of the logs.

The Hu values were also calculated for shear and Young’s

moduli logs (Table 2) for their original and filtered values from

EBH-2. Comparing Hu values of these logs with the corresponding

values of DL and GR shows that the Hu values are close to each

other for shear modulus-GR and Young’s modulus-DL pairs. This

indicates, in the physical sense, that shear modulus is mostly

related to the clay or shale content in a formation and thus to

gamma ray readings, whereas the elastic modulus of a formation

is more related to its bulk density.

Power spectrum analysis is the next technique to be used in

investigating the fractal character and ordering of the data

sequence. Power spectral density is a statistical measure based

on the Fourier transform of the original discrete sequence (e.g. a

well log). In this analysis, the logarithm of the square of amplitude

array (spectral density) is plotted against the logarithm of angular

frequency. The resulting early slope is denoted as b (spectrum

power), which is related to global scaling coefficient, H, by bfGn 1⁄4 

-ð2HfGn - 1Þ for fGn-type series. These values are given in Table 3

for GR and DL of EBH-2.

Spectrum power and global scaling coefficients (H) calculated

using the PS approach are given in Table 3. The table shows that

the global scaling coefficient, HfGn, which is allowed to change

between +1 and -1, is somewhat larger than Hu obtained from R/S

analysis. However, they are close to the values of unfiltered

‘‘original’’ log data and these and the fractal dimensions still show

a pairing between GR-shear modulus and Young’s modulus-DL.
1.00E+00 

1.00E+01 

1.00E+02 

1.00E+00 

R/S 

Gamma-Original

Gamma-Filtered 

R/S 

1.00E+01 1.00E+02 1.00E+03 

1.00E+00 

1.00E+01 

1.00E+02 

Density-Original 

Density-Filtered 

1.00E+00 1.00E+01 1.00E+02 1.00E+03 

Lag Lag 

Fig. 10. R/S plots for original and filtered sequences for GR and DL from EBH-2. The straight lines are the best linear fits to the long-lag regions of the data series.

Table 3

Results of

DL, shear

power spectrum (PS)

modulus and Young’s

analysis for original

modulus (EBH-2).

and filtered data series of GR,

bfGn (original) HfGn (original) bfGn (filtered) HfGn (filtered)

PS analysis on EBH-2 logs

Gamma -0.126

Density -0.525

Shear moduli -0.240

Young’s moduli -0.517

0.563

0.763

0.620

0.759

-0.125

-0.724

-0.150

-0.654

0.562

0.862

0.575

0.827

Fig. 11 shows the power spectrum of the filtered GR and DL

data of EBH-2. Since the power spectrum is obtained via Fourier

transform, the power spectrum of filtered (Fourier transformed)

data is related to the power spectrum of the original sequence by

the convolution theorem [16]. This results in a power spectrum

plot of filtered function versus frequency that is constant up to

some frequency value, ~f. This value is the ‘‘roll-off frequency,’’

which is the value at which a filter adds a constant to within some

percentage to the spectral density of the original sequence from

-N ~ ~ to f. Above f, the amount added rapidly decreases [16].

The results presented in this section show similarities between

the univariate and fractal statistical measures of the original and

filtered data. In other words, the filtering technique applied to the

original data honors the properties of the original data sequence

without critical distortions.

5. Development of a radial basis function network for prediction

of Young and shear moduli from gamma and density logs

A neural network simulates a highly interconnected, parallel

computational structure of the brain with many relatively simple

individual processing elements, called neurons [18] to approximate nonlinear or complex data. Neurons are networked

(topology) in a number of ways depending on problem type and

complexity. Within the hidden layer, the inputs are summed and

processed by a nonlinear function, called a transfer function or

axon. The nonlinear nature of the function plays an important role

in the performance of a neural network.

The process of finding a suitable set of weights for the network

is called network optimization or ‘‘training’’. The most common

way of training the networks is via supervised training algorithms,

which require repeated showings (epochs) of both input vectors

and the expected outputs of the training set to the network to

allow it to learn the relations. However, in some specific cases

(Kohonen self-organizing maps, radial basis functions), unsupervised networks are also used with success.

The neural network computes its output at each epoch and

compares it with the expected output (target) of each input vector

in order to calculate the error. An error is defined for a given

pattern (input–output) and summed over all output neurons over

the entire epoch. The error is summed over all neurons to give an

average mean squared error (MSE) [18]. Minimizing this error is

the goal of the training process. During error minimization, it is

preferable to find the global optimum rather than the local

optimum. Once the training phase is complete, the performance of

the network needs to be validated on an independent data set.

Cross validation is a model evaluation method that indicates the

performance of the ANN when it is confronted with data it has not

yet seen. It is important that validation data not be used as part of

training. More detail is found in [19,20].

One of the most widely used topologies is the feed forward (FF)

network because it can be applied in almost every kind of

modeling, general classification, and regression. Radial basis

function is the most widely used network after FF. The key to

successful implementation of these networks is to find suitable

centers for the Gaussian functions. Radial basis functions are

typically used to build up function approximations of the form

N

yðxÞ1⁄4 X wifðkx - cikÞ (9)

i1⁄41

where the approximating function y(x) is represented as a sum of

N radial basis functions, each associated with a different center ci,

and weighted by an appropriate coefficient wi. In this equation, f

is the radial basis function. Approximation schemes of this kind

have been particularly used in time series prediction and in

control of nonlinear systems exhibiting sufficiently simple chaotic

behavior.

The sum (Eq. (9)) can also be interpreted as a rather simple

single-layer type of artificial neural network called a radial basis

function network, with the radial basis functions taking on the

role of the activation functions of the network. It can be shown

that any continuous function on a compact interval can in

principle be interpolated with arbitrary accuracy by a sum of this

form if a sufficiently large number, N, of radial basis functions is

used. Radial basis function networks have a static Gaussian

function as the nonlinearity for the hidden layer processing

elements. One of the most commonly used types of radial basis

functions is Gaussian function:

fðrÞ1⁄4 e-br2; b40 (10)

The Gaussian function responds only to a small region of the

input space where the Gaussian is centered. Suitable centers can

be found using either supervised or unsupervised learning
approaches, the latter which may in some cases produce better

results. In order to take advantage of both learning methods, a

hybrid supervised-unsupervised topology was applied in this

study to build the RBFs.

1.00E-13

1.00E-11

1.00E-09

1.00E-07

1.00E-05

1.00E-03

1.00E-01

1.00E+01

1.00E+03

1.00E-03

Angular Frequency

Spectral Density

Gamma-Original

Gamma-Filtered

1.00E-23

1.00E-20

1.00E-17

1.00E-14

1.00E-11

1.00E-08

1.00E-05

1.00E-02

1.00E+01

1.00E+04

1.00E-03

Angular Frequency

Spectral Density

Density-Original

Density-Filtered

1.00E-02 1.00E-01 1.00E+00 1.00E-02 1.00E-01 1.00E+00

Fig. 11. Power spectrum plots of original and filtered GR and DL (EBH-2). The straight solid lines are the linear fits to the short-lag regions of the data series. Plots also show

the frequency roll-off when applying FT to filtered data prior to PS analysis.

In this approach, the simulation started with the training of an

unsupervised layer. This layer was trained based on ‘‘FullConscience’’ competition. Its function derives the Gaussian

centers and the widths from the input data. These centers are

encoded within the weights of the unsupervised layer using

competitive learning. During unsupervised learning, the widths of

the Gaussians are computed based on the centers of their

neighbors. The output of this layer is derived from the input data

weighted by a Gaussian mixture [21]. By using the ‘‘FullConscience’’ algorithm in this process, the winning rate of each

neuron is tracked and the weights are adjusted so that there will

not be a single neuron dominating in the whole process but,

rather, each neuron will have a chance in the learning process. At

the conclusion of training, each neuron’s weight represents

reference data, which is later used for classification.

Once the unsupervised layer has completed its training, the

supervised segment then sets the centers of Gaussian functions

(based on the weights of the unsupervised layer) and determines

the width (standard deviation) of each Gaussian. A 2-layer

multilayer perceptron (MLP) was used for processing the

weighted input at this stage.

For developing an RBF-based network, filtered GR and DL data

were used as inputs instead of the original data, since they

included less noise and were capable of representing the original

logs. An RBF network having Gaussian inner functions and trained

on a full-conscience unsupervised learning algorithm was used

and the output was fed into a 2-step multilayer perceptron-based

network, which was trained using a back-propagation algorithm.

The graphical representation of the network and its major

elements and parameters are shown in Fig. 12.

∑

∑

Unsupervised: RBF 

24 Euclidian distance 

Gaussian 

Full-Conscience competition 

100 iterations 

Supervised: MLP 

21 – 10 Neurons for Layer 1 – Layer 

2, respectively 

0.7 momentum 

Hyperbolic tangent activation 

600 iterations

Gamma 

(cps) 

Density 

(g/cc)

Shear M. 

(GPa)

Young M. 

(GPa)

Fig. 12. Architecture of RBF-MLP network developed to calculate shear and Young’s moduli from GR and DL readings.

5.1. Optimization of the network

In optimization or training of the RBF network, filtered log data

from GR, DL, shear and Young’s moduli from EBH-2 were used.

Each of these logs had 1231 evenly spaced data points. First, these

data were randomized to eliminate the biases. Of the randomized

input series (1231 log data points), 180 were tagged for cross

validation, 180 were tagged for testing of the network response

and the rest were allocated for training purposes. In optimization

of the network, it was trained six times in such a way that the 180

optimization inputs were shifted in every training step within the

whole training set in order to change the order of the inputs and

thus to minimize the error caused by ordering of the data in the

randomized set.
The optimization and the testing of the network were

completed with mean squared errors of 0.015 for the optimization

stage, and 0.089 and 0.290 for the testing stage when predicting

shear and Young’s moduli, respectively. The results showed that

the expected and predicted values were close to each other

indicating that the network achieved a reasonable predictive

performance. The performance indicators of the testing phase are

given in Table 4.

Table 4

Performance of network in testing phase for the prediction of shear and Young’s

moduli.

Performance measure Shear-modulus Young-modulus

MSE

Nominal MSE

Mean abs error (GPa)

Min abs error (GPa)

Max abs error (GPa)

R

0.0894

0.3089

0.2209

0.0001

1.0228

0.8337

0.2906

0.1592

0.4213

0.0100

2.0058

0.9183

The moduli data used in the testing phase of the network

(calculated from sonic log) and the RBF predicted values based on

GR and DL are also plotted as a function of depth (Fig. 13).

200

205

210

215

220

225

230

235

240

245

0

Modulus (GPa)

Depth (m)

Shear-Sonic

Shear-RBF Network

Young-Sonic

Young-RBF Network

5 10

Fig. 13. Values calculated from sonic log and the predicted values based on gamma

and density logs in EBH-2.

This

graph shows that the modulus values calculated from the sonic

log are in good agreement with the modulus values predicted by

the network using gamma and density logs, encouraging its use as

a predictive tool in other wellbores where gamma and density

logs are present but sonic log is not available.

5.2. Application of RBF network for predicting shear and Young’s

moduli for EBH-1 and EBH-2 boreholes

200

205

210

215

220

225

230

235

240

245

1

Shear M (GPa)

Depth (m)

200

205

210

215

220

225

230

235

240

245

2

Young M (GPa).

Depth (m)

Sonic log 

RBF 

network 

Sonic log 

RBF 

network 

234 4 6 8 10

Fig. 14. Sonic-log derived shear and Young moduli versus network predictions

using gamma and density logs (EBH-2).

The verification of the network was performed independently

from the network optimization and testing protocol presented in

the previous section. For this process, the arithmetic averages of

GR and DL values within different stratigraphic intervals shown in

Fig. 2 were calculated. Shear and Young’s moduli values were also

calculated using the same method of averaging. To test and verify

the network response on this set of data, the average GR and DL

values were given to the network as inputs and shear and Young’s
ARTICLE IN PRESS

moduli were predicted as outputs. The data from this prediction

process data shown in Fig. 14. The network produced reasonable

results for both moduli parameters with an accuracy of 0.92.

For further verification, the bulk modulus and Poisson’s ratio

calculated using sonic velocities were averaged within different

stratigraphic intervals as described in the previous paragraph, and

plotted against the bulk modulus and Poisson’s ratio calculated

using Young’s and shear moduli from RBF output in Eqs. (6) and

(7) as input (Fig. 15).

210 

215 

Bulk M (GPa) Poisson's Ratio 

0 2 4 6 8 10 0.0 0.1 0.2 0.3 0.4 0.5 0.6 

Sonic log 

RBF 

network 

Depth (m) 

200 200 

205 205 

215 

220

Depth (m)

225 225 

210 

Sonic log 

RBF 

network 

220 

230 230 

235 235 

240 240 

245 245 

Fig. 15. Sonic-velocity derived bulk modulus and Poisson’s ratio versus calculated

ones (Eqs. (6)–(7)) using network-predicted Young and shear moduli from gamma

and density logs (EBH-2).

The comparison shows that errors associated

with RBF calculation of Young’s and shear moduli are not

amplified and bulk modulus and Poisson’s ratio can also be

calculated with a reasonable accuracy.

By using a similar technique, the GR and DL from EBH-1, for

which a sonic log was not available, were averaged within the rock

intervals. These average values were given to the network as

inputs to predict shear and Young’s moduli. The results of the

network are presented in Fig. 16.

Shear M (GPa) Young M (GPa) 

0 1 2 3 1 234567 

200 

205 

210 

215

220

225 

Depth (m)

230 

Depth (m)

235 

240 

245 

250 

200 

205 

210 

215 

220 

225 

230 

235 

240 

245 

250 

Fig. 16. Predicted average shear and Young’s moduli for each strata sections in

EBH-1 (Fig. 3) using the RBF network.

255 255 

6. Summary and concluding remarks

In this study, GR, DL and sonic log data from two exploration

boreholes from a mining area were analyzed for formation

boundaries, shale contents, in situ porosity measurements and

for rock elastic properties.

An analysis algorithm to process available log data was

developed. According to this procedure, univariate statistical

techniques and fractal statistics were used to study similarity

and ordering of the log data in the depth scale. For this procedure,

the data was smoothed using a Fourier transformation with a lowpass filter. Filtered data was tested using basic statistics and

evaluated with rescaled range and power spectrum analyses to

see if this data represented the information characteristics of the

original data set.

Based on the fractal characteristics of the original and filtered

data sequences, filtered data were used in a new prediction

technique for shear and elastic moduli from DR and DL using

radial basis function networks. This technique accurately predicted shear and elastic moduli, showing that it can calculate

shear and Young moduli of the formations when sonic logs are

either not available or too costly to obtain.

Appendix A. Fractional Brownian motion

Fractional Brownian motion is a non-stationary self-affine

random process, defined by Mandelbrot and Van Ness [22] as

1 0

BHuðtÞ1⁄4 " Z Hu10:5

G fjt 1 lj

Hu10:5 1 jl Hu 0:5 j g dBðlÞ ð þ Þ 11

Z t

þ jt 01 lj

Hu10:5 dBðlÞ

]

(A.1)
where G(x) is the gamma function, and BHu(t) an ordinary

Gaussian process with zero mean and unit variance.

In this integral, the parameter Hu is the Hurst coefficient

(0oHuo1). When Hu is 0.5, the repeated integral produces

uncorrelated noise (white noise). The sequence with 0.5oHu o1.0

is a persistent sequence. In other words, if the trend is increasing,

an increasing trend can be expected in the future. If 0oHuo0.5,

then the process is anti-persistent in nature [23,24].

Appendix B. Fractional Gaussian noise

Fractional Gaussian noises are family of random processes with

stationary and Gaussian properties and long statistical dependence [25,26]. The fGn is more varied than fBm since it is fBm’s

derivative, which can be defined by

BHuðx þ dÞ - B Þ 1⁄4 Hu Z ðxÞ ðx lim	 (B.1) d!0 d

It has zero mean and its variance is

Var{BHuðx þ dÞ - BHuðxÞ }

1⁄4 VHud2Hu-2	 (B.2) d

where VHu is

cosðpHu G Þ ð1 - 2HuÞ	 (B.3) pHu

Appendix C. Rescaled range analysis

Hurst coefficient can be evaluated numerically by rescaled

range analysis (R/S analysis), which is a measure of how a

sequence varies as the distance (lag) between data points

increase. This measure has been initially observed and implemented for records of natural data in time. In this analysis, first,

the natural records are transformed into a new variable, X(t, N),

the so-called accumulated departure of the natural record in time

in a given year t (t 1⁄4 1,2,y,N), from the average, x ̄ðtÞ, over a period

of N years. The transformation formula is

X t t

Xðt;NÞ1⁄4 ðxi - x ̄NÞ1⁄4 X xi

i

- tx ̄N (C.1)

i1⁄41 1⁄41

! 

The rescaled range R/S, the range can be defined as

RðNÞ 1⁄4 max1ptpNNðt;NÞ - min1ptpNNðt;NÞ, and the standard deviation by

1=2

1 N

SðNÞ1⁄4 x x ̄ 2 (C.2) N

X 1⁄2 t - N

t

 

1⁄41

!

by dividing the range with its standard deviation, many natural

phenomena can be studied using this power relationship [27,28].

R=S 1⁄4 ðaNÞ	Hu (C.3)

where Hu is the Hurst coefficient. In general, the longer the

interval studied (the lag), the greater the rescaled range. In R/S

analysis, the logarithm of rescaled range is plotted against the

logarithm of lag. If the resulting plot gives a straight line for the

region of long lags, the slope of this line is denoted as H.

R/S is considered as a robust method for investigating the

presence of correlations in random events, or even if the data obey

Gaussian statistics or not [28]. However, some authors suggest

that spectral density technique should be applied first to confirm

fBm and R/S should be applied later to some incremental series of

the data [26].

Appendix D. Power spectrum and spectral density analysis

The spectral density is a statistical measure based on the

Fourier transform of the original data. All of the information in a

discrete sequence created by Fourier transform is retained in

amplitude and phase, which can be converted back to the same

sequence by inverse transform. The square of the amplitude array

is called the ‘‘spectral density’’. In spectral density analysis, one

should plot logarithm of spectral density against logarithm of

frequency. The resulting early slope is denoted as b (spectrum

power), which is related to global scaling coefficient, H, by

bfGn 1⁄4 -ð2HfGn - 1Þ for fractional Gaussian noise; fGn (D.1)

and

bfBm 1⁄4 -ð2HfBm þ 1Þ for fractional Brownian motion; fBm (D.2)

For fGn, b (spectrum power) is between +1 and -1, inclusive. For

fBm, b is between -1 and -3, exclusive. Also, bfGn ffi bfBm - 2 [26].

References

[1] Broch E, Franklin JA. The point-load strength test.	 Int J Rock Mech Min Sci

1972;9:669–97.

[2] Schmidt E. A non-destructive concrete tester. Concrete 1951;59:51–4.

[3] Bemer E, Vincke O, Longuemare P. Geomechanical log deduced from porosity

and mineralogical content. Oil Gas Sci Tech (Rev IFP) 2004;59:405–26.

[4] Palchik V. Use of stress–strain model based on Haldane’s distribution function

for prediction of elastic modulus. Int J Rock Mech Min Sci 2006;44:514–24.

[5] Watt P, Davies GF, O’Connell RJ. The elastic properties of composites material.

Rev Geophys Space Phys 1976;14:541–63.

[6] Sanchez-Palencia E. Non-homogeneous media and vibration theory.	 Berlin:

Springer; 1980.

[7] Sabatakis N, Koukis G, Tsiambos G, Papanakli S. Index properties and strength

variation controlled by microstructure for sedimentary rocks. Eng Geol

2008;97:80–90.

[8] Grima	 MA. Neuro-fuzzy modeling in engineering geology. Rotterdam:

Balkema; 2000.

[9] Verwaal W, Mulder A.	 Estimating rock strength with the equotip hardness

tester. Int J Rock Mech Min Sci Geomech Abstr 1993;30:659–62.

[10] Sonmez	 H, Gokceoglu C, Nefeslioglu A, Kayabasi A. Estimation of rock

modulus: for intact rocks with an artificial neural network and for rock

masses with a new empirical equation. Int J Rock Mech Min Sci

2006;43:224–35.

[11] Maji VB, Sitharam TG. Prediction of elastic modulus of jointed rock mass using

artificial neural networks. Geotech Geolog Eng 2008;26:443–52.

[12] Takahashi	 T, Takeuchi T, Sassa K. ISRM suggested methods for borehole

geophysics in rock engineering. Int J Rock Mech Min Sci 2006;43:337–68.

[13] Dresser Atlas. Log review I. Dresser Atlas Division, Dresser Industries, 1974.

[14] Schlumberger.	 Log interpretation principles/applications. Schlumberger

Wireline and Testing, Sugarland, Texas, 1991.

[15] Karacan CO ̈ . Reservoir rock properties of coal measure strata of the Lower

Monongahela Group, Greene County (Southwestern Pennsylvania), from

methane control and production perspectives. Int J Coal Geol 2009;78:47–64.

[16] Hardy HH,	 Beier RA. Fractals in reservoir engineering. Singapore: World

Scientific Publication; 1994.

[17] Yeten B, Gumrah F. The use of geostatistics and artificial neural networks for

carbonate reservoir characterization. Transp Porous Media 2000;41:173–95.

[18] Eberhart RC, Dobbins RW.	 Neural network PC tools: a practical guide. San

Diego: Academic Press; 1990.

[19] Maier HR, Dandy GC. Neural networks for the prediction and forecasting of

water resources variables: a review of modeling issues and applications.

Environ Modelling Software 2000;15:101–24.

[20] Karacan CO ̈ . Modeling and prediction of ventilation methane emissions of US

longwall mines using supervised artificial neural networks. Int J Coal Geol

2008;73:371–87.

[21] NeuroDimension. NeuroSolutions v. 5.0. Gainesville, FL, 2006.

[22] Mandelbrot BB, Van Ness JW. Fractional Brownian motions, fractional noises

and applications. SIAM Rev 1968;10:422–37.

[23] Arizabalo	 RD, Oleschko K, Korvin G, Ronquillo G, Pardo E. Fractal and

cumulative trace analysis of wireline logs from a well in a naturally fractured

limestone reservoir in the Gulf of Mexico. Geofı ́s Int 2004;43(3):467–73.
[24] Pan S-Y, Hsieh B-Z, Lu M-T, Lin Z-S. Identification of stratigraphic formation

interfaces using wavelet and Fourier transforms. Comput Geosci

2008;34:77–92.

[25] Zeybek AD, Onur M. Conditioning fractal (fBm/fGn) porosity and permeability

fields to multiwell pressure data. Math Geol 2003;35:577–603.

[26] Li F-C. Rescaled-range and power spectrum analysis on well logging data.

Geophys J Int 2003;153:201–12.

[27] Chamoli A, Bansal AR, Dimri VP. Wavelet and rescaled range approach for the

Hurst coefficient for short and long time series. Comput Geosci

2007;33:83–93.

[28] Pallikari F, Boller E. A rescaled range analysis of random events. J Sci Explor

1999;13:25–40.

ARTICLE IN PRESS