JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 102, NO. B2, PAGES 2969-2981, FEBRUARY 10, 1997
Monte Carlo analysis of seismic reflections from Moho
and the W reflector
Klaus Mosegaard
Dept. of Geophysics, Niels Bohr Institute for Astronomy, Physics and Geophysics, Copenhagen
Denmark
Satish Singh and David Snyder
British Institutions Reflection Profiling Syndicate, Bullard Laboratories, University of Cambridge
Cambridge, England
Helle Wagner
Dept. of Geophysics, Niels Bohr Institute for Astronomy, Physics and Geophysics, Copenhagen
Denmark
Abstractø Near-normal-incidence eflections have been used to image the Moho
and the W reflector structure in the lithosphere, offshore northern Scotland. To
determine the impedance variations at these reflectors, we use a Monte Carlo
technique which allows ncorporation of geologically ealistic a priori information as
well as an extensive exploration of the model space, after testing it on a synthetic
data set. The method is based on Bayesian nversion theory. The modeled Moho
consists of a series of layers with a total thickness of •- 2.4 _q- .3 km with an overall
positive mpedance contrast. Inversion of the W reflector results n a model of
five-seven ayers with a total thickness f about 3.7 _+ 0.6 km and mostly nonpositive
impedance contrasts. The implied fine-scale mpedance structure of the Moho is
consistent with the broader velocity structure determined rom previous wide-angle
reflection/refraction rofiles. However, he overall nonpositive mpedance ontrast
at the W reflector requires a structure which is overlain or underlain by a broad
increase n velocity in order to match amplitudes of reflected phases observed at
large offsets nterpreted previously o srcinate at a similar depth.
Introduction
In the last few decades, deep seismic eflection pro-
filing has provided spectacular mages of the continen- tal lithosphere worldwide, of which subhorizontal eflec-
tions from the lower crust, reflections rom the Moho, and reflections from the upper mantle have been par- ticularly notable. Some reflections have been associ- ated with known structures; for example, bright re- flections near the Moho depth have been associated with the crust-mantle transition zone, and dipping re-
flections in the crust have been associated with near
surface aultis nd known ubduction ones Klemperer
and Hobbs, 1991; Clowes et al., 1992; BABEL Work- ing Group, 1993; Zhao et al., 1993; Clowes and Green,
1994]. Many features, mainly subhorizontal nes, e-
main incompletely understood. Although various mod-
els have been proposed or these eatures, one ultimately
requires outcrops, drill holes, or the physical proper-
Copyright 997 by the American Geophysical nion.
Paper number 6JB02566.
0148-0227/97/96JB-02566509.00
ties (density nd seismic elocity) f these eflectors o
choose rom among hese models. Deep seismic eflec- tions, as such, are incapable f providing direct nforma-
tion on physical properties. Instead, one derives some
estimates of the physical properties rom the seismic
record using a modeling strategy. Conventional seismic data processing allows a com-
paratively ast but rough analysis of large volumes of
seismic data. It is essentially based on a linear model of
the seismic race, the so-called onvolutional model, n
which he seismic race s(t) is approximated y a con-
volution of the subsurface eflectivity (t) and a source wavelet w(t):
= ß
According o the convolution heorem, one consequence
of this approximate model s that only hose requencies that are present n the source wavelet can be retrieved
from the reflectivity series of a recorded seismic race.
Seismic sources used in conventional deep seismic ma-
rine experiments ypically have 5 to 80-Hz bandwidths
[Hobbs nd Snyder, 993]. Due o attenuation, seismic
wavelet at a given wo-way ime will have ittle energy
above 60 Hz, and the convolutional model alone cannot
2969
2970 MOSEGAARD ET AL.- MONTE CARLO ANALYSIS OF DEEP REFLECTIONS
WEST
o
EAST
• 40
6o
8O
Figure 1. The DRUM seismic ection howing he Moho (M) and W (W) reflectors nd the
location of the data used n the Monte Carlo nversion. The section was migrated using a two-
dimensional velocity ield derived rom nearby efraction esults Snyder nd Flack, 1990] and
then depth converted sing he same velocities. The Flannan reflector can also be observed,
dipping rom 30-kin depths at the western edge of the profile o 80-kin depths n the east.
provide any nformation bout he reflectivity at spatial frequencies quivalent o _• 60 Hz at that two-way ime.
At frequencies ess han 5 Hz, reflected eismic nergy
is normally absent or obscured by ship noise. Conse-
quently, the convolutional model is unable to predict any low-degree rend or nonzero average value of the
reflectivity.
If one s interested n extracting nformation outside
the limited passband f the source wavelet, t is neces-
sary to incorporate a priori information about the sub-
surface, bandon he convolutional odel and replace t
with the correct relationship between seismic data and
Earth model. In this paper, we propose a flamework
under which variations n impedance nside and outside
the frequency passband f the seismic wavelet can be es-
timated by nonlinear nversion using prior information
on the model in terms of probabilities. A Monte Carlo
inversion echnique Mosegaard nd Tarantola, 995]
provides models which fit the observed data and esti-
mates errors and resolution f these models. We apply
this lamework o shot ecords rom he DRUM (Deep
Reflections rom the Upper Mantle) reflection rofile
from he north of Scotland o estimate mpedance aria-
tions at the Moho and at the W reflector Figures and 2) and to provide estimates f the resolution f these
impedances. The DRUM line was designed o investi- gate the reflectivity of the lower continental ithosphere
and was recorded or 30 s two-way ravel ime (TWT), corresponding o about 110-km depth [McGeary and
Warner, 1985]. Various ipping eflections re observed
in the upper crust 0 to 5 s) and a series f subhorizon- tal reflections n the lower crust (6 to 9 s), the base of
which, at -• 9 s TWT, coincides with the Moho defined
by nearby efraction bservations Barton, 992]. Three
strong reflectors occur in the mantle: a subhorizontal
reflector round 13-15 s TWT (the W reflector), n east-
erly dipping eflector the Flannan eflector) ecorded
from 9 s down to at least 27 s and possibly 30 s, and a
15-kin-long anded one at 23 s (Figure 1) [McGeary
and Warner, 1985]. Derivation f the physical roper-
ties of these deep reflectors rom the seismic ecords s critical to understanding the formation and evolution of the continental lithosphere around the British Isles.
A Probabilistic Formulation of the
Problem
Analysis of the resolution of subsurface structures
from seismic data requires a probabilistic formulation
of the inverse problem. Due to the often strongly non-
MOSEGAARD ET AL.' MONTE CARLO ANALYSIS OF DEEP REFLECTIONS 2971
(a)
13
Trace Number
10 20 30 40 50 60
15
Figure 2. The data used or inversion. a) The shot gather rom 13 to 15 s at shot point 3564
of the DRUM profile, a TWT zone for the W reflector. These 60 traces have increasing offset
from eft (209 m) to right (3210 m). These races were summed n groups f six to produce 0
traces with enhanced ignal-to-noise atios for the inversion. b) The summed races over the
time range of 8-10 s that contains he Moho reflection t about 8.9 s. (c) The summed races
over he the time range of 13-15 s that contains he W reflection t about 13.8 s. (d) Amplitude
spectrum f data with Moho reflection. e) Amplitude spectrum f data with W reflection. f)
Estimated marginal noise distribution or a single data sample. Moho data). The solid curve s
a Gaussian istribution. g) Estimated oise utocorrelation Moho data).
linear relationship between the subsurface mpedance
model and seismic data, the distribution of errors n the
observed data is mapped into the model space as a com-
plex error distribution. Among the important patholo-
gies of this relationship s an inherent nonuniqueness f
models hat fit the data. Further complexity s added to
the model distribution when we, on geological grounds,
must introduce complex, data-independent a priori in-
formation to further weigh models based on their geo- logical feasibility.
We use a Bayesian formulation of the inverse prob-
lem. In this formulation, the state of information about
the subsurface after incorporation of both a priori infor- mation and data information is completely described by
the a postertort robability ensity (m) over he model
space [Tarantola and Valette, 1982]. From c(m) it is
possible o calculate the probability that the true model
belongs o a given class 4 of models:
P(m elongs o 4)-/A (m)dm,
The a postertort probability density s the complete so-
lution to the inverse roblem [Tarantola nd Valette,
1982]. It contains ll the available rior information
such as the approximate sizes of reflection coefficients
and layer thicknesses, nd all the data information, as
"seen hrough the glasses" of the theoretical relation-
ship between model and data in the form of the wave equation. In the Bayesian analysis, all the input infor- mation is preserved and no uncontrolled subjective bias
is introduced.
The a postertort probability distribution for the in- verse problem is given by
a(m) - •'M mldm2"' (m) (m) (1)
[Tarantola nd Valette, 982]. The a postertort roba-
bility density a(m) equals except or a normalization
constant) he a priori probability ensity p(m) times
a likelihood unction L(m) measuring he fit between
observed data and synthetic data calculated from the
model m.
The likelihood function is typically of the form
L(m) = C exp -$(m)]
where C is a constant nd $(m) is a misfit function.
$(m) measures he difference etween he observed
2972 MOSEGAARD ET AL.' MONTE CARLO ANALYSIS OF DEEP REFLECTIONS
(f)
2soJ
200 -
150 -
100 -
50-
(b) (c)
8.0 8.5 9.0 9.5
10.0
13.0
13.5
14.0
14.5 15.0
(d)
(e)
_- _
I I I I I I
10 20 30 40 50 60
Frequency Hz]
i I I
Noise mplitude Arbitrary nit]
0.8 0.6 0.4 0.2 -0.2 -0.4
-0.2
(g)
1
0 10 20 30 40 50 60
Frequency Hz]
,
i i i
-o.1 o o.1
Time lag [s]
Figure 2. (continued)
0.2
MOSEGAARD T AL.: MONTE CARLO ANALYSIS OF DEEP REFLECTIONS
2973
data and synthetic ata calculated rom he model m.
If, for nstance, he observational oise onsists f nde-
pendent aussian rrors, S(m) is proportional o the
sum of squared ifferences etween bserved nd cal-
culated ata values the L• misfit), and he likelihood
function s Gaussian. f, instead, he noise consists f
independent aplace istributed rrors, $(m) is the
Lx misfit where quares re replaced y absolute al-
ues), and he likelihood s Laplacian.
Monte Carlo Sampling of the A
Posteriori Distribution
Real, quantitative eological priori knowledge an-
not be described y means of simple mathematical x-
pressions or p(m). Such nowledge s often available
as statistical information: histograms giving the oc-
currence requency f, for instance, ertain ithologies or physical ock parameters, bserved n outcrops r
nearby wells.
Mosegaard nd Tarantola 1995] rovide method
for Bayesian onte Carlo nversion hat overcomes his
problem. This method has wo major advantages, s
compared o previously ublished ethods e.g., toffa
and Sen, 1991]. First of all, the method s exact n the
sense hat it will provably ample he posterior roba-
bility density. Second, t allows ncorporation f arbi-
trarily complex tatistical priori nformation nto the
inversion.
The algorithm onsists ssentially f two nteracting
parts. The irst part s an a priori model enerator hat
is able o produce andom ubsurface odels, onsistent
with the available priori knowledge. onsistency ere
means hat the generated odels ave exactly he same
statistical roperties s those we have obtained rom
observations n the real Earth. The models produced
by the a priori model enerator re ed nto he second
part, an algorithm hat decides f the a priori model an
pass a data fitting test.
One iteration of the algorithm uns as follows: First,
given he current model, new model s chosen y
the a priori model generator by perturbing he current
model), nd he probability or a given ew model o be
chosen s proportional o its a priori probability. he
new model mn•,v s now accepted r rejected ccording
to the following rule:
1. If the value of the likelihood (mn•,v) of the new
model s arger han or equal o the ikelihood (n•u•)
of the current model, he model mn•w s accepted ith
probability 1.
2. If the value of the likelihood (mn•,v) s smaller
than the likelihood (mcu•), the model mn• is ac-
cepted as he next "current odel") nly with proba-
bility L(mn•)
Paccept --- (mcu•)
If the new model is rejected, he current model also
becomes the next current model.
The series of "current models" produced by this two-
part algorithm re, asymptotically, amples rom he a posterJori istribution r(m) [Mosegaard nd Tarantola,
1995]. After a large umber f terations, he number f
times a given model m occurs n the collection f "cur-
rent models" s approximately roportional o or(m).
A large collection f such amples rom or(m) provides
the raw material from which various characteristics of
the models can be obtained [Mosegaard nd Taran-
tola, 1995]. A set of statistically ndependent amples
of the accepted models allows structures n the sub-
surface hat are well-resolved o be distinguished rom
those that are ill-resolved. A well-resolved structure
will appear n most of the accepted models, whereas
an ill-resolved tructure will appear n only a few mod-
els. The probability hat a certain structure exists s
roughly roportional o its frequency f occurrence n
the set of a posteriori samples. Therefore approximate
a posteriori robability distributions or model parame-
ters can be represented y normalized istograms f the
parameters alues. The peak of the histogram ill be at
the most probable model parameter, nd the deviation
from this peak will provide a measure f uncertainty.
Monte Carlo Sampling of Lithosphere
Reflectivity
We use he Bayesian Monte Carlo algorithm escribed
above o generate eflectivity models or selected arts
of the DRUM reflection rofile. We have concentrated
our efforts nly at a location Figure ) where he Moho
and W reflectors re brightest nd subhorizontal s he
algorithm equires arge computation ime. We have
analyzed total of 60 unprocessed races n the nterval
8.0 s - 10.0 s TWT bracketing he Moho and he nterval
13.0 15.0 s covering he "W reflector." We chose he 60
traces rom one shot gather 3564) with the best signal
to noise atio among 0 neighboring hot ecords. he
receiver roup pacing as 50 m along 3000-m-long
streamer. The raw shot gather shows oherent nergy
between 3.7 s and 14.2 s TWT for the W reflector Fig-
ure 2a). To ncrease he signal o noise atio, we applied
a correction or dip move out and stacked ix adjacent
traces, ielding 0 traces or the nversion Figures b
and 2c). The normalized requency pectra how om-
inant nergy etween and 30 Hz (Figure d and 2e).
Since he Moho nd W reflectors re Subhorizontal nd
are at large depths, he incident eismic aves re es-
sentially ertically raveling lane waves. t is therefore
possible o perform rather areful alculation f syn-
thetic seismograms sing fast one-dimensional ropa-
gator matrix method Haskell, 953].
Wavelet Estimation
The DRUM profile was shot by GECO Geophysical
Company f Norway) sing n 8536-cubic-inch ir gun
array at 7.5 m below he sea surface McGeary nd
Warner, 985]. A simulated ar-field ource ignature

Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.