- Email: [email protected]
Use of the kepstrum in signal analysis
Geoexploration, 16 (1978) 55-73 o Elsevier Scientific Publishing Company, Amsterdam - Printed in The Netherlands
55
USE OF THE KEPSTRUM IN SIGNAL ANALYSIS
MANUEL T. SILVIA’
and ENDERS
1Naval Underwater ‘Northeastern
A. ROBINSON2
Systems Center, Newport, R.I. 02840 University, Boston, Mass. 02115 (U.S.A.)
(Received September
(U.S.A.)
3,1977)
ABSTRACT Silvia, M.T. and Robinson, E.A., 1978. Use of the kepstrum in signal analysis. In: C.H. Chen (Editor), Computer-Aided Seismic Analysis and Discrimination. Geoexploration, 16: 55-73. The kepstrum has evolved as an important mathematical construct in the area of digital signal processing. In particular, kepstral analysis has been used for the deconvolution of reflection seismograms generated by the exploration for oil and natural gas. Here, we review the development of the normal-incidence reflection seismogram from a physical point of view and show how the kepstrum relates to the fine-grain structure of the physical model. Moreover, we develop symmetry properties of the kepstrum, introduce the concept of minimum-advance, and discuss the relationship between spectral factorization and the kepstrum. In the context of deconvolution, we discuss the methods of predictive deconvolution and homomorphic deconvolution and their dependence on the physical model of the reflection seismic process.
INTRODUCTION
The idea of the kepstrum appears in the classical work of Poisson (1823), Schwarz (1872), Szego (1915), and Kolmogorov (1939), and has been applied to geophysical problems by Robinson (1954), Bogert et al. (1963), Schafer (1969), Oppenheim and Schafer (1975), Tribolet (1977), and others. In this paper we bring all this work together and introduce additional properties of the kepstrum so as to bring out its essential features and symmetries. As a result, the kepstrum evolves as an important construct in the theory of signal analysis and time series, and takes the place among such established functions as the autocorrelation and spectrum. For further discussion see Beth (1968, 1974) and Kulhanek (1975). In their classical works on potential theory, Poisson (1823) and Schwarz (1872) addressed the problem of determining a potential function whose real part was an assigned value on the unit circle. As is well known in potential theory (see Bateman, 1944), Schwarz’s classical expression is:
56
. ,
H(z) = ,‘n
1
(:
‘~~~~,)
Re[H(w’)]dw’,
lz I < 1
(1)
with w’ as the integratioh variable and z = rewiw . This equation gives the potential function H(z) whose real part on the unit circle is Re [H(o)]. The Schwarz expression follows from the famous result of Poisson on potential theory. An interesting result arises when we consider a logarithmic potential log H(z), which can be written: log H(z) = log I H(z) I + i arg [H(z)] The real part of the logarithmic potential evaluated on the unit circle z = eAiw is then log 1H(z = eszW) I = log I H(w) I. In terms of the logarithmic potential, Schwarz’s expression becomes: .
log H(z) = $
_i:
(: + ;$)
I
log I H(o’)
I dw’, I z I < 1
where log I H(w) I is the assigned or known value on the unit circle and the function H(z) has no zeroes or poles inside or on the unit circle, i.e., H(z) is minimum-delay. (This definition of minimum-delay is based on the Laplace definition of the z-transform, i.e., the z-transform of a sequence h, is defined by H(z) = 5 h,zn .) As a result, log H(z) has a Taylor-Maclaurin series exn=-m pansion in the unit circle: logH(z)=
cf,z”,lzl
(2)
n=O
The coefficients fn of this Taylor-Maclaurin series are the coefficients of a stable causal system. We call these coefficients the kepstrum of the minimumdelay system H(z). However, as we shall see, the concept of the kepstrum is not limited to minimumdelay systems. In addition to its appearance in the logarithmic potential problem, the kepstrum was used by Szego (1915) and Kolmogorov (1939) to solve the problem of factoring the power spectrum @( w ) of a random process in order to extract a stable causal system H(z). A power spectral density function (a(w ) must satisfy the following conditions in order to be factored so as to yield a causal system: (1) Q(w) is positive on the interval-n < w < 71. (2) @(w )has finite area on the interval --71< o < n. (3) log @((w.)has finite area on the interval --71< o < n. Now the problem of spectral factorization is the problem of extracting a causal system H(z) whose magnitude spectrum is the square root of the power spectrum. In symbols:
57
I H(w) I = [Q(w)1
‘/i
Szego (1915) and Kolmogorov (1939) recognized that a solution to this problem could be obtained by making use of the Schwarz expression (1). Now log I H(w) I is the real part of the function log H(w). Further: Re [log N(w)] = log I H(w)
I = ?h log (a(o)
Thus, the required solution is found by merely substituting ‘/z log (a(w) for Re [H(w)] and substituting log H(z) for H(z) in Schwarz’s expression to obtain: . H(z)
log
=
I
k i (‘,y;fi:.);log
@(w’)dw’,
Iz i < 1
-lT . . 1s.
Hence, the Szego-Kolmogorov
= exp [log H(z)] = exp [&
H(z)
solution to the spectral factorization problem
i (:r;!;:) -77
logrp(,‘)da$
:z I< 1
The desired sequence h,,, hl, h2, . . . can be written in terms of the kepstrum fo, fl, f2, - * * as: h,+hIz+h2z2+...=exp[fO+fIz+f2z2+...]
We call this equation Kolmogorou’s that: ho =
(3) equation.
Putting z = 0 in eq.3 we observe
exp(fo).= exp
Thus, the leading kepstral coefficient fO is the logarithm of the leading impulse response coefficient ho, i.e., fo = log h,
As a result, it is often convenient to normalize an impulse response function that its leading coefficient is unity; thus, the leading kepstral coefficient is zero, i.e., fO = log 1 = 0. With this normalization, Kolmogorov’s equation (eq.3) is rewritten as: h, so
1+;2.+;.*+...= 0
exp (fIz + f2z2 + f3z3 + . . .)
(4)
0
We shall define the power series expansion flz + fizz + f3z3 + . . . appearing in the exponent of Kolmogorov’s equation (4) as the Kolmogorov Equation Power Series (KEPS). Further, since KEPS is the z-transform of the sequence 0, f1, f2, f3, * * -9we can think of fn as being a time response. Thus, we shah de-
58
fine the sequence 0, fI, fi, f3, . . . as the Kolmogorov Equation Power Series Time Response and call this sequence the KEPSTRUM, where we have added the Latin singular ending “urn” to denote one kepstrum and shall use the Latin plural ending “a” to denote more than one kepstrum, i.e., kepstra. Whereas Szego and Kolmogorov were only concerned with the real part of the function log H(z) on the unit circle, Robinson (1954) was concerned with the imaginary part as well; namely, both the real and the imaginary parts have physical meaning. Thus, it is important to consider the logarithm of the spectrum H(w). Let us elaborate. For a minimum-delay system H(z), the kepstrum fn was defined by the relation:
logH(z)=g f,z~,lzl
The complex variable z represents a point within the unit circle (i.e., I z I < 1). However, we can consider the limit of log H(z) as the interior point z approaches a point on the circumference of the unit circle; that is, we consider the limit of log H(z) as z-teTZW. This limit is: log H(w) = 5
f,eminw
n=O
Separating the above into real and imaginary parts, we obtain: logIH(w)I+iOo-(f~+Re[
‘&fneminw])
+i
(Im[~lfne-i”w])(5)
By equating real parts in eq.5, it can be shown that: fit = 2fk for n = 1, 2, 3,. . . f. = f;
for n = 0
f-“, = 2fl, for n = -1, -2, -3, . . .
where : f,‘, = &
i log IH(w) I einw dw for n = 0, +l, +2, . . . -A
(6)
and the superscript * denotes the complex conjugate, i.e., f,” is the complex conjugate of fn. Thus, the kepstrum fn is found in terms of log I H(w) I, thereby allowing us to obtain the minimum-delay system H(z) via (2). Let us now look at the imaginary part of eq.5, namely, by equating the imaginary parts in (5) it can be shown that:
59
e(o)
=zpj cot
(yj
log I H(d)
I dw’
(7)
--n where we have incorporated the result of (6) into (7); the symbol P denotes that the integral has its Cauchy principal value. The kepstrum was used as an intermediate step in deriving (7), which shows that for a minimum-delay system H(z), knowledge of only the log magnitude spectrum is sufficient to determine the phase spectrum 0 (w). In the mathematics literature, eq.7 and various versions of the Schwarz formula (1) are often referred to as Poisson’s formulas. In the context of linear system theory, eq.7 is known as the Hilbert transform, where e(w) and log I H(w) I form a Hilbert transform pair. Bogert et al. (1963) were the first to use the kepstrum as a signal processing operation. Bogert et al. were primarily interested in measuring the time interval between the arrival of certain seismic waves, excited by either natural or man-made disturbances. They defined the term “cepstrum” as the power spectrum of the logarithm of the power spectrum of the observed time series. We see that the “cepstrum” of Bogert et al. involved no phase spectrum information, i.e., it utilized magnitude spectrum and power spectrum information only. Their “cepstrum” was effective in detecting or measuring the time delay of a single echo, even when other methods (such as the autocovariance method) were ineffective. However, it is important to note that Bogert et, al. recognized the limitations of their technique and the difficulties encountered in the multiple echo determination problem. Although first used in geophysical applications, the cepstrum has made its way into other fields of research, for example, speech and sonar signal processing. Schafer (1969) was concerned with the recovery of signals generated by a convolutional process as opposed to the determination of echo arrival times and called his method homomorphic deconvolution. Through a logarithmic (homomorphic) transformation of the z-transform of the observed process, convolution of two signal components is mapped into the addition of their kepstra. By appropriate filtering (sometimes called liftering) in the kepstral domain, the desired signal can be recovered. Schafer used the term “complex cepstrum” to account for the fact that both the magnitude and phase spectra of the observed signal are utilized. To date, Oppenheim and Schafer (1975), Tribolet (1977), and many others have studied the effectiveness of the “complex cepstrum” in deconvolving reflection seismograms. Our word “kepstrum” means the same as their term “complex cepstrum”. Because the kepstrum of a real-time sequence is real, the use of the word “kepstrum” is less confusing than the term “complex cepstrum”. In summary, we have traced the mathematical evolution of the kepstrum, starting with the logarithmic potential problem of the 1800’s to the current deconvolution process in oil exploration. In the next section, we shall review the normal-incidence reflection seismogram so as to link the physics of reflection seismology to the mathematical properties of the kepstrum.
60
THE ARMA MODEL USED IN REFLECTION
SEISMOLOGY
It is important that we review the evolution of the well known normalincidence convolutional model used as a basis for many deconvolution techniques. This model in its final form is simple and linear, and as a result, can be easily evaluated by conventional computational methods. However, without understanding the basic underlying assumptions and physical approximations involved in the model development, one cannot view this simplicity in the proper perspective. If we are going to use the kepstrum as a deconvolution technique, let us begin by examining its relationship to the physics of reflection seismology. Reflection seismology is a method of mapping the subsurface structure of the earth. A source at or near the surface is set off. The source may be a dynamite explosion, a weight-dropping machine, a vibration-inducing machine, or some other means to transmit seismic energy into the earth. The seismic waves travel from the source to subsurface layers within the earth where they are reflected. The returning seismic waves are then recorded by instruments on the surface. If the entire process is assumed to be linear, then the surface source is the excitation, the reflection seismogram is the output or set of observables, and the earth is the governing linear system. Part 1. Transmission of a deep source through the earth Consider an inhomogeneous system (the earth) bounded by (sandwiched in) two homogeneous infinite half-spaces, the air and the basement rock. Fig.1 depicts the situation. Now in earthquake seismology, the source is considered to be embedded deep within the earth’s subsurface. The classical approach to determining the response of the earth to a deep source excitation was to solve the governing partial differential equations, which described the wave propagation through the earth into the air, subject to the initial-value and boundaryair
Velocity Profile Velocity
‘, distributed system
earth
i
,y depth
source
Basement Rock Fig.1. The inhomogeneous earth excited by a deep source and considered as a distributed parameter system. The associated velocity profile of this system is considered as a continuous function of depth.
61
value conditions. Hence, due to the presence of the governing partial differential equations, the earth was considered as a distributed-parameter system. Research in classical seismology is concerned with the analytic solutions of these partial differential equations, and except in the simplest cases, such solutions are mathematically extremely difficult to find. Consider now the idea of modeling the earth by a lumped parameter system. For example, we could quantize the continuous velocity profile (distribution) shown in Fig.1. This quantization is equivalent to considering the earth as a multi-layered system, with the velocity of propagation in each layer and the depth of a layer determined by the quantization procedure. A description of the earth as a lumped-parameter system is given in Fig.2. (velocity profile) (output) air earth
1
Lumped parameter system (N-layers)
J
(input) Basement Rock Fig.2. The multi-layered earth excited by a deep source
and considered
as a lumped-
parametersystem. The continuous velocity profile is now represented as a discontinuous
profile.
Now if the time of signal propagation through a layer is short compared with the duration of the signal, then the lumped-parameter assumption is valid. By choosing the depth of each layer to be very small, i.e., considering many layers, we can satisfy the lumped-parameter conditions. In terms of transmission lines and linear network theory, Fig.1 represents a distributed circuitelement transmission line while Fig.2 represents a lumped circuit-element transmission line. Note that the model for the lumped-parameter system contains a discrete velocity profile distribution, with each discontinuity in the velocity representative of one transmission line section. Thus, the lumpedparameter system is modelled by many sections of transmission lines. Now if we further assume that the system of Fig.2 is linear and timeinvariant, we can expect that the input (en) and output (x,) satisfy the linear difference equation: Xn + UlXn-1
+ . . . UNX,_N
= E,
(8)
where the an’s represent the difference equation operator (i.e., constant parameters associated with our lumpedIparameter system which has N degrees of freedom represented by the N layers). Our reason for considering a lumped-
62
parameter system with the difference eq.8 is that difference equations are easily solved by digital computers. Thus, the difficulty of analytic solutions of wave motion is eliminated by the numerical solution of eq.8. On physical grounds, we know that eq.8 represents a stable system, i.e., for every bounded input we observe a bounded output. Now given the sequences x, and e,, we shall define their respective z-transforms by:
X(z) =
_c
x,zn,E(z)
n=--m
=
c
En9
n =-ca
With the definition in eq.9, the z-transform of each side of eq.8 gives: X(z) [l + a,2 + a*2 + . . . + a#]
= E(z)
or: 1 X(z) = E(z) 1+a,z+a,z* +...+a&v
(10)
Eq.10 represents the transfer function of an “all-pole” linear time-invariant system. Our physical claim that the lumped-parameter model represents a stable system is equivalent to the mathematical condition that the transfer function X(z)/E(z) contains no poles inside the unit circle. Let us now formulate this property in terms of the difference equation operator (1, a,, a,, . . , UN). The z-transform of this operator yields the denominator of eq.10. Thus, the property that the transfer function X(z)/E(z) has no poles inside the unit circle is equivalent to the property that: 1+cIz+c2z2
+...
+cNzN#Ofor(zI
(11)
We can therefore state that the operator (1, a,, a,, . . ., UN) yields a stable difference equation if and only if its z-transform has no zeroes inside the unit circle. Now if we replace the input to our lumped-parameter system by a unit impulse, i.e., let E, = 6, where: l,n=O 6, = i O,n#O then the impulse response, denoted by b,, is a stable time sequence. This condition follows from physical reasoning. But since:
B(z)=
1
1 +a12 +a$* +...+aNzN
(12)
it follows that B(z)has no zeroes or poles inside the unit circle. Now b, is the
63
sequence which is transmitted through the earth into the air when a unit impulse (deep source) excites the earth. From actual observation of field seismic data, b, has the bulk of its energy concentrated at its beginning and is rapidly attenuated due to successive reflections and refractions within the earth. Sequences that possess this “front-loaded” property are called minimum-delay sequences. Now if A (z) = 1 + a,z + a22 + . . . + a~$ is the z-transform of the operator (1, a,, a,, . . ., a~), then from eq.12 we observe that: B(z)A(z)
=1
or: b,*a,
= 6,
which shows that the operator a, and transmission impulse response b, are mutually inverse sequences. But since uR and b, are mutually inverse sequences, then we can conclude that a, is a stable difference equation operator if and only if b, is minimum-delay. That is, from purely physical grounds, we arrive at the important result: stable difference equation operator a,
*
minimum-delay response b,
Because b, represents all the multiple reflections and refractions resulting from an impulsive input 6,, b, also represents the system reverberation; b, = transmission (impulse) response = layered system reverberation wavelet. Thus, aR can also be viewed as the deconvolution operator which removes the layered system reverberation. Part 2. Reflections caused by surface source (internal primary reflections) Next, let us consider the N-layered (lumped-parameter) system of Fig.2 with the “deep” upgoing source replaced by a downgoing unit impulse (source) initiated at the surface. This situation, along with a pictorial definition of an internal primary reflection, is described in Fig.3. Now in Part 1, we claimed that an upgoing unit impulse (source) incident on the lowest interface gave rise to the transmitted (i.e., system. reverberation) waveform 13,. If we consider this situation in the context of linear timeinvariant system theory, then the unit impulse 6, is the cause and the system reverberation b, is the effect. Now if the cause is delayed by m units of time, i.e., 6,_, , then the effect is also delayed by m units of time, i.e., b,_,. Also, if the cause is scaled by a constant C, i.e., C6,, then the effect is also scaled by the constant C, i.e., Cb,. These properties, together with the property of superposition, characterize linear time-invariant systems. Thus, instead of considering the downgoing unit impulse (source), let us replace it by a set of hypothetical (i.e., mathematical) upgoing sources at the reflecting interfaces of our N-layered system. Since we know the “effect” b, due to the upgoing
64
r air
/
1st primary reflection
/
/
r2nd I
downgoing surface
primary reflection rNth /
primary reflection
N basement rock Fig.3. Internal primary reflections caused by a downgoing unit-impulse applied at the surface. For clarity, ray paths are drawn at oblique incidence but wave-motion analysis is for normal incidence.
“cause” 6,) we can then apply the linear time-invariant system properties to deduce the effect due to a downgoing unit impulse at the surface. Fig.4 shows the downgoing source impulse replaced by upgoing sources of strength E,, = r, where r, is the reflection coefficient of interface n. From physical experience gained in seismic prospecting, it is known that generally these reflection coefficients are small in magnitude, that is, I E, ( = I r, I< 0.1. Each of these hypothetical sources gives rise to a reverberation wavelet due to the layering. Further, to a first approximation, each reverberation wavelet is the same shape, namely the system reverberation wavelet b, . Silvia (1977) examines this assumption in some detail and shows the mathematical conditions under which it is a good approximation to a more exact stratified earth model. Taking into account the time delays associated with the source locations (see Fig.4), we reason as follows: transmission impulse-response b, M If an upgoing unit spike 6, (gives rise to)
then e161t--l %&-2
fN&-N
-+ --,
El&-1 EzL.--2
-+
%‘bn--N
Thus, using the superposition property of linear systems, we can add the responses due to all the hypothetical sources to obtain the total response h, Doing so, we get: h, = el b,_l
= +,*bn
+ ezbn+
.
+ . . . + Wbn-N
(13)
65
Eq.13 represents the reflection seismogram resulting from a downgoing unit impulse initiated at the surface. Further, h, is the resultant of all the system reverberations (reverberation wavelets b,) due to the upgoing hypothetical 0
\
\
‘\
El&r \
1
Single impulsive source of unit strength
\
Ez+-
\ =
2
\
N
Multiple impulsive sources of N-Layered strengths system ! El,
\ \
0, n#O
2
\ \
.
1, n=O 6 = n
1
\/
rEN=rN
E2,..EN
J
I
Fig.4.The downgoing unit impulse surface source 6, replaced by hypothetical mathematical) upgoing impulsive sources of strength ek i rk .
(i.e.,
sources E, with strength given by the reflection coefficient sequence r,. The z-transform of eq.13 gives: H(z) = E(z)B(z)
(14)
where: E(z) =
E1Z +
f*2 + . . . + eNzN
(15)
Substitution of eqs.12 and 15 into 14 yields: H(z)
=
El2
+
f22* + . . . + ENzN
1 +alz+a2z2
+.
.+aNzN
(16)
which gives the transfer function of a normal-incidence reflection seismogram and is the ARMA model used in reflection seismology. It is important to note that the same result can be obtained by deriving the exact expression for an N-layered earth model then making the approximation that all the reflection coefficients are small in magnitude, i.e., I r, I < 0.1. Under this approximation, fl G rl, e2 A r2,. . . , EN A rN; that is, the hypothetical sources are exactly the reflection coefficients (Silvia (1977)). Since the layered system is assumed to be both linear and time-invariant, then the reflection seismogram x, due to an arbitrary source pulse s, is: x, = h,*s, = En”b,“S, =
= En*( b,*+)
(17)
% *Wl
where w, = bn*sn is defined as a composite wavelet consisting of the reverberation wavelet b, and source pulse s, . Thus, eq.17 describes the normal-
66
incidence reflection seismogram, where b, represents the autoregressive component and s, *E, the moving-average component of the reflection seismogram. The basic deconvolution problem is to filter X, such that we can best recover the reflection coefficient sequence E, . Before we discuss the role of the kepstrum in deconvolution, let us investigate some of its basic properties. PROPERTIES
OF THE KEPSTRUM
For a minimum-delay system H,h(z), the hepstrum efficients of the Taylor-Maclaurin series: log H&(Z)
= F(z) = 5
fn is defined
f+z” , I2 I < 1
as the co-
(18)
n=O
However, the concept of the kepstrum is not limited to minimum-delay systems. Let us now consider a nonminimum-delay system H(z). If we assume that H(z) has no zeroes or poles in an annular region including the unit circle, then log H(z) is analytic in this annular region. As a result, log H(z) has a Laurent expansion: logH(2)
=
5 f,z” n=-ca
(19)
in this annular region. The coefficients for this Laurent series are the coefficients of a stable noncausal system. We call these coefficients the kepstrum of the nonminimum-delay system H(z). Let us now return to the case of a minimum-delay system:
5
Hmin(Z) =
hn.9
(20)
n=O
From eqs.18 and 20 we obtain: &i&r)
= exp
[JY~)l
or: hnzn --exP
2
[
jjofdn]
(21)
n=O
A formula
for obtaining
the minimum-delay
sequence
h, from the kepstrum
fn can be obtained as follows. Taking the derivative of eq.21 with respect to z gives:
& (jjo hnz") =exP[ j. fnzn]
2 (jjo fnzn)
(22)
6’7
Dividing eq.22 by eq.21 we obtain:
(23) n=O
Expanding eq.23 we get:
5
(n + l)hn+lZn
=
(C h,a”)( 2
(n +
n=O
n=O
n=O
l)f,+lZn)
(24)
The inverse z-transform of eq.24 gives the recursive relation (Oppenheim and Schafer, 1975):
(n + l)bn+l
= 5
hk( 11 + 1 -
k)fn+l-k for n = 0, 1, 2,
.. .
(25)
k=O
which is useful for extracting the minimum-delay sequence ha from the kepstrum fn. For IZ= 0, 1,2, 3, . . . eq.25 gives: n = 0:
h, = hof,
n=l:
h2 =
hof, +
2 fl =Mf,
n = 2: h3 = hOf3+ ihlfz
+$r:,
(
+ th2fl = ho.f3 + flf2
Let us now interpret the kepstral coefficient z = 0 that: ho = exp
1 +-f:6
fo. We have from eq.21 with
(fo) or f,, = log ho
Thus, it is often convenient to normalize the sequence h, so that its leading coefficient is unity, i.e., ho = 1; thus, the leading kepstral coefficient is zero, i.e., f. = 0. Let us now examine some of the properties of the kepstrum. At this point it is convenient to introduce the symbols a, and b, for the time sequences and
68
the corresponding Greek symbols (Yeand & for their respective kepstra. Thus, the kepstrum of a minimum-delay sequence a, with leading coefficient a0 = 1 = 0. If we let L(a,) is the causal sequence al, a2, as, . . . with leading term CY~ denote the z-transform of u, , then the kepstrum can be described by the composition of transformations L, log, and L-’ given by: % = L-l {log
(26)
[Jm,)lI
Hence, to find the kepstrum of a minimum-delay sequence h, whose leading coefficient is unity, we compute the z-transform of the sequence, take the natural logarithm of this transform, and find the inverse z-transform of the quantity log H(z). For convenience, we shall define the kepstral operator K given by the successive transformations L, log and L-’ such that: K(*) = L-l {log [L(m)]}
(27)
where K(h,) = fn and ho = 1 Let us now consider various properties of the kepstrum, i.e., convolution, reverse sequence, inverse sequence, and reverse-inverse sequence. Convolution. Given the minimum-delay sequences a, and b, whose leading coefficients are unity, then:
~@n*bl) = K(%) + K(b)
(28)
where: K(G)
=O,%,~2,%,..*
KC&l)
=
O,Pl,P?,PJ,
* * *
Hence, the kepstrum maps convolution into addition. To see this, we write: L(unfbn)
= A (z)B(z)
Then: log [L (a,*&
)I
= 1% [A WW)l = log A(z) + log B(z)
Taking the inverse z-transform of the above result we get: L-’ {log [L(un*bn)]
} = L-’ [log A(z)]
+ L-l [log B(z)]
or:
Reverse sequence.
Given the minimum-delay sequence b, with leading coeffi-
69
cient bO = 1, we form the reverse sequence b_, by folding b, about theorigin. The resulting sequence is called minimum-advance. Hence: K(b-n)
(29)
= P-n
We can show this result by first noting that: L(b,)
=
5
b,P =B(z)
n=O
and: L(b_,)
b,z_” =&z-l)
= jj n=O
Recall that: logB(2)
=&2+&z’
+P3z3 +. . .
Then: logB(C)
=&z-i
+pzO
+fi3z-3 + . . .
and we observe that given K(bn) = &, K(b_,) origin. In words:
is found by folding
Pn about&
minimum-delay sequence ++ causal kepstrum minimum-advance sequence tf purely noncausa2 kepstrum Inverse sequence. Given the minimum-delay sequence b, with leading coefficient b,, = 1, we define the inverse sequence b;’ by the relation: l,n=O b,*b,’
= 6, = O,n#O
Then : K(b,‘)
= -&,
(30)
Since L(bn*b;‘)
= 1, then B(z)L(b;‘)
or: L(b,‘)
= &-)
=
P(~)l-’
Now,
log [L(b,‘)]
= log [B(z)]-’
= -log B(z)
=1
70
Thus, it follows that: K(b,‘)
= -on
In words: minimum-delay sequence ++ causal kepstrum inverse of minimum-delay sequence ++ negative of causal kepstrum Reverse-inverse sequence. Given the minimum-delay sequence b, with leading coefficient b,, = 1, we define the sequence bZk as the inverse of the reverse sequence. Then: K(bI;)
= -/I+
(31)
To see this, we write: L(b$)
= [&z-l)]-’
Then: log [L(bz;)]
= log [B(P)]-’
= -log B(z_‘)
From the reverse-sequence property we can state: K(bZh) =-/3-n In words: minimum-delay sequence tf causal kepstrum reverse of inverse of minimum-delay sequence ++ negative of purely noncausal kepstrum The above properties bring out the essential features and symmetries of the kepstrum. Thus, given the kepstrum on of a minimum-delay sequence b, with leading coefficient b,, = 1, we can find the kepstrum of the reverse sequence b-,, inverse sequence b;’ and inverse-reverse sequence b1; by simple negations and reflections of the kepstrum on. For additional kepstral properties, see Oppenheim and Schafer (1975) and Tribolet (1977). In the above discussion, we developed the relation (eq.24): 2 (n + l)hn+lzn n=O
= ( jj
h,z”)
n=O
(
5
(n + l)f,+#)
n=O
which can be written as: fj n=O
(n + l)fn+lzn =
&) jj
(n + l)hn+lzn
n-0
where H(z) = 5 n=O
hnzn and ho = 1
(32)
71
The inverse z-transform K(h,)
= fn+l
=
of eq.32 gives:
5
hk+lh$-k
(33)
k=O
for n = 0, 1, 2, . . . and h,_,= 1 which explicitly gives the kepstrum f, in terms of the minimum-delay sequence h,. Let us now relate the mathematical properties of the kepstrum to the physics of the normal-incidence seismogram by next discussing the topic of deconvolution. THE ROLE
OF THE KEPSTRUM
IN DECONVOLUTION
Recall that the reflection seismogram x, can be approximately the output of an ARMA model with transfer function: H(z)
=-
E(z) A (2)
El2
=
1 +
+
U,Z
regarded
as
+Z2 + . . . + fNz N +
UzZZ + . . +aNZN
where E, = the reflection coefficient sequence and is the desired output of a deconvolution filter. Unfortunately, this system is not always minimum-delay, since the numerator polynomial coefficients fl, e2, . . . , EN are the reflection coefficients of an N-layered system and are dependent upon a particular geographic area. Thus, we are not always sure if all the zeroes of E(z) lie outside the unit circle, and in practice, the reflection coefficient sequence is generally mixed-delay, i.e., it contains zeroes both inside and outside the unit circle. However, the reverberation sequence b,, where B(z) = l/A (z) is always minimum-delay. If we excite the earth with an arbitrary mixed-delay source pulse sn, then the reflection seismogram x, is the convolution: x,
=sn
*en*bn
(34)
= En*Wn where w, = s, *b, is the composite mixed-delay wavelet. The method of predictive deconvolution (Robinson (1954)) assumes that the reflection coefficient sequence is a purely random sequence with a flat power spectrum. This assumption is usually met in practice (White and O’Brien, 1974). Such an assumption allows the method of predictive deconvolution to design least-squares deconvolution (noise-whitening) filters. From the convolution property of the last section, the kepstrum K (x, ) of the reflection seismogram is: K(x,)
=
K(E,*s,*b,)
= in +an
(35) +
pn
72
where K(E,) = &, K(s,) = a, and K(&) = &. The method of homomorphic deconvolution is concerned with processing the kepstrum K(x,) so that uII (source pulse kepstrum) and + (reverberation kepstrum) may be eliminated, leaving the desired kepstrum E, . Applying the inverse kepstral operator K-’ to the “filtered” kepstrum & would then yield the desired reflection coefficient sequence. It appears that deconvolution via the kepstrum is less restrictive than the method of predictive deconvolution, where one assumes that the reflection coefficients are an tumor-relatedsequence and that the source pulse be minimumdelay in order to achieve “optimum” results. Further, the kepstrum appears to be a more robust processor in the sense that the technique is not tied to any physical attributes of the model and that a minimum-delay source pulse assumption is not needed. However, both the method of predictive deconvolution and the method of kepstral deconvolution require an assumption of the phase spectrum of the source pulse (Treitel and Robinson, 1977). For example, let the canonical representation of the source pulse be: s, =pn*sp
where pn is a causal all-pass operator, ~$0)is the minimum-delay source wavelet, and s, is the nonminimumdelay source wavelet. The reflection seismogram for a nonminimum-delay source pulse can then be expressed as: xll
= f,
*b,*p,*s$f)
(37)
If we knew the all-pass system pn , i.e., if we had knowledge of the source phase characteristic, then we could operate on the seismogram with the inverse all-pass p;’ to yield: -l*
Pn
Xn = 62 *bn*s,(0)
which is more suitable to the method of predictive deconvolution. Often, the source pulse kepstrum u, is clustered around the kepstral origin (also called “low-time” portion of the kepstrum). See Stoffa et al. (1974) and Ulrych (1971). By a suitable window with the proper cutoff “quefrencies”, the source pulse can be theoretically recovered. In essence, the selection of the cutoff quefrencies determines how well we can estimate the source pulse phase and magnitude characteristics. Thus, both deconvolution techniques require some knowledge of the phase characteristic of the source pulse. CONCLUSIONS
The idea of the kepstrum appeared in the classical works of Poisson and Schwarz concerning the logarithmic potential problem. Szego and Kolmogorov used the kepstrum to factor the power spectrum of a non-deterministic process.
73
Because the normal-incidence reflection seismogram is approximated by a convolutional process, the kepstrum offers a solution to the deconvolution of reflection seismograms. We developed additional properties of the kepstrum which bring out its essential features and symmetries. As a result, the kepstrum evolves as an important mathematical construct in the theory of signal analysis and time series, and takes the place among such established functions as the spectrum and autocorrelation.
REFERENCES Bateman, H., 1944. Partial Differential Equations of Mathematical Physics. Dover, New York. B&h, M., 1968. Mathematical Aspects of Seismology. Elsevier, Amsterdam. BLth, M., 1974. Spectral Analysis in Geophysics. Elsevier, Amsterdam. Bogert, B.P., Healy, M.J.R. and Tukey, J.W., 1963. The quefrency analysis of time series for echoes. In: M. Rosenblatt (Editor), Proc. Symp. Time Series Analysis. Wiley, New York, pp.209-243. Kolmogorov, A.N., 1939. Sur l’interpolation et extrapolation des suites stationnaires. C.R. Acad. Sci. Paris, 208: 2043-2045. Kulhanek, O., 1975. Introduction to Digital Filtering in Geophysics. Elsevier, Amsterdam. Oppenheim, A.V. and Schafer, R.W., 1975. Digital Signal Processing. Prentice-Hall, Englewood Cliffs, N.J. Poisson, S.D., 1823. Sur la distribution de la chaleur dans les corps solides. J. EC. R. Polytech., Ser.1, 19: l-162. Robinson, E.A., 1954. Predictive Decomposition of Time Series with Applications to Seismic Exploration. Thesis, MIT, Cambridge, Mass. (Also in Geophysics, 32: 418-484.) Schafer, R.W., 1969. Echo removal by discrete generalized linear filtering. Res. Lab. Electron. MIT, Tech. Rep., 466. Schwarz, H.A., 1872. Zur Integration der partiellen Differentialgleichung. J. Reine Angewandte Math., pp.218-254. Silvia, M.T., 1977. Deconvolution of Geophysical Time Series. Thesis, Northeastern University, Boston, Mass. Stoffa, P.L., Buhl, P. and Bryan, G.M., 1974. The application of homomorphic deconvolution to shallow-water marine seismology - Part I. Models. Geophysics, 39: 401-416. Szego, G., 1915. Ein Grenzwertsatz iiber die Toeplitzschen Determinanten einer reellen positiven Funktion. Math. Ann., 76: 490-503. Treitel, S. and Robinson, E.A., 1977. Deconvolution -homomorphic or predictive? IEEE Trans. Geosc. Electron., GE-15(l): 11-13. Tribolet, J., 1977. Seismic Applications of Homomorphic Signal Processing. Thesis, Dept. of Electr. Eng. Comput. Sci., MIT, Cambridge, Mass. Ulrych, T.J., 1971. Application of homomorphic deconvolution to seismology. Geophysics, 36: 650-660. White, R.E. and O’Brien, P.N.S., 1974. Estimation of the primary seismic pulse. Geophys. Prosp., 22: 627-651.
55
USE OF THE KEPSTRUM IN SIGNAL ANALYSIS
MANUEL T. SILVIA’
and ENDERS
1Naval Underwater ‘Northeastern
A. ROBINSON2
Systems Center, Newport, R.I. 02840 University, Boston, Mass. 02115 (U.S.A.)
(Received September
(U.S.A.)
3,1977)
ABSTRACT Silvia, M.T. and Robinson, E.A., 1978. Use of the kepstrum in signal analysis. In: C.H. Chen (Editor), Computer-Aided Seismic Analysis and Discrimination. Geoexploration, 16: 55-73. The kepstrum has evolved as an important mathematical construct in the area of digital signal processing. In particular, kepstral analysis has been used for the deconvolution of reflection seismograms generated by the exploration for oil and natural gas. Here, we review the development of the normal-incidence reflection seismogram from a physical point of view and show how the kepstrum relates to the fine-grain structure of the physical model. Moreover, we develop symmetry properties of the kepstrum, introduce the concept of minimum-advance, and discuss the relationship between spectral factorization and the kepstrum. In the context of deconvolution, we discuss the methods of predictive deconvolution and homomorphic deconvolution and their dependence on the physical model of the reflection seismic process.
INTRODUCTION
The idea of the kepstrum appears in the classical work of Poisson (1823), Schwarz (1872), Szego (1915), and Kolmogorov (1939), and has been applied to geophysical problems by Robinson (1954), Bogert et al. (1963), Schafer (1969), Oppenheim and Schafer (1975), Tribolet (1977), and others. In this paper we bring all this work together and introduce additional properties of the kepstrum so as to bring out its essential features and symmetries. As a result, the kepstrum evolves as an important construct in the theory of signal analysis and time series, and takes the place among such established functions as the autocorrelation and spectrum. For further discussion see Beth (1968, 1974) and Kulhanek (1975). In their classical works on potential theory, Poisson (1823) and Schwarz (1872) addressed the problem of determining a potential function whose real part was an assigned value on the unit circle. As is well known in potential theory (see Bateman, 1944), Schwarz’s classical expression is:
56
. ,
H(z) = ,‘n
1
(:
‘~~~~,)
Re[H(w’)]dw’,
lz I < 1
(1)
with w’ as the integratioh variable and z = rewiw . This equation gives the potential function H(z) whose real part on the unit circle is Re [H(o)]. The Schwarz expression follows from the famous result of Poisson on potential theory. An interesting result arises when we consider a logarithmic potential log H(z), which can be written: log H(z) = log I H(z) I + i arg [H(z)] The real part of the logarithmic potential evaluated on the unit circle z = eAiw is then log 1H(z = eszW) I = log I H(w) I. In terms of the logarithmic potential, Schwarz’s expression becomes: .
log H(z) = $
_i:
(: + ;$)
I
log I H(o’)
I dw’, I z I < 1
where log I H(w) I is the assigned or known value on the unit circle and the function H(z) has no zeroes or poles inside or on the unit circle, i.e., H(z) is minimum-delay. (This definition of minimum-delay is based on the Laplace definition of the z-transform, i.e., the z-transform of a sequence h, is defined by H(z) = 5 h,zn .) As a result, log H(z) has a Taylor-Maclaurin series exn=-m pansion in the unit circle: logH(z)=
cf,z”,lzl
(2)
n=O
The coefficients fn of this Taylor-Maclaurin series are the coefficients of a stable causal system. We call these coefficients the kepstrum of the minimumdelay system H(z). However, as we shall see, the concept of the kepstrum is not limited to minimumdelay systems. In addition to its appearance in the logarithmic potential problem, the kepstrum was used by Szego (1915) and Kolmogorov (1939) to solve the problem of factoring the power spectrum @( w ) of a random process in order to extract a stable causal system H(z). A power spectral density function (a(w ) must satisfy the following conditions in order to be factored so as to yield a causal system: (1) Q(w) is positive on the interval-n < w < 71. (2) @(w )has finite area on the interval --71< o < n. (3) log @((w.)has finite area on the interval --71< o < n. Now the problem of spectral factorization is the problem of extracting a causal system H(z) whose magnitude spectrum is the square root of the power spectrum. In symbols:
57
I H(w) I = [Q(w)1
‘/i
Szego (1915) and Kolmogorov (1939) recognized that a solution to this problem could be obtained by making use of the Schwarz expression (1). Now log I H(w) I is the real part of the function log H(w). Further: Re [log N(w)] = log I H(w)
I = ?h log (a(o)
Thus, the required solution is found by merely substituting ‘/z log (a(w) for Re [H(w)] and substituting log H(z) for H(z) in Schwarz’s expression to obtain: . H(z)
log
=
I
k i (‘,y;fi:.);log
@(w’)dw’,
Iz i < 1
-lT . . 1s.
Hence, the Szego-Kolmogorov
= exp [log H(z)] = exp [&
H(z)
solution to the spectral factorization problem
i (:r;!;:) -77
logrp(,‘)da$
:z I< 1
The desired sequence h,,, hl, h2, . . . can be written in terms of the kepstrum fo, fl, f2, - * * as: h,+hIz+h2z2+...=exp[fO+fIz+f2z2+...]
We call this equation Kolmogorou’s that: ho =
(3) equation.
Putting z = 0 in eq.3 we observe
exp(fo).= exp
Thus, the leading kepstral coefficient fO is the logarithm of the leading impulse response coefficient ho, i.e., fo = log h,
As a result, it is often convenient to normalize an impulse response function that its leading coefficient is unity; thus, the leading kepstral coefficient is zero, i.e., fO = log 1 = 0. With this normalization, Kolmogorov’s equation (eq.3) is rewritten as: h, so
1+;2.+;.*+...= 0
exp (fIz + f2z2 + f3z3 + . . .)
(4)
0
We shall define the power series expansion flz + fizz + f3z3 + . . . appearing in the exponent of Kolmogorov’s equation (4) as the Kolmogorov Equation Power Series (KEPS). Further, since KEPS is the z-transform of the sequence 0, f1, f2, f3, * * -9we can think of fn as being a time response. Thus, we shah de-
58
fine the sequence 0, fI, fi, f3, . . . as the Kolmogorov Equation Power Series Time Response and call this sequence the KEPSTRUM, where we have added the Latin singular ending “urn” to denote one kepstrum and shall use the Latin plural ending “a” to denote more than one kepstrum, i.e., kepstra. Whereas Szego and Kolmogorov were only concerned with the real part of the function log H(z) on the unit circle, Robinson (1954) was concerned with the imaginary part as well; namely, both the real and the imaginary parts have physical meaning. Thus, it is important to consider the logarithm of the spectrum H(w). Let us elaborate. For a minimum-delay system H(z), the kepstrum fn was defined by the relation:
logH(z)=g f,z~,lzl
The complex variable z represents a point within the unit circle (i.e., I z I < 1). However, we can consider the limit of log H(z) as the interior point z approaches a point on the circumference of the unit circle; that is, we consider the limit of log H(z) as z-teTZW. This limit is: log H(w) = 5
f,eminw
n=O
Separating the above into real and imaginary parts, we obtain: logIH(w)I+iOo-(f~+Re[
‘&fneminw])
+i
(Im[~lfne-i”w])(5)
By equating real parts in eq.5, it can be shown that: fit = 2fk for n = 1, 2, 3,. . . f. = f;
for n = 0
f-“, = 2fl, for n = -1, -2, -3, . . .
where : f,‘, = &
i log IH(w) I einw dw for n = 0, +l, +2, . . . -A
(6)
and the superscript * denotes the complex conjugate, i.e., f,” is the complex conjugate of fn. Thus, the kepstrum fn is found in terms of log I H(w) I, thereby allowing us to obtain the minimum-delay system H(z) via (2). Let us now look at the imaginary part of eq.5, namely, by equating the imaginary parts in (5) it can be shown that:
59
e(o)
=zpj cot
(yj
log I H(d)
I dw’
(7)
--n where we have incorporated the result of (6) into (7); the symbol P denotes that the integral has its Cauchy principal value. The kepstrum was used as an intermediate step in deriving (7), which shows that for a minimum-delay system H(z), knowledge of only the log magnitude spectrum is sufficient to determine the phase spectrum 0 (w). In the mathematics literature, eq.7 and various versions of the Schwarz formula (1) are often referred to as Poisson’s formulas. In the context of linear system theory, eq.7 is known as the Hilbert transform, where e(w) and log I H(w) I form a Hilbert transform pair. Bogert et al. (1963) were the first to use the kepstrum as a signal processing operation. Bogert et al. were primarily interested in measuring the time interval between the arrival of certain seismic waves, excited by either natural or man-made disturbances. They defined the term “cepstrum” as the power spectrum of the logarithm of the power spectrum of the observed time series. We see that the “cepstrum” of Bogert et al. involved no phase spectrum information, i.e., it utilized magnitude spectrum and power spectrum information only. Their “cepstrum” was effective in detecting or measuring the time delay of a single echo, even when other methods (such as the autocovariance method) were ineffective. However, it is important to note that Bogert et, al. recognized the limitations of their technique and the difficulties encountered in the multiple echo determination problem. Although first used in geophysical applications, the cepstrum has made its way into other fields of research, for example, speech and sonar signal processing. Schafer (1969) was concerned with the recovery of signals generated by a convolutional process as opposed to the determination of echo arrival times and called his method homomorphic deconvolution. Through a logarithmic (homomorphic) transformation of the z-transform of the observed process, convolution of two signal components is mapped into the addition of their kepstra. By appropriate filtering (sometimes called liftering) in the kepstral domain, the desired signal can be recovered. Schafer used the term “complex cepstrum” to account for the fact that both the magnitude and phase spectra of the observed signal are utilized. To date, Oppenheim and Schafer (1975), Tribolet (1977), and many others have studied the effectiveness of the “complex cepstrum” in deconvolving reflection seismograms. Our word “kepstrum” means the same as their term “complex cepstrum”. Because the kepstrum of a real-time sequence is real, the use of the word “kepstrum” is less confusing than the term “complex cepstrum”. In summary, we have traced the mathematical evolution of the kepstrum, starting with the logarithmic potential problem of the 1800’s to the current deconvolution process in oil exploration. In the next section, we shall review the normal-incidence reflection seismogram so as to link the physics of reflection seismology to the mathematical properties of the kepstrum.
60
THE ARMA MODEL USED IN REFLECTION
SEISMOLOGY
It is important that we review the evolution of the well known normalincidence convolutional model used as a basis for many deconvolution techniques. This model in its final form is simple and linear, and as a result, can be easily evaluated by conventional computational methods. However, without understanding the basic underlying assumptions and physical approximations involved in the model development, one cannot view this simplicity in the proper perspective. If we are going to use the kepstrum as a deconvolution technique, let us begin by examining its relationship to the physics of reflection seismology. Reflection seismology is a method of mapping the subsurface structure of the earth. A source at or near the surface is set off. The source may be a dynamite explosion, a weight-dropping machine, a vibration-inducing machine, or some other means to transmit seismic energy into the earth. The seismic waves travel from the source to subsurface layers within the earth where they are reflected. The returning seismic waves are then recorded by instruments on the surface. If the entire process is assumed to be linear, then the surface source is the excitation, the reflection seismogram is the output or set of observables, and the earth is the governing linear system. Part 1. Transmission of a deep source through the earth Consider an inhomogeneous system (the earth) bounded by (sandwiched in) two homogeneous infinite half-spaces, the air and the basement rock. Fig.1 depicts the situation. Now in earthquake seismology, the source is considered to be embedded deep within the earth’s subsurface. The classical approach to determining the response of the earth to a deep source excitation was to solve the governing partial differential equations, which described the wave propagation through the earth into the air, subject to the initial-value and boundaryair
Velocity Profile Velocity
‘, distributed system
earth
i
,y depth
source
Basement Rock Fig.1. The inhomogeneous earth excited by a deep source and considered as a distributed parameter system. The associated velocity profile of this system is considered as a continuous function of depth.
61
value conditions. Hence, due to the presence of the governing partial differential equations, the earth was considered as a distributed-parameter system. Research in classical seismology is concerned with the analytic solutions of these partial differential equations, and except in the simplest cases, such solutions are mathematically extremely difficult to find. Consider now the idea of modeling the earth by a lumped parameter system. For example, we could quantize the continuous velocity profile (distribution) shown in Fig.1. This quantization is equivalent to considering the earth as a multi-layered system, with the velocity of propagation in each layer and the depth of a layer determined by the quantization procedure. A description of the earth as a lumped-parameter system is given in Fig.2. (velocity profile) (output) air earth
1
Lumped parameter system (N-layers)
J
(input) Basement Rock Fig.2. The multi-layered earth excited by a deep source
and considered
as a lumped-
parametersystem. The continuous velocity profile is now represented as a discontinuous
profile.
Now if the time of signal propagation through a layer is short compared with the duration of the signal, then the lumped-parameter assumption is valid. By choosing the depth of each layer to be very small, i.e., considering many layers, we can satisfy the lumped-parameter conditions. In terms of transmission lines and linear network theory, Fig.1 represents a distributed circuitelement transmission line while Fig.2 represents a lumped circuit-element transmission line. Note that the model for the lumped-parameter system contains a discrete velocity profile distribution, with each discontinuity in the velocity representative of one transmission line section. Thus, the lumpedparameter system is modelled by many sections of transmission lines. Now if we further assume that the system of Fig.2 is linear and timeinvariant, we can expect that the input (en) and output (x,) satisfy the linear difference equation: Xn + UlXn-1
+ . . . UNX,_N
= E,
(8)
where the an’s represent the difference equation operator (i.e., constant parameters associated with our lumpedIparameter system which has N degrees of freedom represented by the N layers). Our reason for considering a lumped-
62
parameter system with the difference eq.8 is that difference equations are easily solved by digital computers. Thus, the difficulty of analytic solutions of wave motion is eliminated by the numerical solution of eq.8. On physical grounds, we know that eq.8 represents a stable system, i.e., for every bounded input we observe a bounded output. Now given the sequences x, and e,, we shall define their respective z-transforms by:
X(z) =
_c
x,zn,E(z)
n=--m
=
c
En9
n =-ca
With the definition in eq.9, the z-transform of each side of eq.8 gives: X(z) [l + a,2 + a*2 + . . . + a#]
= E(z)
or: 1 X(z) = E(z) 1+a,z+a,z* +...+a&v
(10)
Eq.10 represents the transfer function of an “all-pole” linear time-invariant system. Our physical claim that the lumped-parameter model represents a stable system is equivalent to the mathematical condition that the transfer function X(z)/E(z) contains no poles inside the unit circle. Let us now formulate this property in terms of the difference equation operator (1, a,, a,, . . , UN). The z-transform of this operator yields the denominator of eq.10. Thus, the property that the transfer function X(z)/E(z) has no poles inside the unit circle is equivalent to the property that: 1+cIz+c2z2
+...
+cNzN#Ofor(zI
(11)
We can therefore state that the operator (1, a,, a,, . . ., UN) yields a stable difference equation if and only if its z-transform has no zeroes inside the unit circle. Now if we replace the input to our lumped-parameter system by a unit impulse, i.e., let E, = 6, where: l,n=O 6, = i O,n#O then the impulse response, denoted by b,, is a stable time sequence. This condition follows from physical reasoning. But since:
B(z)=
1
1 +a12 +a$* +...+aNzN
(12)
it follows that B(z)has no zeroes or poles inside the unit circle. Now b, is the
63
sequence which is transmitted through the earth into the air when a unit impulse (deep source) excites the earth. From actual observation of field seismic data, b, has the bulk of its energy concentrated at its beginning and is rapidly attenuated due to successive reflections and refractions within the earth. Sequences that possess this “front-loaded” property are called minimum-delay sequences. Now if A (z) = 1 + a,z + a22 + . . . + a~$ is the z-transform of the operator (1, a,, a,, . . ., a~), then from eq.12 we observe that: B(z)A(z)
=1
or: b,*a,
= 6,
which shows that the operator a, and transmission impulse response b, are mutually inverse sequences. But since uR and b, are mutually inverse sequences, then we can conclude that a, is a stable difference equation operator if and only if b, is minimum-delay. That is, from purely physical grounds, we arrive at the important result: stable difference equation operator a,
*
minimum-delay response b,
Because b, represents all the multiple reflections and refractions resulting from an impulsive input 6,, b, also represents the system reverberation; b, = transmission (impulse) response = layered system reverberation wavelet. Thus, aR can also be viewed as the deconvolution operator which removes the layered system reverberation. Part 2. Reflections caused by surface source (internal primary reflections) Next, let us consider the N-layered (lumped-parameter) system of Fig.2 with the “deep” upgoing source replaced by a downgoing unit impulse (source) initiated at the surface. This situation, along with a pictorial definition of an internal primary reflection, is described in Fig.3. Now in Part 1, we claimed that an upgoing unit impulse (source) incident on the lowest interface gave rise to the transmitted (i.e., system. reverberation) waveform 13,. If we consider this situation in the context of linear timeinvariant system theory, then the unit impulse 6, is the cause and the system reverberation b, is the effect. Now if the cause is delayed by m units of time, i.e., 6,_, , then the effect is also delayed by m units of time, i.e., b,_,. Also, if the cause is scaled by a constant C, i.e., C6,, then the effect is also scaled by the constant C, i.e., Cb,. These properties, together with the property of superposition, characterize linear time-invariant systems. Thus, instead of considering the downgoing unit impulse (source), let us replace it by a set of hypothetical (i.e., mathematical) upgoing sources at the reflecting interfaces of our N-layered system. Since we know the “effect” b, due to the upgoing
64
r air
/
1st primary reflection
/
/
r2nd I
downgoing surface
primary reflection rNth /
primary reflection
N basement rock Fig.3. Internal primary reflections caused by a downgoing unit-impulse applied at the surface. For clarity, ray paths are drawn at oblique incidence but wave-motion analysis is for normal incidence.
“cause” 6,) we can then apply the linear time-invariant system properties to deduce the effect due to a downgoing unit impulse at the surface. Fig.4 shows the downgoing source impulse replaced by upgoing sources of strength E,, = r, where r, is the reflection coefficient of interface n. From physical experience gained in seismic prospecting, it is known that generally these reflection coefficients are small in magnitude, that is, I E, ( = I r, I< 0.1. Each of these hypothetical sources gives rise to a reverberation wavelet due to the layering. Further, to a first approximation, each reverberation wavelet is the same shape, namely the system reverberation wavelet b, . Silvia (1977) examines this assumption in some detail and shows the mathematical conditions under which it is a good approximation to a more exact stratified earth model. Taking into account the time delays associated with the source locations (see Fig.4), we reason as follows: transmission impulse-response b, M If an upgoing unit spike 6, (gives rise to)
then e161t--l %&-2
fN&-N
-+ --,
El&-1 EzL.--2
-+
%‘bn--N
Thus, using the superposition property of linear systems, we can add the responses due to all the hypothetical sources to obtain the total response h, Doing so, we get: h, = el b,_l
= +,*bn
+ ezbn+
.
+ . . . + Wbn-N
(13)
65
Eq.13 represents the reflection seismogram resulting from a downgoing unit impulse initiated at the surface. Further, h, is the resultant of all the system reverberations (reverberation wavelets b,) due to the upgoing hypothetical 0
\
\
‘\
El&r \
1
Single impulsive source of unit strength
\
Ez+-
\ =
2
\
N
Multiple impulsive sources of N-Layered strengths system ! El,
\ \
0, n#O
2
\ \
.
1, n=O 6 = n
1
\/
rEN=rN
E2,..EN
J
I
Fig.4.The downgoing unit impulse surface source 6, replaced by hypothetical mathematical) upgoing impulsive sources of strength ek i rk .
(i.e.,
sources E, with strength given by the reflection coefficient sequence r,. The z-transform of eq.13 gives: H(z) = E(z)B(z)
(14)
where: E(z) =
E1Z +
f*2 + . . . + eNzN
(15)
Substitution of eqs.12 and 15 into 14 yields: H(z)
=
El2
+
f22* + . . . + ENzN
1 +alz+a2z2
+.
.+aNzN
(16)
which gives the transfer function of a normal-incidence reflection seismogram and is the ARMA model used in reflection seismology. It is important to note that the same result can be obtained by deriving the exact expression for an N-layered earth model then making the approximation that all the reflection coefficients are small in magnitude, i.e., I r, I < 0.1. Under this approximation, fl G rl, e2 A r2,. . . , EN A rN; that is, the hypothetical sources are exactly the reflection coefficients (Silvia (1977)). Since the layered system is assumed to be both linear and time-invariant, then the reflection seismogram x, due to an arbitrary source pulse s, is: x, = h,*s, = En”b,“S, =
= En*( b,*+)
(17)
% *Wl
where w, = bn*sn is defined as a composite wavelet consisting of the reverberation wavelet b, and source pulse s, . Thus, eq.17 describes the normal-
66
incidence reflection seismogram, where b, represents the autoregressive component and s, *E, the moving-average component of the reflection seismogram. The basic deconvolution problem is to filter X, such that we can best recover the reflection coefficient sequence E, . Before we discuss the role of the kepstrum in deconvolution, let us investigate some of its basic properties. PROPERTIES
OF THE KEPSTRUM
For a minimum-delay system H,h(z), the hepstrum efficients of the Taylor-Maclaurin series: log H&(Z)
= F(z) = 5
fn is defined
f+z” , I2 I < 1
as the co-
(18)
n=O
However, the concept of the kepstrum is not limited to minimum-delay systems. Let us now consider a nonminimum-delay system H(z). If we assume that H(z) has no zeroes or poles in an annular region including the unit circle, then log H(z) is analytic in this annular region. As a result, log H(z) has a Laurent expansion: logH(2)
=
5 f,z” n=-ca
(19)
in this annular region. The coefficients for this Laurent series are the coefficients of a stable noncausal system. We call these coefficients the kepstrum of the nonminimum-delay system H(z). Let us now return to the case of a minimum-delay system:
5
Hmin(Z) =
hn.9
(20)
n=O
From eqs.18 and 20 we obtain: &i&r)
= exp
[JY~)l
or: hnzn --exP
2
[
jjofdn]
(21)
n=O
A formula
for obtaining
the minimum-delay
sequence
h, from the kepstrum
fn can be obtained as follows. Taking the derivative of eq.21 with respect to z gives:
& (jjo hnz") =exP[ j. fnzn]
2 (jjo fnzn)
(22)
6’7
Dividing eq.22 by eq.21 we obtain:
(23) n=O
Expanding eq.23 we get:
5
(n + l)hn+lZn
=
(C h,a”)( 2
(n +
n=O
n=O
n=O
l)f,+lZn)
(24)
The inverse z-transform of eq.24 gives the recursive relation (Oppenheim and Schafer, 1975):
(n + l)bn+l
= 5
hk( 11 + 1 -
k)fn+l-k for n = 0, 1, 2,
.. .
(25)
k=O
which is useful for extracting the minimum-delay sequence ha from the kepstrum fn. For IZ= 0, 1,2, 3, . . . eq.25 gives: n = 0:
h, = hof,
n=l:
h2 =
hof, +
2 fl =Mf,
n = 2: h3 = hOf3+ ihlfz
+$r:,
(
+ th2fl = ho.f3 + flf2
Let us now interpret the kepstral coefficient z = 0 that: ho = exp
1 +-f:6
fo. We have from eq.21 with
(fo) or f,, = log ho
Thus, it is often convenient to normalize the sequence h, so that its leading coefficient is unity, i.e., ho = 1; thus, the leading kepstral coefficient is zero, i.e., f. = 0. Let us now examine some of the properties of the kepstrum. At this point it is convenient to introduce the symbols a, and b, for the time sequences and
68
the corresponding Greek symbols (Yeand & for their respective kepstra. Thus, the kepstrum of a minimum-delay sequence a, with leading coefficient a0 = 1 = 0. If we let L(a,) is the causal sequence al, a2, as, . . . with leading term CY~ denote the z-transform of u, , then the kepstrum can be described by the composition of transformations L, log, and L-’ given by: % = L-l {log
(26)
[Jm,)lI
Hence, to find the kepstrum of a minimum-delay sequence h, whose leading coefficient is unity, we compute the z-transform of the sequence, take the natural logarithm of this transform, and find the inverse z-transform of the quantity log H(z). For convenience, we shall define the kepstral operator K given by the successive transformations L, log and L-’ such that: K(*) = L-l {log [L(m)]}
(27)
where K(h,) = fn and ho = 1 Let us now consider various properties of the kepstrum, i.e., convolution, reverse sequence, inverse sequence, and reverse-inverse sequence. Convolution. Given the minimum-delay sequences a, and b, whose leading coefficients are unity, then:
~@n*bl) = K(%) + K(b)
(28)
where: K(G)
=O,%,~2,%,..*
KC&l)
=
O,Pl,P?,PJ,
* * *
Hence, the kepstrum maps convolution into addition. To see this, we write: L(unfbn)
= A (z)B(z)
Then: log [L (a,*&
)I
= 1% [A WW)l = log A(z) + log B(z)
Taking the inverse z-transform of the above result we get: L-’ {log [L(un*bn)]
} = L-’ [log A(z)]
+ L-l [log B(z)]
or:
Reverse sequence.
Given the minimum-delay sequence b, with leading coeffi-
69
cient bO = 1, we form the reverse sequence b_, by folding b, about theorigin. The resulting sequence is called minimum-advance. Hence: K(b-n)
(29)
= P-n
We can show this result by first noting that: L(b,)
=
5
b,P =B(z)
n=O
and: L(b_,)
b,z_” =&z-l)
= jj n=O
Recall that: logB(2)
=&2+&z’
+P3z3 +. . .
Then: logB(C)
=&z-i
+pzO
+fi3z-3 + . . .
and we observe that given K(bn) = &, K(b_,) origin. In words:
is found by folding
Pn about&
minimum-delay sequence ++ causal kepstrum minimum-advance sequence tf purely noncausa2 kepstrum Inverse sequence. Given the minimum-delay sequence b, with leading coefficient b,, = 1, we define the inverse sequence b;’ by the relation: l,n=O b,*b,’
= 6, = O,n#O
Then : K(b,‘)
= -&,
(30)
Since L(bn*b;‘)
= 1, then B(z)L(b;‘)
or: L(b,‘)
= &-)
=
P(~)l-’
Now,
log [L(b,‘)]
= log [B(z)]-’
= -log B(z)
=1
70
Thus, it follows that: K(b,‘)
= -on
In words: minimum-delay sequence ++ causal kepstrum inverse of minimum-delay sequence ++ negative of causal kepstrum Reverse-inverse sequence. Given the minimum-delay sequence b, with leading coefficient b,, = 1, we define the sequence bZk as the inverse of the reverse sequence. Then: K(bI;)
= -/I+
(31)
To see this, we write: L(b$)
= [&z-l)]-’
Then: log [L(bz;)]
= log [B(P)]-’
= -log B(z_‘)
From the reverse-sequence property we can state: K(bZh) =-/3-n In words: minimum-delay sequence tf causal kepstrum reverse of inverse of minimum-delay sequence ++ negative of purely noncausal kepstrum The above properties bring out the essential features and symmetries of the kepstrum. Thus, given the kepstrum on of a minimum-delay sequence b, with leading coefficient b,, = 1, we can find the kepstrum of the reverse sequence b-,, inverse sequence b;’ and inverse-reverse sequence b1; by simple negations and reflections of the kepstrum on. For additional kepstral properties, see Oppenheim and Schafer (1975) and Tribolet (1977). In the above discussion, we developed the relation (eq.24): 2 (n + l)hn+lzn n=O
= ( jj
h,z”)
n=O
(
5
(n + l)f,+#)
n=O
which can be written as: fj n=O
(n + l)fn+lzn =
&) jj
(n + l)hn+lzn
n-0
where H(z) = 5 n=O
hnzn and ho = 1
(32)
71
The inverse z-transform K(h,)
= fn+l
=
of eq.32 gives:
5
hk+lh$-k
(33)
k=O
for n = 0, 1, 2, . . . and h,_,= 1 which explicitly gives the kepstrum f, in terms of the minimum-delay sequence h,. Let us now relate the mathematical properties of the kepstrum to the physics of the normal-incidence seismogram by next discussing the topic of deconvolution. THE ROLE
OF THE KEPSTRUM
IN DECONVOLUTION
Recall that the reflection seismogram x, can be approximately the output of an ARMA model with transfer function: H(z)
=-
E(z) A (2)
El2
=
1 +
+
U,Z
regarded
as
+Z2 + . . . + fNz N +
UzZZ + . . +aNZN
where E, = the reflection coefficient sequence and is the desired output of a deconvolution filter. Unfortunately, this system is not always minimum-delay, since the numerator polynomial coefficients fl, e2, . . . , EN are the reflection coefficients of an N-layered system and are dependent upon a particular geographic area. Thus, we are not always sure if all the zeroes of E(z) lie outside the unit circle, and in practice, the reflection coefficient sequence is generally mixed-delay, i.e., it contains zeroes both inside and outside the unit circle. However, the reverberation sequence b,, where B(z) = l/A (z) is always minimum-delay. If we excite the earth with an arbitrary mixed-delay source pulse sn, then the reflection seismogram x, is the convolution: x,
=sn
*en*bn
(34)
= En*Wn where w, = s, *b, is the composite mixed-delay wavelet. The method of predictive deconvolution (Robinson (1954)) assumes that the reflection coefficient sequence is a purely random sequence with a flat power spectrum. This assumption is usually met in practice (White and O’Brien, 1974). Such an assumption allows the method of predictive deconvolution to design least-squares deconvolution (noise-whitening) filters. From the convolution property of the last section, the kepstrum K (x, ) of the reflection seismogram is: K(x,)
=
K(E,*s,*b,)
= in +an
(35) +
pn
72
where K(E,) = &, K(s,) = a, and K(&) = &. The method of homomorphic deconvolution is concerned with processing the kepstrum K(x,) so that uII (source pulse kepstrum) and + (reverberation kepstrum) may be eliminated, leaving the desired kepstrum E, . Applying the inverse kepstral operator K-’ to the “filtered” kepstrum & would then yield the desired reflection coefficient sequence. It appears that deconvolution via the kepstrum is less restrictive than the method of predictive deconvolution, where one assumes that the reflection coefficients are an tumor-relatedsequence and that the source pulse be minimumdelay in order to achieve “optimum” results. Further, the kepstrum appears to be a more robust processor in the sense that the technique is not tied to any physical attributes of the model and that a minimum-delay source pulse assumption is not needed. However, both the method of predictive deconvolution and the method of kepstral deconvolution require an assumption of the phase spectrum of the source pulse (Treitel and Robinson, 1977). For example, let the canonical representation of the source pulse be: s, =pn*sp
where pn is a causal all-pass operator, ~$0)is the minimum-delay source wavelet, and s, is the nonminimumdelay source wavelet. The reflection seismogram for a nonminimum-delay source pulse can then be expressed as: xll
= f,
*b,*p,*s$f)
(37)
If we knew the all-pass system pn , i.e., if we had knowledge of the source phase characteristic, then we could operate on the seismogram with the inverse all-pass p;’ to yield: -l*
Pn
Xn = 62 *bn*s,(0)
which is more suitable to the method of predictive deconvolution. Often, the source pulse kepstrum u, is clustered around the kepstral origin (also called “low-time” portion of the kepstrum). See Stoffa et al. (1974) and Ulrych (1971). By a suitable window with the proper cutoff “quefrencies”, the source pulse can be theoretically recovered. In essence, the selection of the cutoff quefrencies determines how well we can estimate the source pulse phase and magnitude characteristics. Thus, both deconvolution techniques require some knowledge of the phase characteristic of the source pulse. CONCLUSIONS
The idea of the kepstrum appeared in the classical works of Poisson and Schwarz concerning the logarithmic potential problem. Szego and Kolmogorov used the kepstrum to factor the power spectrum of a non-deterministic process.
73
Because the normal-incidence reflection seismogram is approximated by a convolutional process, the kepstrum offers a solution to the deconvolution of reflection seismograms. We developed additional properties of the kepstrum which bring out its essential features and symmetries. As a result, the kepstrum evolves as an important mathematical construct in the theory of signal analysis and time series, and takes the place among such established functions as the spectrum and autocorrelation.
REFERENCES Bateman, H., 1944. Partial Differential Equations of Mathematical Physics. Dover, New York. B&h, M., 1968. Mathematical Aspects of Seismology. Elsevier, Amsterdam. BLth, M., 1974. Spectral Analysis in Geophysics. Elsevier, Amsterdam. Bogert, B.P., Healy, M.J.R. and Tukey, J.W., 1963. The quefrency analysis of time series for echoes. In: M. Rosenblatt (Editor), Proc. Symp. Time Series Analysis. Wiley, New York, pp.209-243. Kolmogorov, A.N., 1939. Sur l’interpolation et extrapolation des suites stationnaires. C.R. Acad. Sci. Paris, 208: 2043-2045. Kulhanek, O., 1975. Introduction to Digital Filtering in Geophysics. Elsevier, Amsterdam. Oppenheim, A.V. and Schafer, R.W., 1975. Digital Signal Processing. Prentice-Hall, Englewood Cliffs, N.J. Poisson, S.D., 1823. Sur la distribution de la chaleur dans les corps solides. J. EC. R. Polytech., Ser.1, 19: l-162. Robinson, E.A., 1954. Predictive Decomposition of Time Series with Applications to Seismic Exploration. Thesis, MIT, Cambridge, Mass. (Also in Geophysics, 32: 418-484.) Schafer, R.W., 1969. Echo removal by discrete generalized linear filtering. Res. Lab. Electron. MIT, Tech. Rep., 466. Schwarz, H.A., 1872. Zur Integration der partiellen Differentialgleichung. J. Reine Angewandte Math., pp.218-254. Silvia, M.T., 1977. Deconvolution of Geophysical Time Series. Thesis, Northeastern University, Boston, Mass. Stoffa, P.L., Buhl, P. and Bryan, G.M., 1974. The application of homomorphic deconvolution to shallow-water marine seismology - Part I. Models. Geophysics, 39: 401-416. Szego, G., 1915. Ein Grenzwertsatz iiber die Toeplitzschen Determinanten einer reellen positiven Funktion. Math. Ann., 76: 490-503. Treitel, S. and Robinson, E.A., 1977. Deconvolution -homomorphic or predictive? IEEE Trans. Geosc. Electron., GE-15(l): 11-13. Tribolet, J., 1977. Seismic Applications of Homomorphic Signal Processing. Thesis, Dept. of Electr. Eng. Comput. Sci., MIT, Cambridge, Mass. Ulrych, T.J., 1971. Application of homomorphic deconvolution to seismology. Geophysics, 36: 650-660. White, R.E. and O’Brien, P.N.S., 1974. Estimation of the primary seismic pulse. Geophys. Prosp., 22: 627-651.