`United States Patent
`5,986,973
`[11] Patent Number:
`
`[45] Date of Patent: Nov. 16, 1999
`Jericéevié et al.
`
`US005986973A
`
`[54] ENERGY MINIMIZATION IN SURFACE
`MULTIPLE ATTENUATION
`
`(57]
`
`ABSTRACT
`
`[75]
`
`Inventors: Zeljko Jerivevié, William H.
`Dragoset, Jr., both of Houston, ‘Tex.
`
`[73] Assigncc: Western Atlas International, Inc.,
`Houston, Tex.
`
`[21] Appl. No.: 08/923,150
`
`[22]
`
`Filed:
`
`Sep. 4, 1997
`
`Tmt, C00 ieee cccsseecsssssssssessessssnseeeesssssneees GO1V 1/38
`[SD]
`[52] U.S. CI. ees csseeseeeneeeeee 367/24; 367/45; 367/24;
`367/53
`[58] Field of Search 0... 367/21, 24, 38,
`367/45, 53
`
`[56]
`
`References Cited
`U.S. PATENT DOCUMENTS
`
`5,138,583
`5,448,531
`5,587,965
`5,719,821
`5,729,506
`
`8/1992 Wasonet al. on eee eceeecne rere 367/21
`9/1995 Dragosetet al. cece 367/21
`12/1996 Dragoset, Jr. et al. eee. 367/24
`2/1998 Sallas et al. cee eeeece reese 367/21
`3/1998 Dragosetet al...eee 367/21
`
`Primary Examiner—Christine K. Oda
`Assistant Examiner—Anthony Jolly
`
`A marine seismic signal is transformed from time domain
`into frequency domain and represented by matrix D. The
`marine data signal is truncated in time, transformed into the
`frequency domain and represented by matrix D; Eigenvalue
`decomposition Dj_5,,"° of matrix D; is computed. Matrix
`product D S is computed and saved in memory. Conjugate
`transpose [D S]* is computed and saved in memory. Matrix
`product [D S]*[D S] is computed and saved in memory.
`Matrix inverse S~* is computed and saved in memory.
`Conjugate transpose (S~*)* is computed and saved in
`memory. Matrix product S"1(S~*)* is computed and saved in
`memory. An initial estimate of the source wavelet wis made.
`Source wavelet w is optimized by iterating the steps of
`computing the diagonal matrix [I-w7*A], computing matrix
`inverse [I-w-‘A]*, computing conjugate transpose [(I-w~
`1A)"'}*, retrieving matrix products [D S]*[D S] and S-'(S~
`1)* from memory, and minimizing the total energy in trace
`of matrix product
`
`S“(st)>[(-wtA)4[D SID S[I-wAP
`
`Primary matrix P representing the wavefield free of surface
`multiples is computed by inserting computed value for w
`into the expression [D S][I-wtA] ‘S~". Primary matrix P is
`inverse transformed from frequency domain into time
`domain.
`
`17 Claims, 4 Drawing Sheets
`
`100.
`
`108
`
`
`
`0X
`110 -
`LC]
`
`168
`
`166
`
`
`
`120
`=
`KX
`XL,
`SS
`SS
`SS
`SS
`KR
`SX
`SS
`DY
`RE 102
`SEK
`CER
`RR
`x
`S
`RR
`yA
`KR 114
`SARL
`<
`SK
`
`
`BI
`DIOS
`
`CARROLL
`WOONONA
`
`170
`
`112
`
`Page 1 of 13
`
`SAMSUNG EXHIBIT 1029
`
`Page 1 of 13
`
`SAMSUNG EXHIBIT 1029
`
`
`
`U.S. Patent
`
`Nov.16, 1999
`
`Sheet 1 of 4
`
`5,986,973
`
`N9
`
`ms
`
`“
`
`i
`SINSIESISISISISISSYSS
`CMCCCCAR
`IOSOIIILIIIINININ
`
`
`»©»oh6wAMPHPwrKCLLEGEETEOKRCROOKEX
`
`VWIYIN
`SSPDAHHHAOOPOPHS>}.HYmWM9»
`
`SPOIL}SWAYSDYDSYYAYASS
`
`VONRARRRRRRNA
`
`
`SSISNRELRANANENANANASSOONSOsRRORRRPROOIRTRRRORS
`~IN—INXSNKXKAKARAR
`SPAHLYYOOA
`
`
`
`
`
`“
`x
`CRY
`Te
`Xe
`4KK
`AGZ,
`ZRNS
`x
`XX
`
`“
`
`\
`
`291
`
`991
`
`~~
`
`PEOld
`
`Page 2 of 13
`
`Page 2 of 13
`
`
`
`
`
`U.S. Patent
`
`Nov. 16, 1999
`
`Sheet 2 of 4
`
`5,986,973
`
`200.
`
`INPUT DATA
`
`202
`
`206
`
`FOURIER TRANSFORM DATA
`
`212
`
`204
`
`TRUNCATE TRACES IN
`TIME DOMAIN
`
`208
`FOURIER TRANSFORM
`TRUNCATED TRACES
`
`210
`
`214
`
`
`
`
`
`
`
`
`MODIFY TRACES IN
`FREQUENCY DOMAIN
`
`
`
`CONSTRUCT MODIFIED
`CONSTRUCT DATA
`
`DATA CUBE Dm(s,r,f)
`CUBE D(s,r,f)
`
`
`218
`216
`
`EXTRACT MATRIX Dw(s.r)
`EXTRACT MATRIX D(s,r)
`
`FROM DATA CUBE Dx(s,r,f)
`FROM DATA CUBE U(s,r,f)
`
`
`
`COMPUTE PRIMARY MATRIX P(s,r)
`
`222
`INSERT MATRIX P(s,r)
`
`INTO DATA CUBE P(s,r,f)
`
`224
`
`INVERSE FOURIER TRANSFORM P
`
`FIG, 2
`
`Page 3 of 13
`
`Page 3 of 13
`
`
`
`U.S. Patent
`
`Nov. 16, 1999
`
`Sheet 3 of 4
`
`5,986,973
`
`300 /
`
`302
`COMPUTE EIGENVALUE DECOMPOSITION Du= S A S7!
`
`304
`COMPUTE PRODUCT [0-S], ADJOINT [D-S]*
`AND PRODUCT [D:S]*-[D-S]
`
`306
`STORE PRODUCT [D S]*-[D-S] IN MEMORY
`
`308
`COMPUTE INVERSE S7~', ADJOINT (S7')*
`AND PRODUCT (S~') - (S7')*
`
`310
`STORE PRODUCT (S7')- (S7')* IN MEMORY
`
`CONTINUE TO STEP 402 (FIG. 4)
`
`FIG. 3
`
`Page 4 of 13
`
`Page 4 of 13
`
`
`
`U.S. Patent
`
`Nov. 16, 1999
`
`Sheet 4 of 4
`
`5,986,973
`
`#00.
`
`CONTINUE FROM STEP 312 (FIG.3)
`
`402
`
`ESTIMATE WAVELET w
`
`404
`
`COMPUTE DIAGONAL MATRIX [
`
`406
`I-w7!A]
`
`408
`COMPUTE DIAGONAL MATRIX INVERSE [ I-w7!A]7!
`
`410
`COMPUTE ADJOINT [(l-w ~'A)7']*
`
`412
`RETRIEVE [D-S]*-[D-S] AND
`(S~') - (S7'1)* FROM MEMORY
`
`414
`
`COMPUTE TRACE OF MATRIX PRODUCT
`;
`
`TR $S7'-(S71)* [Gews7ta)71]* [D-s]*- (D-S]
`l-w7ta}7'}
`[
`
`
`
`
`
`MINIMIZE TOTAL ENERGY IN TRACE TO SOLVE
`FOR w USING AN ITERATIVE LOOP
`
`416
`
`
` INSERT FINAL w TO COMPUTE
`P=(D-S]
`[
`I-w7'aA]-? Su)
`
`FIG. 4
`
`Page 5 of 13
`
`Page 5 of 13
`
`
`
`5,986,973
`
`1
`ENERGY MINIMIZATION IN SURFACE
`MULTIPLE ATTENUATION
`BACKGROUNDOF THLE INVENTION
`1. Field of the Invention
`
`The present invention relates generally to marine seismic
`surveying and, more particularly, to a method for attenuating
`the effect of surface multiples in a marine seismic signal.
`2. Description of the Related Art
`Seismic surveying is a method for determining the struc-
`ture of subterranean formations in the earth. Seismic sur-
`veying typically utilizes seismic energy sources which gen-
`erate seismic waves and seismic receivers which detect
`
`seismic waves. The seismic waves propagate into the for-
`mations in the earth, where a portion of the wavesreflects
`from interfaces between subterranean formations. The
`amplitude and polarity of the reflected waves are determined
`by the differences in acoustic impedance between the rock
`layers comprising the subterranean formations. The acoustic
`impedance of a rock layer is the product of the acoustic 5,
`propagation velocity within the layer and the density of the
`layer. The seismic receivers detect
`the reflected seismic
`waves and convert the reflected waves into representative
`electrical signals. The signals are typically transmitted by
`electrical, optical, radio or other means to devices which
`record the signals. Through analysis of the recorded signals,
`the shape, position and composition of the subterranean
`formations can be determined.
`
`Marine seismic surveying is a method for determining the
`structure of subterranean formations underlying bodies of 3
`water. Marine seismic surveying typically utilizes seismic
`energy sources and seismic receivers located in the water
`which are either towed behind a vessel or positioned on the
`water bottom from a vessel. The energy source is typically
`an explosive device or compressed air system which gen-
`erates seismic energy, which then propagates as seismic
`wavesthrough the body of water and into the earth forma-
`tions below the bottom of the water. As the seismic waves
`strike interfaces between subterranean formations, a portion
`of the seismic waves reflects back through the earth and
`water to the seismic receivers, to be detected, transmitted,
`and recorded. The seismic receivers typically used in marine
`seismic surveying are pressure sensors, such as hydro-
`phones. Additionally, though, motion sensors, such as geo-
`phonesor accelerometers may be used. Both the sources and
`receivers may bestrategically repositioned to cover the
`survey area.
`Seismic waves, however, reflect from interfaces other
`than just those between subterranean formations, as would
`be desired. Seismic wavesalso reflect from the water bottom
`
`ey}wn
`
`40
`
`and the water surface, and the resulting reflected waves
`themselves continueto reflect. Waves whichreflect multiple
`times are called “multiples”. Waves which reflect multiple
`times in the water layer between the water surface above and
`the water bottom beloware called “water-bottom multiples”.
`Water-bottom multiples have long been recognized as a
`problem in marine seismic processing and interpretation, so
`multiple attenuation methods based on the wave equation
`have been developed to handle water-bottom multiples.
`However, a larger set of multiples containing water-hottom
`multiples as a subset can be defined. The larger set includes
`multiples with upward reflections from interfaces between
`subtcrrancan formations in addition to upward reficctions
`from the water bottom. The multiples in the larger set have
`in commontheir downwardreflections at the water surface
`and thus are called “surface multiples”. FIG. 1, discussed
`below, provides examples of different types of reflections.
`
`60
`
`Page 6 of 13
`
`2
`FIG. 1 shows a diagrammatic view of marine seismic
`surveying. The procedure is designated generally as 100.
`Subterranean formations to be explored, such as 102 and
`104, lic below a body of water 106. Seismic energy sources
`108 and seismic receivers 110 are positioned in the body of
`water 106, typically by one or more seismic vessels (not
`shown). A seismic source 108, such as an air gun, creates
`seismic wavesin the body of water 106 and a portion of the
`seismic waves travels downward through the water toward
`the subtcrrancan formations 102 and 104 bencath the body
`of water 106. When the seismic waves reach a seismic
`reflector, a portion of the seismic wavesreflects upward and
`a portion of the seismic waves continues downward. The
`seismic reflector can be the water bottom 112 or one of the
`interfaces between subterranean formation, suchasinterface
`114 between formations 102 and 104. Whenthe reflected
`wavestravelling upward reach the water/air interface at the
`water surface 116, a majority portion of the wavesreflects
`downward again. Continuing in this fashion, seismic waves
`can reflect multiple times between upwardreflectors, such as
`the water bottom 112 or formation interfaces below, and the
`downward reflector at
`the water surface 116 above, as
`described more fully below. Each time the reflected waves
`propagate past the position of a seismic receiver 110, the
`receiver 110 senses the reflected waves and generates rep-
`resentative signals.
`Primary reflections are those seismic waves which have
`reflected only once,
`from the water bottom 112 or an
`interface between subterranean formations, before being
`detected by a seismic receiver 110. An exampleof a primary
`reflection is shown in FIG. 1 by raypaths 120 and 122.
`Primaryreflections contain the desired information about the
`subterranean formations whichis the goal of marine seismic
`surveying. Surface multiples are those waves which have
`reflected multiple times between the water surface 116 and
`any upward reflectors, such as the water bottom 112 or
`formation interfaces, before being sensed by a receiver 110.
`An example of a surface multiple which is specifically a
`water bottom multiple is shown by raypaths 130, 132, 134
`and 136. The surface multiple starting at raypath 130 is a
`multiple of order one, since the multiple contains one
`reflection from the water surface 116. Two examples of
`general surface multiples with upwardreflections trom both
`the water bottom 112 and formation interfaces are shown by
`raypaths 140, 142, 144, 146, 148 and 150 and by raypaths
`160, 162, 164, 166, 168 and 170. Both of these latter two
`examples of surface multiples are multiples of order two,
`since the multiples contain two reflections from the water
`surface 116. In general, a surface multiple is of order 1 if the
`multiple contains i reflections from the water surface 116.
`Surface multiples are extraneous noise which obscures the
`desired primary reflection signal.
`Surface multiple attenuation is a prestack inversion of a
`recorded waveficld which removesall orders of all surface
`multiples present within the marine seismic signal. Unlike
`some wave-equation-based multiple-attenuation algorithms,
`surface multiple attenuation does not require any modeling
`of or assumptions regarding the positions, shapes andreflec-
`tion cocfficicnts of the multiple-causing reflectors. Instead,
`surface multiple attenuation relies on the internal physical
`consistency between primary and multiple events that must
`exist in any properly recorded marine data set. The infor-
`mation neededfor the surface multiple attenuation processis
`already contained within the scismic data.
`In the following discussion,
`let
`the original seismic
`wavefields,
`the corresponding recorded data sets, or the
`corresponding data cubes or matrices be represented bythe
`
`Page 6 of 13
`
`
`
`5,986,973
`
`4
`can be found, then the primary wavefield P, free of surface
`multiples, can be computed directly from the recorded
`wavelield D. The operator O being suitable means that the
`operator O must be both geophysically and mathematically
`plausible. The operator O being geophysically plausible
`meansthat the operator O satisfies Eq. (4). The operator O
`being mathematically plausible mcansfirstly that the fac-
`torizations in Eq. (5) are valid and secondly that the inverse
`of the bracketed expression in Eq.(7) exists and thus Eq.(8)
`is valid. The Kirchhoff integral, a mathematical statement of
`Huygens’ principle, honors the wave equation and can be
`used as the basis of a gcophysically suitable operator O.
`Use of the Kirchhoff integral provides the appropriate
`two- or three-dimensional generalization of the inverse of
`the recursive forward modeling equation for surface mul-
`tiple attenuation, as given by Eq. (8). The Kirchhoff integral
`must be made compatible with Eqs. (4) through (8). First,
`the recorded marine seismic data are Fourier transformed
`from the time domain to the frequency domain. Let indi-
`vidual seismic traces or events within the wavefields or data
`
`sets be represented by lower-case letters. Thus m,, is a
`multiple cvent of order i within a scismic trace d in the
`wavefield D. Let p and m represent single-frequency com-
`ponents of Fourier-transformed seismic traces. For example,
`m(s,r) is one frequency component of the seismic trace
`whose source and receiver were at positions s and 1,
`respectively, and which contains only surface multiples of
`order i. Let my; represent m_, after being modified to
`include the scale and phase corrections and the obliquity
`factor required by the Kirchhoff integral. The Kirchhoff
`modification is given by
`
`muy i-1@, = (1 - pVe/ 4a cosfsinUV /wa1, 4),
`
`(9)
`
`where
`
`x=inline coordinate,
`j=(-1)!2
`=angular frequency,
`k,=x-component of wavenumber vector, and
`V=speed of sound in water.
`Because of k,, the modification of m,_, is dip-dependent. In
`the frequency domain, the Kirchhoff integral can be written
`as
`
`J w
`
`m(sr)=— Jp(s, X)ayia (xr)dx.
`
`(10)
`
`The minussign is due to the negative reflection coefficient
`of the water surface.
`In practice, recorded wavefields are not continuous in x,
`so the Kirchhoff integral in Eq. (10) hasto be replaced by the
`following discrete summation over x
`
`MAS,P)=-ZpP(SX)May1G).
`
`(11)
`
`(11) is the formula for
`Except for the minus sign, Eq.
`computing an elementof the product of two matrices. Thus,
`define matrix M,_,, as the matrix whose columns are the
`common-receiver records, m,_,(x,r), define matrix M,,,_, as
`the matrix whose columns are the Kirchhoff-modified
`common-receiver records, My,,_,(x,r), and define matrix P
`as the matrix whose rows are the common-shot records,
`p(s,x). Then Eq. (11) becomes
`
`3
`same upper-case letters. Thus let D represent a marine
`seismic data set corresponding to a wavefield D. The wave-
`field D can be divided into twoparts,
`
`D=P+M
`
`@)
`
`‘lhe primary wavefield,P, represents that portion of D which
`contains no surface multiples. The surface multiples
`wavefield, M, represents that portion of D which contains
`surface multiples of any order. Surface multiple attenuation
`is a processing method for removing the multiples wavefield
`M trom the recorded wavefield D to yield the desired
`primary wavefield P.
`Tor each i from 1 to 2%, let M; represent that portion of M
`containing surface multiples of order i. Then the surface
`multiple wavefield M can be further decomposed into an
`infinite sum of different orders,
`
`M=M,+Mo+...4+Mp...
`
`(2)
`
`Recorded data sets have a finite duration, so only a finite
`number of terms from Eq. (2) are needed to represent the
`corresponding wavefield. Substituting an appropriately trun-
`cated Eq. (2) into Eq. (1) yields
`
`i}wn
`
`D=P+M,+M>+ ...+M,,
`for some value n.
`
`(3)
`
`The process of surface multiple attenuation assumes that
`surface multiple events M, of order i can be predicted from
`knowledgeof both the surface multiple events M,_, of order
`i+1 and the primary wavefield P. This assumption meansthat
`there exists some mathematical operator O such that
`
`M=POM,...
`
`(4)
`
`Inserting Eq. (4) into Eq. (3) and factoring out first P and
`then O yields
`
`ey}wn
`
`40
`
`D=P+POP+ POM, +...+ POM,1
`
`= P[L+O(P+ My, +...+ My-1)].
`
`Define a truncated version D; of D by
`
`Dp =P+M,+...4+M,-1
`
`=D-M,,
`
`(5)
`
`(6)
`
`In practice, as will be discussed later, D, would be approxi-
`mated by truncating the seismic traces in D in time rather
`than actually constructing and subtracting M,, from D.
`Inserting Eq.
`(6) into Eq.
`(5) yiclds the compact form
`D=P[1+0D,]+tm(7)
`
`Eq. (7) is a formula for recursive forward modeling of
`surface multiples. Eq. (7) represents adding the events of
`order n to the wavefield containing all events up to and
`including order n-1. If the bracketed expression in Eq. (7)
`has an inverse, then Eq. (7) can be inverted to yield
`
`60
`
`P=D[1+OD,}'.
`
`(8)
`
`Eq. (8) is the inverse of the recursive forward modeling
`equation, Eq. (7). Eq. (8) states that if a suitable operator O
`
`M=PMy,;-1-
`
`(12)
`
`Page 7 of 13
`
`Page 7 of 13
`
`
`
`5
`Since the matrix indices are the shot and receiver
`coordinates, the zero-offset seismic traces lie along the main
`diagonal of each matrix.If the operator O in Eq. (4) is matrix
`multiplication and the quantities in uppercase are matrices,
`then Eq. (4) becomes Eq. (12) and so Eq. (8) becomes
`
`P=D[I-Dy|".
`
`(13)
`
`is the identity matrix and the “-1” superscript
`where I
`indicates matrix inversion.
`Eq. (13) provides a simple algorithm for two-dimensional
`surface multiple attenuation for idcal data. By “idcal” is
`meant that the wavefield is recorded broadband, contains no
`noise, has all wavelet effects, including source and receiver
`ghosts, removed, and has a trace-offset range that begins at
`zero offset. Furthermore, each individual sample within the
`idcal data sct D must have a truc relative amplitude with
`respect to every other sample within D. For non-ideal data,
`wavelet effects have not been removed and are the major
`factor which must be taken into consideration. In theory, this
`problem is easily fixed by redefining operator O to include
`a convolution by the inverse of the source wavelet w. Eq.
`(13) becomes
`
`P=D[I-w"D,J",
`
`(14)
`
`where w'' is the wavelet inverse. Since Eq. (14) is in the
`frequency domain, convolution is accomplished by multi-
`plication.
`the source wavelet w is initially 3
`In practice, however,
`unknown. The source wavelet w can be found by minimiz-
`ing the total energy in P in Eq. (14). Thus surface multiple
`attenuation becomes an L,-norm minimization problem,
`which has standard methodsoffinding the solution, such as
`the conjugate gradient technique. One could also minimize
`other measures of the surface multiple energy in P.
`The computation of P as described above requires com-
`puting the source wavelet w by minimizing the total energy
`in Pin aniterative loop. This requires computing the inverse
`of the matrix quantity in brackets in Eq. (14) during each
`iteration for a new value of w. Computing the inverse of the
`matrix may be accomplished by a series expansion or an
`exact matrix inversion. The series expansion approach is
`approximate and can be slowly converging. Previous exact
`approaches to surface multiple attenuation require a large
`number of inversions of large square matrices. Square
`matrix inversion is computationally expensive.
`US. Pat. No. 5,587,965 to Dragoset & Jericevic, titled
`“Surface Multiple Attenuation Via Eigenvalue Decomposi-
`tion” describes a method using eigenvalue decomposition in
`an exact approach to solving Eq. (14) which is specially
`designed to make computing the matrix inverse in each
`iteration muchless costly. The eigenvalue decomposition of
`matrix D,, is computed, yielding
`
`ey}wn
`
`40
`
`5,986,973
`
`P=D § [-wtapts,
`
`(17)
`
`Now the expression in brackets in Eq. (17) is a diagonal
`matrix, since both I and A are diagonal matrices and w~ is
`a complex scalar. ‘hus Eq. (17) can also be written as
`
`10
`
`P=D S diag [1/U-wA)]s™*.
`
`The expression in brackets in Eq. (17) can be inverted many
`times at low cost.
`This derivation is described in U.S. Patent No. 5,587,965
`to Dragoset & Jericevic. In their method, the matrix product
`D S and the inverse S~' are computed only once and then
`saved in memory. After that initial cost, the incremental cost
`of recomputing Eq. (17) is that of a matrix multiplication
`and a diagonal matrix inversion, whereas the incremental
`cost of recomputing Eq. (14) is that of a matrix multiplica-
`tion and a square matrix inversion. If the dimension of the
`square matrices involved is N, then recomputing Eq. (14)
`takes approximately 2N® order of magnitude floating point
`operations for each iteration, while recomputing Eq. (17)
`takes only approximately (N°+N)order of magnitude opera-
`tions for each iteration but the first.
`
`SUMMARYOF THE INVENTION
`
`The present invention is a method for substantiallyattenu-
`ating surface multiples from a marine seismic signal. ‘The
`marine seismic signal is truncated in time. Both the marine
`seismic signal and the truncated signal are transformed from
`time domain into frequency domain and represented by
`matrices D and D,, respectively. Next, eigenvalue decom-
`position D;=SAS™ is computed. Here A is the diagonal
`matrix whose elements are the eigenvalues of D;, and S is
`the square matrix whose columnsare the corresponding
`eigenvectors of D;. S7'
`is the matrix inverse of S. Matrix
`product D S is computed and its conjugate transpose [D S]*
`is computed. Matrix product [D S]*[D S] is computed and
`saved in memory. Matrix inverse S~* is computed andits
`conjugate transpose (S")* is computed. Matrix product
`(S-*\(S"1)* is computed and saved in memory. Aqinitial
`estimate for the source wavelet w is made. Diagonal matrix
`[I-w-tA] is computed, matrix inverse [I-w*A]" is
`computed, and its conjugate transpose [(I-w7'A)~*]* is
`computed. The matrix products [D S]*[D S] and (S~')(S~1)*
`are retrieved from memory and the matrix product
`
`SS4*[(-wtay*}(D SLD Slewtay*
`
`is computed. The source wavelet w is optimized by mini-
`mizing the total energy in the matrix trace of the matrix
`product
`
`Dy=SES1.
`
`(15)
`
`{S(S)*[(-wAyD SID ST-wAT},
`
`Ilere A is the diagonal matrix whose diagonal elements are
`the eigenvalues ., of D,,, that is, A=diag [A,] . S is the square
`matrix whose columns are the corresponding eigenvectors of
`D,, and S~ is the matrix inverseof S. Inserting Eq. (15) into
`Eq. (14) yields
`
`60
`
`using an iterative loop. Then primary matrix P is computed
`by inserting the optimized value for the source wavelet w
`into the expression P=[D S][I-w-tA]'‘S~1. Finally, matrix P
`is inverse transformed from the frequency domain into the
`time domain.
`
`BRIEF DESCRIPTION OF THE DRAWINGS
`
`P=D{I-w *S AS *] +.
`
`(16)
`
`Substituting I=S S~* and factoring S and S~* outto the left
`and right, respectively, allows one to rewrite Eq. (16) as
`
`A better understanding of the benefits and advantages of
`the present invention may be obtained from the appended
`detailed description and drawing figures, wherein:
`
`Page 8 of 13
`
`Page 8 of 13
`
`
`
`5,986,973
`
`8
`
`E=ai{psu-wtay'sipsa —wilaytsty}
`
`(20)
`
`=i{DSU—w AySse)" [U-wlay] sy}
`
`7
`FIG. 1 is a diagrammatic view of marine seismic
`surveying, showing the source of surface multiples;
`FIG. 2 is a schematic diagram of the general method of
`surface multiple attenuation;
`5
`= m{soty|a—wtay](sypsu-wt ay}
`FIG. 3 is a schematic diagram ofthe preliminary steps of
`the method of the present invention for surface multiple
`attenuation; and
`FIG. 4 is a schematic diagram of the remaining steps
`comprising the iterative loop of the method of the present
`invention for surface multiple attenuation.
`DESCRIPTION OF THE PREFERRED
`EMBODIMENTS
`
`Substituting the matrices A and B defined by
`
`10
`
`and
`
`Asst)
`
`B=(D 8)*D S=S*D*D S,
`
`the energy expression given by Eq. (20) can be rewritten as
`
`Comparison of computational paths for prior algorithm,
`Eq. (14), and eigenvalue decomposition algorithm, Eq. (17),
`reveals that for all but the first iteration in the optimization
`of w7', the N? process involved in the construction and the
`inversion of matrix (I-w''D,,) in the prior algorithm is
`traded for N multiplications and divisions in the eigenvalue
`decomposition algorithm. However, both algorithms still
`contain the full matrix multiplication operation in each
`iteration, which is again an N°?process.
`‘lhe reason forthis latter computational expenseis that the
`total energy E in primary matrix P is computed by multi-
`plying each element in P by its complex conjugate and
`summing, all the products
`isn
`jeN
`'
`= » d-wlaySS’a,Ja-wtay|e,
`2
`J N
`=I
`Fl
`(18)
`E=
`P:, gPh i
`iz
`
`E=tr{A [(-w *A) *]*BU-w +A) *}.
`
`(21)
`
`Define the elements of matrices A and B by A={a,,} and
`B={b,;}. Then, in explicit form, Eq. (21) is given by
`
`E=
`
`i=N j=N
`f=1 j=l
`
`$
`a, JQ -wlay|b,- way
`
`(22)
`
`where the dagger t indicates the complex conjugate of the
`matrix element P,,. The problem with using Eq. (18) as the
`energy formula is that P, being a function of w~*, must be
`explicitly evaluated in each iteration. Thus each iteration
`requires a full matrix multiplication with N* multiplications.
`If, however,
`the energy minimization problem is
`expressed in linear algebra terms using properties of a trace
`of a matrix, designated tr{X} for any matrix X, then the
`computations can be considerably simplified. The trace of a
`matrix product is defined as
`
`i=N CN
`mX VYp=a Sxjpn
`1 jel
`
`for any matrices X={x,,} and Y={y,;}. Let the asterisk *
`denote conjugate transpose, or adjoint, of a matrix. Then the s
`energy expression given by Eq. (18) can be rewritten as
`
`R=tr{P P*}.
`
`(19)
`
`The conjugate transpose of a matrix product obeys the
`following product rule
`
`The above expressions in Eq. (22) reveal that for all but the
`first iteration, the number of operations in evaluating the
`energy E can be reduced from N* to 3N*. The number of
`operations can be reduced even further to 2N?,if the product
`of the terms a,; and b;; are precomputed once andstored.
`Alternatively, the optimization problem can be solved by
`finding the root of the following equation
`
`40
`
`OE
`=0
`an 7
`
`In the case in which matrix rows are tapered to remove
`edge effects, the approach remains the same, except that
`matrix B is changed. Let T be a tapering matrix, which is real
`and diagonal. Using the definition of the primaries matrix P
`given in Eq. (17), the energy expression of Eq. (19) now
`becomes
`
`E=0TPTPY}
`
`(23)
`
`= o{{rpsu- wt aytsirpsu— waysty}
`
`= {TDs -witay's-lis-ly[u-wtay"| (sy7}
`
`= assisty[a= wtay| wsy reps wl ay}.
`
`(X Y)"=¥* X*
`
`Substituting the matrix A from above and the matrix B,
`defined by
`
`60
`
`for any matrices X and Y. The trace of a matrix product
`obeys the following circularity rule
`
`B=(D $)*f T D S=S*D*T TD S,
`
`tr{X ¥ Zlate{¥ Z X}=tr{Z XY}
`
`then the cnergy expression given by Eq. (23) can be written
`as
`
`for any matrices X, Y and Z. Using these two rules, the
`energy expression given by Eq. (19) can be rewritten as
`
`E=tr{A[U-wtA)1Bewtay}.
`
`Page 9 of 13
`
`Page 9 of 13
`
`
`
`5,986,973
`
`10
`
`15
`
`9
`This is the same as Eq. (21) above, except that matrix B is
`replaced by matrix B,.
`In the case in which the matrix columnsare tapered to
`remove edge effects, the approach remains the same, but
`nowmatrix A is changedinstead of matrix B. Define matrix
`A, by
`
`A-ST 1(S*.
`
`Then the energy expression of Eq. (19) becomes
`
`E=tr{PT(PT)}
`
`=in{ps—w ‘ays ‘T{psu-w tay's '1P}
`= ilpsa-w Ayts Ts! ny
`fu-wtay'Tiwsyt}
`
`{ {
`
`= rlsrnsty[uwtay] sy psu wt Ay}.
`1s
`=a d-wtay] BU-wtay}.
`{
`
`Again, this is the same as Eq. (21) above, except that now
`matrix A is replaced by matrix A,.
`All of the seismic traces in a data set are Fourier-
`transformed and inserted into a data cube, D(s,1,f) . Here s
`is source position, r is receiver position, and f is frequency.
`Next, the original seismic traces are also truncatedin time,
`Fouricr-transformed, Kirchhoff-modificd, and inserted into
`another data cube, D,,s..f). For each frequency f, matrices
`T(s,r) and D,,(s,r) are extracted from the data cubes N(s,1,f)
`and D,/(s,r,f), respectively. An initial estimate for the value 3
`of the source wavelet w component for a particular fre-
`quency is madeandinserted into the energy expression Eq.
`(20) along with D¢s,r) and D,,fs,r). Eq.
`(20) gives an
`expression for the total energy in the primary matrix P(s,r)
`as a function of the complex scalar w. A value for the
`wavelet w is computed by minimizing the total energy in the
`expression for P(s,r). This minimizationis done byiterating
`through values of wavelet w in Eq.
`(20). Inserting the
`computed value for w into Eq. (14) yields matrix P(s,r),
`which is inserted into an output data cube, P(s,r,f). Eq. (20)
`has to be solved for all frequencies independently. Finally,
`each seismic trace in data cube P(s,r,f) is inverse Fourier
`transformed and reorganized into gathers.
`FIG. 2 shows a schematic diagram of the general method
`for attenuating surface multiples from a marine seismic
`signal. The method is designated generally as 200. A two-
`dimensional marine data set is recorded as a seismic signal,
`as shownin Block 202. The seismic traces which comprise
`the data sel are truncated in lime, as shown in Block 204.
`Both the data set from Block 202 and the truncated seismic
`traces from Block 204 are Fourier transformed from the time
`
`40
`
`10
`matrix P(s,r) is computed from matrices D(s,r) and D,,(s,n),
`as shown in Block 220. The computation of the matrix P in
`Block 220 involves the solving of Eq. (8) as given by Eq.
`(20) and will be discussed in greater detail below.
`In
`particular, the present invention is an improved method for
`computing P. For each frequency f,
`the matrix P(s,r) is
`inserted into data cube P(s,r,f), as shown in Block 222, thus
`constructing the primary data cube P. Finally, the data cube
`P is inverse Fourier transformed from the frequency domain
`into the time domain, as shown in Block 224.
`FIGS.3 and 4 show schematic diagrams of the method of
`the present invention for computing primary matrix P(s,r), as
`designated in Block 220 of FIG. 2, for non-ideal data. The
`preliminary steps of the method are shownin FIG.3 and are
`designated generally as 300. Eigenvalue decomposition
`Dy=S A S$“of the modified matrix D,(s,r) is computed, as
`shown in Block 302. Here the columnsof matrix S(s,r) are
`constructed from the eigenvectors of modified matrix D,,
`and matrix A(s,r) is a diagonal matrix whose diagonal
`elements are the corresponding eigenvalues of modified
`matrix D,,. Matrix product D S is computed once, conjugate
`transpose [D S]* is computed oncc, and matrix product [D
`S}*[D S]is computed once, as shownin Block 304. Thenthe
`matrix product [D S]*[D S] is saved in memory, as shown
`in Block 306. Matrix inverse S~* of the matrix S is also
`computed once, conjugate transpose [D S] is computed
`once, and matrix product S-"(S~)* is computed once, as
`shown in Block 308. The matrix product S-*(S~1)* is also
`saved in memory, as shownin Block 310. The point at which
`the schematic diagram of FIG. 3 continues onto FIG. 4 is
`shown by Block 312.
`The remaining steps comprising the iterative loop in the
`methodof the present invention for computing matrix P(s,r)
`are shownin FIG. 4 and are designated generally as 400. The
`5 point at which the schematic diagram of FIG. 4 continues
`from Block 320 of FIG. 3 is shown by Block 402. An initial
`estimate for the value of the complex scalar representing the
`source wavelet w is made, as shown by Block 404. Diagonal
`matrix [I-w-*A] is computed, as shown in Block 406. Here
`Lis the identity matrix. The matrix inverse [I-w*A]? of the
`diagonal matrix [I-w~‘A] is computed, as shown in Block
`408. Conjugate transpose [I-w-'A]+ of the matrix inverse
`[d-w-1A)-"]* is computed, as shown in Block 410. The
`matrix product [D S]*[D S] and the matrix product S-'(S~
`1)* are retrieved from memory, as shown in Block 412. An
`expression for the energy in the primary matrix P(s,r) is
`computed from the matrix trace of the matrix product given
`by Eq. (20),
`
`
`
`E=ir{S(S)*[(-wAy?]}*[D SD 8] U-wtay},
`
`domain into the frequency domain, as shown in Blocks 206
`and 208, respectively. The transformed truncated seismic
`traces are modified in the frequency domain to include the
`obliquity, scale and phase factors required by the Kirchhoff
`integral and given by Eq. (9), as shown in Block 210. The
`transformed data are used to construct a data cube D(s,r,f) ,
`as shown in Block 212. Here s is source location,
`r is
`receiver location, and f is frequency. The modified seismic
`traces are used to construct another data cube, the modified
`data cube D,(s,rf), as shown in Block 214. Next, the series
`of steps in Blocks 216, 218, 220 and 222 are repeated for
`each value of frequency f in data cube D(s,r,f). Thus, for
`each frequency f, a matrix D(s,r) is extracted from data cube
`D(