BACKGROUND

Field of the invention

[0001] The present invention concerns determining the position of a drill bit in an underground formation during drilling.

Prior and related art

[0002] In the present application, a formation comprises several layers of different materials, e.g. solid rocks, sands and salt deposits. The layers are folded and displaced along faults and fractures, and in some regions form structures of commercial interest, hereinafter reservoirs. For example, traps or pockets in the formation may contain hydrocarbons or hot water for the energy industry, and aquifers may be suitable for depositing flue gas and/or C0 ₂ . Regardless of application, the reservoirs are accessed through drilled wells.

[0003] Before drilling, a geophysical model of the formation is built using information from various sources. An important source is seismic data acquisition and subsequent inversion. However, inversion does not produce a unique solution - several possible formations can explain the observed data. Hence, an analyst must select the most likely interpretation based on the geology of the area, experience etc. Still, there will be uncertainties that can only be resolved during drilling.

[0004] Drilling a borehole in these applications involves a rotating drill bit on the end of a hollow drill string. The drill string may comprise threaded joints about 9 m (30 feet) each or 'coiled tubing', i.e. a continuous tube with smaller diameter. A drilling mud injected through the drill string cools the drill bit and removes cuttings through an annulus formed between the drill string and the formation. Decades ago, a rotating jointed string with a fixed drill bit was the only option. This method is still widely used for vertical and near vertical wells.

[0005] Today, horizontal wells are typically drilled with the use of a downhole assembly (DHA), which typically comprises a mud motor driven by the drilling mud to rotate the drill bit. The DHA also typically comprise sensors and a controller to adjust the drill bit's angle of attack relative to the DHA housing. This permits geosteering, also called geonavigation. At least 25 companies provide geosteering software and services.

[0006] Geosteering and related applications such as logging- while-drilling (LWD) may comprise resistivity and/or gamma-ray measurements, analysis of cuttings at the surface etc. These measurements are not part of the invention and need no further explanation herein. [0007] Transmitting data to the surface and/or control signals to the DHA is a challenge. Drill strings with embedded copper or fibre lines are expensive and impractical, especially at threaded connections. Electrical signals conveyed directly through the drill string is severely distorted by capacitances and reflections at every threaded joint, so this alternative is also impractical. Accordingly, mud-pulsing is a widely used alternative to transmit information. However, the transmission rate of mud-pulsing is limited to about 8 b/s or 1 B/s, so much processing is performed in the DHA. This requires accurate data about the formation, which may not be available in a complex geological structure with multiple faults and fractures. Thus, an improved transmission rate would improve flexibility and accuracy in geosteering.

[0008] Eidsvik & Hokstad (2006) [1] and associated patent application WO2006/106337 ('337) propose using a seafloor array of seismic receivers for locating a drill bit, in particular by a hyperbolic and a Bayesian approximation. The latter is linearized in a Markov-model and implemented with a Kalman filter (KF). The article [1] also mentions a non-linear filter (Extended KF - EKF), and '337 contains approximate matrices. Embodiments with the drill bit or DHA as source/receiver for complementary receivers and sources on a surface or seafloor are disclosed.

[0009] It is generally acknowledged that EKF can be difficult to tune, mainly because EKF depend on partial derivatives (in a Jacobian), which are approximated by finite differences. Hence, a derivative free unscented KF (UKF), is usually preferred over EKF for estimating the state of a dynamic or non-linear system. In particular, a square-root variety, SR-UKF, is preferred because it guarantees positive definite matrices and eliminates numerical instabilities in very small matrix elements.

[0010] From a signal processing perspective, detecting a signal in a noisy environment is crucial to an improved transmission rate through a drill string and to improve the method of Eidsvik & Hokstad. A few decades ago, advanced sampling was limited to Gibbs sampling and/or the Metropolis-Hastings algorithm. Similarly, noise reduction was limited to filtering, possibly in a transform domain such as FK (2D Fourier transform from time-space TX), tau-p or a wavelet domain. NO20170488 entitled "Method for denoising seismic data from a seafloor array" and assigned to the present assignee provides an example of filtering that exploits the tau-p transform's ability to distinguish acoustic sources.

[0011] Notable contributors to advances in signal and image processing since the 1980s include Ingrid Daubechies, David Donoho, Emmanuel Candes and their co-workers. For example, the JPEG2000 standard uses Daubechies wavelets to transform an image and store the most significant wavelet coefficients. Even with a compression rate 1:43, an inverse transform (synthesis) recover the original image with barely visible degradation. Donoho's work includes development of effective methods in multiscale geometric analysis. In the early 2000s, Donoho and Candes published a series of papers on certain wavelets called curvelets and associated curvelet transforms useful for seismic applications and image denoising.

[0012] The current field of compressed sensing (CS aka sparse sampling) may be traced back to a seminal paper by Candes et al. (2004) [2], which disclosed statistical conditions for recovery of a signal from highly incomplete frequency data, i.e. a sparse representation without an explicit synthesis formula such as an inverse wavelet transform. A comprehensive list of Candes' papers is available at https://statweb.stanford.edu/~candes/publications.html. A google search reveals contributions from Daubechies, Donoho and others to CS.

[0013] In general, advances in electronics, computer technology, information technology, statistics etc. provide an increasingly large toolbox for building applications. For example, computer networks and statistical inference with a dash of Bayesian statistics enable generally known applications such as Google and Facebook to provide a user with customized advertisements and/or search results in near real time based on the user's search history. In addition to these 'big data' applications, machine-learning algorithms using support vectors, forests and neural networks use statistical techniques. Essentially, statistical techniques may extract a small amount of information from a large 'noisy' dataset based on similarities in the required information. In contrast, an explicit search in a structured database or unstructured data requires explicit a priori information, i.e. known or assumed similarities in the required information.

[0014] In the present context, lifting algorithms enable wavelet constructions without explicit expressions for the wavelet, and sparse sampling may replace a need for a priori noise reduction in a full dataset. Specifically, increased randomness and sparsity (up to a limit) increases the probability of recovery. See, for example, Dragotti (tutorial 2015) [3] for examples and further references.

[0015] Riel et al. (2014) [4] describes efficient transient detection using linear sparse estimation techniques, in particular to detect onset and duration of slips of a few mm in geodetic time series over several years in arrays of GPS-stations. The onsets occur randomly, and are not modelled by a Gaussian distribution. Rather, Bayesian sampling provides uncertainties for transient amplitudes and durations. The research article also describes spatial weighting filters for spatially coherent signals.

[0016] US 4460059 (Katz) discloses a method and system for seismic continuous bit position, in which a signal from a rotating drill bit is recorded by several geophones. The geophones are 'focused' on a 3D grid near the drill bit. 'Focusing' means taking account of direction and polarisation of incoming waves, and thereby automatically reduce noise from the drilling rig and other external sources. Traces from the geophones are time shifted and summed or multiplied to cancel incoherent noise and enhance a coherent signal from the drill bit. The resulting signal is then compared to a database of precomputed wave patterns resulting from sources in different cells in the 3D-grid. The drill bit position corresponds to the precomputed pattern with greatest coherence to the observed signal.

[0017] A common feature of prior art methods is that they essentially precompute a wave pattern from a cell in a spatial 3D grid and compare the result to a signal in time and space. Thus, an accurate bit position depends on an accurate velocity model. In other words, if the velocity model is inaccurate, then the assumed bit position corresponding to the maximum coherence signal will be inaccurate. Katz proposes using a measured length of the drill string to improve the estimate of drill bit position.

[0018] A main objective of the present invention is to provide a method for improved drill bit positioning while retaining benefits from prior art.

SUMMARY OF THE INVENTION

[0019] This is achieved by a method according to claim 1. Further features and benefits appear from the dependent claims.

[0020] More particularly, the invention concerns a method for determining a position of a drill bit in a formation during drilling. The method comprises the steps of:

a) For at least two wave modes (P ,S), recording a separate observed pattern of

observations from a seismic array on a land surface or a seafloor;

b) Computing a predicted pattern for each wave mode (P, S);

c) Comparing the observed and predicted patterns for each wave mode (P, S) to obtain first estimates of the drill bit position;

d) Determining an improved estimate from the first estimates.

[0021] The improved estimate is due to different impulse responses for different wave types.

[0022] It should be noted that the drill bit is positioned in real time / near real time. This also gives the application distinctive character and is an important element for the value of the invention.

[0023] Further, the drill bit application can be combined with real-time detection of fracturing of formations (caused by a drilling operation). An example where the combination may be important is in connection with the loss of drilling fluid into a formation. If one measures the drill bit position and at the same time manage to locate where a loss of drilling fluid in the borehole takes place (along with conventional methods for measuring drill bit depth and monitoring return of circulated drilling fluid / drilling fluid pressure) one will be able to make faster and/or better decisions, i.e. to reduce risk.

[0024] In a first embodiment, recording the separate observed pattern comprises the steps of:

Focusing on a small region near the drill bit and obtain a set of traces;

Time shifting the traces; and

On the time shifted traces, using a technique selected from a group consisting of summation, multiplication and autocorrelation.

[0025] Apart from autocorrelation, this corresponds to the method of Katz applied to signals for each wave type. Recall that focusing implies noise reduction. In particular, this applies to cultural noise and/or standing waves in a water column above a seafloor.

[0026] In the first embodiment, computing and comparing the predicted patterns preferably comprises the steps of:

- From each cell in a set of cells near the drill bit, computing a response at a set of points {xi} in Euclidean space R ³ on the land surface or seafloor;

Storing the computed responses as a set of predicted patterns;

Correlating the separate observed pattern with the each predicted pattern and selecting a best correlated pattern; and

- Using the cell position corresponding to the best correlated pattern as first estimate.

[0027] These steps are similar to those of Katz applied to each wave type. The patterns (and accordingly the seismic array) are specified in 3D rather than on a plane as in Katz. The depth dimension has been found to improve the results.

[0028] In a second embodiment, computing and comparing is performed by means of a sigma point method. As generally known, sigma point methods is a subset of linear regression based on the observation that it may be easier to estimate the first few moments, at least mean and variance, of a probability density function (pdf) than to estimate parameters directly. The sigma point method may be an unscented Kalman filter.

[0029] The first and second embodiments may further comprise the step of determining a characteristic frequency spectrum for the drill bit signal. The characteristic frequency spectrum is preferably determined in a downhole assembly near the drill bit and conveyed to a drilling rig through a drillstring.

[0030] Conveying data through the drillstring is optional, and may involves orthogonal frequency division multiplexing (OFDM) similar to techniques used in cellular networks. In particular, the adapted OFDM algorithm comprises the step of dividing the frequency band used for transmission through the drillstring (not the frequency data of the drill bit) into several orthogonal subchannels (C1-C8) and maintaining a symbol in a subchannel (C1-C8) during a timeslot (S 1-S3) sufficiently long to permit summation of several samples. As usual, summation cancels noise (incoherent signals) and enhances a coherent signal.

[0031] All embodiments may benefit from well-established multiscale wavelet analysis wherein observed and predicted patterns are represented by dyadic wavelets. Specifically, the patterns probably has a sparse representation in a wavelet domain. Dyadic wavelets have scales 2 ^; (j an integer), and the finest scale is usually selected at a noise level. Each wavelet coefficient has information on a parameter and its Fourier equivalent, that is (x, k) and/or (t, f). The discrete wavelet transform (DWT) may be regarded as multiple matching filters applied to parameters and their Fourier equivalents at once. For N data points, the DWT completes in time O(N), whereas the Fast Fourier Transform completes in time O(NlogN). Further details are available in textbooks and online.

[0032] Some embodiments may benefit from compressed sensing techniques, in particular recovering a signal from highly incomplete frequency information. Note that older techniques such as minimizing an h norm, e.g. basis pursuit denoising, are computationally demanding.

[0033] The calculated drill bit position can be fed back to a drilling system for adjustment of drilling parameters in order to correct for any deviations from a planned drilling path or to adjust for any obstacles that are detected. The feedback closes the loop in a closed-loop system as known from the art of control systems.

[0034] Denoising of the data by use of various filtering techniques may be performed prior to calculation of drill bit position in order to increase the signal-to-noise-ratio of the signal.

[0035] The method used in the invention include the use of a network of sensors on the seafloor directly connected to a data acquisition and analysis system on the drilling rig. The layout of the network is carefully planned to optimize drill bit detection capabilities and drill bit detection position accuracies. The signal from the drill bit illuminates the formation and provides finer details than, for example, a velocity model obtained from seismic surveying. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] The invention will be explained with reference to exemplary embodiments and the accompanying drawings, in which

Fig. 1 (prior art) illustrates communication over several subchannels using phase shift keying (PSK); Fig. 2 (prior art) illustrates an observed noisy signal in two adjacent timeslots;

Fig. 3 illustrates mean frequencies estimated at the DHA during drilling;

Fig. 4 illustrates steps in a preliminary version of the invention;

Fig. 5 illustrates the position of an oilrig in relation to a selected area and a cell structure; Fig. 6 illustrates a chunk of a formation with orthotropic properties;

Fig. 7 illustrates an estimated drilling path; and

Fig. 8 illustrates the main steps of the invention.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT

[0037] The drawings are schematic and not necessarily to scale. For ease of understanding, numerous details known to the skilled person are omitted from the drawings and following description.

[0038] Fig. 1 (prior art) illustrates an orthogonal frequency multiplexing (OFDM) scheme, e.g. as used in mobile communication. In the present example, a frequency band is divided into eight subchannels C1-C8, and a harmonic or sinusoid signal may have one of four phases P1-P4 during a timeslot Si, each lasting long enough to detect the phase P1-P4 reliably, e.g. -Is. Fig. 1 shows three timeslots S 1-S3. Noise in one channel does not necessarily affect the signal in adjacent channels.

[0039] Each phase P1-P4 represent a 2-bit combination, namely 00, 01, 10 and 11. This is an example of phase shift keying (PSK), in particular 4-PSK. Current mobile communication use 16-PSK and 32-PSK where the 16=2 ⁴ or 32=2 ⁵ phases represent 4, respectively 5, bits.

[0040] Assuming that each timeslot S 1-S3 lasts one second, the scheme in Fig 1 transmits 2 bits times 8 channels = 16 bits per second (b/s). This number includes a payload and parity bits for error correction. Forward error correction codes are well known, and include

Hamming [7, 4] (1950), Reed-Solomon (1960), BCH (1959 - 1960), LDPC (1960) and Turbo (early 1990s). The fundamental patent for turbo coding, US 5,446,747, expired in 2013. Error correction and other aspects of coding theory are outside the scope of the present invention. Here, we observe that the effective transmission rate exceeds 8 b/s whenever the number of parity bits is less than the number of payload bits. Moreover, the number of channels need not be limited to 8, and 8-PSK may replace 4-PSK.

[0041] Applied to the task of transmitting signals between a DHA and equipment on the surface, acoustic signals through the drillstring is a potentially suitable channel as it avoids the impracticalities of separate communication lines and the distortions associated with electrical signals through a jointed string and acoustic signals through a highly non-linear formation. Acoustic signals through the drillstring should be compared to the currently used mud-pulsing.

[0042] The scheme in Fig. 1 depends on the ability to separate subchannels C1-C8. This is conventionally done by quadrature filters, which incidentally led to the development of wavelet transforms and lifting algorithms. To produce the signal, assume one 'hammer' per channel pinging on the drillstring. Disregarding the actual Bessel-signals produced by a hammer, the hammering produces an envelope signal that may be approximated by a sinusoid. To reduce the effect of reflections at the threaded joints, suitable carrier frequencies are in the range 500-600 Hz. These carrier frequencies correspond to wavelengths greater than the length of joints, i.e. > ~9m (30 feet), and depend on the speed of sound -6000 m/s in (stainless) steel. In other words, the wavelength of the envelope signal makes the threaded joints 'invisible'.

[0043] Considering the length of the timeslots S 1-S3 about 1 s, a sampling frequency of 20 Hz permits addition of about 10 samples. The summation enhances a coherent signal and removes incoherent (random) noise, which contributes negatively equally well as positively. The addition of several signals also reduces the need for large amplitudes in the transmitted signal. Thus, depending on the sampling frequency, piezoelectric devices with smaller amplitudes may replace 'hammers', e.g. driven by induction that provide larger amplitudes.

[0044] The scheme in Fig. 1 also depends on the ability to separate phases P1-P4. Fig. 2 illustrates a transient (p(t) shifting from phase P2 to PI at a known point in time, nT, corresponding to a boundary between adjacent timeslots Si, S(i+1) in Fig. 1. A shaded area illustrates one standard deviation from the transient due to ambient Gaussian (white) noise. Sparse sampling may recover the transient signal (p(t) sufficiently well to determine whether the phase in one timeslot is P2 or PI. Noise reduction before signal processing is not strictly required. Similar statistical techniques form the basis of support vectors in machine learning. Applied to Fig.2, such an algorithm would determine whether an observation falls into 'class' P2 or Pl.

[0045] Fig. 3 illustrates mean frequencies f(t) estimated at the DHA during drilling. The observed signal is noisy, and the noise is well approximated by a Gaussian distribution. The dip in frequencies at tl corresponds to the drill bit entering a fault. The raise at t2 is caused by the drill bit entering a harder rock, which is indicated by higher frequencies. The events at tl and t2 cannot be approximated by a simple global distribution, and they do not appear at fixed intervals such as the phase transitions in Fig. 2. However, tl and t2 can be predicted from a geophysical model of the formation, and they can be used to update the model. For example, a fracture may have a different width than expected, and some fractures may be discovered. Such information is crucial to avoid accidents due to leakage from an oil or gas reservoir.

[0046] Although current (signal) processors in the DHA may provide the most significant components in a frequency spectrum, the times tl and t2 etc., there is a limit to the transmission capacity to the surface. Perhaps more important, there is a need for a feedback loop through the formation in order to determine the position of the drill bit and to verify the geophysical model in an iterative process.

[0047] This objective can be realised by precomputing how a wave with a characteristic frequency, e.g. fl in Fig. 3, would propagate through the formation and be received at an array on the surface or the seafloor. More specifically, the wave may be an S-wave as estimated at the DHA, and the computing might yield an expected pattern, i.e. a prediction of what the array should measure based on the assumed bit-position.

[0048] Fig. 4 illustrates a method 100 according to the invention in greater detail.

[0049] Step 110, recording a pattern, involves collecting data from a seismic array on a land surface or a seafloor and arrange the data in a suitable observed pattern. This step may include step 111, denoising the collected data as an explicit step before further processing. Step 111 is optional because denoising may be implicit in CS-techniques.

[0050] Step 120 involves estimating a drill bit position. The estimation may, for example, involve a Kalman filter, in which the a posteriori estimate is an a priori estimate updated with a weighted sum of an observation and a prediction. Step 121 illustrates that the possible positions are limited by physical constraints such as last position, drilling speed and/or inclination of the DHA.

[0051] Step 130, computing a predicted pattern from a drill bit-position, involves propagating a wave from an assumed drill bit position through the geophysical model using known methods. Variational methods may speed up the computation.

[0052] Preferably, step 130 includes a step 131: Selecting a characteristic frequency. The characteristic frequency should be 'typical' for the rock being drilled such as fl in Fig. 3.

[0053] In step 140, comparing observed and predicted patterns, a simple comparison involves subtracting one dataset from the other. Alternatively, any known pattern recognition technique may be employed. Franek et al. (2007) [5] describes a fast algorithm and contains references to other known algorithms.

[0054] Test 150 determines whether the predicted pattern matches the observed pattern with sufficient precision. If not, the process loops back 151 to step 120 for a new estimate of the drill bit-position. This process is repeated until the patterns match. Then, the process loops 152 to step 110 to collect new observations from the array.

[0055] Fig. 5 illustrates the position of an oilrig 200 in relation to a selected area 202 and for providing a focus area divided into a cellstructure around the planned drilling wellpath. A drill bit 204 is connected to the oilrig 200 through a drillstring 201. The drill bit 204 is emitting signals 205 - which is noise arising from the drilling. The signals 205 are picked up by a seafloor sensor network/array 206 placed at the seafloor 203.

[0056] A focus area and a cell structure is selected. Then a curve for theoretical arrival time for events for each cell to each receiver is calculated. A measured arrival time is recorded and the maximum coherence between this value and the precalculated theoretical value is found.

[0057] As noted in the introduction, estimating the bit position from arrival times depends on a velocity model. In some applications such as geosteering, the most relevant position is relative to the formation. This also applies to a planned drill path 207. In these applications, the velocity model need not be very accurate or complicated as long as the same velocity model is used to estimate the drill bit position and locate the geobodies, for example the facies illustrated by different hatchings and shadings in Fig. 5.

[0058] We want to illuminate the formation with the signal from the drill bit, e.g. in order to discover fractures and drilling hazards on a finer scale than that of a velocity model resulting from a seismic survey. Moreover, we have external datasets that do not depend on a specific velocity model, for example magnetic and gravitational data, rock mechanical data and attenuation as function of frequency and P- and S-wave speeds, i.e. Q(f; Vp, Vs). The list is not complete. In order to improve our geomodel, we want to correlate datasets that depend on a velocity model with datasets that do not. For this, we want to estimate an 'absolute' drill bit position, e.g. in conventional Earth bound (N, E, z) coordinates, as opposed to the relative position used for geosteering.

[0059] The absolute drill bit position should be estimated from on-site measurements of particle motions, in particular by 4C sensors, e.g. 210a and b in Fig. 5, with known impulse responses over a known frequency band, e.g. in the range 1 - 100 Hz. Since the attenuation increases with frequency, the upper limit of a practical frequency band decreases with depth. Similarly, the upper relevant frequency decreases with distance. However, at long offsets the drill bit 204 is viewed from slightly different angles by different sensors 210 a, b. Only longwave/ low frequency information is relevant in these traces, similar to a set of low- resolution images. The information may still be interpolated on a high -resolution grid and then deblurred or deconvoluted to obtain a crisp image. Similar superresolution techniques are well-known from image processing, SAR etc.

[0060] For the present application, we assume a non-linear model f(x) of a region with different impulse responses for different wave typess, for simplicity P- and S waves. A key intuition is that S-waves may improve resolution because their phase velocities are about half the phase velocities of P-waves, and also that their attenuation are different from those of P- waves. For example, a region with multiple liquid filled fractures will attenuate S-waves much more than P-waves because S-waves do not propagate through fluids. Opposite, solid rocks are expected to attenuate S-waves less due to the lower phase velocities. As a result, a geomodel with several acoustic contrasts, that is, rock/rock boundaries, fractures and faults, will be non-linear and have different Green's functions (boundary conditions) for P-waves and S-waves. The different Green's functions or impulse responses further depend on temporal and spatial frequencies. For example, small fractures are 'invisible' to long wavelengths and short wave vectors. This is similar to the previous example with acoustic waves along a drillstring.

[0061] Thus, we propose to expand the method of Katz mentioned in the introduction with a secondary comparison of P-wave patterns and S-wave patterns at selected locations for selected temporal and spatial frequencies.

[0062] Before describing alternative signal representations and methods for secondary comparisons, we will briefly describe a family of forward models.

[0063] Since time is not a limiting factor for precomputing patterns, the geomodel producing patterns or velocities may be anisotropic and in general be relatively complex. For example, one may limit the assumption to linear elasticity, which is reasonable for the small amplitudes considered here. The preferred linear elastic model is σ = c-ε where σ is the stress tensor, ε is the strain tensor and c is the elastic stiffness tensor. In orthonormal coordinates and noting symmetries in c, the stress-strain relation can be expressed by 6-element vectors σ and ε and a 6 x 6 matrix c.

[0064] Fig. 6 shows a chunk 300 of a seismic facies with several dips or rotations relative to conventional Earth bound coordinates (N, E, z). The chunk 300 has local orthonormal coordinates (el, e2, e3), and we assume that an uplift in the el-e3 plane has increased the local stress field such that there is an increased risk for fracturing in or near the e2-e3 plane. Similar increases in the stress field due to an uplift in the e2-e3 plane increases the risk for fracturing in or near the el-e3 plane. The uplifts may of course be caused by different events on different timescales, for example tectonic movement and melting of the latest arctic icecap. [0065] From a geophysical point of view, the chunk 300 essentially has two planes of symmetry, and the material within may be considered orthotropic. Then, the stiffness tensor c may be further simplified to a block matrix of 3 x 3 submatrices Cij, where C12 = C21 = O3 (null matrices) and C22 is a diagonal matrix. For isotropic media, the non-zero elements in C22 are equal, for orthotropic materials they are not. In isotropic media, the elements in c depend on two elastic moduli, e.g. λ and μ, and yield one P-wave velocity for longitudinal waves and one for S-waves regardless of propagation direction and polarity. In anisotropic media, the phase velocities of S-waves and 'pseudo S-waves' are different in different directions, and the group and phase velocities have different directions. See, for example. Crampin (1981) [7]. Either way, large chunks and a large number of zeros in each chunk-matrix simplify computation. Moreover, translation, rotation and scaling (affine transformations) with orthonormal bases are easily accomplished by matrix multiplication of augmented 4 x 4 matrices. Further details are available in textbooks and online.

[0066] Models such as the above may be built using a large number of parameters, including lithology, density, elastic moduli, pore pressure and formation geometry. Moreover, well logs, core samples and known rock properties update previous estimates of the geomodel, e.g. from inversion of a seismic survey, whenever additional data become available.

[0067] Further assumptions give simpler models. For example, Thomsen (1986) [6] assumes weak anisotropy and derives the well-known anisotropy parameter δ. Eidsvik & Hokstad [1] modify the phase velocities with δ in several of their examples.

[0068] Next, we select an appropriate non-linear forward model f(x,...) and compute a number of patterns M _g (-) (M for measurement, g for geophone) starting at points x _s (s for source) corresponding to small cells in the underground. Each small cell may, for example, be a cube with sides ~1 - 10 m in (N, E, z) around x _s . The computation from each x _s ends in values at several points x _g , and both x _s and x _g are point in in Euclidean space R ³ . We compute different patterns for different wave types, for simplicity P-waves and S-waves. A real embodiment may account for other wave types such as head waves and surface waves.

[0069] Measurements or observations for each wave type P, S, ... are obtained as in Katz, i.e. by focusing on small cells (here the set {x _s }) to reduce noise from the drill rig etc., then time- shifting and summing (or averaging or multiplying or autocorrelating ... ) traces.

[0070] A first embodiment has two stages. In the first stage, we may store the precomputed patterns in a database and compare them to a pattern observed on-site in real time. For each wave type, the precomputed pattern with best match or correlation to the observed pattern is a first estimate for the drill bit position, x _s . Next, we determine a new and better estimate from the first estimates. The second stage on may involve s weighted sum, for example such that the S-based estimate is attributed more or less weight relative to the P-based estimate depending on the conditions around the drill bit.

[0071] In an alternative embodiment, the first stage involves a sigma point filter as briefly described previously. Specifically, we consider an unscented Kalman filter a good choice. The second stage of determining a best estimate from first P and B-estimates is equal in all embodiments.

[0072] Fig. 7 illustrates an estimated drill path as obtained from P or S-wave based estimates with a UKF. A main benefit over extended Kalman filters is that the UKFs need no estimates for derivatives (a Jacobian), which may be hard to estimate by finite differences in a highly non-linear model. The UKF starts with assumptions of a probability distribution (pdf) and estimates its first statistical moments. The curve 231 is generated from means (first order moment) and the tube 232 is generated from second order moments (variance). UKF algorithms may converge faster if higher order moments can be assumed or inferred. For example, spherical symmetry of the drill bit signal may lead to an assumed third order moment (skew) equal to zero. The fourth order moment (kurtosis) indicates outliers, and is known to be 3.0 for purely Gaussian distributions.

[0073] Fig. 8 illustrates the main steps of a method 300 according to the invention and summarises major options.

[0074] In step 310, separate patterns for several wave types, for simplicity P and S, are recorded by time shifting and summing traces. Real embodiments may take advantage of any wave type, including head waves and surface waves.

[0075] Step 311 illustrates optional denoising, and usually includes focusing the seismic array toward a small region near the drill bit.

[0076] In step 320, patterns for P and S-waves are computed by means of a forward model of arbitrary complexity, e.g. an elastic medium + assumption of orthotropy. The physical constrains in step 321 include known drilling parameters such as heading and inclination of the DHA. The drill bit does not jump around arbitrarily.

[0077] Step 330 ends when first estimates are obtained from P and S-wave based data. Step 331 means that the first estimates can be enhanced by frequency information. For example, the temporal frequency spectrum emitted from the drill bit may contain a dominant frequency that depends on rotation and lithology. Similarly, high frequencies are attenuated in the formation and can be considered noise at large offsets. As illustrated in Fig. 5, the method may be employed at shallow to medium depths, e.g. less than 1 km below the seafloor 203. The drilling depth is, of course, important to determine suitable cut-off frequencies.

[0078] Steps 320 and 330 may I) involve correlating observed data with a set of computed data, or alternatively, II) employ unscented Kalman filters or other sigma point methods.

[0079] In step 340, an improved estimate is obtained from the P and S-based estimates. This may involve a weighted sum (step 341).

[0080] While the invention has been described by way of examples, its scope is defined by the attached claims.

References

[1] Jo Eidsvik and Ketil Hokstad: "Positioning drill-bit and look- ahead events using seismic travel time data", Geophysics 71(4) July- August 2006, pp 79-90

[2] Emmanuel J. Candes, Justin Romberg and Terence Tao; "Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information", IEEE Trans. Inform. Theory, 52 (2004) pp 489-509. [3] Pier Luigi Dragotti: "Sparse Signal Processing Part 2: Sparse Sampling", tutorial 2015. www.commsp.ee.ic.ac.uk/~pld/talks/TutorialSparseSP_Part2.pdf

[4] Bryan Riel, Mark Simons, Piush Agram and Zhongwhen Zhan: "Detecting transient signals in geodetic time series using sparse estimation techniques", J. Geophys. Res. Solid Earth, 2014, 119

[5] Frantisek Franek, Christopher G. Jennings and W.F. Smyth: "A simple fast hybrid pattern- matching algorithm", J. Discrete Algorithms, Vol 5 issue 4, December 2007, pp 682-695 [6] Thomsen, L: "Weak elastic anisotropy", Geophysics, vol 51, 1986, pp 1954-1966

[7] Stuart Crampin, "A review of wave motion in anisotropic and cracked elastic-media", Wave Motion 3 (1981), pp 343 - 391.