Engineering, 2013, 5, 32-37
http://dx.doi.org/10.4236/eng.2013.510B007 Published Online October 2013 (http://www.scirp.org/journal/eng)
Copyright © 2013 SciRes. ENG
Study on Multi-Exponential Inversion Method for NMR
Relaxation Signals with Tikhonov Segularization
Shanshan Chen1, Hongzhi Wang1,2, Xuelong Zhang1,2
1Shanghai Medical Instrumentation College, Shanghai, China
2The Laboratory of Medical Imaging Engineering, University of Shanghai for Science and Technology,
Shanghai, China
Email: shanshan_1987917@163.com
Received September 2012
ABSTRACT
The analysis of NMR data in terms of inversion of relaxation distribution is hampered by the ill-posed nature of the
required solution about the Fredholm integral equation. Naive approaches such as multi-exponential fitting or standard
least-squares algorithms are numerically unstable and often failed. This paper updates the application of Tikhonov re-
gularization to stabilize this numerical inversion problem and demonstrates the method for automatically choosing the
optimal value of the regularization parameter. The approach is computationally efficient and easy to implement using
standard matrix algebra techniques, which is based on mathematical ware MATLAB. Example analyses are presented
using both synthetic data and experimental resu lts of NMR studies on the liqu id samples like as oils and yoghurt.
Keywords: NMR; Relaxometry; Inversion Relaxation; MATLAB
1. Introduction
Low-field (LF) NMR (Nuclear Magnetic Resonance)
relaxometry is a new technology, which has a very big
development potential. It recognizes the function of
qualitative and quantificational for the oil and moisture
content in samples by the use of analysing the hydrogen
proton relaxation [1].
LF-NMR relaxometry in the current technology is un-
paralleled to others, it is a very useful tool for under-
standing various chemical and physical phenomena in
complex multiphase systems [2,3]. It has a very broad
application prospect in agriculture, biological materials,
food security, life science, pharmaceutical industry and
petrochemical industry. Nevertheless, LF-NMR relaxo-
metry technology is a kind of indirect measurement
technology. To obtain parameters that reflected sample
attribute information, an essential step is the inversion
process which is called NMR multi-exponential fit for T1
or T2 distributions. So it is still a gordian technique to
study how to build a fast and practical inversion algo-
rithm with a higher resolution versus relaxation times
and less dependence on the measured SNR [4].
At present, there are many inversion algorithms such as
SVD (singular value decomposition), NNLS (nonnegative
least squares), and SIRT (solid iteration rebuilding tech-
nique) methods which can solve the problem of multi-
exponenti a l i nve rs ion at a di fferent angle.
Some of those ways are not perfect: high requirement
of echo signal-to-noise ratio, take too long, relatively
complicated process of nonnegativity restrictions, and at
the same time the inversion algorithm parameters such as
the number of observation data, the way of stationing and
the signal-to-nois e ratio will affect the result [5-8].
This paper updates the application of Tikhonov regu-
larization to stabilize this numerical inversion problem
and demonstrates the method for automatically choosing
the optimal value of the regularizatio n parameter. By u se
of numerical simulation with MATLAB, we come up
with a LF nuclear magnetic resonance (NMR) inversion
algorithm. Example analyses are presented using both
synthetic data and experimental results of NMR studies
on the liquid samples like as oils and yoghurt. The actual
calculation indicates that the algor ithm has rapid calcula-
tion speed, high stability and credibility.
2. Basic Theory
The NMR relaxation measurements, either longitudinal
(T1) or transverse relaxation (T2), can be represented in
the form of well known Fredholm integral equation of
first kind. The T2 relax ation signal measured by a simple
CPMG spin echo method, for example, can be written as
follows after making it d iscrete:
2
1
( )(t)
i
j
t
mT
ij
j
ytx e
ε
=
= ⋅+
(1)
S. S. CHEN ET AL.
Copyright © 2013 SciRes. ENG
33
where
()
i
yt
is echo signal raw data; xj is the contribu-
tion of the ith relaxation unit to the total signal;
ε
(t) is
random noise.
The algorithms and simulations in this article are
based on the inversion of T2. For T1 spectrum algorithms,
we substitute the basis vector
1
[12exp(/ T)]t
for basis
vector
2
[exp(/ T)]t
The objective of inversion is to extract xj vs. T2j from
(1). Equation (1) is expressed as the form:
Y= AX (2)
where the matrix Yn×1= (y1, y2, …, yn)T is the measured
relaxation signal, the matrix X1 = (x1, x2, …, xm)T is the
amplitude to be inversed with T2j time, and the coeffi-
cient matrix
( )
n miji2j
nm nm
Aaexp Tt
×××


== −
 
. The
T2j pre-assigned relaxation time constants in T2 distribu-
tion range, which are m bins chosen with logarithms
equally spaced between the minimal and maximal relax-
ation times (T2min and T2max). We must take into account
the problem that the abscissa of T2 is more than three
orders of magnitude (T2min = 0.1 ms and T2max = 10,000
ms). The signal will attenuate to 0 over time and the T2
component will reduce too. Thus we choose with loga-
rithms equally spaced between the min and max times.
Equation (2) is essentially an ill-conditioned deconvo-
lution in mathematics, which is over-determined, deter-
mined or underdetermined, so an iteration solution has
been employed. Like any other inversion problem case,
those processes have the drawback of non-uniqueness. A
small change in a single measurement can lead to an
enormous change in the estimated models. The same raw
data can satisfy several different model solutions very
nicely. The SVD, NNLS, and SIRT algorithms can solve
the problem of multi-exponential inversion of T2 spec-
trum. Several physically convincing constraints like non-
negativity are also introduced to dampen the unnecessary
oscillations in the solutions. Some of above ways are not
perfect, in the paper we updates the application of Tik-
honov regularization to stabilize this numerical inversion
problem. The implementation steps for the new method
are as follows:
Firstly, the approach uses Tikhonov regularization to
stabilise the solution to (2) by imposing a penalty on the
recovered solution. The nature of this penalty can be
chosen by the user. The regularized solution is obtained
by the following minimisation:
22
2
min( )y AxLx
λ
−+
(3)
The formal solut i on to which c a n be wr i tten as:
T2T1 T
()xAALLAy
λ
λ
= +
(4)
The matrix L allows a priori expectations about the
data to be included such as requiring the solution to be
smooth. Here this matrix is chosen as the identity matrix,
giving preference to solutions with smaller norms. The
second term in (3) serves to impose some constraint or
prior knowledge on the size of the returned distribution,
with the power of this term controlled by the regulariza-
tion parameter λ. The appropriate choice of this parame-
ter is therefore very important [9,10].
Then, we find the solution of (4) when λ = 0. Mean-
while, we obtain the chi-square which can be defined as
follows:
(5)
Here we assume each point measurement error is ran-
dom, independent and standard deviation is all the same.
So, considering the smoothness and the non-negativity of
the result we observe chi-square value within the scope
of the change of λ, the linear change of λ from zero make
the chi-square value will be changed too. Then we plug
the λ which is obtained by interpolation method that can
make the chi square value close to echo number into (4),
out popped the relationship between relaxation time
component and the echo of the corresponding amplitude.
This algorithm has a certain degree of randomness and
subjectivity. Through trial and error we find the most
reasonable way to choose the regularization parameter to
obtain a go od invers i on result.
The process of the algorithm procedure is shown as
Figure 1.
3. Methods
For algorithm analysis, we synthesize data from a known
model of T2 distribution and then try to calculate back the
true model using inversion algorithms. It will provide an
estimate of accuracy of the algorithm through compari-
Figure 1. The flow diagram of the regularization inversion.
S. S. CHEN ET AL.
Copyright © 2013 SciRes. ENG
34
son with the true model. Here we configure a typical
double-peak T2 spectrum as a model, and then we obtain
the noiseless echo string by forward modeling. The algo-
rithm is developed and the codes are written in Matlab,
you can obtain the codes by requesting the first author.
The algorithms can take the users inputs : number of echo
string, number of pre-assigned relaxation time constants,
SNR of the measured data.
Those procedures can use the noiseless echoes as input
signals, and relaxation inversion spectrums are the output
of the algorithms. Compared with the configured spec-
trum, we can test the accuracy of algorithm. Because
noise is one of the important factors affecting the inver-
sion precision, SNR effect is studied by creating the si-
mulated data which is added Gaussian white noise with
SNR level from 10 through 100. Those algorithms are
run on the simulated data to obtain the computed distri-
butions. Meanw hile for experimental validation of the
methods, we collected three groups of spin-echo signals
produced by yoghurt and oils in PQ-001 NMR analyzer
in our laboratory We must take into account the problem
that the abscissa of T2 spectra is more than three orders
of magnitude (T2min = 0.1 ms and T2max = 10000 ms). As
time increases, the signal will attenuate to none and the
T2 component will reduce too. Thus we choose with lo-
garithms equally spaced from 0.1 ms to 10000 ms [11-
13].
4. Results
4.1. Simulated Ideal Model
Simulated data were created for two components T2 dis-
tribution with value of 5 ms and 100 ms. The algorithm
is tested with simulated data, the figure below shows the
inversi o n re s ults.
Comparing with the configured spectrum from Fig ure
2, we can see that inversion spectrum is in conformance
Figure 2. T2 inversion spectrum by the algorithm (circle)
and the confi gu red spectrum (cross).
with the configured one. So the result authenticates that
our method is feasible.
Figure 3 shows the procedure of regularization para-
meter selection. Here the final regularization parameter
value is 0.06 while chi-square value is 1017 approxi-
mately equal to the number of echo strings.
4.2. Impact of Noise
In order to facilitate analysis, through adding a random
white Gaussian noise point by point to the simulated
CPMG echo trains, then input six echo trains into the
three procedures with SNR = , 100, 50, 30, 20, 10. The
inversi o n re s ults are shown in Figure 4.
Inversion spectrum with different SNR is shown in the
Figure 3. We can see that inversion algorithm results
will be ideal only if the SNR is larger th an 30, and SIRT
is 20. With low SNR (SNR = 10) the spectrum results
show that big peak of real spectrum can be inversed
while the small peak is far close to the real spectrum
structure. From the point of the program runtime, in ideal
model, the inversion algorithm takes 3 s.
4.3. Experimental Validation
In order to test the feasibility of the procedure, we col-
lected three groups of spin-echo signals in the PQ-001
NMR analyzer in our laboratory. One is produced by
pure vegetable oil added a few drops of water, one yog-
hurt and another soybean oil. Meanwhile we make the T2
spectrums derived from the Niumag inversion software
as referenced spectrums. Through running the new algo-
rithm in MATLAB with the measured three echo trains
gathered from the PQ-001 NMR analyz er, the results are
as shown in Figures 5-7 respectively.
Figures 5-7 show those three algorithms inversion
spectrum for yoghurt sample compared with the real
spectrum obtained by the PQ-001. The experimental re-
sults obtained are roughly the same as the referenced
ones. The feasibility and the stability of the algorithm
can meet the needs of the MRI research and test.
4.4. Problems
The broadening of computed peak is observed because
that echo trains collected from PQ-001 are corrupt by
noise. It is necessary to deal with raw data before input-
ting them into the algorithm. Here, we use the following
two methods for data processing: Firstly, the front data of
echo train contains short T2 relaxation information, and
its frequency is high while the followed data include long
T2 information with low frequency. We can use the index
filter method and at the same time the window width can
be adjustable because signal is decaying with time as a
simple exponential function of time; Secondly, through
experiments, we can see that the attenuation of relaxation
S. S. CHEN ET AL.
Copyright © 2013 SciRes. ENG
35
Figure 3. The most reasonable way to choose the regularization parameter.
Figure 4. T2 inversion spectrum after adding the white Gaussian noise with different SNR (, 100, 50, 30, 20, 10) by the new
algorithm.
Figure 5. T2 inversion spectrum (point) and the referenced spectrum (circle) obtained by pure vegetable oil added a few drops
of water.
S. S. CHEN ET AL.
Copyright © 2013 SciRes. ENG
36
Figure 6. T2 inversion spectrum (point) and the referenced spectrum (circle) by yoghurt.
Figure 7. T2 inversion spectrum (point) and the referenced spectrum (circle) obtained by soybean oil.
signals terminated in about 0.5 s, after that the signal
amplitude became very small. So we can average the
value of the signal after the time 0.5 s, then the noise
signal can almost be cancelled and the rest of relaxation
signal should be the baseline information. In order to
increase SNR of the original echo train, all of the meas-
ured signal data should minus the baseline value. Al-
though a lot of efforts are being spent on improving those
methods, the efficient and effective inversion algorithm
has yet to be developed.
5. Discussion
In this paper, the application of Tikhonov regularization
is updated to stabilize this numerical inversion problem.
We demonstrate the method for automatically choosing
the optimal value of the regularization parameter. Exam-
ple analyses are presented using both synthetic data and
experimental results of NMR studies on the liquid sam-
ples like as oils and yoghurt. Through simulated models
and experiment validation, the feasibility, resolution and
antinoise ability of the new algorithm are compared. It
enabled us to quantify the accuracy of the computed dis-
tributions and help us reduce the uncertainty in the esti-
mated distribution. Meanwhile, this study also provides
guidance to get a reliable NMR T2 inversion distribution
in the field of MRI research and test.
S. S. CHEN ET AL.
Copyright © 2013 SciRes. ENG
37
6. Acknowledgements
We wish to thank the engineer Yang Peiqiang for discus-
sions and help. The authors would like to acknowledge
the support from Shanghai Niumag Corporation.
This work is supported by the young college teachers
training subsidy scheme fund of Shanghai (2012-
slg12034), and the scientific research initial funding of
SMIC (A25001302-14).
REFERENCES
[1] L. z. Xiao, “The Principle and Application of NMR Im-
aging Logging and Rock NMR,” Science Press, Beijing,
1998.
[2] G. X. Xiong and L. B. Li, “The Principle of Magnetic
Resonance Imaging,” Science Press, Beijing, 2007.
[3] H. Z. Wang and X. L. Zhan, “Experimentation of Mag-
netic Resonance Imaging,” Science Press, Beijing, 2008.
[4] L. Z. Xiao, Z. D. Wang and T. Y. Liu, “Application of
Multi-Exponential Inversion Method to NMR Measure-
ments,Petroleum Science, Vol. 1, No. 1, 2004, pp. 20-
22.
[5] H. Wang and G. Y. Li, “Combination of Inversion and
Fitting as an Effective Method for the Analysis of NMR
Relaxation Data,” Acta Physica Sinica, Vol. 54, No. 3,
2005, pp. 1431-1436.
[6] A. H. Weng and Z. B. Li, “On High Rresolution Inversion
of NMR Logging Data,” Well Logging Technology, Vol.
26, No. 6, 2002, pp. 455-459.
[7] W. M. Wang, P. Li and C. H. Ye, “Multi-Exponential
Inversion of NMR Relaxation Signal,” Science of China
(A), Vol. 38, No. 8, 2001, pp. 730-736.
[8] R. D. Jiang, Y. P. Yao and S. Miao, “Improved Algorithm
for Singular Alue Decomposition Inversion of T2 Spec-
trum in Nuclear Magnetic Resonance,” Acta Petrolei Si-
nica, Vol. 26, No. 6, 2005, pp. 57-59.
[9] K. P. Whittall and A. L. Mac Kay, “Quantitative Interpre-
tation of NMR Relaxation Data,” Journal of Magnetic
Resonance, Vol. 84, 1989, pp. 134-152.
[10] G. C. Borgia, R. J. S. Brown and P. Fantazzini, “Uniform
Penalty Invertion of Multi-Exponential Decay Data (II),”
Journal of Magnetic Resonance, Vol. 147, 2000, pp. 273-
285. http://dx.doi.org/10.1006/jmre.2000.2197
[11] Y. Q. Song and L. P. Robert, “Assigning Uncertainties in
the Inversion of NMR Relaxation Data,” Journal of Mag-
netic Resonance, Vol. 174, 2005, pp. 314-324.
http://dx.doi.org/10.1016/j.jmr.2005.03.002
[12] S. Ghodh, M. K Ke nvin and Y. Pan, “A Simulation Based
Method to Assess Inversion Algorithms for Transverse
Relaxation Data,” Journal of Magnetic Resonance, Vol.
191, 2008, pp. 226-230.
http://dx.doi.org/10.1016/j.jmr.2007.12.021
[13] F. Lin, Z. W. Wang and J. Y. Li, “Study on Algorithms of
Low SNR Inversion of T2 Spectrum in NMR,” Applied
Geophysics, Vol. 3, 2011, pp. 233-238.