METHOD OF ADAPTIVE FILTERING OF MULTIPLE SEISMIC REFLECTIONS

Abstract

The invention is a method for processing seismic data to eliminate the coherent noise arising from multiple seismic reflections. A decomposition procedure is applied to decompose along N directions of decomposition seismic data into a set of N components. At least one model of multiple reflections is decomposed according to the same decomposition procedure and along the same N directions. For each direction of decomposition, a relative concentration between the component of the model of multiples and the component of the data in the direction concerned is calculated. Next a procedure of adaptive filtering of the multiple reflections is applied to each of the seismic components. The filtered seismic components are recombined with the recombination being weighted by a weighting dependent on the relative concentrations calculated for each direction of decomposition.

Claims

1-15. (canceled)

16. A method for constructing a filtered seismic image of multiple reflections, based on a recording of seismic data comprising primary reflections and multiple reflections, based on at least one model of the multiple reflections, comprising: a) applying a decomposition procedure to the seismic data to decompose the seismic data along N directions of decomposition into a set of N components of the seismic data; b) applying the decomposition procedure to at least one of the models of the multiple reflections to decompose at least one of the models along the N directions of decomposition, into a set of N components of the at least one model; c) calculating for each direction of decomposition and for at least one of the models of the multiple reflections a relative concentration between the component of the seismic data in each direction and the component of the model in each direction; d) based on at least one of the models of the multiple reflections, adaptively filtering the multiple reflections is of each of the N components of the seismic data, to obtain a set of N components of the filtered seismic data; and e) calculating a weighted recombination of the components of the filtered seismic data, based on a weighting calculated for each direction of decomposition as a function of the relative concentrations calculated in each direction of decomposition.

17. The method as claimed in claim 16, wherein the decomposition procedure use a wavelet transform.

18. The method as claimed in claim 16, wherein the decomposition procedure use a dual-tree M-band wavelet transform.

19. The method as claimed in claim 17, wherein the decomposition procedure use a dual-tree M-band wavelet transform.

20. The method as claimed in claim 16, wherein the adaptive filtering comprises a unary Wiener filter using a complex wavelet frame.

21. The method as claimed in claim 17, wherein the adaptive filtering comprises a unary Wiener filter using a complex wavelet frame.

22. The method as claimed in claim 18, wherein the adaptive filtering comprises a unary Wiener filter using a complex wavelet frame.

23. The method as claimed in claim 19, wherein the adaptive filtering comprises a unary Wiener filter using a complex wavelet frame.

24. The method as claimed in claim 16, wherein the relative concentration between a two-dimensional signal S and a two-dimensional signal S′ is calculated according to the following formula: $\mathrm{CR}\left(S\left(s_1,s_2\right),S^{\prime }\left(s_1,s_2\right)\right)=\frac{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S\left(s_1,s_2\right).Math.}^2}}}{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S^{\prime }\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S^{\prime }\left(s_1,s_2\right).Math.}^2}}}$ where s1 is the number of samples in one direction in space and s2 is the number of samples in the other direction in space.

25. The method as claimed in claim 17, wherein the relative concentration between a two-dimensional signal S and a two-dimensional signal S′ is calculated according to the following formula: $\mathrm{CR}\left(S\left(s_1,s_2\right),S^{\prime }\left(s_1,s_2\right)\right)=\frac{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S\left(s_1,s_2\right).Math.}^2}}}{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S^{\prime }\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S^{\prime }\left(s_1,s_2\right).Math.}^2}}}$ where s1 is the number of samples in one direction in space and s2 is the number of samples in the other direction in space.

26. The method as claimed in claim 18, wherein the relative concentration between a two-dimensional signal S and a two-dimensional signal S′ is calculated according to the following formula: $\mathrm{CR}\left(S\left(s_1,s_2\right),S^{\prime }\left(s_1,s_2\right)\right)=\frac{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S\left(s_1,s_2\right).Math.}^2}}}{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S^{\prime }\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S^{\prime }\left(s_1,s_2\right).Math.}^2}}}$ where s1 is the number of samples in one direction in space and s2 is the number of samples in the other direction in space.

27. The method as claimed in claim 20, wherein the relative concentration between a two-dimensional signal S and a two-dimensional signal S′ is calculated according to the following formula: $\mathrm{CR}\left(S\left(s_1,s_2\right),S^{\prime }\left(s_1,s_2\right)\right)=\frac{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S\left(s_1,s_2\right).Math.}^2}}}{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S^{\prime }\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S^{\prime }\left(s_1,s_2\right).Math.}^2}}}$ where s1 is the number of samples in one direction in space and s2 is the number of samples in the other direction in space.

28. The method as claimed in claim 16, wherein the weighted recombination is calculated as follows: $\mathrm{SR}=\alpha .Math.\underset{n=1,N}{.Math.}.Math.{\mathrm{.Math.}}_n.\widehat{D^{\prime }}\left(0,\mathrm{.Math.}.Math.,{\mathrm{SDF}}_n,\mathrm{.Math.}.Math.,0\right)$ where SDF.sub.n is the component in the direction n (n=1, N) of the filtered seismic data, α is a constant, ε.sub.n is the weighting for the direction of decomposition n, and {circumflex over (D)}′ is an inverse of the decomposition procedure.

29. The method as claimed in claim 17, wherein the weighted recombination is calculated as follows: $\mathrm{SR}=\alpha .Math.\underset{n=1,N}{.Math.}.Math.{\mathrm{.Math.}}_n.\widehat{D^{\prime }}\left(0,\mathrm{.Math.}.Math.,{\mathrm{SDF}}_n,\mathrm{.Math.}.Math.,0\right)$ where SDF.sub.n is the component in the direction n (n=1, N) of the filtered seismic data, α is a constant, ε.sub.n is the weighting for the direction of decomposition n, and {circumflex over (D)}′ is an inverse of the decomposition procedure.

30. The method as claimed in claim 18, wherein the weighted recombination is calculated as follows: $\mathrm{SR}=\alpha .Math.\underset{n=1,N}{.Math.}.Math.{\mathrm{.Math.}}_n.\widehat{D^{\prime }}\left(0,\mathrm{.Math.}.Math.,{\mathrm{SDF}}_n,\mathrm{.Math.}.Math.,0\right)$ where SDF.sub.n is the component in the direction n (n=1, N) of the filtered seismic data, α is a constant, ε.sub.n is the weighting for the direction of decomposition n, and {circumflex over (D)}′ is an inverse of the decomposition procedure.

31. The method as claimed in claim 24, wherein the weighted recombination is calculated as follows: $\mathrm{SR}=\alpha .Math.\underset{n=1,N}{.Math.}.Math.{\mathrm{.Math.}}_n.\widehat{D^{\prime }}\left(0,\mathrm{.Math.}.Math.,{\mathrm{SDF}}_n,\mathrm{.Math.}.Math.,0\right)$ where SDF.sub.n is the component in the direction n (n=1, N) of the filtered seismic data, α is a constant, ε.sub.n is the weighting for the direction of decomposition n, and {circumflex over (D)}′ is an inverse of the decomposition procedure.

32. The method as claimed in claim 16, wherein the weighted recombination is calculated as follows: $\mathrm{SR}=\alpha .Math.\underset{n=1,N}{.Math.}.Math.\widehat{D^{\prime }}\left(0,\mathrm{.Math.}.Math.,{\mathrm{.Math.}}_n,{\mathrm{SDF}}_n,\mathrm{.Math.}.Math.,0\right)$ where SDF.sub.n, is the component in the direction n (n=1, N) of the filtered seismic data, α is a constant, ε.sub.n is the weighting for the direction of decomposition n, and {circumflex over (D)}′ is an inverse of the decomposition procedure.

33. The method as claimed in claim 28, wherein a residual corresponding to a difference between the seismic data and a result of an inverse of a decomposition procedure is added to the components of the seismic data.

34. The method as claimed in claim 28, wherein the weighting for the direction of decomposition n (n=1, N) is calculated as follows: ${\mathrm{.Math.}}_n=\frac{\underset{i=1,I}{.Math.}.Math..Math.{\mathrm{.Math.}}_n^i}{I}$ where ε.sub.n.sup.i is a weighting for the direction n and a model of multiple reflections i (i=1, I).

35. The method as claimed in claim 28, wherein the weighting for the direction of decomposition n (n=1, N) is calculated as follows: ${\mathrm{.Math.}}_n=\underset{i=1,I}{\max }.Math.\left({\mathrm{.Math.}}_n^i\right)$ where ε.sub.n.sup.i is a weighting for the direction n and a model of multiple reflections i (i=1, I).

36. The method as claimed in claim 34, wherein the weighting for the direction of decomposition n and the model of multiple reflections i (with i=1, I and n=1, N) is calculated as follows: ${\mathrm{.Math.}}_n^i=W\left(x\right),\mathrm{with}.Math..Math.x=\min \left({\mathrm{CR}}_n^i,\frac{1}{{\mathrm{CR}}_n^i}\right),$ where W is an increasing function, CR.sub.n.sub.i is the relative concentration between the component of the seismic data in the direction n and the component of the model i in the direction n.

37. The method as claimed in claim 36, wherein W is defined such that for a given x, W(x)=x.sup.p with p≧0.

38. The method as claimed in claim 34, wherein p equals 4.

39. A method for exploiting an underground formation using a method for constructing a filtered image of multiple reflections based on recording of seismic data comprising primary reflections and multiple reflections based on at least one model of the multiple reflections including a) applying a decomposition procedure to the seismic data is applied to decompose, along N directions of decomposition, the seismic data into a set of N components of the seismic data; b) applying the decomposition procedure to at least one of the models of the multiple reflections, to decompose, along the N directions of decomposition, at least one of the models of multiple reflections into a set of N components of the at least one model; c) calculating for each direction of decomposition and for at least one of the models of the multiple reflections, a relative concentration between the component of the seismic data in the direction and the component of the model in the direction; d) based on at least one of the models of the multiple reflections, adaptively filtering the multiple reflections is applied to each of the N components of the seismic data, to obtain a set of N components of the filtered seismic data; and e) calculating a weighted recombination of the components of the filtered seismic data, based on the basis of a weighting calculated for each direction of decomposition as a function of the relative concentrations calculated in the direction of decomposition comprising the steps: constructing a geological model representative of the formation being studied based on at least a determined seismic image; determining an optimal exploitation scheme for the reservoir is determined based on the determined geological model; and exploiting the reservoir by implementing the optimal exploitation scheme.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Other characteristics and advantages of the method according to the invention will become apparent on reading the description hereinafter of nonlimiting exemplary embodiments, while referring to the appended figures described hereinafter.

[0046] FIG. 1 illustrates a device for acquiring seismic data as well as examples of a trajectory of primary reflections and of multiple reflections generated by this device.

[0047] FIGS. 2A and 2B show an example of seismic data before summation and the corresponding model of multiple reflections.

[0048] FIG. 3 illustrates a set of different wavelets with each wavelet representing by a particular range of frequencies and a particular range of orientations in space.

[0049] FIGS. 4A, 4B, 4C, 4D illustrate the decomposition coefficients for the data presented in FIG. 2A, using 4 different oriented wavelets, which makes possible obtaining preferential directions on four components SD.sub.1, SD.sub.2, SD.sub.3, and SD.sub.4.

[0050] FIGS. 4E, 4F, 4G, 4H illustrate the decomposition coefficients for the model of multiples presented in FIG. 2B, using 4 different oriented wavelets, which makes possible obtaining preferential directions on four components MD.sub.1, MD.sub.2, MD.sub.3 and MD.sub.4.

[0051] FIG. 5A shows the seismic data already presented in FIG. 2A, before random-noise filtering; FIG. 5B shows the result obtained after applying the procedure described in Ventosa et al. (2012), and FIG. 5C shows the result obtained after applying the method according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

[0052] The following definitions are used in the course of the description of the invention: [0053] model of multiples: This model is an approximate model of the multiple reflections contained in seismic data. Several methods for obtaining these models of multiple reflections exist. One method obtains an approximate version of the primary reflection (for example by filtering), and then convolves it with itself or else with the initial seismic trace (see for example Pica et al. (2005)). A model of multiples can also be obtained by solving the wave equation in the medium being considered. In general, models of multiples are satisfactory approximations of multiple reflections. They can however be shifted on the vertical axis (time or depth axis), have amplitudes and/or a frequency spectrum that differ with respect to the actual multiple reflections. Because of this lack of precision, the attenuation processing of the multiple reflections can apply to an adaptation of the models of multiples to the seismic data, by adaptive filtering, also called registration. It should be noted that a model of multiples may have limited relevance, being for example representative of a multiple reflection for a limited range of offsets, or else for a given order of multiples. It is then possible to apply several models of multiples to simulate, completely, a multiple reflection recorded in seismic data. [0054] wavelet decomposition procedure: This procedure makes it possible to decompose seismic data into various components represented by wavelets. A formula for obtaining such wavelets is given hereinafter:

[00007]$g\left(x,y\right)=\exp \left(2.Math..Math.j.Math..Math.\pi .Math.\left(u_0.Math.x+v_0.Math.y\right)+\phi \right)\times \exp \left(-\left(\frac{{\left(x-x_0\right)}^2}{{\sigma }_x}+\frac{{\left(y-y_0\right)}^2}{{\sigma }_y}\right)\right)$

where the coordinates x and y correspond to the position of the sample considered in the seismic data, the pair (u0,v0) defines the coordinates of a vector in a given direction, the norm of this vector defining the frequency, and the terms σ.sub.x and σ.sub.x correspond to the envelope widths in the directions x and y. Each choice of parameters provides a particular directional wavelet. FIG. 3 presents various directional wavelets, each being characterized by an envelope, a main frequency and a particular orientation in space. [0055] energy of a signal: In the case of a two-dimensional signal S characterized by s1 samples in one direction and s2 samples in the other direction, the energy of a signal is defined by the following formula:

[00008]$\underset{s_1,s_2}{.Math.}.Math.{.Math.S\left(s_1,s_2\right).Math.}^2.$

[0056] The subject of the present invention is a method for constructing a filtered seismic image of the multiple reflections present in seismic recordings, on the basis of one or more models of multiples. This invention uses a decomposition both of the seismic data and of the model of multiples according to preferential directions of decomposition, and a weighted recombining of each of the filtered seismic components according to the corresponding components of the model of multiples. The seismic data may have been recorded by two-dimensional or three-dimensional acquisition devices and be organized according to any type of collection (that is as a common-shot collection, common-receiver collection, etc). The seismic data may correspond equally to pre-summation or post-summation collections, pre-migration or post-migration collections. The present invention may be applied at any stage of the seismic processing, and preferentially before the step of interpreting the seismic image resulting from the application of the seismic processing as a whole. The invention requires that at least one model of multiple reflections be available.

[0057] The present invention comprises at least the following steps:

[0058] a) applying a decomposition procedure to decompose, along N directions of decomposition, the seismic data into a set of N components of the seismic data;

[0059] b) applying for at least one of the models of the multiple reflections, a decomposition procedure to decompose, along the N directions of decomposition, the model of multiple reflections into a set of N components of the model;

[0060] c) calculating for each direction of decomposition and for at least one of the models of the multiple reflections, a relative concentration between the component of the seismic data in the direction and the component of the model in the direction;

[0061] d) based on at least one of the models of the multiple reflections, applying adaptive filtering of the multiple reflections to each of the N components of the seismic data, to obtain a set of N components of the filtered seismic data;

[0062] e) calculating a weighted recombination of the components of the filtered seismic data, based on a weighting calculated for each direction of decomposition as a function of the relative concentrations calculated in the direction of decomposition.

[0063] The main steps of the present invention are described hereinafter in the case of two-dimensional post-summation seismic data S, but the method can equally well be applied to three-dimensional seismic data. The main steps of the present invention are specifically described hereinafter in the case of a single model of multiples M. The case of several models of multiples is particularized in a variant of the present invention.

[0064] a) Decomposition of the Seismic Data into N Components

[0065] This step entails decomposing the seismic data S into N components according to a decomposition procedure {circumflex over (D)}. This step better distinguishes certain characteristic features of the seismic data. A series of components SD.sub.n with n=1, N associated with N different directions of decomposition are thus obtained.

[0066] According to one embodiment of the present invention, a decomposition procedure {circumflex over (D)} is chosen which makes it possible to decompose seismic data according to various frequency bands and various orientation ranges in space. In this case, a direction of decomposition is a pair made up of a frequency band and an orientation range in space.

[0067] According to one embodiment of the present invention, a dual-tree M-band wavelet decomposition procedure, such as described in Chaux et al. (2006), is used. This entails particular directional wavelets, obtained by choosing a set of filters, termed a primal bank of filters with each filter being defined for a specific frequency band, and with the set of filters making it possible to cover a predefined frequency range. Another set of filters, which are deduced from the previous ones and termed a dual bank of filters, is calculated thereafter. This set is obtained by shifting each filter by a half-coefficient by conventional interpolation techniques or by calculation in the Fourier domain. These filters are thereafter applied separately to the horizontal and vertical directions, and combined to provide diagonal directions. In practice, the P frequency bands are defined and the dual bank of filters is constructed as described previously so as to obtain the 4P decompositions according to the various ranges of orientations in space and frequency bands. According to a preferred embodiment of the present, P lies between four and eight. It should be noted that for wavelets of this type, the envelopes are defined by the choice of the frequency ranges and of the orientations in space.

[0068] According to one embodiment of the invention, a decomposition procedure is used for which at least one inverse {circumflex over (D)}′ is known exactly, that is to say such that {circumflex over (D)}′={circumflex over (D)}.sup.−1. An inverse of the decomposition procedure, also termed a recombining procedure, makes it possible to recombine the various components SD.sub.n, with n=1, N, so as to produce the initial seismic data S.

[0069] According to one embodiment of the invention, a decomposition procedure is used for which an inverse {circumflex over (D)}′ is known in an approximate manner, which is also termed a pseudo-inverse {circumflex over (D)}′ or else an approximate recombining procedure. In this case, in order to produce the initial seismic data S, the residual corresponding to the difference between the recorded seismic data S and the result of the chosen pseudo-inverse {circumflex over (D)}′ (SD.sub.n) is added beforehand to the various components SD.sub.n, with n=1, N.

[0070] The result of the decomposition procedure described in Chaux et al. (2006), is applied to the seismic data presented in FIG. 2A, which is shown in FIGS. 4A to 4D. Thus, the seismic data S have been decomposed using 4 different directional wavelets, making it possible to obtain four components SD.sub.1, SD.sub.2, SD.sub.3, and SD.sub.4 oriented along 4 different directions in space, which more precisely are a low-frequency isotropic component (FIG. 4A), and three directional components, approximately covering the following angular spans: −23° to +23° (FIG. 4B), 22° to 68° (FIG. 4C), 67° to 113° (FIG. 4D).

[0071] b) Decomposition of the Model of Multiples into N Components

[0072] This step decomposes the model of multiples M into N components, with the same decomposition procedure {circumflex over (D)} and along the same N directions as those defined in step a) of decomposing the seismic data.

[0073] A series of components MD.sub.n with n=1, N is thus obtained.

[0074] The result of the decomposition procedure described in Chaux et al. (2006), applied, along the same directions as those used to produce FIGS. 4A to 4D, to the model of multiple reflections presented in FIG. 2B, is shown in FIGS. 4E to 4H. Thus, the model of multiple reflections M has been decomposed into four components MD.sub.1, MD.sub.2, MD.sub.3, and MD.sub.4, oriented along the same directions as those defined to undertake the decomposition of the seismic data.

[0075] c) Calculation of the Relative Concentration of the N Components

[0076] This step calculates, for each direction of decomposition n (with n=1, N) defined during the previous steps, the relative concentration between the component SD.sub.n arising from the seismic data in the direction being considered and the component MD.sub.n arising from the model of multiples in the same direction.

[0077] The relative concentration between two signals corresponds to the ratio of the concentration of each of the signals.

[0078] According to one embodiment of the invention, the concentration C of a two-dimensional signal S which is defined by s1 samples in one direction and s2 samples in the other direction, is defined in the following manner:

[00009]$C\left(S\left(s_1,s_2\right)\right)=\frac{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S\left(s_1,s_2\right).Math.}^2}}}{\sqrt{s_1.Math.s_2}-1}$

[0079] The defined concentration therefore lies between 0 and 1. More precisely, a very concentrated signal, for example such that a single value of this signal is non-zero, has a concentration C equal to 1. On the other hand, a very unconcentrated signal, for example for which all the values of the signal are equal to one another, will have a concentration C equal to 0. The concentration makes it possible in particular to rate the quality of a processing aimed at concentrating the energy of a given signal in an optimal manner. For example, a sinusoidal signal is converted by Fourier transformation into a well-localized frequency spike. Its concentration index is then maximal.

[0080] According to one embodiment of the invention, the relative concentration CR between a signal S and a signal S′ which are both two-dimensional and characterized by s1 samples in one direction and s2 samples in the other direction may be then written as:

[00010]$\mathrm{CR}\left(S\left(s_1,s_2\right),S^{\prime }\left(s_1,s_2\right)\right)=\frac{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S\left(s_1,s_2\right).Math.}^2}}}{\sqrt{s_1.Math.s_2}-\frac{\underset{s_1,s_2}{.Math.}.Math..Math.S^{\prime }\left(s_1,s_2\right).Math.}{\sqrt{\underset{s_1,s_2}{.Math.}.Math.{.Math.S^{\prime }\left(s_1,s_2\right).Math.}^2}}}$

[0081] The relative concentration CR.sub.n for a given direction n (with n=1, N) between a component SD.sub.n arising from the seismic data in this direction n and the component MD.sub.n arising from the model of multiples in the same direction indicates the capacity of a decomposition procedure to concentrate in an equivalent manner the seismic data and the model of multiples for this direction n. Thus, if the relative concentration CR.sub.n for a given direction n is greater than 1, the seismic signal is more concentrated than the model of multiples for the direction being considered. And conversely, if the relative concentration is less than 1, then the model of multiples is more concentrated than the seismic signal for this same direction.

[0082] d) Adaptive Filtering of the N Seismic Components

[0083] This step uses a procedure {circumflex over (F)} of adaptive filtering of the primary reflections and the multiple reflections, on the basis of a model of multiples. More precisely, in the course of this step, an adaptive filtering procedure {circumflex over (F)} is applied to each of the components SD.sub.n arising from the seismic data and calculated in step a). Thus, N filtered seismic components SDF.sub.n of the multiple reflections are obtained in the following manner:

SDF.sub.n={circumflex over (F)}(SD.sub.n,MD.sub.n) with n=1,N.

[0084] Thus, the adaptive filtering is not applied to the seismic data themselves but instead to the seismic components arising from the decomposition of the seismic data according to favored frequencies and orientations in space. In this manner, the seismic components are simpler and/or more homogeneous than the original seismic data, the parameters of the adaptive filtering are simpler to estimate, and as a consequence, the adaptive filtering is more efficient. More precisely, the goal of applying an adaptive filter to each component arising from the seismic data is to better eliminate the multiple reflections and to better safeguard the frequency content, component by component.

[0085] According to one embodiment of the invention, the adaptive filtering procedure {circumflex over (F)} is chosen for its capacity to safeguard the frequency content, and especially the low and high frequencies of the primary reflections.

[0086] According to a preferred embodiment of the invention, a procedure for adaptive filtering of multiple reflections by unary Wiener filters using a single-dimensional complex wavelet frame such as described in Ventosa et al. (2012) is used. This procedure makes it possible in fact to best safeguard the frequency content, especially the low and the high frequencies of the seismic data.

[0087] e) Weighted Recombining of the N Filtered Seismic Components

[0088] The objective of this step is to recombine the N filtered seismic components SDF.sub.n of the multiple reflections obtained during step d), the recombination being weighted by taking into account the relative concentrations CR.sub.n (with n=1, N) calculated in step c). Thus, this step recombine the various filtered seismic components, while account for of the efficiency of the decomposition procedure applied to the seismic data and to the model of multiples.

[0089] According to one embodiment of the present invention, to carry out the recombination of the various filtered seismic components, recourse is had to an inverse {circumflex over (D)}′ of the decomposition procedure {circumflex over (D)} such as defined in step a) is used.

[0090] According to one embodiment of the present invention, to carry out the recombination of the various filtered seismic components, recourse is had to a pseudo-inverse {circumflex over (D)}′, or approximate inverse, of the procedure {circumflex over (D)} such as defined in step a).

[0091] According to one embodiment of the invention, the concentrated components are favored in equivalent ways, by using a maximum weighting when the relative concentration in the direction considered is close to 1, and a minimum weighting in the converse case.

[0092] According to a preferred embodiment of the invention, an increasing function W is considered and a weighting is defined for a given direction n by the value of this function at a point x defined by the minimum between the relative concentration calculated for this direction and its inverse, that is ε.sub.n=W(x) with

[00011]$x=\min .Math..Math.\left({\mathrm{CR}}_n,\frac{1}{{\mathrm{CR}}_n}\right).$

In this case, the point of evaluation x of the function W necessarily lies between 0 and 1. Moreover, with W being chosen to increase, more significant weight is thus given to the values of x close to 1, that is to say to the equitably concentrated components.

[0093] According to one embodiment of the invention, a function is chosen such that for a given value x, W(x)=x.sup.p with p≧0.

[0094] According to a preferred embodiment of the invention, a function is chosen such that for a given value x, W(x)=x.sup.p with p=4.

[0095] According to one embodiment of the invention, an estimation SR of the filtered seismic data of the multiple reflections, resulting from the best combinations of the N seismic components SDF.sub.n filtered by the procedure {circumflex over (F)} and complying with the relative concentration of the various components arising from the decomposition procedure {circumflex over (D)}, is calculated according to the following formula:

[00012]$x=\min \left({\mathrm{CR}}_n,\frac{1}{{\mathrm{CR}}_n}\right).$

where α is a constant. According to one embodiment of the invention, α is chosen so that if zero values are assigned to the model of multiples M, SR is characterized by the same energy as the seismic data S at input. According to a particular embodiment of the invention, α is constant and equal to 1. In this way, the filtered components SDF.sub.n of the multiple reflections are recombined according to an inverse of the procedure {circumflex over (D)}, by zeroing all except one of them, and by assigning the calculated weighting to the recombination of the component in question.

[0096] According to another embodiment, an estimation SR of the filtered seismic data of the multiple reflections is calculated in the following manner:

[00013]$\mathrm{SR}=\alpha .Math.\underset{n=1,N}{.Math.}.Math.\widehat{D^{\prime }}\left(0,\mathrm{.Math.}.Math.,{\mathrm{.Math.}}_n,{\mathrm{SDF}}_n,\mathrm{.Math.}.Math.,0\right)$

where α is a constant. According to one embodiment of the invention, α is chosen in such a way that if zero values are assigned to the model of multiples M, SR is characterized by the same energy as the seismic data S at input. According to a particular embodiment of the invention, α is constant and equal to 1.

[0097] In this way, a filtered seismic image of the multiple reflections is obtained, by virtue of an adaptive filtering applied to seismic data decomposed along preferential directions, followed by a recombining of the filtered components taking account of the efficiency of the decomposition procedure.

Variants

[0098] Case of Several Models of Multiples

[0099] In the case where several models of multiples M′ with i=1, I are available, the following are repeated for each model: [0100] step b), described hereinbelow, of decomposing the model of multiples: a series of decompositions MD.sub.n with i=1, I and n=1, N is thus obtained; [0101] step c) of calculating the relative concentration: the relative concentrations CR.sub.n.sup.i are thus calculated for each direction n and for each model of multiples i, with i=1, I and n=1, N.

[0102] Next, the filtering procedure {circumflex over (F)} is applied, for component n (with n=1, N), while accounting for the various models of multiples. More precisely, the N filtered seismic components are calculated in the following manner:

SDF.sub.n={circumflex over (F)}(SD.sub.n,MD.sub.n.sup.1, . . . ,MD.sub.n.sup.i, . . . ,MD.sub.n.sup.I) with n=1,N.

[0103] Next, the weighting calculation described in step e) of weighted recombining of the filtered seismic components is repeated for each model of multiples. A series of weightings En is thus obtained for all the values of i and of n, with i=1, I and n=1, N.

[0104] Next, for each direction n (with n=1, N), a composite weighting ε.sub.n is calculated f based on the set of weightings ε.sub.n.sup.i (with i=1, I) calculated for each model of multiples.

[0105] According to one embodiment of the present invention, a composite weighting ε.sub.n is determined by calculating the average of the individual weightings ε.sub.n.sup.i i.e.

[00014]${\mathrm{.Math.}}_n=\frac{\underset{i}{.Math.}.Math..Math.{\mathrm{.Math.}}_n^i}{I}.$

[0106] According to a preferred embodiment of the present invention, in order to favor, for a given direction n (with n=1, N), the relative concentration between the seismic data S and at least one of the models of multiples M′, the maximum value among the series of weightings ε.sub.n.sup.i (with i=1, I) calculated for each model of multiples is allotted to the composite weighting ε.sub.n.

[0107] Next, after the calculation of the composite weighting en, an estimation SR of the filtered seismic data of the multiple reflections is calculated by weighted recombination of the filtered seismic components SDF.sub.n, as described previously in step e).

[0108] Exploitation of an Underground Formation

[0109] Furthermore, the invention relates to a method for exploiting an underground formation, in which the following steps are carried out: [0110] a filtered seismic image of the multiple reflections is constructed by use of the method as described above; [0111] a geological model representative of the formation being studied is constructed on a basis of at least the determining seismic image; [0112] an optimal exploitation scheme for the reservoir is determined based on the determined geological model and [0113] the reservoir is exploited by implementing the optimal exploitation scheme.

[0114] With the help of a geological model constructed based on at least on the information arising from a seismic image obtained during the previous steps, several exploitation schemes can be determined which correspond to various possible configurations of exploitation of the underground reservoir: siting of the producer and/or injector wells, target values for the flowrates per well and/or for the reservoir, the type of tools used, the fluids used, injected and/or recovered, etc. For each of these schemes, it is appropriate to determine their production forecasts. These probabilistic production forecasts can be obtained by use of flow simulation software as well as by use of the fitted numerical model of the reservoir. Reservoir simulation is a technique making it possible to simulate the flows of fluids within a reservoir by means of software called a flow simulator. For example, the PumaFlow® (IFP Énergies nouvelles, France) software is a flow simulator.

[0115] One or more possible exploitation schemes suitable for the geological model studied is or are defined. For each of these schemes, the responses are determined by simulation.

[0116] On the basis of the probabilistic production forecasts defined for each exploitation scheme, the exploitation scheme can be chosen by comparison which are to them the most relevant. For example: [0117] by comparing the maximum of the oil volume recovered, it is possible to determine the production scheme liable to provide the maximum recovery or to be the most profitable; and [0118] by comparing the standard deviation of the oil volume recovered, it is possible to determine the least risky production scheme.

[0119] The reservoir is then exploited according to the exploitation scheme defined for example by drilling new (producer or injector) wells, by modifying the tools used, by modifying the flowrates and/or the nature of fluids injected, etc.

[0120] Computer Program Product

[0121] The invention relates, moreover, to a computer program product storing a non transient recording downloadable from a communication network and/or recorded on a medium readable by computer and/or executable by a processor. This program comprises program code instructions for the implementation of the method such as described hereinabove, when the program is executed on a computer.

Examples of Application

[0122] The method according to the invention is applied to a case of two-dimensional real seismic data after summation in the case of the presence of random noise. It is observed in FIG. 5A that the seismic data in question correspond to a superposition of primary seismic reflections and of multiple reflections and the multiple reflections have a very high amplitude relative to the primary reflections. It may also be observed that primary reflections and multiple reflections are strongly disturbed by significant random noise. FIG. 5B presents the result of the application of an adaptive filtering according to the prior art (described in Ventosa et al. (2012)) to the data of FIG. 5A. It is observed that the adaptive filtering according to the prior art correctly reveals the primary reflection of interest in a large part of the zone considered, but that the continuity of this primary reflection is lost in the leftmost part of FIG. 5B. However, this is the part where the reflector that generated the reflection of interest slopes the most. FIG. 5C presents the result of the method according to the invention applied to the seismic data presented in FIG. 5A. The method according to the invention has been implemented by using the technique described in Chaux et al. (2006) for the decomposition procedure and the technique described in Ventosa et al. (2012) for the adaptive filtering. It is very clearly observed that the method according to the invention makes it possible to extract in a sharper manner the primary reflection of interest, in particular in the left part of FIG. 5C. This advantage is in particular obtained through the fact that the method according to the invention applies an adaptive filtering to seismic data decomposed according to favored orientations in space, and that the weighting by the relative concentration strengthens the resemblances between the seismic data and the model of multiples, to the detriment of the undesired noise or disturbances, producing a more precise seismic image of the primary reflections.

[0123] Thus, even in the case of seismic data exhibiting strong random noise, the present invention makes it possible to improve the adaptive filtering of multiple reflections contained in seismic data, by operating in a selective manner according to favored frequency bands and ranges of orientations in space.