Note: Descriptions are shown in the official language in which they were submitted.
S ~ ~
S P E C I F I C A T I O N
Title: S~ISMOGRAPHIC METHOD AND APPARATUS
D~SCRIPTION
This invention relates to a method of determining the
location in the earth of sub-surface boundaries and/or
the acoustic properties of sub-surface layers in the
earth and to apparatus for this purpose.
~ = .
. . .
There is well-kno~ a seismic reflection technique
~ which employs a sound source at or near the earth's
; surface to emit an impulsive sound wave at a known time.
As this.sound-wave:passes through the~-earth it en~ounters
~ .
boundaries between:~the-.di-fferent sub-surface layers~:~At-;~
each boundary some:of.the-sound-is-.transmitted-and some.
i reflectsd. A re-ceiver at or-near--the-sur~aGe-close
to-the-source detects the reflected waves which arrive
at later and-later~times.
A record.(or-seismogram)-made.o~ th~.receiuer response.:--.
is then processed to determine the amplitudes and
ar~i~a~ ;times- o~ t~-e individual reflections. ~hes~.ma~.. ..
.'
. . .. .
s~
\
then be used to determine the locations of the xock
boundaries within the earth and/or the acoustic properties
of the rock either side of each boundary.
The accuracy of such analysis depends on the ability
of the processing technique to separate the individual
reflections from one another. One of the reason that
this is a non-trivial task is that it is extremely
difficult to generate a purely impulsive sound wave~
The sound wave generated by most seismic sound sources
has a duration which is longer than the smallest
separation time interval of which the recording system
is capable. In other words, the series of reflections
which arrives at the receiver is not the series of
impulses one desires (the reflectivity series); it is
a series of over-lapping wavelets. The processing step
which is used to remove the effect of the source from
the recorded signal in an attempt to recover the reflectivity
series is usually known as deconvolution.
The conventional description of a seismic signal regards
the propagation of seismic waves as a linear elastic
process in which the signal xl(t) is obtained as the con-
volution of the impulse response of the earth g(t) with a
far field source wavelet s(t). Usually some additive noise
is also present so that
xl(t) = s(t) * g(t) + nl(t) (1)
- 2 ~
( ~)
S~
where the asterisk * denotes convolution. One wishes to extract
g(t) uncontaminated by either s(t) or n1(t). However, n1~t) is
not normally known, and often s(t) cannot be measured or predicted
and must also be regarded as unknown.
Since s(t), g(t) and n1(t) are all unknowns, the problem of
finding g(t) from the measurable quantity x1(t) is basically
that of solving one equation containing three unknowns. It cannot
be done, of course. ~ven when the noise can be ignored the essen-
tial difficulty remains; that of deconvolving s(-t) and g(t). Unless
s(t) is known g(t) cannot be found without .a lot of assumptions.
The noise term ~ nl(t) will be ~ll compared wlth the signal term
s(t) * g(t) provided there is enough signal energy~ In order to
achieve an adequate signal-to-noise ratio it is sometimes necessary
to repeat the experiment a number of times in the same place, up to
a total of, say, p time~, using the same sound source, or identical
sound sources. ~he series..of-received seismic signals x1(t), x2(t),.
...xp(t~ is summed to form-a composite-signal-x(t), where .
x(t) = ~ Xi t)
i = 1 . '
~quation(1 )then becomes
x(t) = p.s(t) * g(t) + n(t), (1a)
where n(t) is a composite noise signal, given by n(t) = ~ ni(t)~~'i
~his summation i~ known as a "vertical stack" and p i8 the "fold of
stack"and i9 an integer greater than or equal to 1. ~he same, or
similar,-~result may-some.times be obtained by generat~.ng p identical- -
.... . . . . . . ............... - . . = . . .......... - - .
.. ~ .. ~ .. . .
(4)
S~
impulsive sound waves simultaneously. If the sources do not inter-
act, each one will generate an identical far field wavelet s(t),
and the received signal x(t) will be described by equatioh (1a).
~or more than twenty years much ingenuity has been devoted to de-
vising methods for solving equation (1) or (1a) using assumptions
which are as realistic as possible. But the fact remains that
these assumptions are made purely for mathematical convenience.
They are not substitutes for hard information.
The best known.example of such-a method is the least - squares
time - domain inverse filtering method used throughout the industry.
For this method to be valid it is required that
(1) g(t) be a stationar~, white, random sequence of impulses;-..
(2) ~(t) be minimum-phase.and have the same shape throughout the
seismigram;
(3) there be no absorPtion.
All these as~umptions-æ e--very strong, and they must all be correct
simultaneousl~^.if:the.,method-is to:work~ This condition-iæ very-- -
difficult to-satisfy,..especially-since,the assumptions are not.,=~
mutually reinforcing. ~or example..in attempting to æatisfy,,th~ta-..,''
tionarity assumption,: some sort of spherical divergence--correction
must-first be applied.-~his has the effect of.distorting s(t~ un--
evenly down the seismogram whi,ch immediately invalidates the-
assumption that-,the shape of s(t) remains constant;-it:also intro~
duces a tendency for s(t) to-be non-minimum-phase-in-the early part
o~ the seismogram,.
, ,
~ : J
(5
~16~S3
According to a first aspect of the present invention there is
provided a method of determining the location in the earth of
sub-surface boundaries and/or the acoustic propertie~ of sub-
surface features in the earth which method comprises employing
one or more first and one or more second point sound sources to
produce respectively first and second sound waves, the energy
of the elastic radiation of the or each first source differing
by a known factor from the energy of the elastic radiation of the
or each second source, detecting reflections of said first and
second sound waves from within the earth and generating therefrom
respective first and second seismic signals and subjecting these
two seismic signals to analysis and comparison.
According to a second aspect of the present invention there is
provided apparatus for determining the locati.on in the earth of
sub-surface boundaries and/or the acoustic properties of sub-
surface features in the earth which apparatu~ comprises one or
more first point sound sources and one or more second point sound
sources adapted respectively to-produce first and-second sound - I
_....... !
waves in the earth, the energy of the elastic radiation of the. .
or each first source differing-by a known factor from the energy -. .
of the elastic radiation of the or each second source3 receiver.
means for detecting reflections of said first and se.cond sound
waves from within the earth and generating therefrom respective
first and second seismic signals, and means for analysing and
comparing said first and second seismic signals.
In the present inventicn the first and second sound sources may
be individual point sound sources or there may..be employed a
:
--- plurality of identical non-interacting Point sound sources which
=_
will produce a seismic signal having a greater signal-to-noise ratio-O~-l
(5a)
Alternatively the receiver means may be adapted to sum a series
of identical seismic signals obtained by repeated production of
identical sound waves by one or more iden-
tical sound sources.
The term "point source"is employed throughout the specification
to denote a source whose maximum dimension is small compared with
the shortest wavelength of the useful radiation it generates. If
this source is buried in an homogenous isotropic elastic medium
it will generate spherically symmetric radiation at distances
greater than about a wavelength. This is the far-field region in
which any aspherical distortions of thewavefield from this point
source will occur only at high frequencies outside the useful
bandwidth, ~
. /
_ .. .. ..
. / .
/ -- - . ' .
!
-6- ~6~S3~
The method Or the present invention is suitable for buriea
point sources, on land ana at sea. It requires ~one of
the ~ssumptions demanded by known ~et}ods. In pa~ticular,
nothing is assumed about differences in the a~plitude or phase
spectra of s(t) and g(t). The present invention is based
upon the fact that the wavelet obeys a scaling law of the t~pe:
s1~ 2) = c~ s (Cc~ ) (2)
In this equation,~1 and ~2 are both very nearly equal to
~=t-r~ whe~e t is time ~easured from the shot instant, r is
the distance from the sound source to a point in the far
field, and c is the speed of sound in the medium; s1 ~ 2) is
the far field wavelet of a source similar to that which
~enerates s(rj), but which contains o~3 times as much ener~y.
~iure 1 sho~s diacrz~m2tically how this source scaling 12w ,~
affects the far field ~avelet. ~;
There is excellent experi~e~tal evidence- for the existence ¦
of such a scaling law for a variety of point sources and -
this law can readily-be derived for explosives,-for exa~ple~ ~
if the following assumptions-are ~ade:
i. that the elastic radiation from the source possesses---
spherical symmetry; thus it will be applicable to ~ost
marine sources such as a single air gunj a single water-gun~
a marine explosive such as that available under--the Trade ~rk
"~laxipulse", a marine source emplo~ing high pressure steam to
cause an implosion such as that available under the Trade ~ark
"Vaporchoc"l ol a sparker and to;explosives buried on land but
probably not to surface sources because~ --
(7)
39
their radiation is not spherically symmetric;
ii. that the fraction of the total available energy
stored in the explosive which is converted into elastic
radiation is a cbns~t for a given type of explosive
and a given medium;
iii. that the volume of the explosive may be neglected
relative to the ~olume-of the sphere--of anelastic defor-
mation produced by the explosion,
iv. that the elastic radiation produced by the explo-
sion could be obtained by replacing the sphere of anelastic
-~:' deformation by a cavity at the interior of which
there is appl.ied a time-dependent pressl~e function P(t)
and that 2(t) is independent of the mass of the e~plosive
and is',constant for explosives of.the same chemical
composition in the:-same---medium; and ~
._ .
~ v.: that ~r1 -.-.fo~ an--explosi~n-of a first mass and- '
,., ~ .
~',for-an explosion-.:of.. a second-ma&s ~an be,taken-to-be: ,
approximately equal and equal to ~. Thi-s-is sufficientl~ --
accurate: if the time-interval ~^lr ~etween~ and ~ is
unobservable within:the frequency band--of~interest7-i.e.~
~ r =should--be less than about one sam~le interval. -This--
approximation will suffice for values of C<up to about
5 or soc
.. . .
.... . .
. . .
.' .
.
LS~
To exploit the scaling law a seismic signal x(t) as described
by equation (la) is generatedO The experiment is then re-
peated in the same place using a source of the same type but
containing a3 times as much energy. This will generate a
seismogram:
xl(t) = qsl(t) * g(t) + nl(t) (32
where sl(t) is the far field wavelet of the source and is
defined in equation (.2).; g(t) is the same as in equation (1)
because it is the response of the earth to an impulse in the
same placei th,ë. noise nl(.t) may be different from n(t) in
equation (la),; q is a known integer greater than or equal to
1 and which may be different from p in the equation (la~
.
Let us consider these equations together for the case where
the noise is negligibly small:
x(.t) = ps~t) * g'(t) (4~
x (t) = ~s (.t), * g(.t) (5)
s (t), = ~s(t/~) (2)
In these three independent equations there a~e three unknowns:
s(,t).~ sl(,t). and g(,t) r Therefore, in principle when the
noise is negligi~ly small~ we can solve for all three exactly
without maki~ng further assumptions.
.
By ta~ing Fourier transforms and by manipulation we obtain
- the equation
qS(af) = P2 S(f), R(f) (.62
a
where R(,fl is defined as for equation 13 hereinafter.
, .
~ 8 -
s~g
Equation (6) suggests a recu~sive algorithm of the form:
9~5(C~r\4D ) '- ~ 5 ~ ) R (~
n ~ J ~
where ~ is dicated by the highest frequency of interest
and the process must be initiated with a g~less at ~0. If
b~ > 1 equation (7) enables us to worX up the spectrlLm
calculating values at ~(fo , ~2 fo , ,..,.~. ~ fO
star~ing with a guess at fO.
To co~pute values at frequencies less than fO, equation
(6) can be rearranged: /
p S (~ R (f )
~ such that we obtain the recursion:
3 p ~
-~ where M is dictated by--the lowest frequency of interest.--This
now enables the values-:at frequencies ~0 /~ f ~/~7 ' _
to be computed.
. . .
, : " ,
q
(10)
S~
~hus from the recursion scheme of equations ~7) and (8)
we can obtain values at frequencies ~0/~ 7 f~
~/~ o , ~ f o ~ o( f o ~
We can now use an interpolation routine to find a value
at another specified frequency, sa~ f1 and use the
recursion to calculate-values-at-~f1, ~2f1 etc. ~his
procedure is repeated until sufficient values have been
computed. Once S(f~ has been calculated, s(t) is obtained
by taking the inverse Fourier transform.
It should be noted that the quantities involved in the
al~orith~ are complex. One can operate either with the
~odulus (amplitude) and arg~ment (phase) 3 or with the
real and imaginary parts. The~eal and-imaginary parts have
been used in the-example as.-these are considered to be
the ~ost 'basic' components o~-the comp-lex.n~bers.in-a~-
computer~ whereas ampl-itude:-and phase-are---ad~i~tures~~~o~~~
these quantities9
The Initial Guess
The.algorithm.is initiated -with a guess. If this:.guess.is
wrong, the final result will be wrong. The guess at fO ~s
a complex number which, in all probabilitg, will not be
the true value at f . In fact, the-guess SG(.~ -is-related -
~,
~lS~9
to the true value S(f ) in the following way :
S ~ ' S ~
. (9)
i where r e is.the u~known complex error factor. ~f this
error-is not taken.into account tnere will be generated
the values:
~ 5~ 2
~ ~ l
l '
2)
(10)
., which,--.with sufrisient interpolation-yiela~the functio~--
~ SG:-(fO)~forh ~ ~ ...the range :ca~-:be~.extended-t`o~
: the origin by defining~ (o)~ D which:is-compatible.with.
a time-series-8G (t)-with-z-er-o mean.
~he ~ect--of.the-initial error ca~ be seen~by substituti~g
:equation 9 into equation 10; thus:
~S4(~ ~o)=~
~"~2 S(~ 0 ) `
~` ~ '' ' . (11)
"~
~ , . . .
~ 39
It is evident that the error factor is constant fpr ali the
values deduced.from the algorithm. Thus far the algorithm has
allowed computation of the function:
SG, ~ e, S~f), (O ~
(12)
where ~0 has been assumed to be positive.
~wo problems now exist. ~irst-the transform must be co~pleted
by generating values of SG(f) at ne~ative frequenciesO
Secondly the error factor must be found to obtain S(f) from
equ~tion (12). Both these problems can be solved by consider-
ation-of the physical-properties of s(t~, which impose con-
strain~s on the properties-of S(f)...`.
It is known th~t S(t) is real, and therefore-the estimated
wavelet should be real. This constraint imposes Hermitian
symmetry on Stf). That-is, the real and imaginar~-parts.of
S(f) must be eve~ ~nd odd functions, respectivelg. Thus if
S(f) is known-for positive frequencies, S(f) can easil~ be
computed for negative frequencies using this condition.
Howe~er, only SG(f) is known, which is in error by a phase
shift 0 and a scale factor r. The scale factor is unimportant
because it has no effect on the shape of s(t), and conse-
quently cannot affect our estimate of the sampe of g(t). It
can therefore be ignored.
However, the phase error ~ cannot be ignored, because this
will make SG(tl non-causal~ and it is known that s(t? is
causal. That is, s(,t) is zero for times t less than zero.
In the frequency domain causality imposes the condition that
the odd and eyen parts of the Fourier transform are a Hilbert
transform pair. It can be shown that t~is causal relation-
ship is destroyed unless the phase error ~ is zero.
This consideration sug~ests a trial-and-error procedure for
improving the estim,ate of s(t), This is as follows:
lo Compute SG(f) from an initial guess at fO as
described above, notin~ that SG(f~ and S(f)
are related as in ~12)~
2. Multiply S~(,f) by a correction factor e i~G
where~ Q~ is a guess.
3. Impose Hermitian symmetry.
4, Check for causality. If the recoYered wavelet is
non-causal, return to step 2 and repeat, using a
different HG. This procedure is repeated until
the causality condition is met.
mab~c
(14)
~ 6~ 5~9
Thus the equations may be solved in the.frequency domain
using the algorithm described above and applying the con-
straints which follow from two physical properties of s(t):
it is real and causal~ ~he final estimate of s(t) will be
in error only by a scaling factor r, which is tri~ial.
Having obtained~a satisfactory estimate of s(t),g(t) can be
obtained using e~uation (1), by standard methods~
The algorithm described depends on -a complex division in
the frequency doLlain. :~here are two problems associated with
this. First, the ratio will become unstable~-at any freouency
at which the amplitude of the denominator is too s~all.
.~.condly, if the denominator contains non- minil~u~-pha.se
components ~lhich are not conta~ned in the numerator then
the quotient becomes-unstable in the sense that it is non-
realisable.
To solve the ~irst problem-it is usu-al--to-add--a s~all threshold :.
of white noise to the denominator to negate-the possibilit~
of zero -or near zero division. An alternative but ~ore time--
consuming method is ~o.search for -low values in the denominator -
~nd to replace them with small-positive -values~ -
Finding the inverse of non-minimum phase wavelets i9 a well-
kn`own ~roblem. However, the problem can~be?a~oided~-simply b~
(15)
5~
applying an exponential taper of the form e ~t to both x(t) and
x (t). By choosing ~ large enoughthe quotient R(f) can be forced
to be stable, but then the estimates of s(t), s1(t) and g~t) will
be distorted~ In practice the distortion may'be removed simply by
applying the inverse taper e ~t to these functions.
In the presence of noise the problem is to obtain a reliable
estimate of the ratio spectrum R(f), for then the scaling law and
recursive algorithm can be used -to find S(f) as described above.
From equation (6) we define R(f) in the absence of noise as:
x1(f) qS (f)
R(f) = = - ~ ' (1~)
X (~) pS (f)
It follows that
s1(t) = r (-t) * s (t) (14)
.
where r(t) is the inverse Fourier transform of R(f) and, since
s(t) and s1(t) are both real ~nd causal, r(t) mu~t also be real.
However,~~'r(t) will not be causal unless s(t) is minimum-phase.
Both s(t')-and s1(t) must~~e forced to be minimum-phase by applying
the exponential taper to x(t) and x1(t) as described above. Under
these conditions r(t) will be real and causal.
.. ,; . . , . .- , -
, .. ~ .. ..... . ~
(16) ~l ~l 5 39
In the noise-free case it is also true that
,
x1(t) = r(t) * x(t) (15)
and it will be seen that r(t) is simply a one-sided filter which
shapes x(t) into x1(t), provided the correct exponential taper
has been applied. When noise is present the estimate of r(t) must
be stabilized and this can easily be done using a least-square~
approach (N. ~evin~on, in N. Wiener, 1947; ~xtrapolation, Inter-
polation and smo~thing of S-tationary ~ime Series, Wiley~ New York).
A ~hat is, a filter r1(-t) is found which,-for an ~ x(t)-wlll
l'npu~
give an ~u~t which is the best fit in a least-squares sense to
x1(t). This filter r1(t) will be the best estimate of r(t).
In other words, in the presence of noise r(t) can be calculated
in the time domain using standard progra~s, and then its ~ourier
transform taken, whence s(t) etc, can be found as described above. -
It will be understood:that although~-the problem-has been discussed ~
.
: in terms of the scaled-energies of the elastic radiatlon-of the-=-
sources, normally
. . .
(17)
the particle veloc~ies or the sound pressures generated by
the source ma~ be detected and recorded using respectivel~
a geophone or a hydrophone as conventionally employed.
It will further be understood that the individual elements of
the apparatus Or this invention ma~ .be chosen at will to be
suitable for the particular purpose for ~hich they are required
thus air guns, water guns,!IIaxipulse ~Vaporchoc Jsparkers etc.
~i ~ag be employed as the sources. Similarly any suitable
¦ analysers, receivers etc may be em~loyed as necessary.
¦- ~ It is believed that ~may have a value of.from 1.1 bo
~ore pre~erablg fromi1.5 to 3
I
~X~ ~LE
Applying the above-~avelet deconvolution-scheme---to:a -~-
.. s~nthetic-example................................................. -
: Two independent-synthetic far field source-wavelets as:sh~Q ~-
: in ~i~ure 2 were generated. Each wavelet was-calculated:usi~g
a ~odel as-describea in the Geoph~sical Journal of the~Royal
hstronomical Soci-ety 21, 137-161, for the signal generated
. by a~ air gun in water. The ~odel is based on the nonlinear
. : oscillations of.a spherical.bubble-in wateriand itakes into - -.
~- .
. ,`''~ .;
~ S 3~
account nonlinear elastic effects close to the bubble. The
model predicts waveforms which very closely match measure-
ments.
The top wavelet s(t) of Figure 2 was computed for a 10 cubic
inch gun at a depth of 30 feet/ a firing pressure of 2000
p.s.i. and a range of 500 feet from the gun. No sea'f sur-
face reflection has been included. The bottom wavelet sl(t)
was computed using the same computer program for a 80 cukic
inch gun, at the same depth, firing pressure, and range.
In other words, only the volume was changed.
Secondly each of these wavelets was convolved with a syn-
thetic reflectivity series g(t) shown in Figure 3. The
result of performing these convolutions is shown in Figure
4. The top trace x(t) represents the convolution of g(t)
with the upper wavelet s(t) of ~igure l; the bottom trace
x (tl represents the convolution of g(t) with the lower wave-
let sl(t) of Figure-l. Thus these two traces, x(t) and
xl(t), were constructed entirely independently without any
use of the scaling law.
,
It was then assumed that these two traces had been obtained
knowing ly-that they were from the same place and that
the top one was made using a 10 cubic inch gun, while the
bottom was made using an 80 cubic inch gun at the same depth
and pressure.
- 18 -
mak/J t
( 1 9 )
1~15~9
Since only the gun volume was changed, the scaling law can be
invoked. In this case 3 = 8; therefore = 2.
Solving for s(t) and g(t) as described above using the set of
simultaneous equations (4), (5) and (2) with p = q - 1, sub
stituting ~ = 2, the recovered wavelet and reflectivity series
are shown in Figure 5; they compare very well with the top
wavelet of Figure 2 and the original reflectivity series of
~igure 3. ~he small difference between the recovered series
and the original are attributed primarily to computer round-off
error.
This example shows that the method is valid in principle.
~hus by means of the present invention in the absence of noise
the impulse response of the earth can be obtained exactly. In
the presence of noise a stable approximation to this impulse
response can be obtainedj- the accuracy of which approximation^
is dependent on the noise level present.
.
..
.
