WO2023027593A1

WO2023027593A1 - Method of simulating seismic data

Info

Publication number: WO2023027593A1
Application number: PCT/NO2022/050203
Authority: WO
Inventors: Marianne HOUBIERS; Fredrik Hansteen
Original assignee: Equinor Energy As
Priority date: 2021-08-25
Filing date: 2022-08-24
Publication date: 2023-03-02
Also published as: GB202112171D0; NO20240186A1; GB2610201B; GB2610201A

Abstract

A method of using a source of substantially continuous noise to obtain seismic data simulating data obtained from an impulsive source located. The method comprises a) obtaining, for each of a plurality of seismic sensors, a seismic data stream, b) pre- processing and cross-correlating the seismic streams to obtain a plurality of cross-correlation data sets, c) estimating the location of the noise source by searching a grid of locations with respective semblances, where the semblance is obtained by moveout correcting each of said cross-correlation data sets using a velocity model, and semblance stacking the moveout corrected data, d) processing the data sets to enhance components associated with said noise source, and e) generating the simulated data by stacking each of said cross-correlation data sets, performing moveout correction, and reintroducing absolute source-receiver travel times.

Description

Method of Simulating Seismic Data

Technical Field

The present invention relates to the field of simulating seismic data and more particularly to the field of using a substantially continuous source of noise, located within a subsurface formation of the Earth, to obtain seismic data simulating data that would be obtained from an impulsive source located at substantially the same location as the continuous noise source.

Background

It is known that the acoustic signal emitted by a drill bit deep in the subsurface and recorded on seismic sensors disposed on the surface or seabed above the drill bit, or in an adjacent well, can be used for both accurately determining the drill bit location and for seismic imaging while drilling. One such proposal is described in WO2021145778. Related techniques are also discussed in: “Walk- away VSP using drill noise as a source”, Haldorsen et al, Geophysics (1995) 60, 978-997, and “Seismic while drilling: Fundamentals of drill-bit seismic for exploration”, Poletto and Miranda (2004), Elsevier 35, and US4922362 by Miller and Haldorsen. An example acquisition geometry, where all the sensors are located at the seabed, is illustrated in Figure 1 . In practice, hundreds or even thousands of seismic sensors may be employed.

The main challenge when using drill bit noise as a source for seismic imaging is the long random source signature which is both unknown and lacks time synchronization with the seismic receivers. This is very different from a commonly used “impulsive” noise source, such as an explosion or impact, which is strictly time synchronized with the recording instruments.

The physics of the wave propagation is independent of the source duration. The continuous acoustic wave train emitted by a drill bit therefore reflects, refracts, and diffracts as it encounters varying acoustic or elastic properties on its way towards the receivers, in the same way as a seismic impulse. The signal recorded on a sensor can thus be considered a convolution of the source signature with an impulse response function of the subsurface, specific to the given source and receiver locations. In a conventional seismic record, acquired with an impulsive source, arrivals corresponding to the reflections, diffractions, converted waves, etc. are clearly visible and can be interpreted or processed with seismic migration algorithms to produce 3D seismic images of the subsurface. For drill bit seismic to be interpretable or suitable for conventional seismic imaging, the source signature must be deconvolved from the recorded data.

Figure 2 illustrates an exemplary receiver geometry for a Permanent Reservoir Monitoring (PRM) system on the sea bed. A total of 3458 sensors are deployed on cables with 50m inline spacing and 300m cross-line spacing. Receivers up to 1600m lateral distance away from a drill bit location (highlighted) may be used in the analysis of the drill bit signal. Figure 3 shows a raw seismic record having a duration of one second recorded by the selected sensors of the PRM system of Figure 2 during a formation drilling operation. The displayed seismic traces are sorted by cable and by receiver along the cable. The most conspicuous signal, with parabolic moveout, is due to a nearby vessel producing a continuous noise signal. The much weaker drill bit noise is not visually identifiable. Separation and interpretation of the various seismic arrivals coming from the continuous drill bit source, while minimizing the disturbance from other noise sources (e.g. the vessel), is among the key objectives of drill-bit seismic methods.

The problem of separating the source signature from the impulse response function, when these are both unknown, is generally referred to as “blind deconvolution”. It is an ill-posed problem that arises in many situations, from the design of noise cancelling microphones, to ultrasonic testing and under water acoustics. It is a much-studied problem, and a wide range of approaches have been described in literature. As an example from the seismic domain, "Focused Blind Deconvolution", Bharadwaj, Demanet and Fournier, IEEE T ransactions on Signal Processing, (2019) vol. 67, no. 12, pp. 3168- 3180, doi: 10.1109/TSP.2019.2908911 , teaches a method that imposes statistical constraints on the inversion to ensure that the impulse response functions are “maximally white” and “maximally front-loaded”

When the source signature is known, e.g. from measurements at or near the source, deconvolution methods can be directly employed to separate the signature from the impulse response function of the subsurface. However, with current technology, the bandwidth to transfer such measurements, made down-hole near a drill bit, up to the surface - while drilling - is not widely available. Highly stable and synchronized clocks would also be required to preserve information about travel times from the drill bit to surface seismic sensors. An alternative to down-hole pilot recordings is to place a sensor on the well-head such as is taught by “The use of drill-bit energy as a downhole seismic source”, Rector et aL, Geophysics (1991 ), 56(5), 628. This approach records the pilot signal after it has propagated from the bit, along the drill string, to the surface. It is however not a direct representation of the bit signature, but rather a convolution of the source signature with the transfer function of the string that may include also converted modes and multiple reflections from the elements that make up the drill string and bottom hole assembly (BHA). Following careful data processing, such drill string pilot signals have been reported to work well as a substitute for down-hole recordings near the source. Travel time information is, however, still challenging to preserve accurately, as the transit time of the pilot signal along the drill string varies with drilling parameters and with depth as drill pipe is added.

When no dedicated pilot recording is available, the data from a wide aperture seismic sensor array can be used to estimate the pilot signal. One common approach is to stack data from multiple receivers after correction for the moveout of the direct arrival. The concept is known as beam steering, or "focused stack", and emphasizes energy from a particular direction by delaying the successive channels (“Seismic while drilling: Fundamentals of drill-bit seismic for exploration”, Poletto and Miranda (2004), Elsevier 35). In order to obtain the required moveout correction times, the method relies on a coherence analysis of the received signal, often performed iteratively and in combination with deconvolution (“Walk- away VSP using drill noise as a source”, Haldorsen et al, Geophysics (1995) 60, 978-997). The source location and the absolute travel time from source to receiver are initially not available with this approach, but can be estimated by parameterizing travel times in terms of depth and average velocity, and fitting of the moveout correction times (US4922362). If a detailed velocity model is available, estimates of source location and absolute travel times can be improved. The resulting pilot stack can be used for deconvolution or as a correlation template to produce data resembling that from an impulsive source (“Seismic while drilling: Fundamentals of drillbit seismic for exploration”, Poletto and Miranda (2004), Elsevier 35).

Summary According to a first aspect of the present invention there is provided a method of using a source of substantially continuous noise located within or above a subsurface formation of the Earth to obtain seismic data simulating data that would be obtained from an impulsive source located substantially at the continuous noise source location. The method comprises a) obtaining, for each of a plurality of seismic sensors at respective locations above or in said formation, a seismic data stream recorded whilst operating the source of substantially continuous noise; b) for each seismic sensor, pre-processing and cross-correlating the associated seismic stream with the seismic streams of each of the other seismic sensors, to obtain a plurality of cross-correlation data sets, one for each seismic sensor; c) estimating the location of the continuous noise source by searching a grid of locations with respective semblances, where the semblance for each grid point is obtained by

- moveout correcting each of said cross-correlation data sets according to gridpoint and sensor pair dependent travel times computed using a velocity model, and

- semblance stacking the moveout corrected data; d) either before or after moveout correction of the cross-correlation data sets, processing the data sets to enhance components associated with said continuous noise source; and e) generating seismic data simulating data that would be obtained from an impulsive source by stacking each of said cross-correlation data sets, moveout corrected for the estimated location of the continuous noise source, and reintroducing absolute sourcereceiver travel times consistent with the identified source location and velocity model.

Step d) may comprise manually or automatically analysing the cross-correlation data sets to determine properties of the signal and (unwanted) noises, and designing targeted filters.

Step c) may comprise searching said grid to identify a grid point or location intermediate the grid points having a maximum semblance. The processing of step d) may comprise one or more of; applying weighting factors, selection of specific receivers or receiver pairs, applying targeted filters to mute other sources of noise.

The processing of step d) may be performed iteratively with step c) until a suitably accurate location is obtained. The said seismic data streams may be collected over a period of five minutes or less.

The source of substantially continuous noise may be a drill bit in operation within said formation.

According to a second aspect of the present invention there is provided a method of obtaining an image of a subsurface formation of the Earth, the method comprising: using the method of the above first aspect of the invention to obtain seismic data simulating data that would be obtained from an impulsive source located substantially at the continuous noise source location; and performing an inversion or seismic imaging process on the data to obtain said image.

According to a third aspect of the present invention there is provided a non-transitory computer-readable medium storing a program including instructions that, when executed by a processor, cause the processor to implement the method of one of the preceding aspects of the invention.

Brief Description of the Drawings

Figure 1 illustrates schematically a drill bit emitting continuous noise that is recorded on seismic receivers at the seabed;

Figure 2 shows an exemplary PRM array deployed on the sea bed;

Figure 3 presents as an example a one-second long record of raw data from a set of seismic receivers of the PRM array of Figure 2, acquired during a drilling operation;

Figure 4 illustrates in the upper panel the direct wave travel times to each receiver on the seafloor from an assumed subsurface source location, computed by ray tracing through a known velocity model, whilst the lower panel shows a vertical section from a 3D P-wave velocity model;

Figure 5 illustrates interferometric signals

for a selection of receivers, computed by cross-correlating all pairs of receivers where i #= j, with the data moveout corrected with the direct wave travel time differences from source location to each receiver i and j, so that the correlation event between direct arrivals appears flat at zero correlation time lag;

Figure 6 illustrates three perpendicular cross-sections through the semblance volume used for locating the drill bit source by interferometric imaging; Figure 7 illustrates an interferometric panel dj₇(-r) shown for one receiver i = 360 after moveout correction with direct wave travel time differences (one vertical slice of Figure 5);

Figure 8 shows an estimated common shot gather from an impulsive source at the drill bit location, produced by stacking the focused interferometric signals di₇(r) over correlated receivers j and reintroducing the moveout associated with the direct wave from the found source location to each receiver i;

Figure 9 Leftmost panel: The same moveout corrected interferometric traces as shown in Figure 7, here sorted by source-receiver j offset. Middle panel: Variable density plot showing the cumulative stack over offsets from near to far. Spurious events with moveout different from the direct wave are suppressed. The rightmost panel illustrates the effect of cumulative stacking by plotting the stack trace for every 25^th offset trace included. The black trace in bold represents the stack over all offsets; and Figure 10 is a flow diagram illustrating an embodiment.

Detailed Description

The following disclosure relates to a method whereby a source location and a consistent time reference for impulsive data is obtained using a detailed velocity model, available for example from a full waveform inversion (FWI) method or other velocity model building method. Figure 4 shows (in the lower panel) a vertical section taken from such a 3D velocity model and (in the upper panel) an illustration of the travel times obtained by ray tracing from a source in the subsurface to a set of receivers on the seabed above.

For depth imaging of drill bit seismic to be consistent relative to existing 3D seismic, the source location must be precisely determined. The location of the drill bit, commonly derived from magnetic and/or gyroscopic measurements while drilling (MWD) or on wireline after drilling, is typically not known with a lateral accuracy better than several tens of meters. The error associated with MWD measurements accumulates along the well, and can become relatively large, especially for long horizontal wells. As a result, it is not unusual for long horizontal wells to have a lateral positional uncertainty of more than ± 60m at 6000m measured depth (MD). The situation is similar for measurements made using wireline tools after drilling of a well section. With a velocity model and a wide aperture sensor array, bit and well positioning accuracy can be improved and the error accumulation of MWD avoided. Very high lateral accuracy can be achieved relative to other wells, and relative to conventional active source 3D seismic acquired with the same (permanent) sensors and migrated with the same velocity model. Different methods exist to process drill bit signals for source localization. WO2021 145778 describes an approach based on time domain semblance stacking of moveout corrected data, followed by a grid search over candidate subsurface locations. In this case, it is the moveout of the direct wave, modelled by ray tracing through an anisotropic 3D velocity model, that allows discrimination of source locations. Other related approaches include coherence analysis in the Fourier domain (“Walk- away VSP using drill noise as a source”, Haldorsen et al, Geophysics (1995) 60, 978-997) and interferometric imaging of the source (“Interferometric/daylight seismic imaging”, Schuster et al, Geophys. J. Int., (2004), 157, 838-852).

It is worth noting that the drill bit signal is usually much weaker than that from conventional seismic sources, and usually also weaker than the ambient noise from vessels, rigs, storms, and ocean swell, as is illustrated in Figure 2. Methods to enhance signal and suppress noise are therefore essential for successful real-world applications of drill bit noise for seismic imaging.

The disclosed method uses interferometric signals, also known as cross-correlograms, computed from recordings on an array of seismic sensors to a) estimate the drill bit source location via a velocity model, and b) estimate the seismic records (with consistent absolute travel times) that would have been obtained at the receivers were an impulsive (band limited) source to be located at the drill bit location. A key feature of the method is that all processing, from source localization to filtering and stacking to enhance signals and to suppress noise, both from other sources and from spurious events introduced by the cross-correlation, takes place in the interferometric domain.

In essence, what is described is an energy focusing approach to estimating the impulsive source data from continuous source data, with no data input other than seismic array recordings and a velocity model. The method is similar to the focused pilot method. However, by performing cross-correlations before focusing and stacking, the presented method offers a view into the interferometric (correlogram) domain where properties of signal and noise can be assessed, and various filters applied in order to suppress noise and enhance the drill bit signal. By using a velocity model for locating the source, consistent absolute travel times can be recovered for the interferometric shot records.

A signal d_i t') recorded by a seismic receiver i can be described as the source signature s(t) convolved with the impulse response function g^t), also denoted “Green’s function”, that describes the wave propagation from the source to the receiver, including the direct wave arrival, all reflections and multiple reflections, diffractions, dispersion, converted waves, and other wave phenomena.

We define t’=0 as the arbitrary start of a recorded time series and reserve the unprimed variable t to later refer to time relative to our synthesized impulsive source. For simplicity, noise from other uncorrelated sources has been neglected in the following discussion which primarily aims to explain how the source is located and why focusing and stacking of the interferometric signals yields an estimate of the impulse response function.

To prepare the continuous sensor data for use in a conventional seismic imaging process, we aim to extract an estimate of g_t from our measurements

Our approach assumes that a velocity model is available (analytic or gridded, simple or detailed, isotropic or anisotropic), sufficient to compute by ray tracing, Eikonal solvers or other methods, the travel time of a direct wave from any given subsurface location to all the available receivers. A velocity model is illustrated in Figure 3, with an array of receivers installed in the seabed (in this example at about 120m depth). A drill bit source at an unknown location, or a location approximately known from MWD, is emitting an unknown seismic signal s(t) with continuous random character.

The travel time information is used to estimate the source location by interferometric imaging, and subsequently for aligning, or focusing, the interferometric data in such a way that stacking over receiver correlations produces an estimate d₍(t) of the impulse response function g^t) convolved with the source autocorrelation function. For a band limited white noise source, this would represent an estimate of the seismic data that one would record at receiver i if one had an impulsive source at the drill-bit location instead of a continuous source. We compute the interferometric signals, defined here as the cross-correlation of data from all pairs of receivers i and j, i #= j. Figure 5 presents an illustration of the resulting three-dimensional data set (after time aligning traces as discussed later).

Here ® denotes cross-correlation, s_a is the source autocorrelation function and g^ = gi®gj is the interferometric Green’s function of the receiver pair (i,j) for a given source location. The time lag T represents the travel time difference between an event recorded at both receivers i and j. An event with sufficiently large spectral bandwidth will be visible as a peak at a discrete T. The benefit of the cross-correlation is that it reduces long time series of continuous noise data to much shorter correlation series with the event travel time differences between receivers preserved, and with impulsive source character in the cases that the source spectrum is sufficiently white. The main disadvantage is the introduction of unphysical spurious events caused by cross-terms. It is the purpose of the presently described method to efficiently suppress these, as well as other noises.

In practice, the step of cross-correlating long time-windows of sensor data to obtain

will often be preceded by one or more standard signal pre-processing steps, such as rotation of the geophone components to a common reference frame, bandpass filtering to the expected band of the bit noise, and shaping or whitening of the amplitude spectrum. The length of the time series used in the cross-correlation will depend on the signal to noise ratio of the data, the desired time resolution, and the speed of the drilling. With typical penetration rates from drilling of the order of 20 m/hour, the bit moves about 0.3m per minute, and its’ location may safely be considered stationary for time windows of a few minutes. For comparison, the shortest seismic wavelength of a 100Hz signal in a formation with velocity 2000m/s is 20m.

Imaging of the source distribution in passive seismic data via interferometric signals is a concept known from (“Interferometric/daylight seismic imaging”, Schuster et al, Geophys. J. I nt. , (2004), 157, 838-852). To locate the drill bit source, we implement this as a grid search for the maximum semblance, or stack amplitude, over a set of subsurface points covering a volume within which we expect the source to be found. The direct wave travel time from each grid point is computed based on a velocity model of the subsurface, using an Eikonal solver. For real-time applications, travel time tables can be pre-computed and stored on disk or in memory. Specifically, at each grid point, the interferometric signals

are moveout corrected according to the direct wave travel time difference between receivers i and j for that particular grid location, and stacked over receivers i and j within a specified offset range where the signal from the drill bit is expected to be received. For each grid location, the semblance is computed from the stack (this is a known computation). The precise location of the drill bit is associated with the maximum semblance grid point or refined to sub-grid resolution by fitting a function to the semblance values obtained on the compute grid.

Figure 6 shows perpendicular sections through a 3D semblance volume, each centred at the interpolated maximum. The semblance at any location (x,y,z) in the subsurface is formed by stacking over the entire cube of Figure 5 after the moveout correction with the corresponding travel time differences of direct arrivals on receivers i and j for that particular location, and subsequently computing the semblance of the stack. Note that the drill bit location is more sharply defined in the horizontal x-y plane than the vertical (depth) planes due to all seismic sensors being located at a similar depth, at the seabed. Combining the surface sensor geometry with, e.g., down-hole DAS or geophone sensors can improve the depth resolution.

To further illustrate the method, we can neglect dispersion and write the generally complex valued g_t as a sum of delta functions with amplitudes m_ik and peak times t_ik, representing the arrivals of the direct wave (k = 0), reflections, multiples, and other events.

Note that the result of cross-correlating two delta functions <5(t - t_a) and 8(t - t_b), peaking at times t_a and t_b, respectively, is a delta function 8(T - t) peaking at the time difference 8t = t_a - t_b. Assuming that the first term in each sum (k = 0, I = 0) represents the direct wave from source to receivers i and j respectively, the term representing their cross-correlation becomes

and peaks at the direct arrival time difference between i and j,

= t_i0 - t_]0. Thus, interferometric imaging of the source location can be viewed as the search for the subsurface source location r with associated travel times t_i0(r) and t₇₀(r) that optimally flattens the cross-correlation event between the direct waves at the two receivers after moveout correction with +At_i7. Figures 5 and 7 show examples of the optimally flattened direct wave cross-correlation event using the traveltimes for the correct source location. As direct waves commonly have the largest amplitude (m_i0 > m_ik, for k > 0), they will dominate in the image or semblance stack. Note that, if the source approaches a very strong reflecting interface, it may become difficult to get an accurate depth estimate as there is very little difference between the direct and reflected wave moveout. A weighting or selection of the receivers with the most moveout discrimination may improve the resolution.

Raw recorded time series d^t') may also be used to estimate the drill bit source location, by direct semblance stacking and grid search, as shown by (WO2021 145778), but is much more efficiently done using the shorter interferometric signals once they have been computed. The cost of computing the cross-correlations can, however, become relatively high and scales as the number of sensors squared. By taking advantage of modern GPU or FPGA hardware, cross-correlation of long time series can be efficiently parallelized to achieve near real-time processing even for large sensor arrays consisting of thousands of receivers. Not all available receivers must be cross-correlated; a selection can be made according to their distance from the source, the character of noise from other sources that should be suppressed, the strength of the direct or reflected signals, as well as their moveout on different receivers. Furthermore, the time lag can be limited to the range of interest.

Once an estimate of the drill bit location has been found, the interferometric signals c?i₇ (T) are moveout corrected with the travel time differences At_i7 = t_i0 - t_j0 associated with the found source location, and can be stacked along dimension j. Then by re-applying the direct arrival travel times t_i0, one generates reverse VSP gathers suitable for seismic imaging - either directly, in the case of a sufficiently white source spectrum, or after deconvolving the source autocorrelation signature s_a.

In the following we describe the effect of this focused stacking

and why it results in an estimate g_t of the impulse response between a source at the drill bit location and receiver i that is kinematically correct and with relative event amplitudes preserved. We also discuss refinements by filtering or weighting of amplitudes.

Observe that, after the source location has been found, the effect of moveout correcting g_ij(r') with the direct wave travel time differences At_i7 can be expressed as

When the source location and direct wave travel times are known, time shifting by At_i7 aligns all these (g_t vs direct wave at receiver j) correlation events. The result is the desired impulse response at receiver i multiplied by the amplitude m₇₀ of the direct wave on receiver j. Hence, for each correlated receiver j we have available a differently scaled version of the desired g_t (albeit time shifted by t_i0), but time aligned for all j. Note that the time between events in g_t and their relative amplitudes are preserved.

Figure 7 shows an example interferometric panel d_i7(T + At ₇), also called a “correlogram”, for the receiver i = 360 after travel time correction (focusing) with direct wave travel time differences Ati₇. This corresponds to one section from the cube in Figure 5. The panel represents receiver trace d₃₆₀ cross-correlated with traces d₇ from all other receivers j. The applied moveout correction is based on the drill bit source position and aligns the correlation events between the direct wave on receivers j and all events on receiver i originating from the same assumed source, so that these appear as flat in the panel. The direct P-arrival, free surface multiples and a direct shear wave arrival are among the most clearly visible. Other contributions to the cross-correlation (e.g. from other sources or generated by cross-correlation with arrivals other than the direct wave on receivers j) are generally not flattened and can be suppressed by stacking over j, or by filtering prior to stacking. Examples of this are the indicated spurious events at negative lag times that are due to the cross-correlation of the first two orders of free surface multiples on receiver j with the direct wave on receiver i. Noise from sources other than the drill bit, such as nearby vessels and rigs, can be seen as events with different moveout and source signature passing through the array.

Flat events at positive lag represent the events of gt correlated with the direct wave on receiver j. The wavelet s_a appears to resemble a sine function in this example, consistent with a band limited white noise source. Stacking over all, or only selected j, produces an estimate, g_t, of g_t in which noises and spurious events that do not conform to the moveout of the direct wave from the found source location are suppressed. This is apparent by observing that the remaining (Z > 0) terms of g_tj in Equation (7), when subject to the same time shift At_i7, become

and are generally not flat after the direct wave focusing. These spurious events do not represent physical events, but are artifacts of the cross-correlation that appear as time delayed replicas of gt for each event I > 0 in gj (i.e. reflections, multiples, etc.). In addition to having different move out than the direct wave, these are scaled by the factors rriji that are smaller than the direct wave I = 0 amplitudes. Spurious events are known to be difficult to suppress in passive source data (“Interferometric/daylight seismic imaging”, Schuster et al, Geophys. J. Int., (2004), 157, 838-852), and strong reflectors near the source will pose the largest challenge. In our approach, the suppression can be optimized by stacking over selected receivers j that have moveouts sufficiently different from the direct wave. The assessment of which correlated receivers to stack is done on correlograms like the one shown in Figure 7.

When stacking the aligned g^ over selected j, we obtain our estimated Green’s function for receiver i, composed of contributions from the direct wave correlation and all other events in j (where the former is expected to dominate and the latter to stack out).

Time shifting the resulting stack by t_i0 recovers absolute travel times of the Green’s functions estimates

Similarly, in terms of interferometric data di₇(r)

That is, we estimate a common shot gather {<3j(t)} from an impulsive source with signature s_a at the drill bit, from the interferometric signal d_tJ, by moveout correcting with respect to direct wave travel time differences for the found source location, summing over j, and finally re-applying the direct arrival travel times t_i0.

Figure 8 shows a resulting common shot gather after stacking over j and recovery of the absolute travel times. Four orders of free surface multiples can be identified, as well as a direct shear wave arrival. Reflections are also present but harder to distinguish in this display. Spurious events ahead of the direct P-wave (seen at negative lags in Figure 7) are not harmful and can easily be “muted” away. These events originate predominantly from the correlation of the first order free surface multiple on receiver j with the direct wave on receiver i. [Spurious events at positive lag times would be far more harmful and may leak into seismic images but are largely suppressed in this example.]

The noise suppressing effect of stacking focused interferometric traces is further illustrated in Figure 9. Starting from near offsets, increasingly more traces are stacked to achieve a better suppression. Note however, that a more regular selection of offsets may give a better result, and that the optimal selection generally will depend on receiver geometry. Weighting factors w₇ can be introduced to the terms of the sum in equation (11 ) to compensate for the different level of direct wave amplitudes at sensors j, caused by spherical spreading, absorption, scattering, or source directivity. Such weighting, prior to stacking, constitutes a simplistic way to filter the data to improve the signal to noise ratio. For example, a spherical spreading correction due to amplitude decay could be introduced as w₇ =

where p_k is the 3-dimensional distance from the source point to sensor k. Alternatively, a data driven scaling to weight all offsets equally in the interferometric stack may take the form w₇ = where m₇₀ is an estimate of the direct

wave amplitude measured on data in the interferometric domain, e.g. the peak amplitude or a windowed average. Radon or dip filters, to target hyperbolic or linear noise, are other examples of common filters that can be applied in the interferometric domain before stacking to obtain the estimated impulsive source data.

With reference to Figure 6, note in particular how other sources of noise may be filtered to improve the stack, and how a wide aperture array is needed to capture data with sufficient moveout variation to suppress the spurious events. The addition of down-hole sensors, conventional or DAS, may offer scope to record data with better moveout discrimination and thus to further improve the suppression of spurious events.

We note that the commonly used focused pilot method, whereby a source signature is estimated for an assumed source location, or based on coherency analysis, may produce results equivalent to the present method. However, this happens only under the following conditions: the same direct wave focusing time shifts are applied; the same set of receivers {/} enter into the pilot stack; no filtering is applied in the interferometric domain. This follows from the fact that the order of summation over traces and crosscorrelation can be swapped. The advantage of first computing the interferometric signals dtj is the insight they provide into which receivers contain valid arrivals, and which suffer from noise that should be targeted with filters. “Source shaping” by filtering and selective stacking can enhance the signal. A similar assessment of which traces to include in the pilot stack could not easily be conducted on the raw geophone traces with the long-acting source signature.

Embodiments of the invention may be implemented using a suitable program executed on a processor or processors of a computer system. Figure 10 is a flow diagram further illustrating a method of using a source of substantially continuous noise to obtain seismic data simulating data that would be obtained from an impulsive source located substantially at the continuous noise source location. The illustrated steps are as follows:

51 ) obtaining, for each of a plurality of seismic sensors at respective locations above or in said formation, a seismic data stream recorded whilst operating the source of substantially continuous noise;

52) for each seismic sensor, pre-processing and cross-correlating the associated seismic stream with the seismic streams of each of the other seismic sensors, to obtain a plurality of cross-correlation data sets, one for each seismic sensor;

53) estimating the location of the continuous noise source by searching a grid of locations with respective semblances, where the semblance for each grid point is obtained by

- semblance stacking the moveout corrected data;

54) either before or after moveout correction of the cross-correlation data sets, processing the data sets to enhance components associated with said continuous noise source; and

55) generating seismic data simulating data that would be obtained from an impulsive source by stacking each of said cross-correlation data sets, moveout corrected for the estimated location of the continuous noise source, and reintroducing absolute sourcereceiver travel times consistent with the identified source location and velocity model.

It will be appreciated by the person of skill in the art that various modifications may be made to the above described embodiments without departing from the scope of the present invention.

Claims

1 . A method of using a source of substantially continuous noise located within or above a subsurface formation of the Earth to obtain seismic data simulating data that would be obtained from an impulsive source located substantially at the continuous noise source location, the method comprising: a) obtaining, for each of a plurality of seismic sensors at respective locations above or in said formation, a seismic data stream recorded whilst operating the source of substantially continuous noise; b) for each seismic sensor, pre-processing and cross-correlating the associated seismic stream with the seismic streams of each of the other seismic sensors, to obtain a plurality of cross-correlation data sets, one for each seismic sensor; c) estimating the location of the continuous noise source by searching a grid of locations with respective semblances, where the semblance for each grid point is obtained by

- semblance stacking the moveout corrected data; d) either before or after moveout correction of the cross-correlation data sets, processing the data sets to enhance components associated with said continuous noise source; and e) generating seismic data simulating data that would be obtained from an impulsive source by stacking each of said cross-correlation data sets, moveout corrected for the estimated location of the continuous noise source, and reintroducing absolute source-receiver travel times consistent with the identified source location and velocity model.

2. A method according to claim 1 , wherein step c) comprises searching said grid to identify a grid point or location intermediate the grid points having a maximum semblance.

3. A method according to claim 1 or 2, wherein processing of step d) comprises one or more of; applying weighting factors, selection of specific receivers or receiver pairs, applying targeted filters to mute other sources of noise.

4. A method according to any one of the preceding claims, wherein the processing of step d) is performed iteratively with step c) until a suitably accurate location is obtained.

5. A method according to any one of the preceding claims, wherein the said seismic data streams are collected over a period of five minutes or less.

6. A method according to any one of the preceding claims, wherein said source of substantially continuous noise is a drill bit in operation within said formation.

7. A method of obtaining an image of a subsurface formation of the Earth, the method comprising: using the method of any one of the preceding claims to obtain seismic data simulating data that would be obtained from an impulsive source located substantially at the continuous noise source location; and performing an inversion or seismic imaging process on the data to obtain said image.

8. A non-transitory computer-readable medium storing a program including instructions that, when executed by a processor, cause the processor to implement the method of any one of the preceding claims.