# Patent

## Title
Estimations of nuclear magnetic resonance measurement distributions

## Abstract

A nuclear magnetic resonance (NMR) related distribution is estimated that is consistent with NMR measurements and uses linear functionals directly estimated from the measurement indications by integral transforms as constraints in a cost function. The cost function includes indications of the measurement data, Laplace transform elements and the constraints, and a distribution estimation is made by minimizing the cost function. The distribution estimation may be used to find parameters of the sample. Where the sample is a rock or a formation, the parameters may include parameters such as rock permeability and/or hydrocarbon viscosity, bound and free fluid volumes, among others. The parameters may be used in models, equations, or otherwise to act on the sample, such as in recovering hydrocarbons from the formation.

## Background
This application is related to and hereby incorporates by reference herein in its entirety co-owned U.S. Pat. No. 6,462,542 to Venkataramanan et al., entitled “Nuclear Magnetic Resonance Measurements and Methods of Analyzing Nuclear Magnetic Resonance Data,” and U.S. Ser. No. 13/333,232 to Venkataramanan et al., entitled “Estimation of Petrophysical and Fluid Properties Using Integral Transforms in Nuclear Magnetic Resonance,” filed on Dec. 21, 2011.
BACKGROUND
The statements made herein merely provide information related to the present disclosure and may not constitute prior art, and may describe some embodiments illustrating the invention.
Measured nuclear magnetic resonance (NMR) data resulting from a multi-component sample can be denoted by G(t) which represents a multi-exponential decay, with time constants T2 and amplitudes ƒ(T2)

G(t)=∫0∞Pτ(T2)e−t/T2ƒ(T2)dT2  (1)

where the function Pτ(T2) is referred to as the polarization factor and depends on the pulse sequence of the NMR equipment used to probe and measure the sample.

In certain samples containing hydrocarbons and water, e.g., geological formations, the transverse relaxation time T2 is the characteristic de-phasing time for protons in hydrocarbons or water present in pores of a rock or in the bulk fluid. As aforementioned, the function Pτ(T2) of equation (1) depends on the pulse sequence of the NMR equipment used to probe and measure the sample. For example,




  
    
      
        
          
            P
            τ
          
          ⁡
          
            (
            
              T
              2
            
            )
          
        
        =
        
          {
          
            
              
                1
              
              
                
                  CPMG
                  ⁢
                  
                    
                  
                  ⁢
                  pulse
                  ⁢
                  
                    
                  
                  ⁢
                  sequence
                  ⁢
                  
                    
                  
                  ⁢
                  with
                  ⁢
                  
                    
                  
                  ⁢
                  full
                  ⁢
                  
                    
                  
                  ⁢
                  polarization
                
              
            
            
              
                
                  1
                  -
                  
                    2
                    ⁢
                    
                      e
                      
                        
                          -
                          τ
                        
                        /
                        
                          T
                          2
                        
                      
                    
                  
                
              
              
                
                  
                    inversion
                    ⁢
                    
                      
                    
                    ⁢
                    recovery
                  
                  -
                  
                    CPMG
                    ⁢
                    
                      
                    
                    ⁢
                    pulse
                    ⁢
                    
                      
                    
                    ⁢
                    sequence
                  
                
              
            
            
              
                
                  1
                  -
                  
                    e
                    
                      
                        -
                        τ
                      
                      /
                      
                        T
                        2
                      
                    
                  
                
              
              
                
                  
                    saturation
                    ⁢
                    
                      
                    
                    ⁢
                    recovery
                  
                  -
                  
                    CPMG
                    ⁢
                    
                      
                    
                    ⁢
                    pulse
                    ⁢
                    
                      
                    
                    ⁢
                    sequence
                  
                
              
            
          
        
      
    
    
      
        (
        2
        )
      
    
  




where τ is a function of pre-polarization time and longitudinal relaxation, and CPMG refers to the well-known Carr-Purcell-Meiboom-Gill sequence. In downhole applications, the function Pτ(T2) is a complex function of tool geometry (such as length of magnet and design of RF antenna), operational constraints (such as cable speed) as well as the pulse sequence.

In equation (1) the T2 distribution denoted by ƒ(T2) is estimated from indications of the measured data G(t). The conventional approach to estimating ƒ(T2) utilizes an inverse Laplace transform (ILT), See L. Venkataramanan et al., “Solving Freholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions”, IEEE Transactions on Signal Processing, 50:1017-1026, 2002. More particularly, it is assumed that the indications of data G(t) are measured with additive, white, Gaussian noice. Traditionally, assuming Pτ(T2) is known, an inversion algorithm is used to estimate the distribution of relaxation times ƒ(T2) in equation (1) from the measured data indications G(t). Next, linear functionals of the estimated ƒ(T2) can be used to estimate the petro-physical or fluid properties. For example, the area under the T2 distribution can be interpreted as the porosity of the rock. Often, based on lithology, a threshold T2 is chosen as the cut-off characteristic time separating fluid bound to the pore surface and fluid that is not bound to the pore surface and flow more easily. For example, in sandstones, the area under the T2 distribution corresponding to relaxation times smaller than 33 msec has been empirically related to bound fluid volume. The remaining area under the T2 distribution corresponds to free fluid volume. The mean of the distribution lnT2, is empirically related to either rock permeability and/or to hydrocarbon viscosity. The width of the distribution, σln(T2), provides a qualitative range of the distribution of pore sizes in the rock. Moments of relaxation time or diffusion are often related to rock permeability and/or hydrocarbon viscosity. Similar answer products can also be obtained from linear functionals computed from two-dimensional diffusion-relaxation data or T1-T2 relaxation data.
It is well known in the literature that estimation of ƒ(T2) is an ill-conditioned and non-linear problem. Small changes in the indications of measured data G(t) due to noise can result in widely different ƒ(T2). In theory, there are infinitely many ƒ(T2) that fit the data; i.e., there are non-unique solutions. In the literature, this problem is addressed by introducing a choosing the “smoothest” solution ƒ(T2) that fits the data. More particularly, a regularization functional is introduced, and the solution is estimated through minimization of a cost function Q with respect to the underlying distribution f,

Q=∥G−Lƒ∥2+α∥ƒ∥2  (3)

where G is a vector representing the measured data, L is the matrix of the discretized function Pτ(T2)e−t/T2, and f is the discretized version of the underlying density function ƒ(T2). The first term in the cost function is the least squares error between the data and the fit from the model in equation (1). The second term is referred to as the regularization and incorporates smoothness in the relaxation amplitudes into the problem formulation. The parameter α denotes the compromise between the fit to the data and an a priori expectation of the distribution. In equation (3), α is the weight given to the regularization and can be chosen by a number of different methods.

Previously incorporated U.S. Ser. No. 13/333,232 describes a method based on integral transforms that allows direct computation of a linear functional of ƒ(x) (where the variable x may be the transverse relaxation time (T2), the longitudinal relaxation time (T1), Diffusion (D), any of various two-dimensional data sets (such as D-T2 or D-T1 or T1-T2), multi-dimensional data sets (such as D-T2-T1, D-T2-T1-azimuth, D-T2-T1-different depth of investigations) without involving computation of ƒ(x). The central idea of this approach is to compute the integral transform of G(t) using a function k(t) such that the Laplace transform of k(t) denoted by K(T2) has the desired properties in the x domain (hereinafter for simplicity purposes described with respect to the T2 domain). Let the integral transform of the data indications G(t) be denoted by I{G(t)}=A, such that

I{G(t)}=A=∫0∞k(t)G(t)dt  (4)

From equations (1) and (4),

A=∫0∞K(T2)Pτ(T2)ƒ(T2)dT2  (5)

where the functions k(t) and K(T2) form a Laplace-transform pair, with

K(T2)≡∫0∞k(t)e−t/T2dt  (6)

From the right hand side of equation (5), for a desired linear transformation in the T2 domain, the integral transform approach constructs a function k(t) in the time-domain, so that the scalar product of the function with the measured data provides A, the parameter of interest.
Using this approach, direct computation of linear functionals of the T2 distribution may be made without use of the ILT method. For example, as described in previously incorporated U.S. Ser. No. 13/333,232, the Mellin, Exponential Haar, and Fourier-Mellin transform of magnetization data can be used to estimate moments, tapered areas, and the porosity of the distribution directly from the measured data.
More particularly, consider measured CPMG data G(t), its corresponding and unknown distribution ƒ(T2), and a tapered transition K(T2) for computing the desired tapered area of the T2 distribution. The function k(t) corresponding to this tapered function K(T2) may be shown. This function k(t) can be found using multiple methods: (1) using either an analytical form for k(t) (examples of this is shown in Table A), (2) using a numerical approach, (3) using a method of successive approximations, and/or (4) using convolution analysis. All these methods are described below. Once k(t) is estimated, a scalar product of k(t) with G(t) directly provides the tapered area A. Since this approach does not involve solving for ƒ(T2) and then estimating ∫0∞K(T2)ƒ(T2)dT2, it is more straight-forward and not susceptible to the subjectivity of traditional algorithms that involve inversion of an ill-conditioned and non-linear problem.
As discussed below, given a desired K(T2), the function k(t) can be computed in four ways: often, it can be computed analytically or numerically. It can also be computed using the method of successive approximations to K(T2). An alternate method of computing it involves taking advantage of the convolution-multiplication equivalence between the time and T2 domain. All four methods are described below.
For a desired K(T2), when the function k(t) exists analytically or can be computed numerically, the parameter A is obtained from equation (5). However, the function k(t) may not exist ∀K(T2). When it exists, it may also have infinite energy which can be related to infinite uncertainty in the estimated parameter A leading to instability in computing the parameter. Thus, the integral transform approach can provide insight into what linear functionals of the T2 distribution can be directly applied to the data G(t) and are stable in the context of providing low uncertainty in A.
The uncertainty in A can be quantified as a function of the signal-to-noise ratio (SNR) in the measured data. Let σε denote the standard deviation of the additive white Gaussian noise in the data. Equation (5) can be computed in the discrete-time domain as

A=tEΣn=0Nk(ntE)G(ntE)  (6a)

where tE denotes the sampling time or echo spacing. Therefore,

σA2=σε2tE[Σn=0Nk2(ntE)tE]  (6b)

Equation 6b shows that when the function k(t) is square integrable, i.e, k(t) has finite energy E, where E=∫0∞k2(t)dt, then the uncertainty in A is always finite and directly related to the uncertainty in the measurement.
In the sub-sections below, tables are described of integral transforms developed for different polarization factors PτT2 encountered in oilfield NMR applications.
A. Tapered Areas from Fully Polarized Data
Consider NMR data that have been fully polarized, with Pτ(T2)=1∀T2. In this sub-section, a few integral transforms are described where K(T2) corresponds to tapered and sharp Heaviside functions.
Let Tc denote the T2 relaxation time at which the desired cut-off of the tapered Heaviside function is 0.5. The parameter Tc is user-specified and may come from laboratory study of rock and fluid properties or may correspond to a value of T2 expected to separate two fluids in the T2 domain. For example, in sandstones and carbonates, the area under the T2 distribution corresponding to relaxation times smaller than 33 and 100 msec, respectively, has been empirically related to bound fluid volume. Thus, given a value of Tc, a tapered or sharp Heaviside function K(T2, Tc) is sought such that the tapered area can be computed in the time-domain using the corresponding function k(t, Tc). The integral transform for computing tapered and sharp transitions should satisfy the following properties:
1. The function k(t, Tc) should exist ∀t and K(T2, Tc) should exist ∀T2.
2. Based on the underlying petrophysics, it is desirable that K(T2, Tc) be monotonic between 0 and 1 (on the y-axis), with

K(T2,Tc)|T2−0=0  (6c)

limT2→∞K(T2,Tc)=1  (6d)

K(T2,Tc)|T2−Tc=0.5  (6e)

3. It should be possible to adjust the slope m in the log(T2) space of the transition region, with





  
    
      
        
          m
          ≡
          
            (
            
              
                d
                
                  K
                  ⁡
                  
                    (
                    
                      
                        T
                        2
                      
                      ,
                      
                        T
                        c
                      
                    
                    )
                  
                
              
              
                
                  d
                  log
                
                ⁢
                
                  
                
                ⁢
                
                  T
                  2
                
              
            
            )
          
        
        ⁢
        
          |
          
            
              T
              2
            
            =
            
              T
              c
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          f
        
        )
      
    
  




In most oilfield applications, the desired slope varies from m=0.4 for gradual tapered cut-offs to m=4 for sharp cut-offs.
Using analytical means a set of integral transforms are developed that satisfy these properties, and they are summarized in Table A. For ease of reference, suggested names for the transforms are provided based on the function k(t). The energy for some of the transforms is infinity, implying infinite uncertainty in the estimated area. This energy can be decreased by several methods. One such method involves multiplication of the function k(t) by an exponential decaying signal in the time domain. A second method involves restricting the integral transform to a finite time-period. Both methods decrease the energy considerably while also reducing the slope in the transition region. For e.g., as shown in Table A, the Haar transform (HT) (row 3) has infinite energy. On the other hand, the energy of an exponential Haar transform (EHT) (row 5) is finite.
From equation (7), a desired uncertainty in the estimated area σA can be translated to a desired and finite energy in the function. This finite energy can be achieved by suitable choice of parameters of the transform satisfying both the energy criteria as well as properties 1-3 described above. For example, when the desired energy for the EHT is




  
    2
    
      π
      ⁢
      
        
      
      ⁢
      
        T
        c
      
    
  
  ,




then me parameters C, α and β take values provided in the following table.








TABLE A







Integral transforms that give rise to tapered transitions in the log(T2) domain

















Name of


K(T2, Tc)
Parameters
k(t, Tc)|
E
m
Transform



















  
    2
    π
  
  ⁢
  
    
      tan
      ⁢
      
        
      
    
    
      -
      1
    
  
  ⁢
  
    (
    
      αT
      2
    
    )
  








  α
  =
  
    1
    
      T
      c
    
  








  
    2
    π
  
  ⁢
  
    
      sin
      ⁡
      
        (
        αt
        )
      
    
    t
  








  2
  
    π
    ⁢
    
      
    
    ⁢
    T
  




0.32
Sinc


 






  
    
      T
      2
    
    α
  
  ⁢
  
    tanh
    ⁡
    
      (
      
        α
        
          T
          2
        
      
      )
    
  








  α
  =
  
    
      T
      c
    
    0.52219
  







  
    
      1
      α
    
    ⁢
    
      
        (
        
          -
          1
        
        )
      
      n
    
  




  
    
      2
      ⁢
      n
      ⁢
      
        
      
      ⁢
      α
    
    <
    t
    <
    
      2
      ⁢
      
        (
        
          n
          +
          1
        
        )
      
      ⁢
      α
    
  



∞
0.42
Haar


 






  
    α
    2
  
  
    
      α
      2
    
    +
    
      1
      
        T
        2
        2
      
    
  








  α
  =
  
    1
    
      T
      c
    
  




αsin(αt)
∞
0.5
Sine


 






  
    C
    
      (
      
        
          1
          
            T
            2
          
        
        +
        β
      
      )
    
  
  ⁢
  
    tanh
    ⁡
    
      [
      
        α
        ⁡
        
          (
          
            
              1
              
                T
                2
              
            
            +
            β
          
          )
        
      
      ]
    
  







  
    C
    =
    
      0.7213
      
        T
        c
      
    
  




  
    α
    =
    
      
        (
        1.572
        )
      
      ⁢
      
        T
        c
      
    
  




  
    β
    =
    
      0.4087
      
        T
        c
      
    
  



C(−1)ne−βt 2nα < t < 2(n + 1)α




  2
  
    π
    ⁢
    
      
    
    ⁢
    
      T
      c
    
  




0.35
Exponential Haar


 






  
    
      α
      2
    
    +
    
      β
      2
    
  
  
    
      α
      2
    
    +
    
      
        (
        
          β
          +
          
            1
            
              T
              2
            
          
        
        )
      
      2
    
  







  
    α
    =
    
      
        
          4
          ⁢
          E
          ⁢
          
            
          
          ⁢
          β
        
        -
        
          β
          2
        
      
    
  




  
    β
    =
    
      1
      
        
          T
          c
          2
        
        ⁡
        
          (
          
            
              4
              ⁢
              E
            
            -
            
              2
              
                T
                c
              
            
          
          )
        
      
    
  







  
    
      
        α
        2
      
      +
      
        β
        2
      
    
    α
  
  ⁢
  
    e
    
      
        -
        β
      
      ⁢
      
        
      
      ⁢
      t
    
  
  ⁢
  
    sin
    ⁡
    
      (
      
        α
        ⁢
        
          
        
        ⁢
        t
      
      )
    
  








  2
  
    π
    ⁢
    
      
    
    ⁢
    
      T
      c
    
  




0.3
Exponential Sine


 





  
    
      
        g
        0
      
      +
      
        
          ∑
          
            n
            =
            1
          
          ∞
        
        ⁢
        
          
            a
            n
          
          ⁢
          
            
              g
              n
            
            ⁡
            
              (
              x
              )
            
          
        
      
    
    ,
    
      x
      =
      
        
          T
          c
        
        
          T
          2
        
      
    
  




  
    
      
        g
        n
      
      ⁡
      
        (
        x
        )
      
    
    =
    
      
        
          x
          
            
              2
              ⁢
              n
            
            -
            2
          
        
        -
        
          x
          
            2
            ⁢
            n
          
        
      
      
        
          (
          
            1
            +
            
              x
              2
            
          
          )
        
        
          
            2
            ⁢
            n
          
          -
          1
        
      
    
  



ak computed recursively
See Appendix A
∞
0.5 ≦ mn ≦ ∞ Variable Slope
Series Expansion 









B. Transforms on Imperfectly Polarized Data

Consider imperfectly polarized data with

Pτ(T2)=1−e−τ/T2  (6g)

where





  
    τ
    =
    
      
        T
        w
      
      /
      r
    
  
  ,
  
    r
    =
    
      〈
      
        
          T
          1
        
        
          T
          2
        
      
      〉
    
  




and Tw is the wait time. This polarization factor plays an important role in saturation-recovery-CPMG pulse sequence and in enhanced precision mode (EPM) used in downhole applications. We show below that the integral transform approach that was developed on fully polarized data with Pτ(T2)=1∀T2 can be applied to imperfectly polarized data as well. From eqns. (1) and (6g)





  
    
      
        
          G
          ⁡
          
            (
            t
            )
          
        
        =
        
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                e
                
                  
                    -
                    t
                  
                  /
                  
                    T
                    2
                  
                
              
              ⁢
              
                f
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                d
                
                  T
                  2
                
              
            
          
          -
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                e
                
                  
                    -
                    
                      (
                      
                        t
                        +
                        τ
                      
                      )
                    
                  
                  /
                  
                    T
                    2
                  
                
              
              ⁢
              
                f
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                d
                
                  T
                  2
                
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          h
        
        )
      
    
  




Let the fully polarized data be denoted by M(t), where





  
    
      
        
          M
          ⁡
          
            (
            t
            )
          
        
        =
        
          
            ∫
            0
            ∞
          
          ⁢
          
            
              e
              
                
                  -
                  t
                
                /
                
                  T
                  2
                
              
            
            ⁢
            
              f
              ⁡
              
                (
                
                  T
                  2
                
                )
              
            
            ⁢
            
              
                d
                
                  T
                  2
                
              
              .
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          i
        
        )
      
    
  




Equation (6h) can then be cast as follows

G(t)=M(t)−M(t+τ)  (6j)

For a finite time t we have





  
    
      
        
          
            
              lim
              
                N
                ->
                ∞
              
            
            ⁢
            
              M
              ⁡
              
                (
                
                  t
                  +
                  
                    N
                    ⁢
                    
                      
                    
                    ⁢
                    τ
                  
                
                )
              
            
          
          =
          0.
        
        ⁢
        
          

        
        ⁢
        Therefore
      
    
    
      
        (
        
          6
          ⁢
          l
        
        )
      
    
  
  
    
      
        
          
            ∑
            
              n
              =
              0
            
            
              N
              -
              1
            
          
          ⁢
          
            G
            ⁡
            
              (
              
                t
                +
                
                  n
                  ⁢
                  
                    
                  
                  ⁢
                  τ
                
              
              )
            
          
        
        =
        
          
            M
            ⁡
            
              (
              t
              )
            
          
          -
          
            M
            ⁡
            
              (
              
                t
                -
                
                  N
                  ⁢
                  
                    
                  
                  ⁢
                  τ
                
              
              )
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          m
        
        )
      
    
  




For a finite time t and in the limit N→∞, we get





  
    
      
        
          M
          ⁡
          
            (
            t
            )
          
        
        =
        
          
            ∑
            
              n
              =
              0
            
            ∞
          
          ⁢
          
            
              G
              ⁡
              
                (
                
                  t
                  +
                  
                    n
                    ⁢
                    
                      
                    
                    ⁢
                    τ
                  
                
                )
              
            
            .
          
        
      
    
    
      
        (
        
          6
          ⁢
          n
        
        )
      
    
  




Therefore, if τ is either known or estimated, the fully polarized data M(t) can be reconstructed using equation (6n) from the measured data G(t). Integral transforms can be applied to M(t) to directly estimate linear functionals of ƒ(T2).
Effect of Noise
In practice we have a limited number of noisy samples of the form

{tilde over (G)}(itE)=G(itE)+n(itE), i=1, . . . , N.  (6o)

Using equation (6n) will result on a modified noise





  
    
      
        
          
            n
            s
          
          ⁡
          
            (
            
              it
              E
            
            )
          
        
        =
        
          
            ∑
            
              j
              =
              0
            
            N
          
          ⁢
          
            n
            ⁡
            
              (
              
                
                  it
                  E
                
                +
                jτ
              
              )
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          p
        
        )
      
    
  




with variance

σM2=Nσε2.  (6q)

For noisy data, it may be desirable to perform a denoising procedure before applying equation (6n).

C. Transforms on Data from Logging Tools

The function k(t) can also be found using a combination of numerical and analytical methods as follows. This method is illustrated with an example in logging tools where the polarization factor Pτ(T2) tends to be more complex and depends on a number of parameters including hardware design such as length of permanent magnet, cable speed, etc. For simplicity, we have ignored the subscript τ in the polarization term in this sub-section. In logging applications, we have found that over a range of factors, P(T2) is well represented by




  
    
      
        
          P
          ⁡
          
            (
            
              T
              2
            
            )
          
        
        =
        
          1
          
            
              ∑
              
                k
                =
                0
              
              ∞
            
            ⁢
            
              
                a
                k
              
              ⁢
              
                e
                
                  
                    -
                    bk
                  
                  /
                  
                    T
                    2
                  
                
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          r
        
        )
      
    
  




Equation (6r) is a good fitting function ƒ or a wide range of parameters. For example, the imperfectly polarized data in the last sub-section is obtained from equation (6r) with b=τ and ak=1∀k. Similarly, at a range of cable speeds, the fit from equation (6r) matches P(T2) reasonably well. The polarization factor in logging tools is a complex function of tool geometry, operational constraints such as logging speed and and pulse sequence. In a number of circumstances, at logging speeds varying from 800-2000 ft/hour, the fit (solid line) from equation (6r) fits the polarization factor very well.

The fully polarized data M(t) can be reconstructed from G(t). From eqns. (6i) and (6n), we get




  
    
      
        
          
            
              
                M
                ⁡
                
                  (
                  t
                  )
                
              
              =
              
              ⁢
              
                
                  ∫
                  0
                  ∞
                
                ⁢
                
                  
                    e
                    
                      
                        -
                        t
                      
                      /
                      
                        T
                        2
                      
                    
                  
                  ⁢
                  
                    f
                    ⁡
                    
                      (
                      
                        T
                        2
                      
                      )
                    
                  
                  ⁢
                  
                    
                      P
                      ⁡
                      
                        (
                        
                          T
                          2
                        
                        )
                      
                    
                    
                      P
                      ⁡
                      
                        (
                        
                          T
                          2
                        
                        )
                      
                    
                  
                  ⁢
                  
                    d
                    
                      T
                      2
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  ∫
                  0
                  ∞
                
                ⁢
                
                  
                    e
                    
                      
                        -
                        t
                      
                      /
                      
                        T
                        2
                      
                    
                  
                  ⁢
                  
                    f
                    ⁡
                    
                      (
                      
                        T
                        2
                      
                      )
                    
                  
                  ⁢
                  
                    
                      P
                      ⁡
                      
                        (
                        
                          T
                          2
                        
                        )
                      
                    
                    [
                    
                      
                        ∑
                        k
                      
                      ⁢
                      
                        
                          a
                          k
                        
                        ⁢
                        
                          e
                          
                            
                              -
                              bk
                            
                            /
                            
                              T
                              2
                            
                          
                        
                      
                    
                    ]
                  
                  ⁢
                  
                    d
                    
                      T
                      2
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  ∑
                  k
                
                ⁢
                
                  
                    a
                    k
                  
                  ⁢
                  
                    
                      ∫
                      0
                      ∞
                    
                    ⁢
                    
                      
                        e
                        
                          
                            -
                            
                              (
                              
                                t
                                +
                                bk
                              
                              )
                            
                          
                          /
                          
                            T
                            2
                          
                        
                      
                      ⁢
                      
                        P
                        ⁡
                        
                          (
                          
                            T
                            2
                          
                          )
                        
                      
                      ⁢
                      
                        f
                        ⁡
                        
                          (
                          
                            T
                            2
                          
                          )
                        
                      
                      ⁢
                      
                        d
                        
                          T
                          2
                        
                      
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  ∑
                  k
                
                ⁢
                
                  
                    a
                    k
                  
                  ⁢
                  
                    
                      G
                      ⁡
                      
                        (
                        
                          t
                          +
                          bk
                        
                        )
                      
                    
                    .
                  
                
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          s
        
        )
      
    
  




Integral transforms can be applied to M(t) to directly estimate linear functionals of ƒ(T2).
Method 2: The function k(t) can also be found numerically as follows. For example, it is possible that either Pτ(T2) or K(T2) is not well approximated by a closed form expression or an analytical k(t) does not exist for a specified K(T2). In this case, k(t) can be computed numerically as follows. For example, consider the case where the data are fully polarized with Pτ(T2)=1∀T2. The desired K(T2, Tc) can be shown as a trace. A numerical least squares approximation to k(t, T) can be obtained using singular value decomposition (SVD), with

{tilde over (k)}(t)≈VnΣn−1UnTK(T2)  (6t)

Here matrices U, Σ and V are obtained by SVD of function e−t/T2 and n refers to the number of significant singular values.
In another embodiment of the method, the function k(t) can be found such that its Laplace transform minimizes the error with respect to the desired K(T2, Tc) and has a minimal energy,

mink(t)∥∫0∞k(t)et/T2dt−K(T2,Tc)∥2 

such that

∥k∥2<E 

Method 3: The function k(t) can also be found using the equivalence of the convolution-multiplication operation between the time and T2 domain. This is further described below. We show that the product of two functions in the T2 domain corresponds to convolution in the time-domain. This property implies that the integral transforms described in this memo can also be combined in the time domain to estimate other parameters. For example, the moments of a specified region of the T2 distribution can be computed by using a function computed as the convolution of the Mellin operator and the Exponential Haar transform. Consider two different integral transforms of the measured data, where function k(t) in equation (5) is represented by k1(t) and k2(t) respectively,




  
    
      
        
          A
          1
        
        =
        
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                
                  k
                  1
                
                ⁡
                
                  (
                  t
                  )
                
              
              ⁢
              
                G
                ⁡
                
                  (
                  t
                  )
                
              
              ⁢
              
                d
                t
              
            
          
          =
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                
                  K
                  1
                
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                
                  P
                  τ
                
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                f
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                d
                
                  T
                  2
                
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          u
        
        )
      
    
  
  
    
      
        
          A
          2
        
        =
        
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                
                  k
                  2
                
                ⁡
                
                  (
                  t
                  )
                
              
              ⁢
              
                G
                ⁡
                
                  (
                  t
                  )
                
              
              ⁢
              
                d
                t
              
            
          
          =
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                
                  K
                  2
                
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                
                  P
                  τ
                
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                f
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                
                  d
                  
                    T
                    2
                  
                
                .
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          v
        
        )
      
    
  




Here, the functions K1(T2) and K2(T2) correspond to different linear functionals. Our interest is in evaluation of A3, where




  
    
      
        
          
            A
            3
          
          =
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                
                  K
                  3
                
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                
                  P
                  τ
                
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                f
                ⁡
                
                  (
                  
                    T
                    2
                  
                  )
                
              
              ⁢
              
                d
                
                  T
                  2
                
              
            
          
        
        ,
        
          

        
        ⁢
        and
      
    
    
      
        (
        
          6
          ⁢
          w
        
        )
      
    
  
  
    
      
        
          
            K
            3
          
          ⁡
          
            (
            
              T
              2
            
            )
          
        
        ≡
        
          
            
              K
              1
            
            ⁡
            
              (
              
                T
                2
              
              )
            
          
          ⁢
          
            
              K
              2
            
            ⁡
            
              (
              
                T
                2
              
              )
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          x
        
        )
      
    
  
  
    
      
        =
        
          
            ∫
            0
            ∞
          
          ⁢
          
            
              
                k
                1
              
              ⁡
              
                (
                τ
                )
              
            
            ⁢
            
              e
              
                
                  -
                  τ
                
                
                  T
                  2
                
              
            
            ⁢
            
              d
              τ
            
            ⁢
            
              
                ∫
                0
                ∞
              
              ⁢
              
                
                  
                    k
                    2
                  
                  ⁡
                  
                    (
                    
                      t
                      2
                    
                    )
                  
                
                ⁢
                
                  e
                  
                    
                      -
                      
                        t
                        2
                      
                    
                    
                      T
                      2
                    
                  
                
                ⁢
                
                  d
                  
                    t
                    2
                  
                
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          y
        
        )
      
    
  




From equation (6)





  
    
      
        
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                
                  k
                  3
                
                ⁡
                
                  (
                  t
                  )
                
              
              ⁢
              
                e
                
                  
                    -
                    t
                  
                  /
                  
                    T
                    2
                  
                
              
              ⁢
              
                d
                t
              
            
          
          =
          
            
              ∫
              0
              ∞
            
            ⁢
            
              
                ∫
                0
                ∞
              
              ⁢
              
                
                  
                    k
                    1
                  
                  ⁡
                  
                    (
                    τ
                    )
                  
                
                ⁢
                
                  
                    k
                    2
                  
                  ⁡
                  
                    (
                    
                      t
                      2
                    
                    )
                  
                
                ⁢
                
                  e
                  
                    
                      -
                      
                        (
                        
                          τ
                          +
                          
                            t
                            2
                          
                        
                        )
                      
                    
                    
                      T
                      2
                    
                  
                
                ⁢
                
                  d
                  τ
                
                ⁢
                
                  d
                  
                    t
                    2
                  
                
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              Let
              ⁢
              
                
              
              ⁢
              τ
            
            +
            
              t
              2
            
          
          =
          
            t
            .
            
              
            
            ⁢
            Thus
          
        
      
    
    
      
        (
        
          6
          ⁢
          z
        
        )
      
    
  
  
    
      
        
          
            ∫
            0
            ∞
          
          ⁢
          
            
              
                k
                3
              
              ⁡
              
                (
                t
                )
              
            
            ⁢
            
              e
              
                
                  -
                  t
                
                /
                
                  T
                  2
                
              
            
            ⁢
            
              d
              t
            
          
        
        =
        
          
            ∫
            0
            ∞
          
          ⁢
          
            
              [
              
                
                  ∫
                  0
                  t
                
                ⁢
                
                  
                    
                      k
                      1
                    
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                  ⁢
                  
                    
                      k
                      2
                    
                    ⁡
                    
                      (
                      
                        t
                        -
                        τ
                      
                      )
                    
                  
                  ⁢
                  
                    d
                    τ
                  
                
              
              ]
            
            ⁢
            
              e
              
                
                  -
                  t
                
                /
                
                  T
                  2
                
              
            
            ⁢
            
              d
              t
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          aa
        
        )
      
    
  




Thus, the parameter A3 can be computed as,





  
    
      
        
          A
          3
        
        =
        
          
            ∫
            0
            ∞
          
          ⁢
          
            
              
                k
                3
              
              ⁡
              
                (
                t
                )
              
            
            ⁢
            
              G
              ⁡
              
                (
                t
                )
              
            
            ⁢
            
              d
              t
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          bb
        
        )
      
    
  




where k3(t) is obtained as a convolution of k1(t) and k2(t),





  
    
      
        
          
            k
            3
          
          ⁡
          
            (
            t
            )
          
        
        =
        
          
            ∫
            0
            t
          
          ⁢
          
            
              
                k
                1
              
              ⁡
              
                (
                τ
                )
              
            
            ⁢
            
              
                k
                2
              
              ⁡
              
                (
                
                  t
                  -
                  τ
                
                )
              
            
            ⁢
            
              
                d
                τ
              
              .
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          cc
        
        )
      
    
  




Hence, the product of two functions in the T2 domain corresponds to convolution in the time-domain. This property implies that the integral transforms described in this manuscript can also be combined to estimate other parameters. For example, this property implies that the moments of a specified region of the T2 distribution can be computed by integral transforms of the measured data, using a function obtained as a convolution of the Mellin operator and the Exponential Haar transform.
Method 4: We describe below a method for computing k(t) using method of successive approximations. We illustrate this with an example. Consider K(T2) is an arbitrarily sharp transition in the T2 domain. Let




  x
  =
  
    
      
        T
        c
      
      
        T
        2
      
    
    .
  




Let the Heaviside function H(x) be defined as follows.





  
    
      
        
          H
          ⁡
          
            (
            x
            )
          
        
        =
        
          
            1
            ⁢
            
              
            
            ⁢
            for
            ⁢
            
              
            
            ⁢
            x
          
          <
          1
        
      
    
    
      
        (
        
          6
          ⁢
          ee
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          =
          
            
              0.5
              ⁢
              
                
              
              ⁢
              for
              ⁢
              
                
              
              ⁢
              x
            
            =
            1
          
        
      
    
    
      
        (
        
          6
          ⁢
          ff
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          =
          
            0
            ⁢
            
              
            
            ⁢
            
              otherwise
              .
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          gg
        
        )
      
    
  




In this method, we define g0(x) to be a generating function if it is monotonic and takes values between 0 and 1 and satisfies the following property,




  
    
      
        g
        0
      
      ⁡
      
        (
        x
        )
      
    
    +
    
      
        g
        0
      
      ⁡
      
        (
        
          1
          x
        
        )
      
    
  
  =
  1.




Examples of generating functions that resemble a Heaviside function and satisfy the above property are




  
    
      g
      0
    
    ⁡
    
      (
      x
      )
    
  
  =
  
    
      
        2
        π
      
      ⁢
      a
      ⁢
      
        
      
      ⁢
      
        tan
        ⁡
        
          (
          
            1
            x
          
          )
        
      
      ⁢
      
        
      
      ⁢
      and
      ⁢
      
        
      
      ⁢
      
        
          g
          0
        
        ⁡
        
          (
          x
          )
        
      
    
    =
    
      1
      
        1
        +
        
          x
          2
        
      
    
  




We seek a series of coefficients an and functions g(x), n=1, . . . , ∞ such that





  
    
      
        
          H
          ⁡
          
            (
            x
            )
          
        
        =
        
          
            
              g
              0
            
            ⁡
            
              (
              x
              )
            
          
          +
          
            
              ∑
              
                n
                =
                1
              
              ∞
            
            ⁢
            
              
            
            ⁢
            
              
                a
                n
              
              ⁢
              
                
                  g
                  n
                
                ⁡
                
                  (
                  x
                  )
                
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          ii
        
        )
      
    
  




For n≧1, we seek functions gn(x) such that the functions satisfy the following properties:
1. gn(x) should have unique inverse Laplace transform in closed-form and should exist for all x.
2. gn(x) should be anti-symmetric in log−x, with




  
    
      
        
          
            g
            n
          
          ⁡
          
            (
            x
            )
          
        
        =
        
          -
          
            
              
                g
                n
              
              ⁡
              
                (
                
                  1
                  x
                
                )
              
            
            .
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          jj
        
        )
      
    
  




3. When x is small, gn(x): xn.

4. When x is large, gn(x): x−n.

Properties (1), (3) and (4) are self-explanatory. Property (2) follows from the Heaviside function H(x) and generating function g0(x) satisfying eq. (6hh). At the first iteration, the approximate Heaviside function is

H1(x)=g0(x)+a1g1(x)  (6kk)

Therefore,




  
    
      
        
          
            H
            1
          
          ⁡
          
            (
            
              1
              x
            
            )
          
        
        =
        
          
            
              g
              0
            
            ⁡
            
              (
              
                1
                x
              
              )
            
          
          +
          
            
              a
              1
            
            ⁢
            
              
                g
                1
              
              ⁡
              
                (
                
                  1
                  x
                
                )
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          ll
        
        )
      
    
  




From our construction, H1(x) and g0(x) satisfy eq. (6hh). Therefore,




  
    
      
        
          1
          -
          
            
              H
              1
            
            ⁡
            
              (
              x
              )
            
          
        
        =
        
          1
          -
          
            
              g
              0
            
            ⁡
            
              (
              x
              )
            
          
          +
          
            
              a
              1
            
            ⁢
            
              
                g
                1
              
              ⁡
              
                (
                
                  1
                  x
                
                )
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          mm
        
        )
      
    
  




Thus, from equations (6kk) and (6mm), we have




  
    
      g
      1
    
    ⁡
    
      (
      x
      )
    
  
  =
  
    
      -
      
        
          g
          1
        
        ⁡
        
          (
          
            1
            x
          
          )
        
      
    
    ·
    A
  




similar proof follows for any gn(x), n≧1.

Case 1:





  
    Let
    ⁢
    
      
    
    ⁢
    
      
        g
        0
      
      ⁡
      
        (
        x
        )
      
    
  
  =
  
    
      2
      π
    
    ⁢
    a
    ⁢
    
      
    
    ⁢
    
      
        tan
        ⁡
        
          (
          
            1
            x
          
          )
        
      
      .
    
  




Its Taylor-series expansion is





  
    
      
        
          
            2
            π
          
          ⁢
          a
          ⁢
          
            
          
          ⁢
          
            tan
            ⁡
            
              (
              
                1
                x
              
              )
            
          
        
        =
        
          
            
              H
              ⁡
              
                (
                x
                )
              
            
            -
            
              
                
                  2
                  π
                
                [
                
                  x
                  -
                  
                    
                      x
                      3
                    
                    3
                  
                  +
                  
                    
                      
                        x
                        5
                      
                      5
                    
                    ⁢
                    
                      
                    
                    ⁢
                    ...
                  
                
                ⁢
                
                  
                
                ]
              
              ⁢
              
                
              
              ⁢
              for
              ⁢
              
                
              
              ⁢
              
                
                
                  
                
                
              
            
          
          <
          1.
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          nn
        
        )
      
    
  




If we subtract from g0(x) the terms in the Taylor-series proportional to xn(n≧1), written as a function of gn(x), then, we will obtain a function that converges to 1 for |x|<1 and converges to 0 for |x|>1. Since the Taylor-series expansion has only odd-powers of x, we consider




  
    
      
        
          
            g
            n
          
          ⁡
          
            (
            x
            )
          
        
        =
        
          
            
              -
              
                x
                
                  
                    2
                    ⁢
                    
                      
                    
                    ⁢
                    n
                  
                  -
                  1
                
              
            
            +
            
              x
              
                
                  2
                  ⁢
                  
                    
                  
                  ⁢
                  n
                
                +
                1
              
            
          
          
            
              (
              
                1
                +
                
                  x
                  2
                
              
              )
            
            
              2
              ⁢
              
                
              
              ⁢
              n
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          oo
        
        )
      
    
  




These functions satisfy properties (1)-(4). In addition,





  
    
      
        
          
            ∑
            
              k
              =
              1
            
            ∞
          
          ⁢
          
            
          
          ⁢
          
            
              a
              k
            
            ⁢
            
              
                g
                k
              
              ⁡
              
                (
                x
                )
              
            
          
        
        =
        
          
            ∑
            
              k
              =
              1
            
            ∞
          
          ⁢
          
            
          
          ⁢
          
            
              a
              k
            
            ⁢
            
              
                
                  -
                  
                    x
                    
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        k
                      
                      -
                      1
                    
                  
                
                +
                
                  x
                  
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      k
                    
                    +
                    1
                  
                
              
              
                
                  (
                  
                    1
                    +
                    
                      x
                      2
                    
                  
                  )
                
                
                  2
                  ⁢
                  
                    
                  
                  ⁢
                  k
                
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          pp
        
        )
      
    
  




Using the series expansion for





  1
  
    
      (
      
        1
        +
        
          x
          2
        
      
      )
    
    
      2
      ⁢
      
        
      
      ⁢
      k
    
  




around x=0, we get





  
    
      
        
          
            ∑
            
              k
              =
              1
            
            ∞
          
          ⁢
          
            
          
          ⁢
          
            
              a
              k
            
            ⁢
            
              
                g
                k
              
              ⁡
              
                (
                x
                )
              
            
          
        
        =
        
          
            ∑
            
              k
              =
              1
            
            ∞
          
          ⁢
          
            
          
          ⁢
          
            
              
                a
                k
              
              ⁡
              
                (
                
                  
                    x
                    
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        k
                      
                      +
                      1
                    
                  
                  -
                  
                    x
                    
                      
                        2
                        ⁢
                        k
                      
                      -
                      1
                    
                  
                
                )
              
            
            ⁡
            
              [
              
                
                  ∑
                  
                    m
                    =
                    0
                  
                  ∞
                
                ⁢
                
                  
                
                ⁢
                
                  
                    
                      (
                      
                        -
                        1
                      
                      )
                    
                    m
                  
                  ⁢
                  
                    (
                    
                      
                        
                          
                            
                              2
                              ⁢
                              
                                
                              
                              ⁢
                              k
                            
                            +
                            m
                            -
                            1
                          
                        
                      
                      
                        
                          m
                        
                      
                    
                    )
                  
                  ⁢
                  
                    x
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      m
                    
                  
                
              
              ]
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          qq
        
        )
      
    
  




Matching the coefficients for x2n−1 in eqns. (43) and (46) yields





  
    
      
        
          a
          n
        
        =
        
          
            
              
                (
                
                  -
                  1
                
                )
              
              
                n
                -
                1
              
            
            
              
                2
                ⁢
                
                  
                
                ⁢
                n
              
              -
              1
            
          
          +
          
            
              ∑
              
                k
                =
                1
              
              
                n
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  (
                  
                    -
                    1
                  
                  )
                
                
                  n
                  -
                  k
                
              
              ⁢
              
                
                  
                    a
                    k
                  
                  ⁡
                  
                    [
                    
                      
                        (
                        
                          
                            
                              
                                n
                                +
                                k
                                -
                                2
                              
                            
                          
                          
                            
                              
                                n
                                -
                                k
                                -
                                1
                              
                            
                          
                        
                        )
                      
                      +
                      
                        (
                        
                          
                            
                              
                                n
                                +
                                k
                                -
                                1
                              
                            
                          
                          
                            
                              
                                n
                                -
                                k
                              
                            
                          
                        
                        )
                      
                    
                    ]
                  
                
                .
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          rr
        
        )
      
    
  




The first three coefficients are




  
    
      a
      1
    
    =
    
      2
      π
    
  
  ,
  
    
      a
      2
    
    =
    
      
        
          8
          3
        
        ⁢
        
          a
          1
        
        ⁢
        
          
        
        ⁢
        and
        ⁢
        
          
        
        ⁢
        
          a
          3
        
      
      =
      
        
          
            128
            15
          
          ⁢
          
            
              a
              1
            
            .
            
              
            
            ⁢
            Let
          
          ⁢
          
            
          
          ⁢
          τ
        
        =
        
          
            t
            
              T
              c
            
          
          .
        
      
    
  




The inverse Laplace transforms of the first three terms in the series expansion are




  
    
      
    
    ⁢
    
      
        
          L
          
            -
            1
          
        
        ⁡
        
          [
          
            
              x
              -
              
                x
                3
              
            
            
              
                (
                
                  1
                  +
                  
                    x
                    2
                  
                
                )
              
              2
            
          
          ]
        
      
      =
      
        
          1
          
            T
            c
          
        
        ⁡
        
          [
          
            
              τ
              ⁢
              
                
              
              ⁢
              
                sin
                ⁡
                
                  (
                  τ
                  )
                
              
            
            -
            
              cos
              ⁡
              
                (
                τ
                )
              
            
          
          ]
        
      
    
  




  
    
      
    
    ⁢
    
      
        
          L
          
            -
            1
          
        
        ⁡
        
          [
          
            
              
                x
                3
              
              -
              
                x
                5
              
            
            
              
                (
                
                  1
                  +
                  
                    x
                    2
                  
                
                )
              
              4
            
          
          ]
        
      
      =
      
        
          1
          
            T
            c
          
        
        ⁡
        
          [
          
            
              1
              24
            
            ⁢
            
              (
              
                
                  
                    -
                    6
                  
                  ⁢
                  
                    τ
                    2
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                -
                
                  6
                  ⁢
                  τ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                +
                
                  
                    τ
                    3
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
              
              )
            
          
          ]
        
      
    
  




  
    
      
        L
        
          -
          1
        
      
      ⁡
      
        [
        
          
            
              x
              5
            
            -
            
              x
              7
            
          
          
            
              (
              
                1
                +
                
                  x
                  2
                
              
              )
            
            6
          
        
        ]
      
    
    =
    
      
        
          1
          
            T
            c
          
        
        ⁡
        
          [
          
            
              1
              1920
            
            ⁢
            
              (
              
                
                  30
                  ⁢
                  
                    
                  
                  ⁢
                  
                    τ
                    2
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                -
                
                  15
                  ⁢
                  
                    τ
                    4
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                -
                
                  30
                  ⁢
                  τ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                -
                
                  55
                  ⁢
                  
                    τ
                    3
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                +
                
                  
                    τ
                    5
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
              
              )
            
          
          ]
        
      
      .
    
  



A general expression for the inverse Laplace transform for gn(x) can be obtained from the following:




  
    
      
        
          
            L
            
              -
              1
            
          
          ⁡
          
            [
            
              
                
                  x
                  r
                
                -
                
                  x
                  
                    r
                    +
                    2
                  
                
              
              
                
                  (
                  
                    1
                    +
                    
                      x
                      2
                    
                  
                  )
                
                
                  r
                  +
                  1
                
              
            
            ]
          
        
        =
        
          
            
              τ
              
                r
                -
                1
              
            
            
              
                T
                c
              
              ⁢
              
                Γ
                ⁡
                
                  (
                  r
                  )
                
              
              ⁢
              
                Γ
                ⁡
                
                  (
                  
                    r
                    +
                    2
                  
                  )
                
              
            
          
          ⁡
          
            [
            
              
                
                  τ
                  2
                
                ⁢
                
                  
                    Γ
                    ⁡
                    
                      (
                      r
                      )
                    
                  
                  1
                
                ⁢
                
                  
                    F
                    2
                  
                  ⁡
                  
                    (
                    
                      
                        
                          r
                          +
                          1
                        
                        ;
                        
                          
                            r
                            2
                          
                          +
                          1
                        
                      
                      ,
                      
                        
                          
                            r
                            2
                          
                          +
                          
                            3
                            2
                          
                        
                        ;
                        
                          -
                          
                            
                              τ
                              2
                            
                            4
                          
                        
                      
                    
                    )
                  
                
              
              -
              
                
                  
                    Γ
                    ⁡
                    
                      (
                      
                        r
                        +
                        2
                      
                      )
                    
                  
                  1
                
                ⁢
                
                  
                    F
                    2
                  
                  ⁡
                  
                    (
                    
                      
                        
                          r
                          +
                          1
                        
                        ;
                        
                          
                            r
                            2
                          
                          +
                          
                            1
                            2
                          
                        
                      
                      ,
                      
                        
                          r
                          2
                        
                        ;
                        
                          -
                          
                            
                              τ
                              2
                            
                            4
                          
                        
                      
                    
                    )
                  
                
              
            
            ]
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          ss
        
        )
      
    
  




where 1F2 refers to the generalized hypergeometric function.

Case 2:





  
    Let
    ⁢
    
      
    
    ⁢
    
      
        g
        0
      
      ⁡
      
        (
        x
        )
      
    
  
  =
  
    
      1
      
        1
        +
        
          x
          2
        
      
    
    .
  




Its Taylor-series expansion (around x=0) is





  
    
      
        
          1
          
            1
            +
            
              x
              2
            
          
        
        =
        
          
            H
            ⁡
            
              (
              x
              )
            
          
          -
          
            x
            2
          
          +
          
            x
            4
          
          -
          
            
              x
              6
            
            ⁢
            
              
            
            ⁢
            ...
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          tt
        
        )
      
    
  




Since the series expansion has only even powers of x, we consider gn(x) of the form,





  
    
      
        
          
            
              g
              n
            
            ⁡
            
              (
              x
              )
            
          
          =
          
            
              
                x
                
                  
                    2
                    ⁢
                    
                      
                    
                    ⁢
                    n
                  
                  -
                  2
                
              
              -
              
                x
                
                  2
                  ⁢
                  
                    
                  
                  ⁢
                  n
                
              
            
            
              
                (
                
                  1
                  +
                  
                    x
                    2
                  
                
                )
              
              
                
                  2
                  ⁢
                  
                    
                  
                  ⁢
                  n
                
                -
                1
              
            
          
        
        ,
        
          n
          =
          1
        
        ,
        2
        ,
        ...
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          uu
        
        )
      
    
  




These functions satisfy properties (1)-(4). Using series expansion of





  1
  
    
      (
      
        1
        +
        
          x
          2
        
      
      )
    
    
      
        2
        ⁢
        
          
        
        ⁢
        n
      
      -
      1
    
  




around x=0, we get,





  
    
      
        
          
        
        ⁢
        
          
            
              ∑
              
                k
                =
                1
              
              ∞
            
            ⁢
            
              
            
            ⁢
            
              
                a
                k
              
              ⁢
              
                
                  g
                  k
                
                ⁡
                
                  (
                  x
                  )
                
              
            
          
          =
          
            
              
                ∑
                
                  k
                  =
                  1
                
                ∞
              
              ⁢
              
                
              
              ⁢
              
                
                  a
                  k
                
                ⁢
                
                  
                    
                      -
                      
                        x
                        
                          
                            2
                            ⁢
                            
                              
                            
                            ⁢
                            k
                          
                          -
                          2
                        
                      
                    
                    +
                    
                      x
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        k
                      
                    
                  
                  
                    
                      (
                      
                        1
                        +
                        
                          x
                          2
                        
                      
                      )
                    
                    
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        k
                      
                      -
                      1
                    
                  
                
              
            
            =
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          vv
        
        )
      
    
  
  
    
      
        
          ∑
          
            k
            =
            1
          
          ∞
        
        ⁢
        
          
        
        ⁢
        
          
            
              a
              k
            
            ⁡
            
              (
              
                
                  x
                  
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      k
                    
                    -
                    2
                  
                
                -
                
                  x
                  
                    2
                    ⁢
                    
                      
                    
                    ⁢
                    k
                  
                
              
              )
            
          
          ⁡
          
            [
            
              
                ∑
                
                  m
                  =
                  0
                
                ∞
              
              ⁢
              
                
              
              ⁢
              
                
                  
                    (
                    
                      -
                      1
                    
                    )
                  
                  m
                
                ⁢
                
                  (
                  
                    
                      
                        
                          
                            2
                            ⁢
                            
                              
                            
                            ⁢
                            k
                          
                          +
                          m
                          -
                          2
                        
                      
                    
                    
                      
                        m
                      
                    
                  
                  )
                
                ⁢
                
                  x
                  
                    2
                    ⁢
                    
                      
                    
                    ⁢
                    m
                  
                
              
            
            ]
          
        
      
    
    
      
        (
        
          6
          ⁢
          ww
        
        )
      
    
  




Matching the coefficients for x2n−2 in eqns. (49) and (52) yields





  
    
      
        
          a
          n
        
        =
        
          
            
              (
              
                -
                1
              
              )
            
            
              n
              -
              1
            
          
          -
          
            
              ∑
              
                k
                =
                1
              
              
                n
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  (
                  
                    -
                    1
                  
                  )
                
                
                  n
                  -
                  k
                
              
              ⁢
              
                
                  
                    a
                    k
                  
                  ⁡
                  
                    [
                    
                      
                        (
                        
                          
                            
                              
                                n
                                +
                                k
                                -
                                2
                              
                            
                          
                          
                            
                              
                                n
                                -
                                k
                              
                            
                          
                        
                        )
                      
                      +
                      
                        (
                        
                          
                            
                              
                                n
                                +
                                k
                                -
                                3
                              
                            
                          
                          
                            
                              
                                n
                                -
                                k
                                -
                                1
                              
                            
                          
                        
                        )
                      
                    
                    ]
                  
                
                .
              
            
          
        
      
    
    
      
        (
        
          6
          ⁢
          
            
          
          ⁢
          xx
        
        )
      
    
  




The first three coefficients are a1=0, a2=−1 and a3=−3. Let




  τ
  =
  
    
      t
      
        T
        c
      
    
    .
  




The inverse Laplace transforms of the first three terms in the series expansion are




  
    
      
    
    ⁢
    
      
        
          L
          
            -
            1
          
        
        ⁡
        
          [
          
            
              1
              -
              
                x
                2
              
            
            
              (
              
                1
                +
                
                  x
                  2
                
              
              )
            
          
          ]
        
      
      =
      
        
          1
          
            T
            c
          
        
        ⁡
        
          [
          
            
              2
              ⁢
              
                
              
              ⁢
              
                sin
                ⁡
                
                  (
                  τ
                  )
                
              
            
            -
            
              δ
              ⁡
              
                (
                τ
                )
              
            
          
          ]
        
      
    
  




  
    
      
    
    ⁢
    
      
        
          L
          
            -
            1
          
        
        ⁡
        
          [
          
            
              
                x
                2
              
              -
              
                x
                4
              
            
            
              
                (
                
                  1
                  +
                  
                    x
                    2
                  
                
                )
              
              3
            
          
          ]
        
      
      =
      
        
          1
          
            T
            c
          
        
        ⁡
        
          [
          
            
              1
              4
            
            ⁢
            
              (
              
                
                  
                    τ
                    2
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                -
                
                  sin
                  ⁡
                  
                    (
                    τ
                    )
                  
                
                -
                
                  3
                  ⁢
                  τ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
              
              )
            
          
          ]
        
      
    
  




  
    
      
        L
        
          -
          1
        
      
      ⁡
      
        [
        
          
            
              x
              4
            
            -
            
              x
              6
            
          
          
            
              (
              
                1
                +
                
                  x
                  2
                
              
              )
            
            5
          
        
        ]
      
    
    =
    
      
        
          1
          
            T
            c
          
        
        ⁡
        
          [
          
            
              1
              192
            
            ⁢
            
              (
              
                
                  
                    τ
                    4
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                -
                
                  10
                  ⁢
                  
                    τ
                    3
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                -
                
                  21
                  ⁢
                  
                    τ
                    2
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      t
                      )
                    
                  
                
                -
                
                  3
                  ⁢
                  
                    
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
                +
                
                  3
                  ⁢
                  τ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      τ
                      )
                    
                  
                
              
              )
            
          
          ]
        
      
      .
    
  



It can be shown that the Taylor-series expansion of the generating function in terms of anti-symmetric higher-order polynomials systematically leads to convergence of the generating function to a Heaviside function in log(x) space.
As explained in previously incorporated U.S. Ser. No. 13/333,232, the above describes methods for computing linear functionals of the distribution function without first computing the distribution of relaxation times. These methods involve a linear transform of the measured data using integral transforms. Different linear functionals of the distribution function can be obtained by choosing appropriate functions in the integral transforms. This approach can be used such that integral transforms can be computed on data corresponding to longitudinal relaxation time (T1), and such that integral transforms can be computed on data corresponding to diffusion coefficient (D). In addition, this approach can be extended to multiple dimensions.
Using the techniques given in previously incorporated U.S. Ser. No. 13/333,232, it is possible to estimate the uncertainty in the parameters. Let the discretized version of the linear transform be denoted by A=kTG, where k is the discretization of the k(t) in equation (4) and G refers to the discretization of the measured data indications. Let σε2 be the variance of the vector of measurement indications G. Then the variance of A is given by

σA2=σε2∥k∥2  (7)

SUMMARY
This summary is provided to introduce a selection of concepts that are further described below in the detailed description. This summary is not intended to identify key or essential features of the claimed subject matter, nor is it intended to be used as an aid in limiting the scope of the claimed subject matter.
Embodiments are provided for the estimation of an NMR-related distribution that is consistent with NMR measurements and uses linear functionals directly estimated from the measurement indications by integral transforms as constraints in a cost function. This cost function includes indications of the measurement data, Laplace transform elements and the constraints, and an NMR-related distribution estimation is made by minimizing the cost function. The NMR-related distribution estimation may then be used to find parameters of the sample. Where the sample is a rock or a formation, the parameters may include parameters such as rock permeability and/or hydrocarbon viscosity, bound and free fluid volumes, among others. The parameters may then be used, if desired, in models, equations, or otherwise to act on the sample, such as in recovering hydrocarbons from the formation. Examples of NMR-related distributions include transverse relaxation time (T2) distributions, longitudinal relaxation time (T1) distributions, Diffusion (D) distributions, and various multi-dimensional distributions, although embodiments are not limited thereto.



BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 is a high level block diagram showing a method of an embodiment.
FIG. 2 provides a first graph of an assumed T2 distribution and a second graph comparing estimated T2 distributions resulting from ILT processing and processing of an embodiment.
FIGS. 3a-3d provide sets of graphs, each set providing a model of simulated T2 distributions and noise-free echoes used in benchmarking ILT processing and the processing of an embodiment on simulated data.
FIG. 4 is a graph showing values of weighted terms used in the processing of an embodiment.
FIGS. 5a-5b show sets of histograms comparing estimates at different signal to noise ratios obtained by ILT processing and the processing of an embodiment for the model of FIG. 3a. 
FIGS. 6a-6b show sets of histograms comparing estimates at different signal to noise ratios obtained by ILT processing and the processing of an embodiment for the model of FIG. 3b. 
FIG. 7 is a graph showing the fit error as a function of a parameter for the model of FIG. 3a for the ILT processing and for the processing of an embodiment at a given signal to noise ratio.



DETAILED DESCRIPTION OF THE EMBODIMENTS
In one embodiment, constraints obtained using the integral transform approach described previously incorporated U.S. Ser. No. 13/333,232 are applied to a cost function which also includes indications of measured NMR data and Laplace transform elements, and a T2 distribution estimation is made by minimizing the cost function. The T2 distribution estimation may then be used to find parameters of the sample.
In an embodiment seen in FIG. 1, data G(t) 10, indicative of what might have been measured by an NMR tool as a result of having applied a pulse sequence to a sample, is processed at 20 using the integral transform approach described previously incorporated U.S. Ser. No. 13/333,232 to provide a plurality of linear functionals, e.g., T2ω, Tapered areas 30. As indicated, porosity calculations 40 which may be obtained by processing the NMR data G(t), or from other sources can be used as inputs to the integral transform approach processing 20. The linear functionals obtained at 30 are then used as constraints or “priors” in a cost function 50 incorporating them as well as Laplace transform elements and utilizing indications of the NMR data in order to generate a T2 distribution denoted by ƒNEW(T2) 60. The calculated T2 distribution 60 is an answer product that itself may be used for any of many purposes. By way of example only, the T2 distribution 60 may be used at 70 to generate an estimate of one or more parameters or properties of the sample. Where the sample is a rock or a formation, the parameters may include parameters such as rock permeability and/or hydrocarbon viscosity, bound and free fluid volumes, among others. The parameters may then be used at 80, if desired, in models, equations, or otherwise to act on the sample, such as in recovering hydrocarbons from the formation.
It should be appreciated that the data acquired at the NMR tool can result from any of a variety of pulse sequences including, but not limited to CPMG, diffusion editing, inversion-recovery, saturation-recovery, and enhanced precision mode (EPM) where the data are acquired at two different wait times.
According to an embodiment, the linear constraints generated by using the integral transform approach can be included into the problem formulation in a number of ways. While one example is described below with constraints from moments and areas being included in the problem formulation, the problem may be formulated with other constraints.
Constraints from the moments and areas determined using the integral transform approach described previously incorporated U.S. Ser. No. 13/333,232 can be included in the problem formulation as follows. Let G be a vector containing indications of the NMR measurements made by an NMR tool on a sample (e.g., the data itself, the data in compressed form by means of the singular value decomposition or windows sum, or other indications of the data), and let L be a matrix representing the discretization of the function in equation (1). Assume that Nm moments denoted by T2ωi, i=1, . . . , Nm, and Na areas denoted by Bi, i=1, . . . , Na are estimated directly by means of the appropriate integral transforms on G(t) as described in previously incorporated Ser. No. 13/333,232. A cost function with respect to ƒNEW T2ƒNEW(T2) may be minimized,

minƒNEW≧0∥W( G− LƒNEW)∥2+α∥ƒNEW∥2  (8)

where ƒNEW is a discretized vector version of an underlying density function ƒNEW(T2), W, as described in more detail below, is a covariance matrix of uncertainties in the parameters, α is a regularization parameter,





  
    
      
        
          G
          _
        
        =
        
          [
          
            
              
                G
              
            
            
              
                
                  〈
                  
                    T
                    2
                    
                      ω
                      1
                    
                  
                  〉
                
              
            
            
              
                ⋮
              
            
            
              
                
                  〈
                  
                    T
                    2
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                  〉
                
              
            
            
              
                
                  B
                  1
                
              
            
            
              
                ⋮
              
            
            
              
                
                  B
                  
                    N
                    a
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        9
        )
      
    
  




is an extended vector containing the indications of the measurements G as well as the constraints (moments T2ωi and areas Bi) generated by the integral transform approach, and L is the extended matrix,





  
    
      
        
          L
          _
        
        =
        
          [
          
            
              
                
                  
                
              
              
                L
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      min
                    
                    
                      ω
                      1
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      max
                    
                    
                      ω
                      1
                    
                  
                
              
            
            
              
                
                  
                
              
              
                ⋱
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      min
                    
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      max
                    
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                
              
            
            
              
                
                  
                
              
              
                
                  H
                  ⁡
                  
                    (
                    
                      
                        T
                        
                          c
                          1
                        
                      
                      ,
                      
                        T
                        2
                      
                    
                    )
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                ⋮
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  H
                  ⁡
                  
                    (
                    
                      
                        T
                        
                          N
                          a
                        
                      
                      ,
                      
                        T
                        2
                      
                    
                    )
                  
                
              
              
                
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        10
        )
      
    
  




where L is a Laplace transform matrix with components (L)ij=Pτ(T2,j)e−ti/T2,j, φ is the porosity, H(Tc, T2) represent the tapered Heaviside function varying smoothly between 0 and 1 with H(Tc, T2)=0.5 when T2=Tc as described in previously incorporated U.S. Ser. No. 13/333,232, and T2,min and T2,max represent the minimum and maximum value of the discretized T2 vector so that for a given T2 distribution ƒNEW with components ƒNEW(T2,i) the ω-th moment is defined as





  
    
      
        
          〈
          
            T
            2
            ω
          
          〉
        
        =
        
          
            1
            φ
          
          ⁢
          
            
              ∑
              
                i
                =
                min
              
              max
            
            ⁢
            
              
            
            ⁢
            
              
                T
                
                  2
                  ,
                  i
                
                ω
              
              ⁢
              
                
                  f
                  NEW
                
                ⁡
                
                  (
                  
                    T
                    
                      2
                      ,
                      i
                    
                  
                  )
                
              
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




Thus, L contains Laplace transform elements as well as functions of the constraints. As previously mentioned, W is the co-variance matrix of uncertainties in the parameters and may be described by





  
    
      
        W
        =
        
          [
          
            
              
                
                  1
                  
                    σ
                    ε
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  1
                  
                    σ
                    
                      ω
                      1
                    
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  
                
              
              
                ⋱
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  1
                  
                    σ
                    
                      N
                      m
                    
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  1
                  
                    σ
                    
                      B
                      i
                    
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                ⋱
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
              
                
                  1
                  
                    σ
                    
                      N
                      m
                    
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        12
        )
      
    
  




where σε is the standard deviation of noise in the measurements and σωi and σBi are the estimated uncertainties in the moments and areas estimated according to equation (7).

It will be appreciated by those skilled in the art that any of many optimization routines may be utilized to solve the cost function of equation (8) for ƒNEW. By way of example only, MatLab® software sold by MathWorks of Natick, Mass. includes a NNLS (non-negative least squares) routine which may be used to solve equation (8). Likewise, by way of example only, IMSL (International Mathematics and Statistics Library) sold by Rogue Wave Software of Boulder Colorado includes optimization routines which may be used to solve equation (8).
Extensive testing with simulated data shows that the addition of constraints generated by the integral transform approach as part of the data processing tends to eliminate artifacts in the estimated T2 distribution. For example, the bottom graph of FIG. 2 shows an actual T2 distribution, and the top graph of FIG. 2 compares estimated T2 distributions using the standard ILT processing and the processing of the described embodiment (denoted “NEW” in the Figures) for one noisy data set with a noise standard deviation σε=0.2. The general shape of the distribution using the embodiment appears a little closer than the distribution using standard ILT processing to the general shape of the actual distribution. In addition, parameters such as the mean-log of the T2 distribution referred to as T2LM, and the bound fluid volume (BFV) obtained with a sharp cut-off of 33 msec from ILT are considerably closer to the true values for the embodiment (true T2LM=85.8 ms, ILT T2LM=72.3 ms, NEW T2LM=88.5 ms; true BFV=0.17, ILT BFV=0.23, NEW BFV=0.19). For both ILT and the embodiment a nominal value of α=1 was chosen.
An embodiment was benchmarked on simulated data obtained from T2 distributions for four models M1-M4 shown respectively in FIGS. 3a-3d. Data are generated assuming a fixed porosity of φ=1 with additive white Gaussian noise with standard deviation per echo of σε=0.2 (signal-to-noise ratio SNR=50) and σε=0.2 (SNR=5) representing lab (i.e., core) data and log data, respectively. To demonstrate the performance of the algorithm on fully and partially polarized data, it is assumed that data from models M1 and M2 are fully polarized. Data from models M3 and M4 are obtained in the EPM mode at two different wait times, indicated in the FIGS. 3c and 3d by the legend—‘main pulse’ and ‘burst’.
The noisy data are analyzed using the ILT processing and the embodiment. For each realization, the same regularization parameter α was used in ILT and embodiment. The value of α was selected according to the Butler-Reed-Dawson algorithm. See “J. P. Butler, et al., “Estimating solutions of the first kind integral equations with nonnegative constraints and optimal smoothing, SIAM J. Numerical Analysis, 18(3):381-397, 1981. Several moments and areas in equation (8) were tested initially but in the following simulations twelve moments with ω between −0.5 and 1 and three tapered areas with Tc=0.01, 0.1 and 1 sec. are considered. The algorithm is not sensitive to the number of moments and areas. The data was compressed using the truncated singular value decomposition and the number of singular values was selected to achieve a condition number (ratio between first singular value and last singular value used in the calculations) of 1000. FIG. 4 shows the weighted and compressed data (first 10 values), moments (next 12 values), and areas (last 3 values) used in equation (8).
As before, two parameters, T2LM and BFV, are used to compare the performance of ILT and the embodiment. FIGS. 5a and 5b respectively show the estimated T2LM for Model M1 with SNR=50 and SNR=5 analyzed with 500 realizations of noise. As seen in FIGS. 5a and 5b, the embodiment (NEW) estimates show a smaller bias and standard deviation in the estimated T2LM than do the ILT estimates. This behavior was also seen for the other models (see, for example, FIGS. 6a and 6b which compare the performance of ILT and embodiment for Model M2). The performance of ILT and the embodiment for estimating T2LM and BFV are also summarized in Tables 1a, 1b, 2a, and 2b where μ indicates the mean of the estimate, σ indicates the standard deviation, {circumflex over (σ)} represents the estimated uncertainty, and rmse is the root mean square error defined for any estimate Est as

rmse√{square root over (Est−True value)2)}

where · represents the mean, and nmrse represents the normalized rmse defined as





  
    
      
        nrmse
        =
        
          
            rmse
            
              True
              ⁢
              
                
              
              ⁢
              Value
            
          
          ×
          100
        
      
    
    
      1
    
  











TABLE 1a







Comparison of T2LM; SNR = 50














True

ILT






T2LM
ILT
nrmse
NEW
NEW
{circumflex over (σ)}T2LM


Models
(ms)
μ ± σ
(%)
μ ± σ
nrmse (%)
μ
















M1
67.2
64.4 ± 3.6
6.8
67.3 ± 1.1
1.7
2.1


M2
85.8
82.7 ± 3.8
5.7
86.4 ± 2.1
2.5
3


M3
54.8
54.1 ± 0.9
2
54.6 ± 0.3
0.6
0.4


M4
73.8
73.1 ± 1.3
1.9
73.8 ± 0.5
0.6
0.5
















TABLE 1b







Comparison of T2LM; SNR = 5
















ILT







ILT
nrmse
NEW
NEW
{circumflex over (σ)}T2LM


Models
True
μ ± σ
(%)
μ ± σ
nrmse (%)
μ
















M1
67.2
58.7 ± 10.6
20.2
70.2 ± 8.1
12.8
16.1


M2
85.8
77.2 ± 13  
18.2
92.8 ± 9.8
14
21


M3
54.8
50.3 ± 6.1 
13.8
53.3 ± 2.3
5
2.9


M4
73.8
 69 ± 6.2
10.7
71.5 ± 3.2
5.3
4.1
















TABLE 2a







Comparison of BFV; SNR = 50














ILT
ILT
NEW
NEW


Models
True
μ ± σ
rmse
μ ± σ
rmse















M1
0.28
0.25 ± 0.019
0.03
0.25 ± 0.013
0.03


M2
0.17
.014 ± 0.012
0.03
0.13 ± 0.009
0.04


M3
0
  0 ± 0.003
0
  0 ± 0.001
0


M4
0
  0 ± 0.003
0
0 ± 0 
0
















TABLE 2b







Comparison of BFV; SNR = 5














ILT
ILT
NEW
NEW


Models
True
μ ± σ
rmse
μ ± σ
rmse















M1
0.28
0.22 ± 0.07 
0.03
 0.16 ± 0.066
0.14


M2
0.17
0.17 ± 0.058
0.03
 0.1 ± 0.055
0.09


M3
0
0.14 ± 0.066
0
0.09 ± 0.06
0.11


M4
0
0.11 ± 0.043
0
0.06 ± 0.04
0.07









These tables show that the estimates of the T2LM from the embodiment (NEW) have less bias and either smaller or comparable variance than the ILT processing. The results were mixed for estimation of BFV, especially in the case of low SNR. The tables also indicate that the uncertainty estimate of parameters such as T2LM compares reasonably well with the standard deviation as obtained from Monte-Carlo simulations of the data.
It was found that the addition of linear functionals as constraints to the estimation of the T2 distribution reduces the dependence of the error in the fit to the regularization parameter. Let χILT2=∥G−KƒILT∥2 be the error in the fit obtained from using the T2 distribution ƒILT estimated by using ILT and χNEW2=∥G−KƒNEW∥2 the corresponding error in the fit using the embodiment. FIG. 7 shows the fit error as a function of α from each. In this example, the data is obtained from model M1 with SNR=5. It is observed that the fit error is independent of α for a larger range when using the embodiment.
According to another embodiment, instead of minimizing a cost function with respect to a transverse relaxation time (T2) to obtain a T2 distribution estimation ƒ(T2), a cost function is minimized with respect to a longitudinal relaxation time (T1) to obtain a T1 distribution estimation ƒ(T1). More specifically, as in FIG. 1, data G(t) indicative of what might have been measured by an NMR tool as a result of having applied a pulse sequence to a sample, is processed using the integral transform approach described previously incorporated U.S. Ser. No. 13/333,232 to provide a plurality of linear functionals, e.g., T1ω, Tapered areas. Porosity calculations which may be obtained by processing the NMR data G(t), or from other sources can be used as inputs to the integral transform approach processing. The linear functionals obtained are then used as constraints or “priors” in a cost function incorporating them as well as Laplace transform elements and utilizing indications of the NMR data in order to generate a T1 distribution denoted by ƒNEW(T1). The calculated T1 distribution may be used for any of many purposes such as (by way of example only) to generate an estimate of one or more parameters of the sample. Where the sample is a rock or a formation, the parameters may include parameters such as rock permeability and/or hydrocarbon viscosity, bound and free fluid volumes, among others. The parameters may then be used, if desired, in models, equations, or otherwise to act on the sample, such as in recovering hydrocarbons from the formation.
It will be appreciated by those of ordinary skill in the art, that when minimizing a cost function with respect to a longitudinal relaxation time (T1) to obtain a T1 distribution estimation ƒ(T1), equation 1 can be modified to

G=(t)=∫0∞Pτ(T1)e−t/T1ƒ(T1)dT1,  (13)

where the function Pτ(T1) depends on the pulse sequence of the NMR equipment used to probe and measure the sample. Now, if G is a vector containing indications of the NMR measurements made by an NMR tool on a sample (e.g., the data itself, the data in compressed form by means of the singular value decomposition or windows sum, or other indications of the data), and L is a matrix representing the discretization of the function in equation (13), and it is assumed that Nm moments denoted by T1ωi, i=1, . . . , Nm, and Na areas denoted by Bi, i=1, . . . , Na are estimated directly by means of the appropriate integral transforms on G(t) as described in previously incorporated Ser. No. 13/333,232, a cost function with respect to ƒNEW(T1) may be minimized,

minƒNEW≧0∥W( G− LƒNEW)∥2+α∥ƒNEW∥2  (14)

where ƒNEW is a discretized vector version of an underlying density function ƒNEW(T1), W is a covariance matrix of uncertainties in the parameters (as described in (12) above), α is a regularization parameter,





  
    
      
        
          G
          _
        
        =
        
          [
          
            
              
                G
              
            
            
              
                
                  〈
                  
                    T
                    2
                    
                      ω
                      1
                    
                  
                  〉
                
              
            
            
              
                ⋮
              
            
            
              
                
                  〈
                  
                    T
                    1
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                  〉
                
              
            
            
              
                
                  B
                  1
                
              
            
            
              
                ⋮
              
            
            
              
                
                  B
                  
                    N
                    a
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        15
        )
      
    
  




is an extended vector containing the indications of the measurements G as well as the constraints (moments T1ωi and areas Bi) generated by the integral transform approach, and L is the extended matrix,





  
    
      
        
          L
          _
        
        =
        
          [
          
            
              
                
                  
                
              
              
                L
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      min
                    
                    
                      ω
                      1
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      max
                    
                    
                      ω
                      1
                    
                  
                
              
            
            
              
                
                  
                
              
              
                ⋱
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      min
                    
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      max
                    
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                
              
            
            
              
                
                  
                
              
              
                
                  H
                  ⁡
                  
                    (
                    
                      
                        T
                        
                          c
                          1
                        
                      
                      ,
                      
                        T
                        1
                      
                    
                    )
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                ⋮
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  H
                  ⁡
                  
                    (
                    
                      
                        T
                        
                          N
                          a
                        
                      
                      ,
                      
                        T
                        1
                      
                    
                    )
                  
                
              
              
                
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        16
        )
      
    
  




where L is a Laplace transform matrix with components (L)ij=Pτ(T1,j)e−ti/T1,j, φ is the porosity, H(Tc, T1) represent the tapered Heaviside function varying smoothly between 0 and 1 with H(Tc, T1)=0.5 when T1=Tc as described in previously incorporated U.S. Ser. No. 13/333,232, and T1,min and T1,max represent the minimum and maximum value of the discretized T1 vector so that for a given T1 distribution ƒNEW with components ƒNEW(T1,i) the ω-th moment is defined as





  
    
      
        
          〈
          
            T
            1
            ω
          
          〉
        
        =
        
          
            1
            φ
          
          ⁢
          
            
              ∑
              
                i
                =
                min
              
              max
            
            ⁢
            
              
                T
                
                  1
                  ,
                  i
                
                ω
              
              ⁢
              
                
                  f
                  NEW
                
                ⁡
                
                  (
                  
                    T
                    
                      1
                      ,
                      i
                    
                  
                  )
                
              
            
          
        
      
    
    
      
        (
        17
        )
      
    
  




Thus, L contains Laplace transform elements as well as functions of the constraints.

It will be appreciated by those skilled in the art that any of many optimization routines may be utilized to solve the cost function of equation (14) for ƒNEW.
In a similar manner, according to another embodiment, instead of minimizing a cost function with respect to a relaxation time to obtain a relaxation time distribution estimation, a cost function is minimized with respect to NMR diffusion (D) to obtain a D distribution estimation ƒ(D). More specifically, as in FIG. 1, data G(t) indicative of what might have been measured by an NMR tool as a result of having applied a pulse sequence to a sample, is processed using the integral transform approach described previously incorporated U.S. Ser. No. 13/333,232 to provide a plurality of linear functionals, e.g., Dω, Tapered areas. Porosity calculations which may be obtained by processing the NMR data G(t), or from other sources can be used as inputs to the integral transform approach processing. The linear functionals obtained are then used as constraints or “priors” in a cost function incorporating them as well as Laplace transform elements and utilizing indications of the NMR data in order to generate a D distribution denoted by ƒNEW(D). The calculated D distribution may be used for any of many purposes such as (by way of example only) to generate an estimate of one or more parameters of the sample. Where the sample is a rock or a formation, the parameters may include parameters such as rock permeability and/or hydrocarbon viscosity, bound and free fluid volumes, among others. The parameters may then be used, if desired, in models, equations, or otherwise to act on the sample, such as in recovering hydrocarbons from the formation.
It will be appreciated by those of ordinary skill in the art, that when minimizing a cost function with respect to diffusion (D) to obtain a D distribution estimation ƒ(D), equation 1 can be modified to

G(t)=∫0∞Pτ(D)e−tDƒ(D)dD,  (18)

where the function Pτ(D) depends on the pulse sequence of the NMR equipment used to probe and measure the sample. Now, if G is a vector containing indications of the NMR measurements made by an NMR tool on a sample (e.g., the data itself, the data in compressed form by means of the singular value decomposition or windows sum, or other indications of the data), and L is a matrix representing the discretization of the function in equation (18), and it is assumed that Nm moments denoted by Dωi, i=1, . . . , Nm, and Na areas denoted by Bi, i=1, . . . , Na are estimated directly by means of the appropriate integral transforms on G(t) as described in previously incorporated Ser. No. 13/333,232, a cost function with respect to ƒNEW(D) may be minimized,

minƒNEW≧0∥W( G− LƒNEW)∥2+α∥ƒNEW∥2  (19)

where ƒNEW is a discretized vector version of an underlying density function ƒNEW(D), α is a regularization parameter,





  
    
      
        
          G
          _
        
        =
        
          [
          
            
              
                G
              
            
            
              
                
                  〈
                  
                    D
                    
                      ω
                      i
                    
                  
                  〉
                
              
            
            
              
                ⋮
              
            
            
              
                
                  〈
                  
                    D
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                  〉
                
              
            
            
              
                
                  B
                  1
                
              
            
            
              
                ⋮
              
            
            
              
                
                  B
                  
                    N
                    a
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        20
        )
      
    
  




is an extended vector containing the indications of the measurements G as well as the constraints (moments Dωi and areas Bi) generated by the integral transform approach, W is a covariance matrix of G, and L is the extended matrix,





  
    
      
        
          L
          _
        
        =
        
          [
          
            
              
                
                  
                
              
              
                L
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    D
                    min
                    
                      ω
                      1
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    D
                    max
                    
                      ω
                      1
                    
                  
                
              
            
            
              
                
                  
                
              
              
                ⋱
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    D
                    min
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    D
                    max
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                
              
            
            
              
                
                  
                
              
              
                
                  H
                  ⁡
                  
                    (
                    
                      
                        D
                        
                          c
                          1
                        
                      
                      ,
                      D
                    
                    )
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                ⋮
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  H
                  ⁡
                  
                    (
                    
                      
                        D
                        
                          N
                          a
                        
                      
                      ,
                      D
                    
                    )
                  
                
              
              
                
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        21
        )
      
    
  




where L is a Laplace transform matrix with components (L)ij=Pτ(D)e−tD, φ is the porosity, H(Dc,D) represent the tapered Heaviside function varying smoothly between 0 and 1 with H(Dc,D)=0.5 when D=Dc as described in previously incorporated U.S. Ser. No. 13/333,232, and Dmin and Dmax represent the minimum and maximum value of the discretized D vector so that for a given D distribution ƒNEW with components ƒNEW(Di) the ω-th moment is defined as





  
    
      
        
          〈
          
            D
            ω
          
          〉
        
        =
        
          
            1
            φ
          
          ⁢
          
            
              ∑
              
                i
                =
                min
              
              max
            
            ⁢
            
              
                D
                i
                ω
              
              ⁢
              
                
                  f
                  NEW
                
                ⁡
                
                  (
                  
                    D
                    i
                  
                  )
                
              
            
          
        
      
    
    
      
        (
        22
        )
      
    
  




Thus, L contains Laplace transform elements as well as functions of the constraints.

According to other embodiments, the embodiments described above may be extended to extract multidimensional distributions such as (by way of example only) diffusion-T2 distribution functions, diffusion-T1 distribution functions and T1-T2 distribution functions. This may be done by one of ordinary skill in the art after reference to the previous embodiments and by reference to knowledge in the art such as described in M. D. Hurlimann and L. Venkataramanan, “Quantitative Measurement of Two-Dimensional Distribution Functions of Diffusion and Relaxation in Grossly Inhomogenous Fields,” Journal of Magnetic Resonance 157, 31-42 (2002); Y.-Q. Song et al., “T1-T2 Correlation Spectra Obtained Using a Fast Two-Dimensional Laplace Inversion,” Journal of Magnetic Resonance 154, 1-8 (2002); and R. L. Kleinberg et al., “Nuclear Magnetic Resonance of Rocks: T1 vs. T2,” SPE 26470 (1993). Thus, by way of example, for a T1-T2 distribution function, the measured nuclear magnetic resonance (NMR) data resulting from a multi-component sample can be denoted by G(t, τ) which represents a multidimensional, multi-exponential decay, with time constants T1 and T2 and amplitudes ƒ(T1, T2)

G(t,τ)=∫0∞∫0∞Pτ(T2)e−t/T2e−t/T1ƒ(T1,T2)dT1dT2  (23)

where the function Pτ(T2) is referred to as the polarization factor and depends on the pulse sequence of the NMR equipment used to probe and measure the sample. Now, if G is a vector containing indications of the NMR measurements made by an NMR tool on a sample (e.g., the data itself, the data in compressed form by means of the singular value decomposition or windows sum, or other indications of the data), and L is a matrix representing the discretization of the function in equation (23), and it is assumed that Nm moments denoted by T1ωiT2nj, i=1, . . . , Nn and NaNb areas denoted by Bi, i=1, . . . , NaNb are estimated directly by means of the appropriate integral transforms on G(t) as described in previously incorporated Ser. No. 13/333,232, a cost function with respect to ƒNEW may be minimized,

minƒNEW≧0∥W( G− LƒNEW)∥2+α∥ƒNEW∥2  (24)

where ƒNEW is a discretized vector version of an underlying density function ƒNEW(T1, T2), α is a regularization parameter,





  
    
      
        
          G
          _
        
        =
        
          [
          
            
              
                G
              
            
            
              
                
                  〈
                  
                    
                      T
                      1
                      
                        ω
                        1
                      
                    
                    ⁢
                    
                      T
                      2
                      
                        n
                        1
                      
                    
                  
                  〉
                
              
            
            
              
                ⋮
              
            
            
              
                
                  〈
                  
                    
                      T
                      1
                      
                        ω
                        
                          N
                          m
                        
                      
                    
                    ⁢
                    
                      T
                      2
                      
                        n
                        
                          N
                          n
                        
                      
                    
                  
                  〉
                
              
            
            
              
                
                  B
                  1
                
              
            
            
              
                ⋮
              
            
            
              
                
                  B
                  
                    
                      N
                      a
                    
                    ⁢
                    
                      N
                      b
                    
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        25
        )
      
    
  




is an extended vector containing the indications of the measurements G as well as the constraints (moments T1ωiT2nj and areas Bi) generated by the integral transform approach, W is a covariance matrix corresponding to G, and L is the extended matrix,





  
    
      
        
          L
          _
        
        =
        
          [
          
            
              
                
                  
                
              
              
                L
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      min
                    
                    
                      ω
                      1
                    
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      min
                    
                    
                      n
                      1
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      max
                    
                    
                      ω
                      1
                    
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      max
                    
                    
                      n
                      1
                    
                  
                
              
            
            
              
                
                  
                
              
              
                ⋱
              
              
                
                  
                
              
            
            
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      min
                    
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      min
                    
                    
                      n
                      N
                    
                  
                
              
              
                ...
              
              
                
                  
                    1
                    φ
                  
                  ⁢
                  
                    T
                    
                      1
                      ,
                      max
                    
                    
                      ω
                      
                        N
                        m
                      
                    
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      max
                    
                    
                      n
                      N
                    
                  
                
              
            
            
              
                
                  
                
              
              
                
                  
                    H
                    ⁡
                    
                      (
                      
                        
                          T
                          
                            c
                            
                              i
                              ⁢
                              
                                
                              
                              ⁢
                              1
                            
                          
                        
                        ,
                        
                          T
                          1
                        
                      
                      )
                    
                  
                  ⁢
                  
                    H
                    ⁡
                    
                      (
                      
                        
                          T
                          
                            c
                            j
                          
                        
                        ,
                        
                          T
                          2
                        
                      
                      )
                    
                  
                
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                ⋮
              
              
                
                  
                
              
            
            
              
                
                  
                
              
              
                
                  
                    H
                    ⁡
                    
                      (
                      
                        
                          T
                          
                            N
                            a
                          
                        
                        ,
                        
                          T
                          1
                        
                      
                      )
                    
                  
                  ⁢
                  
                    H
                    ⁡
                    
                      (
                      
                        
                          T
                          
                            N
                            b
                          
                        
                        ,
                        
                          T
                          2
                        
                      
                      )
                    
                  
                
              
              
                
                  
                
              
            
          
          ]
        
      
    
    
      
        (
        26
        )
      
    
  




where L is a Laplace transform matrix with components corresponding to Pτ(T1)e−t/T2e−τ/T1, φ is the porosity, H(Tc, T1) and H(Tc, T2) represent the tapered Heaviside function varying smoothly between 0 and 1 with H(Tc, T1) and H(Tc, T2)=0.5 when T1=Tc and T2=Tc as described in previously incorporated U.S. Ser. No. 13/333,232, and T1,min, T2,min and T1,max and T2,max represent the minimum and maximum value of the discretized T1 and T2 vectors so that for a given (T1, T2) distribution ƒNEW with components ƒNEW(T1,i, T2,j) the (ω,n)-th moment is defined as





  
    
      
        
          
            〈
            
              
                T
                1
                ω
              
              ⁢
              
                T
                2
                ω
              
            
            〉
          
          =
          
            
              1
              φ
            
            ⁢
            
              
                ∑
                
                  i
                  =
                  min
                
                max
              
              ⁢
              
                
                  ∑
                  
                    j
                    =
                    min
                  
                  max
                
                ⁢
                
                  
                    T
                    
                      1
                      ,
                      i
                    
                    ω
                  
                  ⁢
                  
                    T
                    
                      2
                      ,
                      j
                    
                    n
                  
                  ⁢
                  
                    
                      f
                      NEW
                    
                    ⁡
                    
                      (
                      
                        
                          T
                          
                            1
                            ,
                            i
                          
                        
                        ,
                        
                          T
                          
                            2
                            ,
                            j
                          
                        
                      
                      )
                    
                  
                
              
            
          
        
        )
      
    
    
      
        (
        27
        )
      
    
  




Thus, L contains Laplace transform elements as well as functions of the constraints.

According to other embodiments, the embodiments described above may be extended to extract multidimensional distributions such as (by way of example only) diffusion-T1-T2 distribution functions, diffusion-T1-T2 functions as a function of depth of investigation as well as diffusion-T1-T2 functions as a function of azimuth. This may be done by one of ordinary skill in the art after reference to the previous embodiments and by reference to knowledge in the art such as described in M.D. Hurlimann and L. Venkataramanan, “Quantitative Measurement of Two-Dimensional Distribution Functions of Diffusion and Relaxation in Grossly Inhomogenous Fields,” Journal of Magnetic Resonance 157, 31-42 (2002); Y.-Q. Song et al., “T1-T2 Correlation Spectra Obtained Using a Fast Two-Dimensional Laplace Inversion,” and L.Venkataramanan, Y. Q. Song and M. D. Hurlimann, “Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions,” vol. 50, No. 5, May 2002 and N. Heaton et al., “4D NMR—Applications of the radial dimension in magnetic resonance logging,” Petrophysics, 49, 2 (2008).
There have been described and illustrated herein several embodiments of methods for estimating distributions of NMR-related measurements. While particular embodiments have been described, it is not intended that the claims be limited thereto, as it is intended that the claims be as broad in scope as the art will allow and that the specification be read likewise. Thus, while particular moments and areas generated by using the integral transform approach were described as being included into the problem formulation, it will be appreciated that the problem may be formulated with other constraints. Also, while certain optimization routines were disclosed as being useful to solve the cost function ƒ or the distributions, it will be appreciated that other routines could be utilized. In addition, while the calculated distributions were described as being useful to find parameters of a rock or a formation such as rock permeability and/or hydrocarbon viscosity, and bound and free fluid volumes, it will be appreciated that the calculated distributions could be used to find other parameters of other types of samples, or other parameters of rocks or formation (e.g., saturates, aromatic, resin and asphaltene analysis of hydrocarbons, see M. D. Hurlimann et al “Hydrocarbon Composition from NMR Diffusion and Relaxation Data,” Petrophysics, 50, no. 2, 116-129 (2009)). Furthermore, while the parameters were described as being used in models, equations, or otherwise to act on the sample, such as in recovering hydrocarbons from a formation, it will be understood that depending upon the parameters and the sample, they could be used in other manners. It will therefore be appreciated by those skilled in the art that yet other modifications could be made without deviating from its spirit and scope of the claims.

## Claims

1. A method of investigating a sample, comprising:
performing nuclear magnetic resonance (NMR) measurements on the sample using an NMR tool to obtain NMR data;
conducting integral transforms on the NMR data to obtain linear functionals associated with said sample;
utilizing said linear functionals as constraints in a cost function which includes said NMR data, said linear functionals, and Laplace transform elements;
minimizing said cost function with respect to an NMR-related distribution to obtain a distribution estimation relating to the sample; and
using said distribution estimation to find values of at least one parameter of the sample.

2. A method according to claim 1, wherein:
said sample is a geological formation, and said at least one parameter includes at least one of permeability, hydrocarbon viscosity, bound fluid volume, and free fluid volume.

3. A method according to claim 2, further comprising:
using said at least one parameter in order to act on said formation to recover hydrocarbons from the formation.

4. A method according to claim 1, wherein:
said NMR-related distribution comprises at least one of T1, T2, and D, and said NMR distribution is one of a T1 distribution estimation ƒ(T1) where T1 is a longitudinal relaxation time for protons in the sample, a T2 distribution estimation ƒ(T2) where T2 is a transverse relaxation time for protons in the sample, and a D distribution estimate ƒ(D) where D is an indication of diffusion.

5. A method according to claim 4, wherein:
said NMR-related distribution is a multidimensional distribution comprising at least two of T1, T2, and D.

6. A method according to claim 4, wherein:
said minimizing said cost function comprises minimizing said cost function with respect to a T2 distribution to obtain a T2 distribution estimation f (T2) which is described by minƒ NEW≧0∥W( G− LƒNEW)∥2+α∥ƒNEW∥2 where W is a covariance matrix of uncertainties in said linear functionals, G is an extended vector containing said NMR data and said linear functionals, L is an extended matrix containing Laplace transform matrix elements, ƒNEW is a discretized vector version of an underlying density function ƒNEW(T2), and α is a regularization parameter.

7. A method according to claim 6, wherein:
said linear functionals comprise moments <T2ωi> and areas Bi.

8. A method according to claim 7, wherein:
said extended matrix is described by





  
    L
    _
  
  =
  
    [
    
      
        
          
            
          
        
        
          L
        
        
          
            
          
        
      
      
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                2
                ,
                min
              
              
                ω
                1
              
            
          
        
        
          ...
        
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                2
                ,
                max
              
              
                ω
                1
              
            
          
        
      
      
        
          
            
          
        
        
          ⋱
        
        
          
            
          
        
      
      
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                2
                ,
                min
              
              
                ω
                
                  N
                  m
                
              
            
          
        
        
          ...
        
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                2
                ,
                max
              
              
                ω
                
                  N
                  m
                
              
            
          
        
      
      
        
          
            
          
        
        
          
            H
            ⁡
            
              (
              
                
                  T
                  
                    c
                    1
                  
                
                ,
                
                  T
                  2
                
              
              )
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          ⋮
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            H
            ⁡
            
              (
              
                
                  T
                  
                    N
                    a
                  
                
                ,
                
                  T
                  2
                
              
              )
            
          
        
        
          
            
          
        
      
    
    ]
  




where L is a Laplace transform matrix with said Laplace transform elements (L)ij=Pτ(T2,j)e−ti/T2,j, φ is a porosity indication, H(Tc, T2) elements represent a tapered Heaviside function and T2,min and T2,max represent minimum and maximum value of said discretized T2 vector so that for a given T2 distribution f with components ƒNEW,i, the ω-th moment is defined as




  
    〈
    
      T
      2
      ω
    
    〉
  
  =
  
    
      1
      φ
    
    ⁢
    
      
        ∑
        
          i
          =
          min
        
        max
      
      ⁢
      
        
          T
          
            2
            ,
            i
          
          ω
        
        ⁢
        
          
            
              f
              NEW
            
            ⁡
            
              (
              
                T
                
                  2
                  ,
                  i
                
              
              )
            
          
          .

9. A method according to claim 8, wherein:
said Heaviside function varies smoothly between 0 and 1 with H(Tc,T2)=0.5 when T2=Tc.

10. A method according to claim 7, wherein:
said co-variance matrix of uncertainties is described by





  W
  =
  
    [
    
      
        
          
            1
            
              σ
              ∈
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            1
            
              σ
              
                ω
                1
              
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          ⋱
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            1
            
              σ
              
                N
                m
              
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            1
            
              σ
              
                B
                i
              
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          ⋱
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            1
            
              σ
              
                N
                m
              
            
          
        
      
    
    ]
  




where σε is the standard deviation of noise in said NMR measurements, and σωi and σBi the estimated uncertainties in said moments and said areas.

11. A method according to claim 4, wherein:
said minimizing said cost function comprises minimizing said cost function with respect to a T1 distribution to obtain a T1 distribution estimation f (T1) which is described by minƒ NEW≧0∥W( G− LƒNEW)∥2+α∥ƒNEW∥2 where W is a covariance matrix of uncertainties in said linear functionals, G is an extended vector containing said NMR data and said linear functionals, L is an extended matrix containing Laplace transform matrix elements, ƒNEW is a discretized vector version of an underlying density function ƒNEW(T1), and α is a regularization parameter.

12. A method according to claim 11, wherein:
said linear functionals comprise moments (T1ωi) and areas Bi.

13. A method according to claim 12, wherein:
said extended matrix is described by





  
    L
    _
  
  =
  
    [
    
      
        
          
            
          
        
        
          L
        
        
          
            
          
        
      
      
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                1
                ,
                min
              
              
                ω
                1
              
            
          
        
        
          ...
        
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                1
                ,
                max
              
              
                ω
                1
              
            
          
        
      
      
        
          
            
          
        
        
          ⋱
        
        
          
            
          
        
      
      
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                1
                ,
                min
              
              
                ω
                
                  N
                  m
                
              
            
          
        
        
          ...
        
        
          
            
              1
              φ
            
            ⁢
            
              T
              
                1
                ,
                max
              
              
                ω
                
                  N
                  m
                
              
            
          
        
      
      
        
          
            
          
        
        
          
            H
            ⁡
            
              (
              
                
                  T
                  
                    c
                    1
                  
                
                ,
                
                  T
                  1
                
              
              )
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          ⋮
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            H
            ⁡
            
              (
              
                
                  T
                  
                    N
                    a
                  
                
                ,
                
                  T
                  1
                
              
              )
            
          
        
        
          
            
          
        
      
    
    ]
  




where L is a Laplace transform matrix with said Laplace transform elements (L)ij=Pτ(T2,j)e−ti/T1,j, φ is a porosity indication, H(Tc, T1) elements represent a tapered Heaviside function and T1,min and T1,max represent minimum and maximum value of said discretized T1 vector so that for a given T1 distribution f with components ƒNEW,i, the ω-th moment is defined as




  
    〈
    
      T
      1
      ω
    
    〉
  
  =
  
    
      1
      φ
    
    ⁢
    
      
        ∑
        
          i
          =
          min
        
        max
      
      ⁢
      
        
          T
          
            1
            ,
            i
          
          ω
        
        ⁢
        
          
            
              f
              NEW
            
            ⁡
            
              (
              
                T
                
                  1
                  ,
                  i
                
              
              )
            
          
          .

14. A method according to claim 13, wherein:
said Heaviside function varies smoothly between 0 and 1 with H(Tc,T1)=0.5 when T1=Tc.

15. A method according to claim 12, wherein:
said co-variance matrix of uncertainties is described by





  W
  =
  
    [
    
      
        
          
            1
            
              σ
              ∈
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            1
            
              σ
              
                ω
                1
              
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          ⋱
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            1
            
              σ
              
                N
                m
              
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            1
            
              σ
              
                B
                i
              
            
          
        
        
          
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          ⋱
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            
          
        
        
          
            1
            
              σ
              
                N
                m
              
            
          
        
      
    
    ]
  




where σε is the standard deviation of noise in said NMR measurements, and σωi and σBi are the estimated uncertainties in said moments and said areas.

16. A method according to claim 4, wherein:
said minimizing said cost function comprises minimizing said cost function with respect to a D distribution to obtain a D distribution estimation f (D) which is described by minƒNEW≧0∥W( G− LƒNEW)∥2+α∥ƒNEW∥2 where W is a covariance matrix of uncertainties in said linear functionals, G is an extended vector containing said NMR data and said linear functionals, L is an extended matrix containing Laplace transform matrix elements, ƒNEW is a discretized vector version of an underlying density function ƒNEW(D), and α is a regularization parameter.

17. A method according to claim 16, wherein:
said linear functionals comprise moments (Dωi) and areas Bi.

18. A method according to claim 17, wherein:
said extended matrix is described by





  
    L
    _
  
  =
  
    [
    
      
        
          
            
          
        
        
          L
        
        
          
            
          
        
      
      
        
          
            
              1
              φ
            
            ⁢
            
              D
              min
              
                ω
                1
              
            
          
        
        
          ...
        
        
          
            
              1
              φ
            
            ⁢
            
              D
              max
              
                ω
                1
              
            
          
        
      
      
        
          
            
          
        
        
          ⋱
        
        
          
            
          
        
      
      
        
          
            
              1
              φ
            
            ⁢
            
              D
              min
              
                ω
                
                  N
                  m
                
              
            
          
        
        
          ...
        
        
          
            
              1
              φ
            
            ⁢
            
              D
              max
              
                ω
                
                  N
                  m
                
              
            
          
        
      
      
        
          
            
          
        
        
          
            H
            ⁡
            
              (
              
                
                  D
                  
                    c
                    1
                  
                
                ,
                D
              
              )
            
          
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          ⋮
        
        
          
            
          
        
      
      
        
          
            
          
        
        
          
            H
            ⁡
            
              (
              
                
                  D
                  
                    N
                    a
                  
                
                ,
                D
              
              )
            
          
        
        
          
            
          
        
      
    
    ]
  




where L is a Laplace transform matrix with said Laplace transform elements (L)ij=Pτ(D)e−tD, φ is a porosity indication, H(Dc,D) elements represent a tapered Heaviside function and Dmin and Dmax represent minimum and maximum value of said discretized D vector so that for a given D distribution f with components ƒNEW,i, the ω-th moment is defined as




  
    〈
    
      D
      ω
    
    〉
  
  =
  
    
      1
      φ
    
    ⁢
    
      
        ∑
        
          i
          =
          min
        
        max
      
      ⁢
      
        
          D
          i
          ω
        
        ⁢
        
          
            
              f
              NEW
            
            ⁡
            
              (
              
                D
                i
              
              )
            
          
          .

19. A method according to claim 18, wherein:
said Heaviside function varies smoothly between 0 and 1 with H(Dc,D)=0.5 when D=Dc.