A new two dimensional inverse Laplace transform algorithm using

By Kristen Spencer,2014-04-17 23:02
12 views 0
A new two dimensional inverse Laplace transform algorithm using

    Header From Top




     1,*234 Åsmund Ukkelberg, Geir Humborstad Sørland, Eddy Walter Hansen, Hege C. Widerøe

     1,3 The University of Oslo (UiO), Boks 1072 Blindern, N-0316 Oslo, Norway 2 Times new Roman 10pt Anvendt Teknologi AS, Munkvollveien 56, N-7022 Trondheim, Norway 4 Capital Each Word Statoil Research Centre, N-7005 Trondheim, Norway


    The problem of fitting a sum of exponentials to second order data is often called the second order (or two dimensional) inverse Laplace transform (2D-ILT). The present work describes ANAHESS, a new inverse Hessian Doc Margins : Left 1 Right 1 algorithm which handles this ill-posed numerical problem. The algorithm uses numerical expressions for the Page size Letter 8.5”*11” gradient and the Hessian matrix, and has performed well on various artificial and experimental data sets. The new algorithm is presented in a Nuclear Magnetic Resonance (NMR) context.

Keywords: Sum of exponentials, TT relaxation, Hessian, Bayesian information criterion, Earth mover’s distance 12

     Times new Roman 1. INTRODUCTION 10pt All CAPS

    2D NMR techniques are Nuclear Magnetic Resonance measurement techniques which produce second order data matrices. Such techniques are a frequently used tool in the characterization of fluid filled porous media, and they have found a wide variety of applications, such as analysis of food products, rock core samples and hydrating cement pastes. In these NMR techniques, the fitting of a sum of exponentials to the measured response surface is a vital tool. Such a fitting is often called the second order Laplace transform (2D-ILT). The goal of 2D-ILT is to transform the data from a sum of exponentially decaying components to a set of parameters characterizing the exponential decay of each component. Depending on the experimental setup, these parameters have a physical interpretation such as spin-lattice relaxation time T, spin-spin relaxation time Tor diffusion coefficient D. 12

    Although such a transform is notoriously ill conditioned numerically [1], algorithms that attempt to perform this transform are in widespread use. Hürlimann et al. [2] investigated the properties of dairy products using second order D-T and TTNMR 212

    experiments. Vital to these experiments was a 2D-ILT algorithm described by Venkataramanan et al. [3] and by

    Song et al. [4]. In food analysis, it is often seen that first order NMR measurements are insufficient in distinguishing between constituents such as fat and water, as their contributions to the total signal overlap too much. If second order experiments are done, however, 2D-ILT calculations and a study of T/T ratios or D/T ratios can be 212

    used to quantify food constituents [5]. In a similar study, McDonald et al. [6] used second order NMR

    measurements and the 2D-ILT algorithm made by Venkataramanan et al. [3] to investigate chemical exchange in

    cement pastes. Clear evidence was found for chemical exchange of water protons between gel and capillary pores in the cement pastes. Toft Pedersen et al. [7] have devised a method for curve resolution of NMR data which is based on the Direct Exponential Curve Resolution Algorithm (DECRA) proposed by Winding and Antalek [8]. This approach is based on numerical tools which are common in chemometrics. Toft Pedersen et al. used their algorithm to determine fat

    content in minced meat, and also pre slaughter stress in processed meat.

    Arns et al. [9] used 2D-ILT software when discussing susceptibility effects in oil reservoir rock core samples, with the objective of developing an NMR response simulation tool for X-ray-CT images, including the simulation of T1

    Trelaxation. It was found that Tand Tvalues depend on bulk vs. surface relaxation, and therefore depend upon 2 1 2 the pore structure within the rock core sample.

    Sørland et al. [10] investigated ways to improve the performance of 2D-ILT on rock cores containing oil and water. They found that by running a series of experiments at varying gradient strength, the different molecule mobilities of

    Footer From Bottom 1


oil and water could be used to separate the NMR signals from oil and from water. Doing this separation prior to 2D-

    ILT, enhanced the ability of 2D-ILT to extract valuable information about the rock core sample. The above applications, with the exception of the method by Toft Pedersen et al. [7] involve algorithms which fit the fractional intensities to a grid of equally spaced values of T and T (This applies to the case of T/T experiments. 1212

    This discussion could equally well be applied to D/T experiments. In the following, explanations will be restricted 2

    to T/T experiments in order to illustrate the concepts). The widely used algorithm implemented by Godefroy et al. 12

    [11], based on work by Venkataramanan et al. [3] is an example of such an approach. The grid of preset T and T 12

    values is an input to the algorithms, and may not reflect the optimal T and T values of the system investigated. 12Thus, without any prior knowledge of the optimal T and T values, those values may not be represented in the 12

    chosen grid. Instead, the algorithms provide an approximation to a continuous distribution in T and T. In contrast, 12

    we propose a new algorithm called ANAHESS which searches for the T and Tvalues that fit the dataset optimally, 12

    and provides a set of separate and distinct components which together constitute a model for the system investigated,

    instead of an estimate of a continuous distribution in T and T. 12


    Assume an experiment which produces a second order matrix R having NSA (number of samples) rows and NSE

    (number of sensors) columns. Assume also that the measured response is the linear sum of NCO different components, each component decaying exponentially with both the row index and the column index. The row-wise

    exponential decay follows an axis g, i=1..NSA, contained in a vector g of size NSA, while the column-wise iexponential decay follows an axis t, j=1..NSE, contained in a vector t of size NSE. Each component is assumed to j

    exhibit a row-wise decay with a characteristic parameter T, p=1..NCO, while the component’s column-wise decay 1,p

    is assumed to have a characteristic parameter T, p=1..NCO. The magnitude of each component is determined by a 2,p

    factor a, p=1..NCO, while any permanent signal offset present is assigned to the variable a. The spacing between p0

    subsequent element in g and t may or may not be constant. Any element R, i =1..NSA, j =1..NSE in R is then i,j

    assumed to follow the function:

    Rf(a,a,T,T);Ewith,012,,ijijijNCO((1/)(1/))TgTt (1) ij12,p,pf(a,a,T,T)a;ae012,0ijpp

with a, T and T being the elements with the index p=1..NCO in the vectors a, T and T, and subject to the p1,p2,p12

    following inequality constraints:


    T(0,p1..NCO (2) 1,p


    a is a baseline offset which may be positive or negative, and E is the contribution from noise. The parameters a, 0i,jpT and T are thus characteristic properties of the component with index p out of all the NCO components. The 1,p2,p

    objective is to perform a least squares fit of the function f to the data in the matrix RThis may be done by .

    minimizing the function: NSANSENSANSENCO((1/)(1/))TgTt1i2j22,,ppSSF(a,a,T,T){f(a,a,T,T)R}{a;aeR} (3) :::::012012,0,resijpijijijp again subject to the constraints given by the inequalities (2). This problem is ill-conditioned, with very minute

    changes in the data leading to widely different solutions.


    The minimization of the function given in eq. (3) can be reformulated as an unconstrained minimization by

    introducing a new set of variables a, , and so that: 0 2a,p1..NCOpp21/T,p1..NCO (4) 1,pp21/T,p1..NCO2,pp


For all real values of a, , and , minimizing F(a, , , ) will yield a solution within the permitted region of a, 00

    T and Twhich is defined by (2). If the variables a, , and are combined into one vector of variables x (of size 12 03NCO + 1) as follows:


    α (5) xβ


    then an inverse Hessian (or Gauss-Newton) algorithm will be based on the equation:

    1 (6) xxH(x)(x)k;1kkk

in which k is the index of the present iteration, x is the present position in the variable space, is the optimal scalar k

    step length, H(x) is the Hessian matrix of second derivatives of F at x, and (x) is the gradient of F at x (Note kkkk

    that the gradient in this context is the steepest ascent of a mathematical function, not the magnetic field gradient!). -1 is found using a suitable univariate minimization algorithm. If the search direction of -H(x)(x) is not a kk

    descending direction at the present location x, a steepest descent search along -(x) is performed. See Press et al. kk

    [12] for a detailed description of inverse Hessian algorithms. Analytical expressions of (x) and H(x) are given in kkan appendix.

3.1. The Univariate Minimization

    Each iteration in an inverse Hessian algoritm requires an univariate minimization along the search direction of --1 H(x)(x). In the present MATLAB [13] implementation of ANAHESS, this is done using the routine fminbnd, kkwhich is based on golden section search and parabolic interpolation [12]. In the present implementation, the optimal -20step length is always somewhere between a minimum of 10 and a maximum of 20. Also, if the solution strays outside of the interval defined by the minima and maxima of g and t, F is given a very high value.

3.2. Choice of the Initial Estimate For One Specific Choice of NCO

    Kaufmann [14] has developed a non-iterative algorithm for the fitting of a sum of exponentials to first order data.

    The goal of the first order sum of exponentials fit is to fit a sum of exponentials to a vector r so that any element

     r, i =1..N in r follows the function:i rf(a,a,b);ewith0iiiNCObv (7) pif(a,a,b)a;ae00ipp

in which vis the element having index i=1..N in the vector v. Applying Kaufmann’s algorithm on each row in R i 2provides an estimate of =-b, p=1..NCO. Likewise, Kaufmann’s algorithm applied on each column provides an pp222 estimate of =-b, p=1..NCO. Once the and estimates are available, estimates of a and a, p=1..NCO may pppp0pthen be found using the method of linear least squares.

    If a 2D-ILT calculation with NCO components is attempted, and the results from a calculation on the same data set

    with NCO-1 components are available, then the NCO-1 calculation results are compared to the initial estimates from Kaufmann’s algorithm, and taken into account when the NCO calculation is started.

3.3. Choice of Number of Components

    A good criterion for choosing the correct number of components NCO depends on the descriptive power of the model and the parsimony of the model. In this case, the descriptive power rests on the model’s ability to minimize the sum of squared residuals, while the parsimony is the model’s dependence on a low number of free parameters. When dealing with a demanding modelling problem, a model’s descriptive power and the parsimony typically pull

    in opposing directions. Improving one property usually diminishes the other.

One criterion for choosing a model is the Bayesian information criterion (BIC), introduced by Schwartz [15]. Let n

    be the number of observed data points, let p be the number of free model parameters, and ss be the sum of squared res

    residuals. Then if the residuals are normally distributed, BIC has the form:


    ssres (8) BICnln();pln(n)n

In this equation, a good model fit gives a low first term while few model parameters give a low second term. When

    comparing a set of models, the model with the minimal BIC value is selected. In the context of ANAHESS, n is

    NSA x NSE, SS is given in equation 3, and p is equal to 3NCO + 1, since a, a, T and T are the free parameters in res012the model. In ANAHESS, BIC is calculated for each successive value of NCO, and the model with the minimal BIC

    is chosen as the final model.

3.4. Summary of ANAHESS

    The elements outlined above may now be combined into the following algorithm:

I. Acquire the necessary input data: R, g and t.

    II. Choose a suitable maximum number of components NCO. max

    III. For NCO=1 to NCOdo (component loop): max

    a) Find an initial estimate x using Kaufmann’s algorithm along the rows and columns of R. 0

    b) If NCO>1 then compare x to the results from the previous calculation with NCO-1 components, and use 0the results to improve upon x if possible. 0

    While no convergence in F and x do (inverse Hessian loop): k

    (A) Calculate the various functions given by (9).

    (B) Calculate (x) according to (10). k(C) Calculate H(x) according to (11). k

    (D) Check H for positive definiteness. Make H positive definite if necessary.

    (E) Do a univariate optimization of in equation (6), and move from x to x. kk+1(F) Check for convergence in F and x. k

    c) End inverse Hessian loop.

    d) Extract a, T and Taccording to (4) and (5), and store NCO, a, a, T and T. 12 012

    e) Calculate BIC for the given NCO.

IV. End component loop.

At the end, this algorithm has produced a set of 2D-ILT solutions with NCO ranging from 1 to NCO, with one max

    solution for each NCO value. The solution with the minimal BIC is finally chosen as the optimal ANAHESS solution.


    Godefroy et al. [11] have developed a Tikhonov regularization MATLAB routine called twodlaplaceinverse.m [16],

    which is based an a strategy given by Lawson and Hanson [17], in which a regularization function or smoothing

    function is added to the data before minimizing the squared sum of errors. A singular value decomposition is then

    used to reduce data size. Instead of computing a specific set of NCO components with accompanying a, T and T 12

    vectors, the results are computed as a second order distribution and presented as a matrix of A values, with row and column indices corresponding to evenly spaced T and Tvalues. The weighing of the regularization function is 12

    determined by a parameter , chosen by the user. In this model, any element R, i =1..NSA, j =1..NSE in the i,jresponse surface R is assumed to follow the function:

    Rf(A,T,T);Ewithi,j12i,ji,jMN((1/T)g(1/T)t) (9) ij12,n,mf(A,T,T)Ae::12i,jm,nmn

    The regularized least squares function to be minimized now takes on the form:

    NSANSEMNMN11?1222;;SS(fR)(2AAA)(2AAA) (10) ::::::~?resijijmnm;nmnmnmn;mn,,,1,1,,,1,12ijmnmn22??


    2 The first term in this equation is a standard least squares term, while the second term headed by 1/αis a

    regularization term which adjusts the curvature of the solution in TTspace. The first term inside the 12

    regularization term covers the second derivatives along the rows of A, while the second term inside the

    regularization term covers the second derivatives along the columns of A. Adjusting α will thus move the solution

    towards the optimum leasts squares solution (large α) or towards smoothing the TTdistribution. The actual 12

    implementation by Godefroy et al. will be denoted the NNLS (nonnegative least squares) algorithm in the following algorithm comparison.

    Relevant criteria when comparing inverse Laplace transform algorithms are a good model fit, parsimony of the model and model robustness to noise in the data. In the present study, these criteria where studied as follows:

     Parsimony: Compare NCO in the ANAHESS model with minimal BIC with the number of elements in the

    A matrix in equation 15 for which A is not equal to zero. m,n

     Good model fit: Compare SS in equation 3 with SS in equation 10. resres

     Robustness to noise: Add noise to data sets and see how much and in what way the ANAHESS and NNLS

    results change.

    The 2D-ILT results from ANAHESS and NNLS are sparse nonnegative matrices or histograms. Such sparse structures are called signatures in image analysis. A commonly used metric when comparing signatures is the Earth mover’s distance (EMD), first introduced by French mathematician Gaspard Monge in 1781. The Earth mover’s

    distance between two matrices is based on the following intuitive idea: Consider each nonzero element in the two matrices as heaps of earth, the size of which is determined by the size of the element. The EMD is then the minimal

    amount of work (earth mass x element distance) needed to transform the first matrix into the second. If the calculated distance between matrix elements is a true metric, then EMD will be a true metric. EMD may be

    calculated using linear programming, see for instance Rubner et al. [18]. A counterintuitive property of EMD is that

    if large elements in one of the matrices cover smaller elements in the other, then their contribution to EMD will be

    zero. Therefore, in the following algorithm comparison, all the matrices to be compared were scaled so that the sum of all elements was unity.

4.1. A Comparison of Exponential Fit Algorithms on Experimental Data

    The NMR data in question were obtained from second order relaxation experiments done on a set of carbonate rock cores. The objective was to determine rock core properties such as pore sizes and porosity from NMR data. One experiment from this series is presented in table 1.

    Table 1. Settings used in a second order NMR relaxation experiment on a rock core.

     Response surface

    Core diameter 1.5 inches

    Temperature during experiment 308 K 1 Instrument provider Resonance2Laboratory location Anvendt teknologi

    Frequency 12 MHz

    g..g0.0010 .. 8.0 s minmax

    NSA 20

    t..t0.000212 .. 1.6415 s minmax

    NSE 8192

    1. Resonance Instruments, now Oxford Instruments , Tubney Woods, Abingdon, Oxfordshire, OX13 5QX, UK. 2. Anvendt Teknologi AS, Munkvollveien 56, N-7022 Trondheim, Norway. All in all, the set of plugs consisted of three dolomite samples and three calcite samples. The NMR measurements from each of the six samples were subjected to 2D-ILT analysis using both ANAHESS and NNLS. In the following, these calculations are labelled parallel 0. Then artificial noise was added to each of the raw data matrices, creating


    10 parallels labelled 1,2,..10 for each rock core sample. The noise had a normal distribution with a variance equal to 1% of the maximum measurement in the response surface. These parallels were then subjected to ANAHESS and NNLS calculations, and the results studied. 111213In the NNLS calculations, an A resolution of 40 x 20 was chosen. Also, a succession of α values of 10, 10, 10 141011121314and 10 or 10, 10, 10, 10 and 10 was applied, since these choices captured the flattening out of SS as a res

    function of α. In the ANAHESS calculations, the model with the lowest BIC was used consistently. Figure 1 shows

    BIC as a function of NCO for the first dolomite sample, parallel 0 (original data) and parallel 8 (the eighth parallel having added noise). The optimal model according to the BIC function is the one with NCO = 8 and NCO = 6,

    respectively. Figure 2 shows the two accompanying ANAHESS 2D-ILT plots. The results from all six samples are given in table 2. The corresponding NNLS calculations on parallel 0 and 8 are shown in figure 3. The NNLS calculations are summarized in table 3. Data for all the parallels, with NNLS and ANAHESS compared, are given in figures 4, 5 and 6.

    Figure 1. BIC as a function of NCO in ANAHESS. The

    sampleis the first dolomite, parallel 0 and parallel 8.


    Figure 2. ANAHESS plots from the first dolomite sample, parallel 0 (NCO = 8) and parallel

    8 (NCO = 6). The EMD between the two plots is 0.1406 as shown. Crosshair centers

    denote TT placement while circle greyscale denotes a value, as shown in sidebar. 12p

    Table 2. ANAHESS calculations on six rock core samples, with 1 original data

    set (parallel 0) and 10 parallels having added noise for each rock core sample.

    The EMD is between parallel 0 and the 10 parallels having added noise.

     9Sample NCO, parallel 0 SS/10, parallel 0 res

    Dolomite 1 8 0.1344

    Dolomite 2 9 0.1005

    Dolomite 3 6 0.1186

    Calcite 1 6 0.2447

    Calcite 2 6 0.1191

    Calcite 3 7 0.1290 9 NCO, excl. parallel 0, SS/10, excl. parallel 0, EMD from parallel 0, res

    Sample 10 parallels 10 parallels 10 parallels

    Lowest/mean/highest Lowest/mean/highest Lowest/mean/highest

    Dolomite 1 6 6.4 8 4.0406 4.0557 4.0768 0.0544 0.1455 0.1779

    Dolomite 2 4 5.8 6 3.6270 3.6377 3.6625 0.0539 0.0802 0.1532

    Dolomite 3 5 5.3 6 3.9528 3.9761 4.0047 0.0282 0.0473 0.0729

    Calcite 1 5 5.6 6 3.5286 3.5446 3.5936 0.0122 0.0505 0.2942

    Calcite 2 4 5.1 6 2.9609 2.9744 2.9900 0.0466 0.0713 0.1059

    Calcite 3 4 6.7 8 3.2239 3.2404 3.2681 0.0174 0.0750 0.1385


    Table 3. NNLS calculations on six rock core samples, with 1 original data set

     (parallel 0) and 10 parallels having added noise for each rock core sample .The EMD is between parallel 0 and the 10 parallels having added noise.

     9Sample α NCO, parallel 0 SS/10, parallel 0 res11 10 73 1.2200 12Dolomite 10 61 0.1864 131 10 45 0.1358 1410 36 0.1343 11 10 57 0.5853 12Dolomite 10 36 0.8705 132 10 28 0.1186 1410 26 0.1042 11 10 51 0.4992 12Dolomite 10 33 0.1358 133 10 22 0.1283 1410 20 0.1283 10 10 63 1.1509 11Calcite 10 47 0.2885 121 10 35 0.2746 1310 24 0.2726 1410 22 0.2726 10 10 67 1.3324 11Calcite 10 50 0.2767 122 10 36 0.1596 1310 30 0.1255 1410 27 0.1255 11 10 37 1.6079 12Calcite 10 29 0.2827 133 10 25 0.1346 1410 24 0.1346 9 NCO, excl. parallel 0, SS/10, excl. parallel 0, EMD from parallel 0, res 10 parallels 10 parallels 10 parallels

    Lowest/mean/highest Lowest/mean/highest Lowest/mean/highest 11 10 67 73.9 90 4.7218 5.2781 5.8704 0.0420 0.0881 0.1791 12Dolomite 10 48 56.3 68 4.1074 4.2683 4.6708 0.0753 0.1704 0.3087 131 10 41 45.1 53 4.0377 4.0696 4.1353 0.0646 0.2548 0.4997 1410 31 38 44 4.0377 4.0533 4.0717 0.1292 0.6085 1.3778 11 10 44 49.3 63 4.0732 4.4197 4.8994 0.0228 0.0540 0.0854 12Dolomite 10 30 34.4 42 3.8410 4.1235 4.6476 0.0265 0.0483 0.0712 132 10 25 29.7 36 3.6220 3.6394 3.6653 0.0260 0.0650 0.1610 1410 26 28.2 31 3.6220 3.6392 3.6649 0.0238 0.1450 0.5094 11 10 36 40.9 47 4.2904 4.4251 4.7601 0.0183 0.0305 0.0557 12Dolomite 10 24 28.7 35 3.9648 3.9982 4.0260 0.0218 0.0424 0.0798 133 10 20 22.8 27 3.9628 3.9866 4.0149 0.0219 0.0476 0.0816 1410 18 21.8 25 3.9628 3.9857 4.0147 0.0216 0.0479 0.0813 10 10 56 72.6 85 4.3539 4.4780 4.6374 0.0077 0.0290 0.0516 11Calcite 10 36 51.4 65 3.5685 3.5846 3.6105 0.0140 0.0410 0.0848 121 10 20 38.9 51 3.5570 3.5689 3.5874 0.0191 0.0573 0.1082 1310 20 29.7 36 3.5570 3.5685 3.5869 0.0344 0.1125 0.3290 1410 17 25.4 30 3.5569 3.5686 3.5869 0.0434 0.3757 1.0262 10 10 48 74.7 89 3.9542 4.1785 4.5342 0.0195 0.0422 0.0883 11Calcite 10 34 52.8 64 3.0497 3.1231 3.1725 0.0312 0.0579 0.1092 122 10 24 38 54 2.9747 2.9950 3.0233 0.0414 0.0957 0.2478 1310 20 31.2 47 2.9677 2.9804 2.9969 0.0572 0.1366 0.3474 1410 19 27.6 36 2.9676 2.9804 2.9969 0.0680 0.3452 0.7977 11 10 32 42.5 61 4.3295 4.8571 5.6074 0.0190 0.0302 0.0583 12Calcite 10 25 32.5 46 3.2501 3.3210 3.4492 0.0161 0.0467 0.0830 133 10 23 27.7 33 3.2298 3.2455 3.2657 0.0267 0.0627 0.1580 1410 23 26.3 32 3.2298 3.2453 3.2654 0.0417 0.1695 0.3860


Figure 3. 2D-ILT plots from the NNLS algorithm, with 4 settings of α, run on parallel 0

and parallel 8 from the first dolomite sample. Note the EMD comparisons plot by plot.


Figure 4. The number of components NCO as a function of α in NNLS, compared to

the results from ANAHESS on all six rock core samples, 11 parallels per sample.


Report this document

For any questions or suggestions please email