NO343160B1 - Total compression of multiple echo trains using principal component analysis and independent component analysis - Google Patents

Abstract

Spin echo signals are obtained downhole by NMR. Principal component analysis is used to represent the signals by a weighted combination of the principal components, and these weights are telemetry to the surface. At the surface, the NMR spin echo signals are restored and inverted to determine formation properties. 1

Description

BACKGROUND OF THE INVENTION

1. The scope of the invention

[0001] The present invention generally relates to the determination of geological properties of underground formations using nuclear magnetic resonance (NMR - Nuclear Magnetic Resonance) methods for logging boreholes, in particular to represent NMR echo trains with a limited number of functional parameters and thus enable efficient transmission of echo trains from the wellbore.

2. Description of Related Art

[0002] US 2005/0206378 A1 relates to prediction of rock properties, categorization and recognition from NMR echo trains using linear and non-linear regression. Partial least squares regression relates spin echo signals with samples that have a known parameter, such as bound water, clay bound water, BVI (Bound Volume Irreducible), porosity and effective porosity. The regression defines a prediction model that is validated and can then be applied to spin echo signals of unknown samples to directly provide an estimate of the parameter of interest. The unknown samples may include underground formations, in which an NMR sensor unit is transported in a borehole. NMR methods are among the most effective non-destructive methods for material analysis. When hydrogen nuclei are placed in an applied static magnetic field, a small majority of spinning nuclei (spins) will align with the applied field in the lower energy state since the lower energy state is more stable than the higher energy state. The individual spinning nuclei precess about the axis of the applied static magnetic field vector with a resonant frequency, also known as the Larmor frequency. This frequency is characteristic of a given core and proportional to the applied static magnetic field. An alternating magnetic field at the resonant frequency in the radio frequency (RF) range, applied by a transmitter antenna to a subject or sample in the static magnetic field, transfers nuclear spin to a coherent superposition of the lower energy state and the higher energy state. In this superimposed state, the magnetization of the spin nuclei precesses about the axis of the static magnetic field vector and therefore induces an oscillating voltage in a receiving antenna even after the emitted field is removed, whose amplitude and rate of return depend on the physiochemical properties of the material under investigation. The applied RF field is designed to perturb the thermal equilibrium of the magnetized nuclear spins, and the time dependence of the emitted energy is determined by how this system of spinning nuclei loses coherence and returns to equilibrium magnetization. The regression is characterized by two parameters: T1, the longitudinal or spin/lattice relaxation time; and T2, the transverse or spin/spin relaxation time.

[0003] Measurements of NMR parameters for fluid that fills the pore spaces in underground formations, such as relaxation times for the spinning hydrogen nuclei, diffusion coefficient and/or hydrogen density, form the basis for NMR-based well logging. NMR-based well logging instruments can be used to determine properties of subsurface formations, including the volume fraction of pore space and the volume fraction of mobile fluid that fills the pore spaces in the subsurface formations.

[0004] One fundamental problem associated with NMR logging or MRI (imaging) is the enormous amount of data that must be analyzed. When well logging with cable-guided instruments, the processing capacity down the hole is limited, and there is also the possibility of sending data up the hole for further analysis since all the data is typically sent up a cable with limited bandwidth. In the so-called measurement-under-drilling methods, the problem is even greater as a result of the harsh environment in which a downhole processor must run and of the very limited telemetry capacity; data is typically transferred at a rate of no more than twenty bits per second.

[0005] A second problem associated with NMR logging and MRI relates to analysis of the data. As will be explained below, the problem of data compression and the problem of data analysis are closely related.

[0006] Methods for using NMR measurements to determine the volume fraction of pore space and the volume fraction of mobile fluid are described, for example, in Spin Echo Magnetic Resonance Logging: Porosity and Free Fluid Index Determination, M. N. Miller, et al., Society of Petroleum Engineers paper no.20561, Richardson, TX, 1990. In porous media, there is significant discrepancy between the T1 and T2 relaxation time spectra of the fluid mixture filling the pore spaces. For example, light hydrocarbons and gas may have a T1 relaxation time of a few seconds, while T2 may be a thousand times less than this. This phenomenon is a consequence of diffusion effects in internal and external static magnetic field gradients. Internal magnetic field gradients arise as a result of differences in magnetic susceptibility between the matrix in the form of formation rocks and the fluid that fills the pores.

[0007] Since oil is found in porous rock formations, the relationship between porous rocks and the fluids that fill the pore spaces in them is extremely complicated and difficult to model. Nuclear magnetic resonance is sensitive to important petrophysical parameters, but is unable to determine these complex relationships. Oil and water generally occur together in reservoir rocks. Since most reservoir rocks are hydrophilic, droplets of oil sit in the center of pores and are unaffected by the pore surface. The water/oil interface will not normally affect relaxation, and therefore the relaxation rate of oil is essentially proportional to its viscosity. However, the oil itself is a very complex mixture of hydrocarbons that can be considered a wide range of relaxation times. In a simple case of pure fluid in a single pore, there are two diffusion regimes that determine the relaxation rate. Rocks normally have a very wide distribution of pore size and fluid properties. It is therefore not surprising that the magnetization reversal for fluid in rock formations is not exponential. The most widespread method for analyzing relaxation data is to calculate a spectrum of relaxation times. The CPMG (Carr-Purcell-Meiboom-Gill) pulse sequence is used to determine the transverse magnetization reversal. The non-exponential magnetization decays are fitted to the multiexponential form:

where M(t) represents the spin echo amplitudes, uniformly distributed in time, and T2 are predetermined time constants, uniformly distributed on a logarithmic scale, typically between 0.3 ms and 3000 ms. The set m is found using a regularized non-linear least square error method. The function m(T2i), traditionally called a T2 distribution, usually maps linearly to a volume-weighted distribution of pore sizes.

[0008] Calibration of this imaging is discussed in several publications. Known methods seek a solution to the problem of mathematical modeling of the received echo signals using several techniques, extensive use of non-linear regression analysis of the measurement signal and non-linear least square error fitting routines. Other known methods include a number of different signal modeling methods, such as polynomial rooting, singular value decomposition (SVD - Singular Value Decomposition) and various further developments of these to find a better approximation of the received signal. A problem with the signal compression in the known technique is that information is lost.

[0009] Inversion methods discussed in the prior art are generally computationally heavy and still end up with a large number of parameters that have to be sent over and over. In particular, no simple methods have been proposed to take advantage of prior knowledge of the structure of the investigated material and the signal-to-noise (SNR) ratio of the received echo signal. Nor have effective solutions been proposed to combine advanced mathematical models with simple signal processing algorithms to improve the accuracy and numerical stability of the parameter estimates. Finally, existing solutions require the use of significant computing power, which makes the practical application of these methods inefficient, and often impossible to implement in real-time applications.

[0010] US patent application 11/845983 to Thern et al. (publication US 2008/0036457) discloses a method which comprises transporting a nuclear magnetic resonance (NMR) measuring device into a borehole, using the NMR measuring device to produce a signal indicative of the property of the subsurface formation, using a predetermined matrix to estimating from the signal a parametric representation of nuclear spin relaxation expressed by at least one basis function, telemetering the parametric representation to a location on the surface and, on the surface, using the telemetered parametric representation to estimate the property of the subsurface formation. The signal can be a spin echo signal and representation of the relaxation of nuclear spin can comprise a distribution of transverse relaxation time (T2). The at least one basis function can be a Gaussian function, and the parametric representation can include mean value/expectation, standard deviation and amplitude of the Gaussian function. Defining the predetermined matrix can be done by performing a regression analysis of artificially generated NMR signals and/or NMR signals measured on samples with known properties. The dependent variable in the regression analysis can be a spin echo signal. The regression analysis can be a partial least square error, a principal component regression, an inverse least square error, a ridge regression, a neural network, a neural network based partial least square error regression and/or a locally weighted regression. The property determined may be BVI (Bound Volume Irreducible), effective porosity, bound water, clay bound water, total porosity, permeability and/or pore size distribution.

[0011] A potential disadvantage of Thern's method is a lack of adaptability; the number of Gaussian functions used to describe the T2 distribution and the matrix is predefined and is not necessarily equally suitable for all types of subsoil formations and all types of pulse sequences used when acquiring the data. These problems are solved by the present invention.

SUMMARY OF THE INVENTION

[0012] The main features of the present invention appear from the independent patent claims. Further features of the invention are indicated in the independent claims. One embodiment of the invention relates to a method for determining a property of an underground formation. The method includes transporting a nuclear magnetic resonance (NMR) measurement device into a borehole, using the NMR measurement device to acquire a signal indicative of the property of the subsurface formation, representing the NMR signals using a set of eigenfunctions, and telemetry a representation of The NMR signals as a combination of the eigenfunctions of a site on the surface.

[0013] Another embodiment of the invention relates to an apparatus for determining a property of an underground formation. The apparatus comprises a nuclear magnetic resonance (NMR) measuring apparatus arranged to be carried into a borehole and obtain a signal indicative of the property of the subsurface formation. The apparatus also comprises a processor down the wellbore arranged to represent the NMR signals using a set of eigenfunctions and telemeter a representation of the NMR signals as a combination of the eigenfunctions to a location on the surface.

[0014] Another embodiment of the invention relates to a computer-readable medium that can be accessed by a processor. The computer-readable medium comprises instructions that enable the processor to represent at least one signal, which is obtained by an NMR measurement apparatus in a borehole and which is representative of a property of a subsurface formation, with a set of eigenfunctions; and telemetry a representation of the at least one signal as a combination of the eigenfunctions of a location on the surface.

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] The present invention will be most easily understood by referring to the attached figures, where like reference numbers refer to like elements and where:

Figure 1 shows a measuring-while-drilling tool suitable for use with the present invention;

Figure 2 (prior art) shows a sensor section of a measurement-while-drilling device suitable for use with the present invention;

Figure 3 shows examples of principal component signals recovered from a sequence of echo trains;

Figure 4 shows (top): an example of echo train, and (bottom) comparison of an NMR-based T2 with a reconstructed distribution;

Figure 5 shows overall compression of data and a train sequence (trainlet); and Figure 6 is a flow diagram showing several details of the implementation of the invention.

DETAILED DESCRIPTION OF THE INVENTION

[0016] Figure 1 shows a schematic diagram of a drilling system 10 with a drill string 20 carrying a drilling unit 90 (also referred to as a bottom hole unit or "BHA") that is carried in a "wellbore" or "borehole" 26 to drill the wellbore . The drilling system 10 comprises a traditional derrick 11 set up on a floor 12 which supports a rotary table 14 which is rotated by a driving force such as an electric motor (not shown) at a desired rotational speed. The drill string 20 comprises a pipe string such as a drill pipe 22 or a coiled pipe that goes from the surface into the borehole 26. The drill string 20 is driven into the wellbore 26 when a drill pipe 22 is used as a pipe string. In coiled tubing applications, a tubing injector, such as an injector (not shown), is used to feed the tubing string from a source, such as a drum (not shown), to the wellbore 26. The drill bit 50 attached to the end of the drill string breaks up the geological formations as it is rotated to drill the borehole 26. If a drill pipe 22 is used, the drill string 20 is connected to a hoist 30 via a rotary pipe joint 21, a swivel 28 and a line 29 through a pulley 23. During drilling operations, the hoist 30 is activated to control the bit pressure, which is an important parameter that affects the drilling speed. The operation of the hoist is well known to the person skilled in the art and is thus not described in detail here. For the purpose of this description, it is necessary to know the axial speed (drilling speed or ROP - Rate of Penetration) of the downhole unit. Depth information and ROP can be communicated downhole from a location on the surface. Alternatively, the method described in US patent 6769 497 to Dubinsky et al., which is assigned to the same as this application and which is incorporated herein as a reference in its entirety, can be used. Dubinsky's method uses axial accelerometers to determine ROP. During drilling operations, a suitable drilling fluid 31 from a mud tank (source) 32 is circulated under pressure through a channel in the drill string 20 by a mud pump 34. The drilling fluid is fed from the mud pump 34 into the drill string 20 through a desurger (not shown), a fluid pipe 38 and a rotary joint 21. The drilling fluid 31 is led out at the bottom of the drill hole 51 through an opening in the drill bit 50. The drilling fluid 31 circulates uphole through the annulus 27 between the drill string 20 and the drill hole 26, and returns to the mud tank 32 via a return pipe 35. The drilling fluid serves to lubricate the drill bit 50 and carry drilling chips or cuttings away from the drill bit 50. A sensor S1, typically arranged in the pipe 38, provides information on the fluid flow rate. A torque sensor S2 on the surface and a sensor S3 connected to the drill string 20 respectively provide information about the torque on and the rotation speed of the drill string. In addition, a sensor (not shown) connected to the line 29 is used to determine the hook load from the drill string 20.

[0017] In one embodiment of the invention, the drill bit 50 is rotated by only rotating the drill pipe 22. In another embodiment of the invention, a wellbore motor 55 (mud motor) is arranged in the drilling unit 90 to rotate the drill bit 50, and the drill pipe 22 is usually rotated to increase the rotational force, if necessary, and to change the drilling direction.

[0018] In the embodiment in Figure 1, the mud motor 55 is connected to the drill bit 50 via a drive shaft (not shown) arranged in a bearing unit 57. The mud motor rotates the drill bit 50 when the drilling fluid 31 passes through the mud motor 55 under pressure. The bearing unit 57 supports the radial and axial forces on the drill bit. A stabilizer 58 connected to the bearing unit 57 serves as a centering unit for the lower part of the mud motor unit.

[0019] In one embodiment of the invention, a drill sensor module 59 is placed near the drill bit 50. The drill sensor module contains sensors, circuits as well as software and algorithms for processing the dynamic drilling parameters. Such parameters typically include drill bit jumping, jerking of the drilling unit, backward rotation, torque, impact, borehole and annulus pressure, acceleration measurements and other measurements of the condition of the drill bit. A suitable telemetry or communication section 72, which for example uses two-way telemetry, is also provided as illustrated in the drilling unit 90. The drilling sensor module processes the information from the sensors and sends it to the control unit 40 on the surface via the telemetry system 72.

[0020] The communication section 72, a power unit 78 and an MWD tool 79 are all connected along the drill string 20. Bending sections, for example, are used to connect the MWD tool 79 in the drilling unit 90. Such sections and tools form the bottom hole unit 90 between the drill string 20 and the drill bit 50. The drilling unit 90 makes various measurements, including pulsed nuclear magnetic resonance measurements, while the borehole 26 is being drilled. The communications section 72 acquires the signals and measurements and transmits the signals, for example using two-way telemetry, for processing on the surface. Alternatively, the signals can be processed using a downhole processor in the drilling unit 90.

[0021] The control unit or processor 40 on the surface also receives signals from other sensors and devices down the wellbore, and signals from the sensors S1-S3 and other sensors used in the system 10, and processes these signals according to programmed instructions provided to the surface control unit 40. The surface control unit 40 displays desired drilling parameters and other information on a display device 42 which is used by an operator to control the drilling operations. The surface control unit 40 typically comprises a computer or a microprocessor-based processing system, memory for storing programs or models and data, a recorder for recording data and other peripheral equipment. The control unit 40 is typically designed to activate alarms 44 when certain unsafe or undesirable operating conditions occur.

[0022] A suitable device for use of the present invention is described in US patent 6,215,304 to Slade, which is incorporated herein by reference in its entirety. It should be noted that the device shown by Slade is only intended as an example, and the method of the present invention can be used with many other NMR logging devices, and can be used for cabled applications as well as MWD applications.

[0023] Now referring to Figure 2, the tool has a drill bit 107 at one end, a sensor section 102 behind the drill bit and electronics 101. The sensor section 102 comprises a magnetic field generating unit to generate a B0 magnetic field (which is substantially time constant over the duration of a measurement) and an RF system to transmit and receive magnetic RF pulses and echoes. The magnetic field generating unit comprises a pair of axially spread main magnets 103, 104 with opposite pole orientations (ie with equal magnetic poles facing each other) and three ferrite elements 109, 110 axially arranged between the magnets 103, 104. The ferrite elements are made of "soft" ferrite, which differs from "hard" ferrite by the shape of the BH curve which affects both intrinsic coercivity (Hj, the intersection with the H axis) and initial permeability (μi, the gradient of the BH curve in the unmagnetized case). Soft ferrite's μi values typically range from 10 to 10000, while hard ferrite has a μi value of about 1. Soft ferrite therefore has a high initial permeability (typically higher than 10, preferably higher than 1000). The RF system comprises a set of RF transmitter antenna and RF receiver antenna coil windings 105 arranged as a center "field forming" solenoid array 113 and a pair of outer "switch control" solenoid arrays 114.

[0024] The tool has a mud pipe 160 with a free center bore 106 and a number of exit openings 161-164 to lead drilling mud to the drill bit 107, and the main body of the tool is formed by a weight pipe 108. Drilling mud is pumped down the mud pipe 160 by a pump 121 and returns around the tool, and the entire tool is rotated by a drive device 120. Coiled tubing or a drill string may be used to connect the drive device to the wellbore assembly.

[0025] The weight tube 108 has a recess 170 for RF transmitter antenna and RF receiver antenna coil windings 105. Gaps in the pockets between the soft ferrite elements are filled with electrically insulating material 131, 135 (eg: ceramic or high temperature plastic), and The RF coils 113, 114 are then wrapped around the soft ferrite elements 109, 110. The assembly of soft ferrite 109, 110 and RF coils 113, 114 is pressure impregnated with a suitable low viscosity high temperature epoxy resin (not shown) to strengthen the system against the impact of vibration, seal against drilling fluid at well pressure and reduce the likelihood of magnetoacoustic oscillations. The RF coils 113, 114 are then covered with wear plates 111, typically of ceramic or other durable electrically insulating material, to protect them from the rock fragments flowing upwards past the tool in the drilling mud.

[0026] As a result of the opposite magnetic pole arrangement, Slade's device has an axisymmetric magnetic field and examination area 112 which is unaffected by tool rotation. The use of ferrite results in a survey area that is close to the borehole. This is not a major problem on an MWD tool because there is little intrusion into the formation of drilling fluids in the borehole before the logging operation. The survey area is within a shell with a radial thickness of approximately 20 mm and an axial length of approximately 50 mm. The gradient within the survey area is less than 2.7 G/cm. It should be noted that these values are for Slade's device, and as indicated above, the method according to the present invention can also be used with other suitable NMR devices.

[0027] The method according to the present invention is based on a representation of the collected echo train in the underground formation as a weighted combination of principal components derived during data collection. This enables compression of the data; instead of a thousand samples being required to produce a single echo train, typically 10 principal components are transferred for each echo train.

The principal components are derived downhole and can be sent uphole when previously derived principal components do not provide a satisfactory reconstruction of echogenic downholes. At the surface, the received data (which may include process noise) is decompressed and inverted to produce a T2 distribution. We will now briefly describe the principal component analysis (PCA – Principal Component Analysis) method for compressing and decompressing the data.

[0028] We represent a sequence of N echo trains, each M echo long, with the matrix:

The echo trains are typically 1000 echoes long. The mean value of the jth echo is named as:

We then define the covariance matrix of the data as:

there:

The covariance matrix C is decomposed into its eigenvalues and eigenvectors

C = V Λ V<-1>(5),

where V is a matrix whose columns are the eigenvectors of C and Λ is the diagonal matrix of eigenvalues:

with

With this arrangement of the eigenvalues, the eigenvectors in V are the principal components.

[0029] Representation of the echo train data is formed by the transformation

and the inverse transformation

can be used to restore the data. Data compression is achieved by truncating the matrix V to the first k rows corresponding to the dominant eigenvalues in Eq. (7). Table 1 shows an example of dominant eigenvalues for an example sequence of echo trains.

Table 1: Variance distribution

[0030] Figure 3 shows the eigenvectors corresponding to the seven largest eigenvalues for the echo train used in the derivation of table 1. They are sorted according to the absolute value of the eigenvalues 301, 303, 305, 307, 309, 311, 313.

[0031] In the lower part of Figure 4, the curve 451 shows the original T2 distribution corresponding to the noise-free echo train. Using curve 451, an initial echo train (contained in 401) is generated. Noise is added to the original echo train to create a denoised echo train, also contained in 401 in the upper portion of Figure 4. Curve 453 is the result of inverting the denoised echo train. Comparison of 451 and 453 shows that even with a low level of noise, the inversion deviates from the correct result. Compression and decompression of the noisy echo train produces a result that is still contained in 401. As will be explained below, the results of the decompression of the eigenvectors are formed and are not multi-exponential. Inverting the results of the decompression gives the curve 455. Such a result is not satisfactory because the inversion algorithm attempts to fit a multiexponential to a curve that is no longer multiexponential. To avoid problems resulting from the decompression results not being exponential, a small amount of noise is added to the decompressed data before the inversion is performed. The results of inversion of these noisy, decompressed data are shown by 457. One can see that there is a good match between the original T2 distribution and the results of inversion using PCA.

[0032] The PCA method can also be used to compress two or more echo trains in a single operation. The upper part of Figure 5 shows a chaining of two echo trains. The early part 501 was obtained with a short waiting time so that a measurement of rapidly relaxed components in the T2 spectrum was obtained, while the later part 503 is a long echo train intended to restore the slower components in the T2 spectrum. The lower part of Figure 5 shows little difference between the actual T2 distribution and the results from applying PCA to the concatenated echotags.

[0033] It should be noted that the PCA method also works for T1 data. It has been found that joint compression of T1 and T2 data is only satisfactory for a fixed value of T1/T2. As this ratio varies down the wellbore, joint compression of T1 and T2 data has limited value.

[0034] Figure 6 shows a flow diagram summarizing the implementation of the method, including further details of the adaptation method described above. NMR data is collected downhole 607. A truncated eigenvector matrix is applied 609 to the collected echo train X. The compressed data is telemetered to the surface 613, and a reconstruction of the echo train is inverted to determine the T2 distribution 615.

[0035] In one embodiment of the invention, the eigenvector matrix is generated on the surface and the truncated matrix is loaded into memory in the downhole processor. We form artificial single exponential data. This is a pure exponential function, 1000 values uniformly distributed with TE = 0.6 ms, with given T2. We generate such a sequence of artificial data values (a single exponential) for each value of T2 to be considered, e.g. for 0.3ms, 0.35ms, ..., 3000ms. This gives 64 sequences of data values for single exponentials. We note that any imaginable measured echo train can be decomposed into a set of these sequences of data values. We therefore apply the PCA method to these data to find out more about their statistical properties. We want to perform a coordinate system rotation (in a 1000-dimensional vector space), and we now apply PCA to find out which basis vectors must be used to express any multiexponential in the new coordinate system as cheaply as possible. Note that although the original data matrix consisted of exponentials, the eigenvectors, after PCA, are not necessarily exponentials. After the PCA method is performed, the matrix is truncated to the number of rows corresponding to the dominant eigenvalues. See equation (9).

[0036] In an alternative embodiment of the invention, the PCA method is carried out down the wellbore. This requires enormous computational power and should be done sparingly in cases where it has been determined that a previously determined set of eigenvectors does not provide a satisfactory representation of the data. This may be relevant if, for example, parameters for the pulse sequence change, or if there is a major change in the lithology and/or fluid content of the formation.

[0037] The reconstruction of properties of interest can cover T2 distribution, volume measure, permeability, echo train and other rock and fluid properties that are based on NMR data. It should also be noted that the method as such is naturally not limited to downhole applications. As suggested in Hamdan, BVI, effective porosity, bound water, clay bound water and total porosity are among the formation properties that can be determined. From the T2 relaxation spectrum, using an inversion method, it is possible to estimate the pore size distribution. Use of a pore-scale geometric model used when inverting NMR spectra is described, for example, in US patent 7363161 to Georgi et al., which is assigned to the same as the present invention. Determination of permeability is discussed in US 6686 736 to Schoen et al., which is transferred to the same as the present invention.

[0038] In an alternative embodiment of the invention, instead of principal component regression or principal component analysis (PCA), a method called independent component analysis (ICA) can be used. In PCA, the basis vectors are found by solving the algebraic eigenvalue problem

R<T>(XX<T>)R= Λ (10) where X is a data matrix whose columns are training samples (with the mean values removed), R is a matrix of eigenvectors and Λ is the associated diagonal matrix of eigenvalues. With such a representation, the projection of data, Cn= Rn<T>X, from the original p-dimensional space to a subspace spanned by n principal eigenvectors is optimal with respect to mean square errors. Specifically, the projection of Cntilbake to the p-dimensional space has minimum reconstruction error. If n is large enough to include all the eigenvectors with an associated eigenvalue different from zero, the projection is lossless. The goal in the PCA method is to minimize the reconstruction error from compressed data.

[0039] In the ICA method, on the other hand, the goal is to minimize the statistical dependence between the basis vectors. Mathematically, this can be written as WX<T>= U, where the ICA method looks for a linear transformation W that minimizes the statistical dependence between the rows in U, given a training set X (as before). Unlike PCA, the basis vectors in ICA are neither orthogonal nor ranked in order. In addition, there is no analytical expression for determining W. Instead, one must use iterative algorithms. See Baek et al., PCA vs. ICA: A comparison on the FERET data set.

[0040] As stated by Baek, global features are better represented by PCA while local structure is better represented by ICA. Based on a comparison of PCA and ICA, Baek concluded that PCA performed better for face recognition problems (which are holistic in nature). Baek further hypothesized that ICA will be able to provide better results for the evaluation of localized recognition tasks, such as recognition of facial expressions.

[0041] First and foremost, NMR measurements give an indication of pore size distribution in a subsurface formation. Second, they give an indication of fluid types. By its nature, the primary pore size distribution in sedimentary rocks reflects the sedimentation energy, which is episodic. Consequently, a significant amount of local structure can be expected in the pore size distribution. To put it another way, for example, one would expect a high correlation between occurrences of pore sizes of 1 μm and 1.01 μm; this would suggest a local structure in the T2 distribution and the T1 distribution. In addition, the presence of heavy oil in a formation will also indicate a local structure in the relaxation time distributions - once heavy oil is formed, it cannot be converted to light oil.

[0042] We will now deal with the execution of ICA and differences from PCA. NMR relaxation in fluids in rocks exhibits a multiexponential behavior, which can be expressed in a discrete model as follows:

Assuming that T2j= 0.2...8192 using an increment of 2<(1/4)>, T2 will have length 64. This can be expressed in matrix notation by sampling t with TE= 0.6 μs and 1000 samples as :

where Ajer is proportional to the proton content in pores having a relaxation time of T2j, E(t) is the resulting echo train in continuous time and E is the discretized version of E(t).

We first map all possible single-exponential regression constant echo trains into a matrix F. Then, using the ICA method, we decompose the matrix F into 2 matrices.

F is a matrix that spans all single-component regressions in the echo train space.

S is a matrix of independent components (latent variables) of the associated type of data collection (created from ICA) using the algorithm fastICA, available with MATLAB, of the matrix F). M is mixture matrix. Both M and S must be estimated. Once S and M are found, the data compression method is as follows:

Comp is called the compression vector. Equation 4 can then be written as:

Now multiply both sides from the right by the inverse of S => S<-1>.

E1x1000x S<-1>1000x64= Comp1x64x S64x1000x S<-1>1000x64

which gives:

E1x1000x S<-1>1000x64= Comp1x64<(16)>

However, the eigenvalue analysis of the covariance of F tells us that after component 6 there will be almost zero percent variance left, as can be seen from table 1.

[0043] Consequently, equation 16 can be reduced to:

E1x1000x S<-1>1000x6= Comp1x6<(17)>

Eq. 17 is used in the downhole tool for compression of echo trains. Equation 15 becomes equation

18 and is used in the surface system to decompress the mud pulse transmitted data:

E1x1000= Comp1x6x S6x1000<(18)>

Eq. 17 tells us that by producing the inverse of a reduced form of the matrix S, we can compress an echo train of length 1000 (and if we have an echo train of length N, we must generate the matrix S of size 6xN) into a 1x6 matrix . Furthermore, equation 18 tells us that we can recreate the echo train by using the same model (independent components) and the corresponding compression.

[0044] The PCA algorithm differs from the ICA algorithm only in the decomposition method. In principal component analysis, we decompose the matrix F into 2 matrices. F64x1000= Scores64x64x Loads64x1000<(19)>

Here, F is a matrix spanning all single-component regressions, Loads is a matrix of eigenvectors from the current type of data collection (created from principal component decomposition of the matrix F) and Scores are the eigenvalues of the matrix F.

It should be noted that Scores form an orthogonal set (Scoresi<T>Scoresj= 0 for i≠j) and Loads form an orthonormal set (Loadsj<T>Loadsj= 0 for i≠j and = 1 for i = j ) => Loads<T>= Loads<-1>. Scores for Scoresji T is a linear combination of F defined by Loads; which means that Scoresier is the projection of F on Loadsi. By inserting the value of F from equation 10 into equation 2, we get:

E1x1000= A1x64x Scores64x64x Loads64x1000<(20)>

Let

Comp1x64= A1x64x Scores64x64<(20a)>

Comp is what we call a compression vector. Equation 20a can then be written as:

E1x1000= Comp1x64x Loads64x1000<(21)>

Now multiplying from the right by the inverse of Loads => Loads<-1>, and applying the fact that Loads<-1>= Loads<T>, we get:

E1x1000x Loads<T>1000x64= Comp1x64x Loads64x1000x Loads<T>64x1000

which gives:

E1x1000x Loads<T>1000x64= Comp1x64<(22)>

[0045] Equation 22 tells us that we can compress the entire echo train from 1000 points to 64 points without loss of information. PCA analysis tells us that after component 5 there will be almost zero percent variance left, as can be seen from the following table:

Table 2: Analysis of the variance distribution in each PCA component

Consequently, equation 22 can be reduced to:

E1x1000x Loads<T>1000x5= Comp1x5<(23)>and equation 21 becomes

E1x1000= Comp1x5x Loads5x1000<(24)>

[0046] Equation 23 tells us that by producing a reduced form of the matrix Loads we can compress an echo train of length 1000 (and if we have an echo train of length N, we must generate the matrix Loads of size 5xN) into a 1x5- matrix. Furthermore, equation 24 tells us that we can reproduce the echo train by applying the same model and the corresponding compression.

[0047] In brief, the ICA algorithm can be used instead of PCA.

[0048] Implied in managing and processing the data is the use of a computer program incorporated on a suitable machine-readable medium which enables the processor to effect said management and processing. The machine-readable medium may include ROM, EPROM, EAROM, flash memory and optical disc storage.

Claims (16)

PATENT CLAIMS 1. Procedure for determining a property of a subsoil formation, where the procedure includes: transporting a nuclear magnetic resonance (NMR) measuring device into a borehole; using the NMR measuring apparatus to obtain at least two signals indicative of the property of the subsurface formation; representing a concatenation of the at least two signals using a set of eigenfunctions; and telemetry a representation of the at least two signals as a combination of the eigenfunctions of a location on the surface. 2. Method according to claim 1, further comprising deriving the eigenfunctions using at least one of: (i) a principal component analysis and (ii) an independent component analysis. 3. Method according to claim 1, further comprising deriving the eigenfunctions at one of: (i) a location on the surface and (ii) a location down in the wellbore. 4. Method according to claim 1, where the at least two signals are selected from the group consisting of: (i) a spin-echo signal that is representative of the distribution of transverse relaxation time (T2), and (ii) a signal that is representative of the distribution of longitudinal relaxation time (T1). 5. Method according to claim 1, further comprising using the telemetry transmitted representation to produce an estimate of the at least two signals and to estimate the property of the underground formation. 6. Method according to claim 1, where the property is selected from the group consisting of: (i) BVI (Bound Volume Irreducible), (ii) effective porosity, (iii) bound water, (iv) clay bound water, (v) total porosity, ( vi) a permeability and (vii) a pore size distribution. 7. Method according to claim 1, further comprising transporting the NMR measurement device into the borehole of a bottomhole unit using a drill pipe. 8. Apparatus designed to determine a property of an underground formation, where the apparatus comprises: a nuclear magnetic resonance (NMR) measuring apparatus adapted to be carried into a borehole and produce at least two signals indicative of the property of the subsurface formation; and at least one processor adapted to: (A) representing the at least two signals using a set of eigenfunctions; and (B) telemetry a representation of the at least two signals as a combination of the eigenfunctions of a location on the surface. 9. Apparatus according to claim 8, further comprising a drill pipe arranged to carry the NMR measurement apparatus into the borehole of a downhole unit. 10. Apparatus according to claim 8, where the at least one processor is further arranged to derive the eigenfunctions using at least one of: (i) a principal component analysis and (ii) an independent component analysis. 11. Apparatus according to claim 8, where the at least one processor is further arranged to derive the eigenfunctions at one of: (i) a location on the surface and (ii) a location down in the wellbore. 12. Apparatus according to claim 8, wherein the NMR measuring apparatus is further arranged to produce the at least two signals selected from the group consisting of: (i) a spin-echo signal which is representative of the distribution of transverse relaxation time (T2), and ( ii) a signal representative of the distribution of longitudinal relaxation time (T1). 13. Apparatus according to claim 8, further comprising a processor on the surface arranged to use the telemetry transmitted representation to produce an estimate of the at least two signals and to estimate the property of the subsurface formation. 14. Apparatus according to claim 8, wherein the NMR measuring apparatus is further arranged to produce the at least two signals indicating a property selected from the group consisting of: (i) BVI (Bound Volume Irreducible), (ii) effective porosity, (iii) bound water, (iv) clay bound water, (v) total porosity, (vi) a permeability and (vii) a pore size distribution. 15. Computer-readable medium comprising instructions which, when read by a processor, enable the processor to carry out the method, the method comprising: representing, by a set of eigenfunctions, a concatenation of at least two signals representative of a property of a subsurface formation obtained by an NMR measurement apparatus in a borehole; and telemetry a representation of the concatenation of the at least two signals as a combination of the eigenfunctions of a location on the surface. 16. Computer readable medium according to claim 15, further comprising at least one of: (i) a ROM, (ii) an EPROM, (iii) an EAROM, (iv) a flash memory and (v) an optical disc storage.