United States Patent
`
`[19]
`
`[11] Patent Number:
`
`5,986,973
`
`Jeriéevié et al.
`[45] Date of Patent: Nov. 16, 1999
`
`
`
`USOOS986973A
`
`ENERGY MINIMIZATION IN SURFACE
`MULTIPLE ATTENUATION
`
`[57]
`
`ABSTRACT
`
`[54]
`
`[75]
`
`[73l
`
`Inventors: Zeljko Jericevié; William H.
`Dragoset, Jr., both of Houston, Tex.
`
`Assigncc: Western Atlas International, Inc.,
`Houston, Tex.
`
`[21]
`
`Appl. No.: 08/923,150
`
`[22]
`
`[51]
`[52]
`
`[58]
`
`[56]
`
`Filed:
`
`Sep. 4, 1997
`
`Int. Cl.5 ....................................................... G01V 1/38
`US. Cl.
`................................. 367/24; 367/45, 367/24;
`367/53
`Field of Search .................................. 367/21, 24, 38,
`367/45, 53
`
`References Cited
`
`U.S. PATENT DOCUMENTS
`
`5,138,583
`5,448,531
`5,587,965
`5,719,821
`5,729,506
`
`............................. 367/21
`8/1992 Wason et a1.
`9/1995 Dragoset et a1.
`......................... 367/21
`12/1996 Dragoset, Jr. et a1.
`367/24
`2/1998 Sallas et a1.
`.............................. 367/21
`3/1998 Dragoset et a1.
`......................... 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 DT Eigenvalue
`decomposition DEEM—1 of matrix Dr 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'1 is computed and saved in memory.
`Conjugate transpose (S'l)* is computed and saved in
`memory. Matrix product S'1(S'1)* is computed and saved in
`memory. An initial estimate of the source wavelet W is made.
`Source wavelet W is optimized by iterating the steps of
`computing the diagonal matrix [I—W_l/\], computing matrix
`inverse [I—W‘lATl, computing conjugate transpose [(I—W‘
`1A)'1]*, retrieving matrix products [D S]*[D S] and S'1(S'
`1)* from memory, and minimizing the total energy in trace
`of matrix product
`
`S’1(S’1)*[(I—w’1A)’1]*[D S][D S][I—w’1A]’1.
`
`Primary matrix P representing the wavefield free of surface
`multiples is computed by inserting computed value for W
`into the expression [D S][I—W'1A]'lS'l. Primary matrix P is
`inverse transformed from frequency domain into time
`domain.
`
`17 Claims, 4 Drawing Sheets
`
`100 \
`108
`108
`108
`
`
`
`
`
`
`:0:
`110
`:0:
`‘
`110
`:
`-
`-
`
`"
`130
`135
`140
`
`132
`134
`142
`
`
`120
`v
`>//
`\\><\
`’\//\
`\ \\
`w/
`>>/>\>
`>>/\\>
`\\/\\
`,</>//
`\
`/ /
`>//\\>/
`
`e
`’
`I
`\\‘\\ \\\\\\\/\
`,
`,
`>//\//\//\//\//\//\//\// // // // //\// / / // // // /\//\//\Z/\// // // // // ///~/// // // // //\ /,
`
`
`
`®®®®®®®®®®®®®®®®®® x\\/>$\/x\\/,\\\/x\\/x\\>x\/x\/>\\/x\/>§A\A\/x\>x\/t
`
`Page 1 of 13
`
`SAMSUNG EXHIBIT 1029
`
`Page 1 of 13
`
`SAMSUNG EXHIBIT 1029
`
`

`

`US. Patent
`
`Nov. 16, 1999
`
`Sheet 1 0f 4
`
`5,986,973
`
`m:\.
`at>—929:
`
`
`mm:3:
`
`
`
`
`
`
`
`
`/////////////.9?9s<\<\<\<\<\<<\<\<\\\\\\\\\\\\\\\\\\>//\////\//\//\//\//\//\//\////\///\///\//////////////////////
`
`
`
`
`
`<\<9E¥§®®®®®®®®®®®®\\>,
`
`
`
`
`
`wvA/vgyvéavéyVyyyyy/v/Vyyav»,///\\///\/\//\\///\//\/\/\/\/\/\/\/\/\/\/\/\/\/\x
`
`X::,\\//\\\//\\\/\\\/\\\/\\\/\\\/\\\/\\\/\\\/\\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\/\\A
`
`
`:hii\\\\\_‘.\\\>
`QW»»&§%§9§»§999§§
`
`
`
`
`
`
`4&)¢&¢¢u.99v\\/\\y/\\/\\/\\/\\/\\/\\/\\\\\\\\\\\\\\\/\/\/\/\/\/\/\/A\/\/\/\/\/\/\/\/\/\/\/\/>
`
`H.UNR
`
`
`
`
`
`I
`
`
`
`Page 2 of 13
`
`Page 2 of 13
`
`
`

`

`US. Patent
`
`Nov. 16,1999
`
`Sheet 2 0f4
`
`5,986,973
`
`200~\\
`
`INPUT DATA
`
`202
`
`205
`
`FOURER TRANSFORM DATA
`
`204
`
`TRUNCATE TRACESIN
`WME DOMAW
`
`208
`
`FOURER TRANSFORM
`TRUNCATED TRACES
`
`
`
`
`
`
`
`
`MODFY TRACESIN
`FREQUENCY DOMAW
`
`
`210
`
`214
`
`212
`
`
`
`216
`
`218
`
`CONSTRUCT MCDWTED
`CONSTRUCT DATA
`
`DATA CUBE DM(srf)
`CUBE D(&rj)
`
`
`
`EXTRACT MAme DMCLG
`EXTRACT MAme D(&r)
`
`FROM DATA CUBE DM(&rj)
`FROM DATA CUBE D(&rj)
`
`
`
`COMPUTE MWMARY MATRm P(&r)
`
`222
`
`INSERT MATWX P(&r)|NTO DATA CUBE P(srf)
`
`224
`
`INVERSE FOUWER TRANSFORM P
`
`FIG. 2
`
`Page 3 of 13
`
`Page 3 of 13
`
`

`

`US. Patent
`
`Nov. 16,1999
`
`Sheet 3 0f4
`
`5,986,973
`
`300 \
`
`COMPUTE EIGENVALUE DECOMPOSITION DII= S /\ 3-1
`
`302
`
`COMPUTE PRODUCT [D-S], ADJOINT [D-S]*
`AND PRODUCT [D-S]*'[D-S]
`
`304
`
`306
`
`STORE PRODUCT [D S]*-[D-S] IN MEMORY
`
`COMPUTE INVERSE 8-1, ADJOINT (S-1)*
`AND PRODUCT (8“) . (s-')*
`
`308
`
`310
`
`STORE PRODUCT (s-‘)- (s-‘)* IN MEMORY
`
`»N|
`
`lL
`
`CONTINUE TO STEP 402 (FIG. 4
`
`FIG. 3
`
`Page 4 of 13
`
`Page 4 of 13
`
`

`

`US. Patent
`
`Nov. 16,1999
`
`Sheet 4 0f4
`
`5,986,973
`
`400\
`
`CONTINUE FROM STEP 312 (FIGJ)
`
`402
`
`ESTIMATE WAVELET w
`
`404
`
`406
`
`COMRUTE DIAGONAL MATRIX [
`
`I—w—IA]
`
`408
`
`COMPUTE DIAGONAL MATRIX INVERSE [ I—w“/\]'1
`
`COMPUTE ADJOINT [(I—w”1/\)"‘]*
`
`470
`
`RETRIEVE [D S]* [0- S] AND
`(S ‘1) (S )* FRO M MEMORY
`
`472
`
`COMPUTE TRACE OF MATRIX PRODUCT
`TR £8“-<ST‘)* [<I—w‘IAI-‘]* [D-SI
`
`474
`
`
`
`
`WAT"?
`*'[D'S] [1—
`
`
`
`
`
`MINIMIZE TOTAL ENERGY IN TRACE TO SOLVE
`FOR w USING AN ITERATIVE LOOP
`
`476
`
`
`INSERT FINAL w To COMPUTE
`P=[D-S]
`[ I—w"/\]‘
`3-1
`
`
`
`
`
`FIG. 4
`
`Page 5 of 13
`
`Page 5 of 13
`
`

`

`5,986,973
`
`1
`ENERGY MINIMIZATION IN SURFACE
`MULTIPLE ATTENUATION
`BACKGROUND OF THE 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
`
`10
`
`seismic waves. The seismic waves propagate into the for-
`mations in the earth, where a portion of the waves reflects
`from interfaces between subterranean formations. The
`amplitude and polarity of the reflected waves are determined
`by the dilferences in acoustic impedance between the rock
`layers comprising the subterranean formations. The acoustic
`impedance of a rock layer is the product of the acoustic
`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
`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
`waves through 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-
`phones or accelerometers may be used. Both the sources and
`receivers may be strategically repositioned to cover the
`survey area.
`Seismic waves, however, reflect from interfaces other
`than just those between subterranean formations, as would
`be desired. Seismic waves also reflect from the water bottom
`
`and the water surface, and the resulting reflected waves
`themselves continue to reflect. Waves which reflect multiple
`times are called “multiples”. Waves which reflect multiple
`times in the water layer between the water surface above and
`the water bottom below are 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-bottom
`multiples as a subset can be defined. The larger set includes
`multiples with upward reflections from interfaces between
`subterranean formations in addition to upward reflections
`from the water bottom. The multiples in the larger set have
`in common their downward reflections at the water surface
`and thus are called “surface multiples”. FIG. I, discussed
`below, provides examples of different types of reflections.
`
`I\)m
`
`LA)LII
`
`40
`
`60
`
`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, lie 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 waves in the body of water 106 and a portion of the
`seismic waves travels downward through the water toward
`the subterranean formations 102 and 104 beneath the body
`of water 106. When the seismic waves reach a seismic
`reflector, a portion of the seismic waves reflects 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, such as interface
`114 between formations 102 and 104. When the reflected
`waves travelling upward reach the water/air interface at the
`water surface 116, a majority portion of the waves reflects
`downward again. Continuing in this fashion, seismic waves
`can reflect multiple times between upward reflectors, 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 example of a primary
`reflection is shown in FIG. 1 by raypaths 120 and 122.
`Primary reflections contain the desired information about the
`subterranean formations which is 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 upward reflections from 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 i 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 wavefield which removes all 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 and reflec—
`tion coefficients 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 needed for the surface multiple attenuation process is
`already contained within the seismic data.
`In the following discussion,
`let
`the original seismic
`wavefields,
`the corresponding recorded data sets, or the
`corresponding data cubes or matrices be represented by the
`
`Page 6 of 13
`
`Page 6 of 13
`
`

`

`5,986,973
`
`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 two parts,
`
`D=P+M
`
`(1)
`
`The 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 from the recorded wavefield D to yield the desired
`primary wavefield P.
`For each i from 1 to 00, let M,- represent that portion of M
`containing surface multiples of order i. Then the surface
`multiple wavefield M can be filrther decomposed into an
`infinite sum of different orders,
`
`10
`
`M=M1+M2+ .
`
`.
`
`. +M,-+ .
`
`.
`
`.
`
`(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\)m
`
`D=P+M1+M7+ .
`for some value n.
`
`+M
`m
`
`.
`
`.
`
`(3)
`
`The process of surface multiple attenuation assumes that
`surface multiple events M,- of order i can be predicted from
`knowledge of both the surface multiple events Mi,l of order
`i—l and the primary wavefield P. This assumption means that
`there exists some mathematical operator 0 such that
`
`Mi=POMi,1.
`
`(4)
`
`Inserting Eq. (4) into Eq. (3) and factoring out first P and
`then 0 yields
`
`LA)LII
`
`40
`
`D=P+POP+ P0M1 +...+P0Mn,1
`
`= P[1+ 0(P+ M1 +
`
`+ M,,,1)].
`
`Define a truncated version DT of D by
`
`DT=P+M1+...+M,,,1
`
`=D—M,,
`
`(5)
`
`(6)
`
`In practice, as will be discussed later, DT would be approxi-
`mated by truncating the seismic traces in D in time rather
`than actually constructing and subtracting Mn from D.
`Inserting Eq.
`(6) into Eq.
`(5) yields the compact form
`D=P[1+ODT]+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 11—1. If the bracketed expression in Eq. (7)
`has an inverse, then Eq. (7) can be inverted to yield
`
`60
`
`P=D[1+0DT]" .
`
`(8)
`
`4
`can be found, then the primary wavefield P, free of surface
`multiples, can be computed directly from the recorded
`wavefield D. The operator 0 being suitable means that the
`operator 0 must be both geophysically and mathematically
`plausible. The operator 0 being geophysically plausible
`means that the operator 0 satisfies Eq. (4). The operator 0
`being mathematically plausible means firstly 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 geophysically suitable operator 0.
`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 ml, is a
`multiple event of order i within a seismic trace d in the
`wavefield D. Let p and m represent single-frequency com-
`ponents of Fourier-transformed seismic traces. For example,
`mi(s,r) is one frequency component of the seismic trace
`whose source and receiver were at positions s and r,
`respectively, and which contains only surface multiples of
`order i. Let mMyl;1 represent m_1 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
`
`
`
`mM1171 (x, r) = (l 7 j)\/a),/47r eos[sin71(kXV/w)]m;,1(x, r),
`
`(9]
`
`where
`
`x=inline coordinate,
`j=(—1)“2
`ou=angular frequency,
`kx=x-component of wavenumber vector, and
`V=speed of sound in water.
`Because of k,” the modification of m,-_l is dip-dependent. In
`the frequency domain, the Kirchhoff integral can be written
`as
`
`mi(51rj=_JP(va)mM,i71 (x,r)dx.
`
`(10)
`
`The minus sign 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) has to be replaced by the
`following discrete summation over x
`
`mi(5'rj=_2p(stxij/I,i71(xvr)
`
`(11)
`
`(11) is the formula for
`Except for the minus sign, Eq.
`computing an element of the product of two matrices. Thus,
`define matrix Mi_1, as the matrix Whose columns are the
`common-receiver records, m,_1(x,r), define matrix MM_,._1 as
`the matrix whose columns are the Kirchhoff-modified
`common-receiver records, mM3,_1(x,r), and define matrix P
`as the matrix whose rows are the common-shot records,
`p(s,x). Then Eq. (11) becomes
`
`Eq. (8) is the inverse of the recursive forward modeling
`equation, Eq. (7). Eq. (8) states that if a suitable operator 0
`
`Mi=PMMY,-,1.
`
`(12)
`
`Page 7 of 13
`
`Page 7 of 13
`
`

`

`5,986,973
`
`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 0 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—DM]’1,
`
`(13)
`
`P=D s [I—w’lATiS’1.
`
`(17)
`
`Now the expression in brackets in Eq. (17) is a diagonal
`matrix, since both I and A are diagonal matrices and w’1 is
`a complex scalar. Thus Eq. (17) can also be written as
`
`is the identity matrix and the “—1” superscript
`where I
`indicates matrix inversion.
`
`10
`
`P=D S diag [1/(I—w’1kl)]S’1.
`
`Eq. (13) provides a simple algorithm for two—dimensional
`surface multiple attenuation for ideal data. By “ideal” 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
`ideal data set D must have a true 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 0 to include
`a convolution by the inverse of the source wavelet w. Eq.
`(13) becomes
`
`P=D[I—w’1DM]’1,
`
`I\)m
`
`(14)
`
`The expression in brackets in Eq. (17) can be inverted many
`times at low cost.
`This derivation is described in US. Patent No. 5,587,965
`to Dragoset & Jericevic. In their method, the matrix product
`D S and the inverse S"1 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 2N3 order of magnitude floating point
`operations for each iteration, while recomputing Eq. (17)
`takes only approximately (N3+N) order of magnitude opera-
`tions for each iteration but the first.
`
`SUMMARY OF TIIE INVENTION
`
`where w"1 is the wavelet inverse. Since Eq. (14) is in the
`frequency domain, convolution is accomplished by multi-
`plication.
`the source wavelet w is initially
`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 LZ-norm minimization problem,
`which has standard methods of finding 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 P in an iterative 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 rising eigenvalue decomposition in
`an exact approach to solving Eq. (14) which is specially
`designed to make computing the matrix inverse in each
`iteration much less costly. The eigenvalue decomposition of
`matrix DM is computed, yielding
`
`LA)LII
`
`40
`
`The present invention is a method for substantially attenti-
`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 DT, respectively. Next, eigenvalue decom-
`position DJ=SAS'1 is computed. Here A is the diagonal
`matrix whose elements are the eigenvalues of DT, and S is
`the square matrix whose columns are the corresponding
`eigenvectors of Dr. S"1 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'1 is computed and its
`conjugate transpose (S1)* is computed. Matrix product
`(S'1)(S'1)* is computed and saved in memory. An initial
`estimate for the source wavelet w is made. Diagonal matrix
`[I—w'lA]
`is computed, matrix inverse [I—w'lA]'1 is
`computed, and its conjugate transpose [(I—w'lA)'1]* is
`computed. The matrix products [D S]*[D S] and (S'1)(S'1)*
`are retrieved from memory and the matrix product
`
`S’1(S’1)*[(I—w’1A)’1]*[D S]*[D S][I—w’1A]’1
`
`is computed. The source wavelet w is optimized by mini-
`mizing the total energy in the matrix trace of the matrix
`product
`
`DM=SES’1.
`
`(15)
`
`"{se(Serra—weArlrtD S]*[D S][I-W’1A]’1},
`
`Here A is the diagonal matrix whose diagonal elements are
`the eigenvalues 7t,- of DM, that is, A=diag [71,-] . S is the square
`matrix whose columns are the corresponding eigenvectors of
`DM and S'1 is the matrix inverse of 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'lA]'lS'l. Finally, matrix P
`is inverse transformed from the frequency domain into the
`time domain.
`
`BRIEF DESCRIPTION OF THE DRAWINGS
`
`P=D[[—w 15 AS 1] 1.
`
`(16)
`
`Substituting I=S S'1 and factoring S and 3'1 out to 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
`
`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;
`FIG. 3 is a schematic diagram of the 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
`
`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 w'l, the N3 process involved in the construction and the
`inversion of matrix (I-W11DM) 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 N3 process.
`The reason for this latter computational expense is 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
`
`E =
`
`;:A' y:N
`
`[:1 i
`
`:1
`
`,
`P:.jPi,j
`
`(13)
`
`where the dagger t indicates the complex conjugate of the
`matrix element PU. The problem with using Eq. (18) as the
`energy formula is that P, being a function of w'l, must be
`explicitly evaluated in each iteration. Thus each iteration
`requires a full matrix multiplication with N3 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
`
`5
`
`10
`
`I\)m
`
`LA)LII
`
`40
`
`for any matrices X={x,-J-} and Y={y]—),-}. Let the asterisk *
`denote conjugate transpose, or adjoint, of a matrix. Then the
`energy expression given by Eq. (18) can be rewritten as
`
`F.=tr{P P*}.
`
`(19)
`
`The conjugate transpose of a matrix product obeys the
`following product rule
`
`8
`
`I: = ”({wu — w’lAJ’IS’IHDSU — w’lAflS’lF}
`
`(20)
`
`= n-{Dsu — WAWS" ts")*[(1— w*‘/\)’1]‘(DS)*}
`
`= zr{S’1(scl)‘[(1 — urlA)’1]*(DSJ*DS(1— urlA)’1}_
`
`Substituting the matrices A and B defined by
`
`A=S’1(S’1)
`
`and
`
`B=(D 33m S=S*D*D s,
`
`the energy expression given by Eq. (20) can be rewritten as
`
`E=tr{A [(I—w 1A) 1]*B(I—w 1A) 1}.
`
`(21)
`
`Define the elements of matrices A and B by A={aij} and
`B={b]-),-}. Then, in explicit form, Eq. (21) is given by
`
`+
`aid-[(1 — w’lxjy’l] bit-(1 — w’l/L-j’l
`
`(22)
`
`cud-[(1 — w’lijfl]lbj,i.
`
`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 N3 to 3N2. The number of
`operations can be reduced even further to 2N2, if the product
`of the terms aw- and bJ-J- are precomputed once and stored.
`Alternatively, the optimization problem can be solved by
`finding the root of the following equation
`
`615
`
`a = 0.
`
`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 = [r{ TP(TP)*}
`
`(23)
`
`[r{{TL)S[l — w’lAJ’IS’l}{wS(I — w’lA)’ls"1}*}
`
`rr{TDS(I — w’lA)"s*1(s*1)*[(I— w’lA)"]*(DS)*T}
`
`”{3’1 :_s"1)*[(1 — w’lA)’l]*(L)S)*r1Dsu— w’lAJ’I}.
`
`(X Y)*=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 S)*TTD S=S*D*T TD S,
`
`u-{x Y Z}=tr{Y z x}=u-{z x Y}
`
`then the energy cxprcssion 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[(1—w’lA)’1]*Bt(I—w’1A)’1j.
`
`Page 9 of 13
`
`Page 9 of 13
`
`

`

`5,986,973
`
`5
`
`10
`
`I\)m
`
`40
`
`10
`matrix P(s,r) is computed from matrices D(s,r) and DM(s,r),
`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 shown in FIG. 3 and are
`designated generally as 300. Eigenvalue decomposition
`DM=S A S'1 of the modified matrix DM(s,r) is computed, as
`shown in Block 302. Here the columns of matrix S(s,r) are
`constructed from the eigenvectors of modified matrix DM
`and matrix A(s,r) is a diagonal matrix whose diago ial
`elements are the corresponding eigenvalues of modified
`matrix DM. Matrix product D S is computed once, conjugate
`transpose [D S]* is computed once, and matrix product [D
`S]*[D S] is computed once, as shown in Block 304. Then he
`matrix product [D S]*[D S] is saved in memory, as shown
`in Block 306. Matrix inverse S"1 of the matrix S is also
`computed once, conjugate transpose [D S] is compu ed
`once, and matrix product S'1(S'1)* is computed once, as
`shown in Block 308. The matrix product S'1(S'1)* is also
`saved in memory, as shown in Block 310. The point at wh'ch
`the schematic diagram of FIG. 3 continues onto FIG. 4 is
`shown by Block 312.
`The remaining steps comprising the iterative loop in he
`method of the present invention for computing matrix P(s,r)
`are shown in FIG. 4 and are designated generally as 400. The
`LA) .11
`. 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'lA] is computed, as shown in Block 406. Here
`I is the identity matrix. The matrix inverse [I—w’1A]’1 of the
`diagonal matrix [I—w'lA] is computed, as shown in Block
`408. Conjugate transpose [I—w'lA]'l of the matrix inverse
`[(I—w'lA)'1]* is computed, as shown in Block 410. The
`matrix product [D S]*[D S] and the matrix product S'1(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=Lr{S’1(S’1)*[(1—W1A)’1]*[D S]*[D S] (Fifi/0’1},
`
`as a function of the wavelet w, as shown in Block 414. A
`value for the source wavelet w is computed by minimizing
`the total energy in the expression for P in w, as shown in
`Block 416. The minimization is accomplished by an iterative
`loop of Blocks 406, 408, 410, 412, 414 and 416. Finally,
`matrix P(s,r) is computed by inserting the final computed
`value for the source wavelet w into the expression given by
`Eq. (17),
`
`60
`
`P=[D S][I—w’1A]’1S’1,
`
`as shown in Block 418.
`An alternative embodiment of the method of the present
`invention for computing the primary matrix P(s,r) uses the
`form of Eq. (22) rather than Eq. (20) for computing the trace
`
`9
`This is the same as Eq. (21) above, except that matrix B is
`replaced by matrix Bt.
`In the case in which the matrix columns are tapered to
`remove edge effects, the approach remains the same, but
`now matrix Ais changed instead of matrix B. Define matrix
`A, by
`
`At=34r T(S’1)*.
`
`Then the energy expression of Eq. (19) becomes
`
`E = tr{PT(PT)*}
`
`{ l {
`
`tr {DSU—w lm’ls lT}{DS(I—w lAflS lT}’}
`
`tr DS(1— w’1/\)’IS’l T(S’1T)*[(I — wrlA)’ll*(DS)*T}
`
`tr s"11'1'(S’1,r[(1 — w"/\)’1]*(L)5)‘DS(I — w’lAJ’l}.
`
`tr{A,[(1— w’lAferU — w’lAfi}.
`
`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,r,f) . Here s
`is source position, r is receiver position, and f is frequency.
`Next, the original seismic traces are also truncated in time,
`Fourier-transformed, Kirchhoff-modified, and inserted into
`another data cube, DM(s,r,f). For each frequency f, matrices
`D(s,r) and DM(s,r) are extracted from the data cubes D(s,r,f)
`and DM(s,r,f), respectively. An initial estimate for the value
`of the source wavelet w component for a particular fre—
`quency is made and inserted into the energy expression Eq.
`(20) along with D(s,r) and DM(s,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 minimization is done by iterating
`through values of