Optimized suppression of coherent noise from seismic data using the Karhunen-Loeve transform

0.29 MB·English

by  Raul Montagne

Optimized suppression of coherent noise from seismic data using the Karhunen-Lo`eve
transform
Rau´l Montagne∗ and Giovani L. Vasconcelos†
Laborat´orio de F´ısica Te´orica e Computacional, Departamento de F´ısica,
Universidade Federal de Pernambuco, 50670-901 Recife, PE, Brazil
Signals obtained in land seismic surveys are usually contaminated with coherent noise, among
which the ground roll (Rayleigh surface waves) is of major concern for it can severely degrade the
quality of the information obtained from the seismic record. Properly suppressing the ground roll
from seismic data is not only of great practical importance but also remains a scientific challenge.
Hereweproposeanoptimized filterbased ontheKarhunen–Lo´evetransform forprocessing seismic
datacontaminatedwithgroundroll. Inourmethod,thecontaminated regionoftheseismic record,
6 to be processed by the filter, is selected in such way so as to correspond to the maximum of a
0 properly defined coherence index. The main advantages of the method are that the ground roll is
0 suppressed with negligible distortion of the remanent reflection signals and that the filtering can
2 beperformed on the computerin a largely unsupervisedmanner. The method has been devised to
filterseismic data, howeverit could also berelevantfor otherapplications where localized coherent
n
structures,embeddedinacomplexspatiotemporaldynamics,needtobeidentifiedinamorerefined
a
J way.

9
PACSnumbers: 93.85.+q,91.30.Dk,43.60.Wy,43.60.Cg

]
S
P I. INTRODUCTION much stronger in amplitude than the reflected signals.
. Groundrollaresurfacewaveswhoseverticalcomponents
n
are Rayleigh-type dispersive waves, with low frequency
i Locating oil reservoirs that are economically viable is
l and low phase and group velocities.
n oneofthemainproblemsinthepetroleumindustry. This
[ taskisprimarilyundertakenthroughseismicexploration, An example of seismic data contaminated by ground
whereexplosivesourcesgenerateseismicwaveswhosere- roll is shown in Fig. 1. This seismic section consists of
2
land–based data with 96 traces (one for each geophone)
v flections atthe different geologicallayersare recordedat
and 1001 samples per trace. A typical trace is shown in
2 the ground or sea level by acoustic sensors (geophones
2 or hydrophones). These seismic signals, which are later Fig. 2 corresponding to geophone 58. The image shown
0 processed to reveal information about possible oil oc- in Fig. 1 was created from the 96 traces using a stan-
7 dardimagingtechnique. Thehorizontalaxisinthisfigure
currences, are often contaminated by noise and properly
0 correspondstotheoffsetdistancebetweensourceandre-
cleaning the data is therefore of paramount importance
5 ceiverandthe verticalaxisrepresentstime, with the ori-
[1]. Inparticular,thedesignofefficientfilterstosuppress
0
ginlocatedattheupper–leftcorner. Themaximumoffset
/ noise that shows coherence in space and time (and often
n is475m(thedistancebetweengeophonesbeing5m)and
appears stronger in magnitude than the desired signal)
i the maximum time is 1000 ms. The gray levels in Fig. 1
l remains a scientific challenge for which novel concepts
n change linearly from black to white as the amplitude of
andmethodsarerequired. Inaddition,thefilteringtools
:
v developed to treat such kind of noise may also find rel- the seismic signal varies from minimum to maximum.
i Owing to its dispersive nature, the ground roll appears
X evant applications in other physical problems where co-
in a seismic image as a characteristic fan-like structure,
herent structures evolving in a complex spatiotemporal
r which is clearly visible in Fig. 1. The data shown in this
a dynamics need to identified properly.
figurewasprovidedbytheBrazilianPetroleumCompany
In land seismic surveys, the seismic sources generate
(PETROBRAS).
varioustypeofsurfacewaveswhichareregardedasnoise
Standard methods for suppressing ground roll include
since they do not contain information from the deeper
one-dimensionalhigh–pass filtering and two-dimensional
subsurface. Thisso-calledcoherentnoiserepresentsase-
f–k filtering [1]. Such “global” filters are based on the
rious hurdle in the processing of the seismic data since
elimination of specific frequencies and have the disad-
it may overwhelm the reflection signal, thus severely de-
vantage that they also affect the uncontaminated part
grading the quality of the information that can be ob-
of the signal. Recently, “local” filters for suppressing
tainedfromthedata. Asource-generatednoiseofpartic-
the groundrollhavebeen proposedusing the Karhunen-
ularconcernisthegroundroll,whichis themaintypeof
Lo`eve transform [2, 3] and the wavelet transform [4, 5].
coherent noise in land seismic records and is commonly
The Wiener-Levinsonalgorithmhasalsobeen appliedto
extract the ground roll [6].
Filters based on the Karhunen–Lo`eve (KL) transform
∗Electronicaddress: [email protected] are particularly interesting because of the adaptativity
†Electronicaddress: [email protected] of the KL expansion, meaning that the original signal
2
of the data autocorrelation matrix. In applying the KL
transformtosuppressthegroundroll,onemustfirstmap
the contaminatedregionofthe seismicrecordintoahor-
izontal rectangular region. This transformed region is
then decomposed with the KL transform and the first
few principalcomponents are removedto extractthe co-
herent noise, after which the filtered data is inversely
mapped back into the original seismic section. The ad-
vantage of this method is that the noise is suppressed
withnegligibledistortionofthereflectionsignals,foronly
the data within the selected region is actually processed
by the filter. Earlier versions of the KL filter [2, 3] have
however one serious drawback,namely, the fact that the
region to be filtered must be picked by hand—a proce-
dure that not only can be labor intensive but also relies
ongoodjudgmentofthe personperformingthe filtering.
Inthis articlewe proposea significantimprovementof
theKLfilteringmethod,inwhichtheregiontobefiltered
is selected automatically as an optimization procedure.
FIG.1: Aspace-timeplotofseismicdata. Thehorizontalaxis We introduceanovelquantity,namely,the coherencein-
represents the offset distance and the vertical axis indicates dex CI, which gives a measure of the amount of energy
time. Theoriginisattheupper-leftcorner,andthemaximum containedinthemostcoherentmodesforagivenselected
offsetandtimeare475mand1000ms,respectively. Thegray
region. The optimal region is then chosen as that that
scale is such that black (white) corresponds to the minimum
gives the maximum CI. We emphasize that introduc-
(maximum) amplitude of the seismic signal. The ground roll
ing a quantitative criterionfor selecting the ‘best’ region
noise appears as downward oblique lines.
to be filtered has the considerable advantage of yielding
a largely unsupervised scheme for demarcating and effi-
ciently suppressing the ground roll.
Although our main motivation here concerns the sup-
pressionof coherentnoise in seismic data, we should like
to remark that our method may be applicable to other
problems where coherentstructures embedded in a com-
plex spatiotemporal dynamics need to be identified or
characterizedinamorerefinedway. Forexample,theKL
transformhas been recently used to identify and extract
spatialfeaturesfromacomplexspatiotemporalevolution
in combustion experiment [7, 8, 9]. A related method—
the so-calledbiorthogonaldecomposition—has also been
applied to characterize spatiotemporal chaos and iden-
tified structures[10, 11] as well as identify changes in
the dynamical complexity, and the spatial coherence of
a multimode laser [12]. We thus envision that our opti-
mizedKLfiltermayfindapplicationsintheseandrelated
problemsofcoherentstructuresincomplexspatiotempo-
FIG. 2: Seismic signal recorded by a single geophone (trace
58). The amplitude is in arbitrary units and time in ms. ral dynamics.
Thearticleis organizedasfollows. InSec.II wedefine
theKarhunen–Lo`evetransform,describeitsmainproper-
is decomposed in a basis that is obtained directly from ties,anddiscussits relationto the singularvalue decom-
theempiricaldata,unlikeFourierandwavelettransforms position of matrices. In Sec. III we present the KL filter
which use prescribed basis functions. The KL transform and a novel optimization procedure to select the noise-
is a mathematical procedure (also known as proper or- contaminatedregiontobeparsedthroughthefilter. The
thogonal decomposition, empirical orthogonal function results of our optimized filter when applied to the data
decomposition,principalcomponentanalysis,andsingu- shown in Fig. 1 are presented in Sec. IV. Our main con-
lar value decomposition) whereby any complicated data clusions are summarized in Sec. V. In Appendixes A
set can be optimally decomposed into a finite, and often and B we briefly discuss, for completeness, the relation
small,numberofmodes(calledproperorthogonalmodes, between the KL transform and two other similar pro-
empirical orthogonal functions, principal components or cedures known as proper orthogonal decomposition (or
eigenimages) which are obtained from the eigenvectors empirical orthogonal function expansion) and principal
3
component analysis. us denote by the χ~ , i=1,...,m, the elements of the ith
i
row of the KL matrix Ψ, that is,
II. THE KARHUNEN–LOE`VE TRANSFORM χ~1
 χ~2 
Ψ= . (7)
...
 
A. Definition and main properties
 χ~ 
m
Consider a multichannel seismic data consisting of m Then (6) can be written as
traces with n samples per trace represented by a m n
× m
matrixA,sothattheelementA ofthedatamatrixcor-
ij
A= ~u χ~ , (8)
respondstotheamplituderegisteredattheithgeophone X i i
at time j. For definiteness, let us assume that m<n, as i=1
is usually the case. We also assume for simplicity that whereitisimpliedmatrixmultiplicationbetweenthecol-
the matrix A has full rank, i.e., r = m, where r denotes umn vector ~u and the row vector χ~ . The eigenvectors
i i
the rank of A. Letting the vectors ~xi and ~yj denote the ~ui are called empirical eigenvectors, proper orthogonal
elements ofthe ithrowandthejthcolumnofA, respec- modes, or KL modes.
tively, we can write As discussedin Appendix A, the totalenergyE ofthe
data can be defined as the sum of all eigenvalues,
~x1
 ~x2  m
A=(~y1 ~y2 ... ~yn)= ... . (1) E =Xλi, (9)
~xm  i=i
so that λ can be interpreted as the energy captured by
With the above notation we have i
the ith empirical eigenvector ~u . We thus define the rel-
i
ative energy E in the ith KL mode as
A =x =y , (2) i
ij ij ji
λ
i
where aij denotes the jth element of the vector~ai. (To Ei = m λ . (10)
avoid risk of confusion matrix elements will always be Pi=i i
denoted by capital letters, so that a small-cap symbol
WenotefurthermorethatsinceΓisacovariance-likema-
with two subscripts indicates vector elements.)
trixitseigenvaluesλ canalsobeinterpretedasthevari-
i
Next consider the following m m symmetric matrix ance of the respective principal component ~u ; see Ap-
× i
pendixBformoredetailsonthisinterpretation. Wethus
Γ AAt, (3) say that the higher λ the more coherent the KL mode
≡ i
~u is. Inthiscontext,themostenergeticmodesareiden-
i
where the superscript t denotes matrix transposition. It tified with the most coherent ones and vice-versa.
is a well known fact from linear algebra that matrices of AnimportantpropertyoftheKLexpansionisthatitis
theform(3),alsocalledcovariancematrices,arepositive ‘optimal’inthefollowingsense: ifweformthematrixΨ
k
definite [19]. Let us then arrange the eigenvalues λi of Γ bykeepingthefirstkrowsofΨandsettingtheremaining
lientn~uonb-aescthenedcionrgreosrpdoenr,dii.neg.,(λn1or≥mλa2liz≥ed..).e≥igλenmve>ct0o,rsa.nd m−k rows to zero, then the matrix Ak given by
i
The Karhunen-Lo`eve (KL) transform of the data ma- A =UΨ (11)
k k
trix A is defined as the m n matrix Ψ given by
×
isthebestapproximationtoAbyamatrixofrankk <m
Ψ=UtA, (4) in the Frobenius norm (the square root of the sum of
the squares of all matrix elements) [13]. This optimal-
where the columns of the matrix U are the eigenvectors ity property of the KL expansion lies at the heart of its
of Γ: applicationsindatacompression[14]anddimensionality
reduction [13], for it allows to approximate the original
U =(~u1 ~u2 ... ~um). (5) data A by a smaller matrix Ak with minimum loss of in-
formation (in the above sense). Another interpretation
The original data can be recovered from the KL trans- ofrelation(11)isthatit givesa low-lassfilter [15], forin
form Ψ by the inverse relation this case only the first k KL modes are retained in the
filtered data A .
k
A=UΨ. (6) On the other hand, if the relevant signal in the appli-
cation at hand is contaminated with coherent noise, as
WerefertothisequationastheKLexpansionofthedata is the case of the ground roll in seismic data, one can
matrix A. To render such an expansionmore explicit let use the KL transform to removeefficiently such noise by
4
constructing a high-pass filter. Indeed, if we form the Offset
matrix Ψ′k by setting to zero the first k rows of Ψ and A
keeping the remaining ones, then the matrix A′ given
k
by First
me original
A′ =UΨ′ (12) Ti sector
k k B A B
...
is a filtered version of A where the first k ‘most coher- D First transformed
ent’ modes have been removed. However, if the noise is sector
localized in space and time it is best to apply the filter Last ...
original
only to the contaminated part of the signal. In previous
sector
versions of the KL filter the choice of the region to be Last transformed
parsedthroughthefilterwasmadea priori, accordingto sector
the best judgement of the person carrying out the filter-
C D C
ing, thus lending a considerable degree of subjectivity to
the process. In the next section, we will show how one
can use the KL expansion to implement an automated
FIG. 3: Schematic diagram for demarcating the ground roll
filter where the undesirable coherent structure can be
onaseismicsectionandthecorrespondingrectangularsectors
‘optimally’ identified and removed.
obtained byapplying a linear map.
Beforegoingintothat,however,weshallbrieflydiscuss
below an important connection between the KL trans-
form and an analogous mathematical procedure known It thus follows that the decomposition in eigenimages
as the singular value decomposition of matrices. Read- seen in (14) is precisely the KL expansion given in (8).
ersalreadyknowledgeableabouttheequivalencebetween Furthermore the approximation A given in (11) can be
k
these two formalisms (or more interested in the specific written in terms of eigenimages as
application of the KL transform to filter coherent noise)
k
may skip the remainder of this section without loss of
A = σ Q . (17)
continuity. k X i i
i=1
Similarly, the filtered data A′ shown in (12) reads in
k
B. Relation to Singular Value Decomposition terms of eigenimages:
m
Werecallthatthesingularvaluedecomposition(SVD) A′ = σ Q . (18)
of any m n matrix A, with m < n, is given by the k X i i
× i=k+1
following expression:
The SVD provides an efficient way to compute the KL
A=UΣVt, (13) transform,andweshallusethismethodinthenumerical
procedures described in the paper.
whereU isasdefinedin(5),Σisam mdiagonalmatrix
×
with elements σ = √λ , the so-called singular values of
i i
A,andV isam nmatrixwhosecolumnscorrespondto III. THE OPTIMIZED KL FILTER
the m eigenvecto×rs ~v of the matrix AtA with nonzero
i
{ }
eigenvalues. The SVD allows us to rewrite the matrix A
As already mentioned, owing to its dispersive nature
as a sum of matrices of unitary rank:
the ground-roll noise appears in a seismic image as a
m m typical fan-like coherent structure. This space-time lo-
A= σ Q = σ ~u~v t. (14) calization of the ground roll allows us to apply a sort of
X i i X i i i
i=1 i=1 ‘surgical procedure’ to suppress the noise, leaving intact
theuncontaminatedregion. Todothat,wefirstpicklines
In the context of image processing the matrices Q are
i to demarcate the startand endof the groundrolland, if
called eigenimages [16].
necessary,intermediatelinestodemarcatedifferentwave-
Now, comparing (6) with (13) we see that the KL
trains,asindicatedschematicallyinFig.3. Inthis figure
transform Ψ is related to the SVD matrices Σ and V
we have for simplicity used straight lines to demarcate
by the following relation
the sectors but more general alignment functions, such
as segmentedstraightlines, can alsobe chosen[2, 3]. To
Ψ=ΣVt, (15)
makeourdiscussionasgeneralaspossible,let us assume
so that the row vectors χ~i of Ψ are given in terms of the thatwehaveasetofN parameters{θi},i=1,...,N,de-
singular values σ and the vectors~v by scribing our alignment functions. For instance, in Fig. 3
i i
the parameters θ wouldcorrespondto the coefficients
i
χ~ =σ~v t. (16) of the straight li{nes}defining each sector.
i i i
5
Once the region contaminated by the ground roll has as the averagecoherence index of all sectors:
been demarcated, we map each sector onto a horizon-
tal rectangular region by shifting and stretching along 1 l
the time axis; see Fig. 3. The data points between the CI(θ1,...,θN)= lXCIk. (20)
top and bottom lines in each sector is mapped into the k=1
corresponding new rectangular domain, with the map-
As the name suggests, the coherence index CI is a mea-
ping being carriedout via a cubic convolutioninterpola-
sure of the amount of ‘coherent energy’ contained in
tion technique [17]. After this alignment procedure the
the chosen demarcated region given by the parameters
groundrolleventswillbecomeapproximatelyhorizontal, θ N . Thus, the higher CI the larger the energy con-
favoring its decomposition in a smaller space. Since any { i}i=1
tainedinthemostcoherentmodesinthatregion. Forthe
given transformed sector has a rectangular shape it can
purpose of filtering coherent noise it is therefore mostly
be representedbyamatrix,whichinturncanbe decom-
favorableto pickthe regionwiththe largestpossibleCI.
posedinempiricalorthogonalmodes(eigenimages)using
We thus propose the following criterion to select the op-
the KL transform. The first few modes, which contain
timalregiontobefiltered: varytheparameters θ over
i
most of the ground roll, are then subtracted to extract some appropriate range and then choose the v{alu}es θ∗
i
the coherent noise. The resulting data for each trans-
that maximize the coherence index CI, that is,
formed sector is finally subjected to the corresponding
inverse mapping to compensate for the original forward CI(θ1∗,...,θN∗ )=max[CI(θ1,...,θN)]. (21)
mapping. This leaves the uncontaminated data (lying {θi}
outside the demarcated sectors) unaffected by the whole
Once we have selected the optimal region, given by the
filtering procedure.
parameters θ∗ N , we then simply apply the KL fil-
{ i}i=1
The KL filter described above has indeed shown good ter to this region as already discussed: we remove the
performance in suppressing source-generated noise from first few eigenimages from each transformed sector and
seismic data [2, 3]. The method has however the draw- inversely map the data back into the original sectors, so
back that the region to be filtered must be picked by as to obtain the final filtered image. In the next section
hand, which renders the analysis somewhat subjective. wewillapplyouroptimizedKLfilteringproceduretothe
Inordertoovercomethis difficulty, it wouldbe desirable seismic data shown in Fig. 1.
tohaveaquantitativecriterionbasedonwhichonecould
decide what is the ‘best choice’ for the parameters θ
i
{ }
describing the alignment functions. In what follows, we IV. RESULTS
proposeanoptimizationprocedurewherebytheregionto
befilteredcanbeselectedautomatically,oncethegeneric Here we illustrate how our optimized KL filter works
form of the alignment functions is prescribed. by applying it to the seismic data shown in Fig. 1. In
Supposewehavechosenlsectorstodemarcatethedif- this case, it suffices to choose only one sector to demar-
ferentwavetrainsinthe contaminatedregionofthe orig- cate the entire region contaminated by the ground roll.
inal data, and let θ1,...,θN be the set of parameters Thismeansthatwehavetoprescribeonlytwoalignment
characterizingther{espectivea}lignmentfunctionsthatde- functions, corresponding to the uppermost and lower-
fine these sectors. Let us denote by A˜ , k = 1,...,l, the most straight lines (lines AB and CD, respectively) in
k
matrix representing the kth transformedsector obtained Fig. 3. To reduce further the number of free parame-
from the linear mapping of the respective original sec- ters in the problem, let us keep the leftmost point of
tor, as discussed above. For each transformed sector A˜ the upper line (point A in Fig. 3) fixed to the origin, so
k
we then compute its KL transform and calculate the co- that the coordinates (iA,jA) of point A are set to (0,0),
herence index CI for this sector, defined as the relative while allowing the point B to move freely up or down
k
energy contained in its first KL mode: within certain range; see below. Similarly, we shall keep
the rightmost point of the lower line (point C in Fig. 3)
pinned at a point (i ,j ), where i =95 and j is cho-
λk C C C C
CI = 1 , (19) senso that the entire groundrollwavetrainis abovethis
k rk λk point. The other endpoint of the lower demarcation line
Pi=1 i
(point D in Fig. 3) is allowed to vary freely. With such
restrictions, we are left with only two free parameters,
where λk are the eigenvalues of the correlation matrix
Γ˜ = A˜iA˜t and r is the rank of A˜ . Such as defined namely, the angles θ1 and θ2 that the upper and lower
k k k k k demarcation lines make with the horizontal axis. So re-
above, CI represents the relative weight of the most
k ducing the dimensionality of our parameter space allows
coherent mode in the KL expansion of the transformed
sector A˜ . (A quantity analogous to our CI is known in us to visualize the coherence index CI(θ1,θ2) as a 2D
k surface. Forthe casein hand,it is moreconvenienthow-
the oceanographyliterature as the similarity index [18].)
everto expressCI notasa function ofthe anglesθ1 and
Next we introduce an overall coherence index θ2 but in terms of two other new parameters introduced
CI(θ1,...,θN) for the entire demarcated region, defined below.
6
Letthe coordinatesofpointB,whichdefinesthe right
endpointoftheupperdemarcationlineinFig.3,begiven
by (i ,j ), where i = 95. In our optimization proce-
B B B
dureweletpointB movealongtherightedgeoftheseis-
mic section by allowing the coordinate j to vary from
B
a minimum value j to a maximum value j , so
Bmin Bmax
that we can write
j =j +k∆ , k =0,1,...,N (22)
B Bmin B B
where N is the number of intermediate sampling
B
points between j and j , and ∆ = (j
Bmin Bmax B Bmax −
j )/N . Similarly, for the coordinates (i ,j ) of
Bmin B D D
point D in Fig. 3, which is the moving endpoint of the
lower straight line, we have i =0 and
D
j =j +l∆ , l=0,1,...,N , (23)
D Dmin D D
where N is the number of sampling points between
D
j and j , and ∆ =(j j )/N .
Dmin Dmax D Dmax − Dmin D
FIG.5: a)Theselectedregioninthenewdomain;b)itsfirst
eigenimage; c) thesecond eigenimage; and d)theresult after
subtracting thefirst eigenimage.
correspondstoabout33%ofthetotalenergyoftheimage
in Fig. 5a, as can be seen in Fig. 6 where we plot the
relative energy E captured by the first 10 eigenimages.
i
Thesecondeigenimage,showninFig.5c,capturesabout
10% of the total energy, with each successively higher
mode contributing successively less to the total energy;
see Fig. 6. In Fig. 5d we give the result of removing
FIG. 4: The coherence index CI as a function of the indices the first KL mode (Fig. 5b) from Fig. 5a. It is clear
k and l that definethedemarcation lines; see text. inFig.5dthatbyremovingonlythe firsteigenimagethe
mainhorizontalevents(correspondingtothegroundroll)
have already been greatly suppressed.
Foreachchoiceofkandlin(22)and(23),weapplythe
procedure described in the previous section and obtain Performingtheinversemappingofthe imageshownin
thecoherenceindexCI(k,l)ofthecorrespondingregion. Fig.5cyieldsthedataseenintheregionbetweenthetwo
In Fig. 4 we show the energy surface CI(k,l), for the whitelinesinFig.7a,whichshowsthefinalfilteredimage
case in which j = 280, j = 600, j = 864, for this case (i.e., after removingthe firstKL mode from
Bmin Bmax C
j = 0, j = 576, and N = N = 64. We the transformed region). We see that the ground roll in-
Dmin Dmax B D
see in this figure that CI possesses a sharp peak, thus side the demarcated regionin Fig. 7a has been consider-
showingthatthis criterionis indeedquite discriminating ably suppressed, while the uncontaminated signal (lying
with respect to the positioning of the lines demarcating outsidethemarkedregion)hasnotbeenaffectedatallby
the region contaminated by the ground roll. The global the filtering procedure. If one wishes to filter further the
maximumofCI inFig.4islocatedatk =42andl =24, ground roll noise one may subtract successively higher
and in Fig. 5a we show the transformed sector obtained modes. For example, in Fig. 7b we show the filtered im-
from the linear mapping of this optimal region. In this age after we also subtract the second eigenimage. One
figure one clearly sees that the ground roll wavetrains sees that there is some minor improvement, but remov-
appearmostlyashorizontalevents. InFig.5bwepresent ing additionalmodes is not recommendedfor it starts to
the first eigenimage of the data shown in Fig. 5a, which degrade relevant signal as well.
7
0.4 optimization procedure whereby the ground roll region
is selected so as to maximize an appropriately defined
0.3
relative energy0.2 cquousheirdeerieannstcihenepinoudpte,txiomnCilzyIa.tthiWoenegpeenrmoepcriehcdaauslirizegenatmshweantetlolfuaursnmcthteieothnnosudmto,brbeeer-
0.1 of eigenimages to be removed from the selected region.
These may vary depending on the specific applicationat
00 1 2 3 4mod5e num6ber7 8 9 10 hand. However,oncethesechoicesaremade,thefiltering
task can proceed in the computer in an automated way.
FIG. 6: The relative energy of the first 10 KL modes of the
Althoughourmainmotivationherehasbeensuppress-
region shown in Fig. 5a.
ingcoherentnoisefromseismicdata,ourmethodisbyno
means restricted to geophysicalapplications. In fact, we
believe that the method may prove useful in other prob-
lemsinphysicsthatrequirelocalizingcoherentstructures
inanautomatedandmorerefinedway. Wearecurrently
exploring further such possibilities.
Acknowledgments
Financial support from the Brazilian agencies CNPq
and FINEP and from the special research program
CTPETRO is acknowledged. We thank L. Lucena for
many useful conversation and for providing us with the
data.
APPENDIX A: RELATION BETWEEN THE KL
TRANSFORM AND PROPER ORTHOGONAL
DECOMPOSITION
In dynamical systems the mathematical procedure
akinto the KLtransformis calledthe properorthogonal
decomposition(POD).Inthiscontext,onemayvieweach
columnvector~y ofthedatamatrixAasasetofmmea-
j
surements(realornumerical)ofagivenphysicalvariable
f(x,t) performed simultaneously at m space locations
and at a certain time t , that is, y =f(x=x ,t=t ),
j jk k j
k =1,...,m. For example, in turbulent flows the vectors
FIG. 7: a) The filtered seismic section after removing the ~y oftenrepresentmeasurementsofthefluidvelocityatm
i
first eigenimage of the select region shown in Fig. 5a. In b) pointsinspaceatagiventimei. ThedatamatrixAthus
we show the result after removing thefirst two eigenimages.
corresponds to an ensemble ~y of n such vectors, rep-
j
{ }
resenting a sequence of m measurements over n instants
of time. In POD one is usually concerned with finding
a low-dimensional approximate description of the high-
V. CONCLUSIONS dimensional dynamical process at hand. This is done by
finding an ‘optimal’ basis in which to expand (and then
An optimized filter based on the Karhunen–Lo´eve truncate) a typical vector ~y of the data ensemble. Such
transform has been constructed for processing seismic a basis is given by the eigenvectors of the time-averaged
data contaminated with coherent noise (ground roll). A autocorrelation matrix R, which is proportional to the
great advantage of the KL filter lies in its local nature, matrix Γ define above:
meaningthatonlythecontaminatedregionoftheseismic
n
1 1
recordisprocessedbythe filter,whichallowstheground R ~y~yt = ~y ~yt = Γ. (A1)
roll to be removed without distorting most of the reflec- ≡(cid:10) (cid:11) nX i i n
i=1
tionsignal. Anotheradvantageisthatitisanadaptative
method in the sense the the signal is decomposed in an Hence the eigenvectors ~u of Γ are also eigenvectors
i
{ }
empirical basis obtained from the data itself. We have of R. In POD parlance the eigenvectors ~u are called
i
{ }
improved considerably the KL filter by introducing an empiricaleigenvectorsorproperorthogonalmodes. Inthe
8
continuous case, the corresponding eigenfunctions u (x) APPENDIX B: RELATION BETWEEN THE KL
i
of the autocorrelation operator are known as empirical TRANSFORM AND PRINCIPAL COMPONENT
orthogonal functions (EOF). ANALYSIS
From (1), (4) and (5), one can easily verify that
In statistical analysis of multivariate data, the KL
Ψij =~ui ~yj. (A2) transform is known as principal component analysis
·
(PCA). In this case, one views the elements of a row
We thus see that the columns of the KL transform Ψ
vector~xi =(xi1,...,xin) of the data matrix A as being n
correspondtothe coordinatesofthe vectors~y inthe em-
realizations of a random variable X , so that the matrix
i
pirical basis: A itself corresponds to n samples of a random vector X~
m withm components: X~ =(X1,...,Xm)t. Inother words,
~yi =XΨki~uk. (A3) thecolumnvectors~yj correspondtothesamplesofX~. If
k=1 therowsofAarecentered,i.e.,thevariablesXihavezero
mean,thenthematrixΓisproportionaltothecovariance
Itisthisexpansionofanymemberoftheensembleinthe
matrix S of X~ [20]:
empirical basis that is called the proper orthogonal de- X
composition or empirical orthogonalfunction expansion.
1 1
It now follows from (A3) that (SX)ij XiXj = x~i x~j = Γij, (B1)
≡h i n · n
n
1 1 1 or alternatively in matrix notation
~y2 = ~y2 = (ΨΨt) = λ , (A4)
(cid:10) (cid:11) nX i nX ii nX i
i=1 i i 1
S X~ X~ t = Γ. (B2)
where in the last equality we used the fact that X ≡D E n
[NotethatthematricesRandS definedrespectivelyin
ΨΨt =UtΓU =Λ, (A5) X
(A1) and (B2) are essentially the same but have differ-
ent interpretations.] In the PCA context, the diagonal
where Λ is the diagonal matrix Λ = diag(λ1,...,λm).
elementsΓ ofthe matrixΓ arethus proportionalto the
Equation (A4) thus suggests that we can interpret the ii
variance of the variables X , whereas the off-diagonalel-
eigenvalue λ as a measure of the energy in the ith em- i
i
ements Γ , i=j, are proportionalto the covariance be-
piricalorthogonalmode. Forexample,inthe caseoftur- ij 6
tween the variables X and X . Furthermore, the eigen-
bulent flows where the vector ~y contains velocity mea- i j
i
vectors ~u of Γ correspond to the principal axis of the
surements at time i, the left hand of (A4) yields twice i
the average kinetic energy per unit mass, so that 1λ covariancematrix SX. The idea behind PCAis to intro-
2 i duce a new set of m variables P , each of which being a
gives the kinetic energy in the ith empirical orthogonal i
linearcombinationofthe originalvariablesX ,suchthat
mode [13]. Similarly, in the case of seismic data the vec- i
these new variables are mutually uncorrelated. This is
tors ~y represent amplitudes of the reflected waves, and
henceithequantity n ~y2 = n λ maybeviewedas accomplished by projecting the vector X~ onto the prin-
Pi=1 i Pi=i i cipaldirectionsofthe covariancematrix. Moreprecisely,
a measureofthe totalenergyofthe data, thus justifying
we define the principal components P , i = 1,...,m, by
the definition given in (9). i
the following relation
The optimality of the KL expansion also has a nice
physical and geometricalinterpretation, as follows. Sup-
n
posewewriteavector~yinanarbitraryorthonormalbasis P =X~ ~u = u X . (B3)
~e m : i · i X ij j
{ i}i=1 j=1
m
~y = a~e , (A6) In other words, the vector of principal components P~ =
Xi=1 i i (P1,...,Pm)t is obtained from a rotation of the original
vector X~:
where a = ~e ~y. If we now wish to approximate ~y by
i i
·
only its first k <m components, P~ =UtX~. (B4)
k ThecovariancematrixS ofthe principalcomponentsis
~yk = a~e , (A7) P
X i i then given by
i=1
1 1
then the optimality of the KL expansion implies that SP =DP~ P~tE=DUtX~ X~ tUE= nUtΓU = nΛ, (B5)
thefirstk properorthogonalmodescapturemoreenergy
(on average) that the first k modes of any other basis. thus showing that
More precisely, the mean square distance ~y ~yk 2 is
(cid:10)| − | (cid:11)
minimum if we use the empirical basis. P P =0, for i=j, (B6)
i j
h i 6
9
as desired. The first principal component P1 then rep- row of the KL transform Ψ correspond to the n samples
resents the particular linear combination of the origi- or scores of the ith principal component. That is, if we
nal variables X (among all possible such combinations denote the sample vector of the ith principalcomponent
i
that yield mutually uncorrelated variables) that has the by ~pi = (pi1,...,pin), then pij = ~ui ~yj = Ψij. For this
·
largest variance, with the second principal component reasonin the PCA context the KL transform Ψ is called
possessing the second largest variance, and so on. the matrix of scores.
From(4)and(B4)oneseesthattheelementsoftheith
[1] O. Yilmaz, Seismic Data Processing (Society of Explo- [12] F.PapoffandG.D’Alessandro,Phys.Rev.A70,063805
ration Geophysicist, Tulsa, 1987). (2004).
[2] X.Liu, Geophysics 64, 564 (1999). [13] P. Holmes, J. L. Lumley, and G. Berkooz, Turbu-
[3] Y.K.Tyapkin,N.Marmalevskyy,andZ.V.Gornyak,in lence,Coherent Structures, Dynamical Systems and Sym-
EAGE66thConference(2004),expandedAbstractD028. metry (Cambridge University Press, Cambridge, 1996).
[4] A. J. Deighan and D. R. Watts, Geophysics 62, 1896 [14] C.Andrews,J.Davies,andG.Schwartz,Proc.IEEE55,
(1997). 267 (1967).
[5] G. Corso, P. Kuhn, L. Lucena, and Z. Thom´e, Physica [15] S. Freire and T. Ulrych,Geophysics 53, 778 (1988).
A 318, 551 (2003). [16] H. Andrews and B. Hunt, Digital Image Restoration
[6] H.Karsli andY.Bayrak,Journal ofAppliedGeophysics (Prentice Hall, 1977).
55, 187 (2004). [17] S.ParkandR.Schowengerdt,Computer,Vision,Graph-
[7] A. Palacios, G. H. Gunaratne, M. Gorman, and K. A. ics & Image Processing 23, 256 (1983).
Robbins,Phys. Rev.E 57, 5958 (1998). [18] H.-J. Kim, J.-K. Chang, H.-T. Jou, G.-T. Park, B.-C.
[8] K. R. M. Gorman, J. Bowers, and R. Brockman, Chaos Suk, and K. Y. Kim, J. Acoust. Soc. Am. 111, 794
14, 467 (2004). (2002).
[9] P. Blomgren, S. Gasner, and A. Palacios, Chaos 15, [19] This is true if r = min{m,n}, as assumed; if r <
013706 (2005). min{m,n}, thematrix Γ has r nonzero eigenvalues with
[10] S.Bouzat,H.S.Wio,andG.B.Mindlin,PhysicaD199, all m−r remaining eigenvalues equal tozero.
185 (2004). [20] Thepropernormalizationforthecovariancematrixisof-
[11] P. D. Mininni, D. O. G´omez, and G. B. Mindlin, Phys. tenchosentobe 1 insteadof 1,butthisisnotrelevant
n−1 n
Rev.Lett. 89, 061101 (2002). for our discussion here.