Wave separation procedure for borehole seismics with walk away acquisitions

The present invention relates to a method for the separation of elastic waves received and recorded after propagation in a medium and, more particularly, the processing of the signals recorded during an acquisition made by emission of elastic waves in several positions of one or more emitters arranged on the surface of the medium, the receptor (s) being fixed (s) in a well drilled in said medium during the entire series of emission, the receiver (s) being then moved (s) to another position in the well for another series of emissions. Such an acquisition is usually referred to as "walkaway", that is to say that the transmission is carried out with increasing offset from one transmission to the next.

In terrestrial and / or marine well seismic, several emitters separated by a constant distance or not and disposed on the surface of the medium to be explored along an acquisition line are used. A receiving tool is lowered into the well whose direction can be vertical or deviated. The tool generally comprises several receivers separated by a constant distance or not which is generally different from the constant distance or not separating two transmitters or consecutive sources. Each receiver comprises three sensors or geophones which are oriented in three different directions generally orthogonal, one of the directions being determined to have a reference axis. The reference axis may be non-vertical or vertical if these sensors are mounted on a gimbal system. In both cases, the reference axis V of the sensors is perpendicular to the plane defined by the two other directions H₁ and H₂ whose positions are not generally known. The axes V, H₁, H₂ constitute a first axis system.

The signals representative of the waves propagating in the medium are recorded as a function of time in the form of traces. The waves that propagate in the medium are mainly rising or falling P waves (direct waves) and S waves.

The traces recorded along the V, H axes₁', H₂ can be grouped into common emitter trace collections, for example.

A current processing of the traces consists in defining a second system of axes R, N, T, in which the axis R is located in the maximum energy plane. Such a transformation is carried out using what is usually called a double rotation, an angle rotation Φ_H being made around the V axis and the second angle rotation Φ_V being carried out around a direction called H_min which is perpendicular to an H vector_max which is the projection of the direct downward wave in the plane defined by H₁ and H₂, the angle Φ_H being the angle that make between them H₁ and H_max.

The angle Φ_V is the angle formed by the axes V and R. This transformation of the first system of axes V, H₁, H₂ in a second axis system R, N, T is carried out in particular using a software well known to specialists and called SEISLINK, marketed by the company WESTERN.

Traces can be produced in the R, N, T axis system by projecting the trace components marked along the V, H axes.₁, H₂ and constitute collections of common emitter traces identified according to the depth z and the time t, starting from at least the traces along the axis R (or traces R).

The collections of traces, for example collections V, R or N, are then processed to separate the different types of waves, especially the descending and rising P waves and the rising and falling S waves.

A technique for separating said waves can be implemented using the SEPAR software of the Compagnie Générale de Géophysique (CGG). However, this requires knowing the apparent speeds of the types of waves sought. Generally, these apparent speeds are determined using a manual dash on the recorded tracks. The drawbacks of this technique lie in the inaccuracy inherent in manual pointing and sometimes in the impossibility of discerning the types of waves directly on the recorded traces. In addition, the cost of treatment is high for inaccurate results.

Another technique advocated by the company SCHLUMBERGER consists of using software that is capable of analyzing the collections of traces with a common emitter and of deducing the apparent speeds of the waves to be separated when they exist. A disadvantage of this technique lies in its instability and especially in the fact that the separation of waves is inefficient because of the small number of sensors used.

These wave separation techniques are described respectively in the following articles:<ul id="ul0001" list-style="none" compact="compact"><li>B. SEEMAN and L. HOROWICZ: "Vertical seismic profiling: Separation of upgoing and downgoing acoustic waves in a stratified medium." GEOPHYSICS, Vol 48, No. 5 (May 1983), pp. 555-568;</li><li>C. ESMERSOY: "Inversion of P and SV waves from multicomponent vertical offset seismic profiles", GEOPHYSICS, Vol. 55, No. 1, (January 1990), p. 39-50.</li></ul>

The subject of the present invention is a novel method for separating elastic waves, in which a particular treatment of trace collections is used, which makes it possible to remedy the drawbacks of the prior methods as briefly described above.

The method according to the invention for separating waves propagating in a medium is of the type consisting of:<ul id="ul0002" list-style="dash" compact="compact"><li>generating waves in the medium successively from several positions of at least one transmitter and located on the surface of said medium,</li><li>receiving said waves on at least one receiver capable of successively occupying at least two positions in a well drilled in the medium, each receiver comprising three sensors oriented in a first axis system (V, H₁, H₂), one of said axes (V) being oriented in a determined direction, each sensor delivering a signal which is recorded as a function of time in the form of a trace which is marked by the positions of the transceiver pair having produced said trace,</li><li>constituting collections C of common emitter traces (z, t) from at least one set of traces V, H₁, H₂,</li></ul> and it is characterized in that it consists, moreover, in<ul id="ul0003" list-style="none" compact="compact"><li>a) extracting the first signal arrivals C of at least one of the collections of traces C,</li><li>b) producing from said collection of traces C a collection of residual traces C in which each residual trace is the difference between the trace C and the estimate C signal,</li><li>c) selecting in the collection of residual traces C a pivot trace located at a depth z,</li><li>d) passing a line of variable slope p at each sampled point of the pivot trace,</li><li>e) summing along said line, for each value of p, the amplitudes of the seismic events situated at the intersections of the traces with said straight line,</li><li>f) postponing amplitude sums for the corresponding times and slope in a domain (τ, p) in which p represents the values ​​of said slope varying between two extreme limits and τ time,</li><li>g) producing a synthetic seismogram from the domain data (τ, p),</li><li>h) comparing the synthetic seismogram with the common emitter collection of residual traces C,</li><li>i) iterating until the synthetic seismogram is acceptable relative to the initial data of said common emitter collection of residual traces C.</li></ul>

According to an advantageous characteristic, the method consists, moreover, in<ul id="ul0004" list-style="dash" compact="compact"><li>define for each pair of transceiver positions a second system of axes (R, N, T) which is calculated, by known means, from the components along the three axes V, H₁, H₂, of the displacement vector corresponding to the first arrival of the wave of interest and the projection H_max of said vector moving in the plane H₁, H₂, the axis R being located in the plane V, H_max, the axis N being perpendicular to the axis R and located in the plane V, H_max and the axis T being perpendicular to the plane R, N,</li><li>produce traces in the system of axes R, N, T by projecting the components of the traces marked along the axes V, H₁, H₂,</li></ul> and to constitute the collections of traces (C) with common emitter (z, t) starting from the traces R,<ul id="ul0005" list-style="dash" compact="compact"><li>extract the first arrivals R signal from the trace collection R,</li><li>realize from the collection of traces R a collection of residual traces R in which each residual trace is the difference between the trace R and the estimate R signal,</li><li>select in the collection of residual traces R a pivot trace located at a depth z,</li><li>perform steps d) to i) on the residual traces R.</li></ul>

According to another characteristic, the method furthermore consists in:<ul id="ul0006" list-style="dash" compact="compact"><li>define a third system of axes (X, Y, Z), fixed geographically,</li><li>determining in the axis system (X, Y, Z) the deviation (Ψ) and the azimuth (θ) of the well,</li><li>align the V axis of the first axis system (V, H₁, H₂) with the direction of the well,</li><li>define a vector h₁ perpendicular to the axes V and Z and contained in the plane H₁, H₂,</li><li>determine the angle δ between the vector h₁ and the H axis₁.</li></ul>

According to another characteristic, the residual components C or R, as well as the components H₁ and H₂ or N and T are projected into the fixed system (X, Y, Z) to form residual components X_r and Y_r, Z_r, on the X, Y, Z axes, the decomposition (τ, p) being made from the components of the X traces_r and Y_r, Z_r.

According to another characteristic, it calculates, by known means, for each transmitter-receiver pair and in the domains (τ, p), an angular polarization attribute ζ and an angular azimuth attribute γ.

When the wave of interest is the P wave, the polarization ζ, azimuth γ, deflection angle ψ and azimuth θ of the well are used as well as the slope p to calculate the speed V attribute._P P waves in the middle.

From the polar polarization and V speed attributes_P, we calculate the apparent slowness attribute p_z with respect to the Z axis of a seismic event corresponding to values ​​τ and p in the domain (τ, p).

P attributes_z and V_P can be used to make a selection (mute) of the events to be kept in the domain (τ, p), said events being in particular those related to the rising P-waves.

Other advantages and characteristics will emerge more clearly on reading the description of the method according to the invention, as well as the appended drawings in which:<ul id="ul0007" list-style="dash" compact="compact"><li>FIG. 1 is a schematic representation of the transmission and reception device in a medium,</li><li>FIG. 2 is a schematic representation of the reception directions V, H₁, H₂ sensors of a receiver of the receiving device,</li><li>FIG. 3 is a schematic representation of the location of the well of FIG. 1 in a geographically fixed X, Y, Z axis system;</li><li>Figures 4 and 5 are representations of angles Φ_H and Φ_V rotations applied to the system V, H₁, H₂ for a transformation in a system of axes R, N, T,</li><li>FIG. 6 schematically represents collections of traces recorded on the receivers along the V, H axes.₁, H₂,</li><li>FIG. 7 is a schematic representation of trace collections along the R, N, T axes,</li><li>FIG. 8 is a schematic representation of the collections of residual traces R, R signal and R,</li><li>FIG. 9 is a schematic representation of the projections of the residual traces collections R, N and T on the geographically fixed axes X, Y, Z,</li><li>FIG. 10 is a schematic representation of the decomposition of a common emitter trace collection (z, t) in a domain (τ, p),</li><li>FIG. 11 is a schematic representation of the decomposition in the domain (τ, p) of the trace collections X, Y, Z of FIG. 9,</li><li>FIG. 12 schematically represents, in the system of axes X, Y, Z, the synthesis of the pivot trace of FIG. 10, representative of the estimate of the rising P-wave field.</li></ul>

According to the method of the present invention, one or more seismic or emitting sources 1 are used (FIG. 1) which are arranged on the surface 2 of a medium 3 in which a well 4 has been drilled. The well 4 may be vertical or deflected as in FIG. In the deflected well 4 is disposed a receiving device 5 which can occupy at least two positions in the well. The reception device may comprise one or more receivers, for example four receivers 6 each consisting of three sensors whose directions are in a right trihedral formed by three axes, V, H₁, H₂, the axis V being directed along the axis 7 of the well 4, the other axes H₁ and H₂ being located in a plane perpendicular to the V axis (Figure 2).

The distance separating two consecutive positions of the transmitter 1 is generally constant and for example equal to 25 m, whereas the distance separating two consecutive receivers 6 is generally constant and for example equal to 10 m.

According to the present invention, the deflection Ψ and the azimuth θ of the axis 7 of the well 4 are determined in a system of orthogonal axes X, Y, Z, fixed geographically, the Z axis being vertical (FIG. 3). The coordinates of the vector V are given as a function of Ψ and θ by the following formulas, in the system of axes X, Y, Z:<maths id="math0001" num=""><img file="EP0909961B1_D0001.tif" /></maths> the angles ψ and θ being respectively in the intervals [0, π / 2] and [0, 2π].

A new vector h₁ is defined to be perpendicular to both the V vector and the Z axis, and whose coordinates are given by the following formulas in the X, Y, Z axis system:<maths id="math0002" num=""><img file="EP0909961B1_D0002.tif" /></maths>

The vector h₁ is then contained in the plane H₁, H₂.

In the same way, we define a vector h₂ located in the plane H₁, H₂ and perpendicular to h₁. The coordinates of the vector h₂ in the system of axes X, Y, Z are given by the following formulas:<maths id="math0003" num=""><img file="EP0909961B1_D0003.tif" /></maths>

The positions of the axes H₁ and H₂, relative to the vectors h₁ and h₂ are given by the following equations:<maths id="math0004" num=""><img file="EP0909961B1_D0004.tif" /></maths> and<maths id="math0005" num=""><img file="EP0909961B1_D0005.tif" /></maths>

In these formulas, δ represents the angle formed by the axis H₁ and the vector h₁ (Figure 4).

The transmitter 1 successively occupies several positions on the surface 2, said positions of the transmitter 1 being located in given geometrical, rectangular or circular concentric configurations around the drilled well 4.

The offset between the transmitter 1 and the receiving device 5 increases from one program to another, so as to perform a walkaway type acquisition. The transmitter 1 generates elastic waves which propagate in the medium 3 and which are received and recorded on each of the sensors of each receiver 6 of the reception device 5, the waves being recorded as a function of the time t in the form of traces in the directions V, H₁ and H₂. Thus, at each transmission, a set of traces is obtained which are indicated by the depth z of the receiver 6 and the position of the transmitter 1. When the recording corresponding to one emission is completed, another transmission is made for another position of the transmitter and so on for all the predetermined positions of the emission points. After that, the receiving device 5 is moved to another depth in the well 4, as shown in FIG. 1, in which the lower receiver 6 occupies the position 6 '.

Another series of transmissions and recordings is then made for this other position of the receiving device 5.

Each trace recorded on each sensor, that is to say in the directions V, H₁, H₂, is identified by the positions of the transmitter pair 1-receiver 6 having produced said trace.

From at least one set of traces recorded in the directions V, H₁, H₂we create collections of traces C with common emitter (z, t). The processing first consists in extracting the first signal arrivals C of at least one of the collections of traces C and in making a collection of residual traces C from said collection of traces C, each residual trace being the difference between the trace C and the estimate C signal. A pivot trace, located at a depth z is selected in the collection of residual traces C, then we pass a line of slope p which is variable at each sampled point of the pivot trace. For each value of p, a summation is made of the amplitudes of the seismic events which are situated at the intersections of the traces with said straight line. The amplitude sums for the corresponding times and slope are plotted in a domain (τ, p) in which τ is the time and p represents the values ​​of said slope varying between two extreme limits. From the data in the domain (τ, p), a synthetic seismogram is made. The synthetic seismogram is then compared to the common emitter collection of residual C traces. Finally, iterations are carried out until the synthetic seismogram is acceptable compared to the initial data of the common emitter collection of residual traces C.

In a preferred embodiment of the invention, the processing consists in defining, for each pair of transceiver positions, a second system of orthogonal axes R, N, T which is calculated from the components along the three axes V , H₁, H₂ of the displacement vector corresponding to the first arrival of the wave of interest, for example the wave P, and the projection H_max of the displacement vector in the plane H₁, H₂. Most of the energy associated with the wave of interest P or direct wave is located on this axis R which corresponds to said displacement vector.

The axis N is perpendicular to the axis R and it is located in the plane V, H_max and the T axis is perpendicular to the R, N plane (Figure 5).

The passage of the axis system V, H₁, H₂ to the system of axes R, N, T is carried out by a rotation of angle Φ_H around the V axis and another angle rotation Φ_V around a vector H_min which is perpendicular to the vector H_max. The angles Φ_H and Φ_V are calculated from the three components V, H₁, H₂ present in a time window located around the direct arrival. The angles Φ_H and Φ_V can be calculated using for example the software SEISLINK (1996 version) marketed by the company WESTERN.

For each pair of transceivers, the initial data V, H₁, H₂ provided by the sensors are projected in the second system of axes R, N, T. The components along the axes R, N, T are expressed using the following equations:<maths id="math0006" num=""><math display="block"><mrow><msub><mrow><mtext>R = V cos Φ</mtext></mrow><mrow><mtext>V</mtext></mrow></msub><msub><mrow><mtext> + (H</mtext></mrow><mrow><mtext>1</mtext></mrow></msub><msub><mrow><mtext>power factor</mtext></mrow><mrow><mtext>H</mtext></mrow></msub><msub><mrow><mtext>+ H</mtext></mrow><mrow><mtext>2</mtext></mrow></msub><msub><mrow><mtext>sinΦ</mtext></mrow><mrow><mtext>H</mtext></mrow></msub><msub><mrow><mtext>) sinΦ</mtext></mrow><mrow><mtext>V</mtext></mrow></msub></mrow></math><img file="EP0909961B1_D0006.tif" /></maths><maths id="math0007" num=""><math display="block"><mrow><msub><mrow><mtext>N = -VcosΦ</mtext></mrow><mrow><mtext>V</mtext></mrow></msub></mrow></math><img file="EP0909961B1_D0007.tif" /></maths><maths id="math0008" num=""><math display="block"><mrow><msub><mrow><mtext>T = -H</mtext></mrow><mrow><mtext>1</mtext></mrow></msub><msub><mrow><mtext>sinΦ</mtext></mrow><mrow><mtext>H</mtext></mrow></msub><msub><mrow><mtext> H</mtext></mrow><mrow><mtext>2</mtext></mrow></msub><msub><mrow><mtext>power factor</mtext></mrow><mrow><mtext>H</mtext></mrow></msub></mrow></math><img file="EP0909961B1_D0008.tif" /></maths>

The angle δ is calculated for each receiving position, using an inversion process that uses angles Φ_H, the relative positions of the transmitter with respect to the receiver, the deflection ψ and the azimuth θ of the well. The value of the angle δ is calculated by minimizing the following function:<maths id="math0009" num=""><img file="EP0909961B1_D0009.tif" /></maths> in which α represents the angle formed by the plane containing the transceiver pair and the X axis.

In another step and after forming common emitter collections of the trace components along the R, N, T axes (FIG. 7), each common emitter collection of traces R is processed to extract the downward P-wave field.

For this operation, carried out independently on each common emitter collection of traces R, the matrix filtering technique (CEPHAG-IFP) is used. After horizontalising the data using the first-arrival dash, this technique is applied over the entire available time or only a portion of the time (this avoids the need to involve another descending event of significant amplitude that would distort the data. estimation - down-wave S for example). The technique is based on a decomposition of interspectral matrices after frequency averaging and / or space; the filtering parameters are chosen so as to allow the extraction of the descending P field under the best conditions. Frequency averaging is performed so as to retain only the first terms of intercorrelation functions, no spatial averaging being performed. For each initial collection R, the present invention consists in providing, at the output, a common emitter collection called R signal, these signals constituting the estimation of the P field falling for the source position considered. A common emitter R residual or residual R collection is obtained by difference between the initial collection R and the estimation R signal (FIG. 8).

The R collections thus estimated can be used later for the calculation of the anti-multiple operators.

The residual R, N and T components are projected in the fixed reference X, Y, Z, using the angles Φ_H, Φ_V and δ as well as the deflection ψ and the azimuth θ of the well. The components thus obtained are called X residue, Y residue and Z residue (Figure 9).

For each collection of common emitter traces, on the one hand, and each component X residue, Y residue and Z residue, on the other hand, we proceed to a decomposition in a domain (τ, p) using an inversion technique . Figure 10 schematically shows a collection of common emitter residual traces. In this collection of residual traces, a pivot trace 10 is selected at a depth z. A line 11 of variable slope p is passed through each sample t of the pivot trace. The amplitudes of the seismic events 12 located at the intersections of the traces 10 and 13 to 16 with said line 11 are summed along said line 11, for each value of p. The amplitudes A are reported for the corresponding times and slopes in the domain (τ, p) in which p represents the values ​​of the slope and τ the time corresponding to t.

From the data of each of the X residues, Y residue, Z residue in the domain (τ, p), a synthetic seismogram is produced which is compared to the corresponding collection X residue, Y residue or Z residue and is carried out. iterations until the synthetic seismogram is acceptable with respect to the data of the corresponding collection X residue, Y residue or Z residue.

To avoid the spreading of events in the domain (τ, p), an L-shaped regularization is applied in the inversion process.₁. FIG. 11 represents the decomposition in the domain (τ, p) of the trace collections X residue, Y residue and Z residue, of FIG. 9. From the decomposition in the domain (τ, p) of the components X residue, Y residue and Z residue, we compute the polarization attribute ζ and the azimuth attribute γ for each value of p and τ on a window sliding, using the proper elements of the variance and covariance matrix included in the aforementioned software SEISLINK. These attributes ζ or polar (τ, p) and γ or azi (τ, p) are used to calculate the attribute V_P (speed of the P wave in the medium considered), using the following formula:<maths id="math0010" num=""><math display="block"><mrow><msub><mrow><mtext>V</mtext></mrow><mrow><mtext>P</mtext></mrow></msub><mtext>(τ, p) = </mtext><mfrac><mrow><mtext>(Cos [polar (τ, p)] cosΨ + sin [polar (τ, p)] sinΨcos [azi (τ, p) -θ]) </mtext></mrow><mrow><mtext>p</mtext></mrow></mfrac></mrow></math><img file="EP0909961B1_D0010.tif" /></maths>

The apparent slowness p_z(τ, p) with respect to the Z axis of the event considered can also be evaluated by the formula:<maths id="math0011" num=""><math display="block"><mrow><msub><mrow><mtext>p</mtext></mrow><mrow><mtext>z</mtext></mrow></msub><mtext>(τ, p) = </mtext><mfrac><mrow><mtext>cos [polar (τ, p)] </mtext></mrow><mrow><msub><mrow><mtext>V</mtext></mrow><mrow><mtext>P</mtext></mrow></msub><mtext>(Τ, p)</mtext></mrow></mfrac></mrow></math><img file="EP0909961B1_D0011.tif" /></maths>

The Z axis being vertical, the slowness p_z(τ, p) is used to discriminate rising events from down events.

Knowledge of V attributes_P(τ, p) and p_z(τ, p) makes it possible to define a selection (mute) of the events to preserve in the domain (τ, p). It may be useful on this occasion to use the knowledge that one can have of the geological velocity V_P the medium in which is positioned the receiving device.

One can for example apply a mute (selection) in the domain (τ, p) using the different attributes described above so as to keep only the events P amounts. These events are then synthesized to construct a seismogram in the domain (X, Y, Z, t) representative of the rising P-wave field. FIG. 12 shows the estimate of the X, Y and Z components of the rising P wave field corresponding to the pivot trace 10 to which the decomposition in the domain (τ, p) has been applied.

Identical processing can be performed to extract the downward and upward S wavefields.

In the case where the axis 7 of the well 4 is vertical (ψ = 0), the value θ is arbitrarily set to zero.

The method according to the invention is more flexible than previous methods of optimal filtering or parametric inversion. It does not require pointing outside of that of the first arrival and it is not necessary to define the number of waves to extract; moreover, it is not necessary to determine the apparent speeds.

The method according to the invention is therefore very easy to use and can be used as such for 3D acquisitions. In addition, it provides an opportunity to analyze the data and control the result of the separation.

Of course, it is possible to proceed with the decomposition (τ, p) of the data represented in the system of axes R, N, T. It is also possible to calculate the attribute polarization ζ and the attribute azimuth γ relative to at this mark R, N, T, as well as the attribute V_P in the domain (τ, p).

The collections of traces R, N, T can be replaced by collections C obtained from components V, H₁, H₂, the other steps of the method according to the invention being performed on said collections C so as to obtain the collections of residual traces C necessary for the decomposition in the domain (τ, p).

Similarly, it is possible to project in the system of axes (X, Y, Z) the collections of residual traces C and their collections of associated traces to obtain the collections of traces X residue, Y residue and Z residue which are then subject to decomposition in the domain (τ, p).

The present invention can be applied for the fine imaging of the reservoirs as well as for the calibration of the AVO measurements (measurement of the reflection coefficient versus angle of incidence). Indeed, from walkaway data, it is possible to measure the incident and reflected fields and by deconvolution, to make a possible estimate of the reflection coefficient, this measurement being all the more exact as the sensor is located near the objective studied. In addition, the polarization analyzes lead to the knowledge of the angles of incidence. As a result, all the elements are used to evaluate the behavior of the reflection coefficient vis-à-vis the angle of incidence without bias.

Finally, the present invention can be used for the measurement of the anisotropy at the line of the receiving device, said anisotropy being obtainable by a joint analysis of the apparent speed and the polarization of the direct arrival.