Study on Multi-Exponential Inversion Method for NMR Relaxation Signals with Tikhonov Segularization

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 regularization 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 results of NMR studies on the liquid samples like as oils and yoghurt.


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 unparalleled to others, it is a very useful tool for understanding 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 relaxometry 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 algorithm 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 technique) methods which can solve the problem of multiexponential inversion at a different 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-noise ratio will affect the result [5][6][7][8].This paper updates the application of Tikhonov regularization to stabilize this numerical inversion problem and demonstrates the method for automatically choosing the optimal value of the regularization parameter.By use 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 algorithm has rapid calculation speed, high stability and credibility.

Basic Theory
The NMR relaxation measurements, either longitudinal (T 1 ) or transverse relaxation (T 2 ), can be represented in the form of well known Fredholm integral equation of first kind.The T 2 relaxation signal measured by a simple CPMG spin echo method, for example, can be written as follows after making it discrete: where ( ) i y t is echo signal raw data; x j is the contribution 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 T 2 .For T 1 spectrum algorithms, we substitute the basis vector t The objective of inversion is to extract x j vs. T 2j from (1).Equation ( 1) is expressed as the form: Y= A•X (2) where the matrix Y n×1 = (y 1 , y 2 , …, y n ) T is the measured relaxation signal, the matrix X m×1 = (x 1 , x 2 , …, x m ) T is the amplitude to be inversed with T 2j time, and the coefficient matrix ( ) The T 2j pre-assigned relaxation time constants in T 2 distribution range, which are m bins chosen with logarithms equally spaced between the minimal and maximal relaxation times (T 2min and T 2max ).We must take into account the problem that the abscissa of T2 is more than three orders of magnitude (T 2min = 0.1 ms and T 2max = 10,000 ms).The signal will attenuate to 0 over time and the T 2 component will reduce too.Thus we choose with logarithms equally spaced between the min and max times.
Equation ( 2) is essentially an ill-conditioned deconvolution in mathematics, which is over-determined, determined 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 T 2 spectrum.Several physically convincing constraints like nonnegativity 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 Tikhonov 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: The formal solution to which can be written as: 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 regularization parameter λ.The appropriate choice of this parameter is therefore very important [9,10].
Then, we find the solution of (4) when λ = 0.Meanwhile, we obtain the chi-square which can be defined as follows: Here we assume each point measurement error is random, 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 good inversion result.
The process of the algorithm procedure is shown as Figure 1.

Methods
For algorithm analysis, we synthesize data from a known model of T 2 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- son with the true model.Here we configure a typical double-peak T 2 spectrum as a model, and then we obtain the noiseless echo string by forward modeling.The algorithm 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 spectrum, we can test the accuracy of algorithm.Because noise is one of the important factors affecting the inversion precision, SNR effect is studied by creating the simulated 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 distributions.Meanwhile 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 T 2 spectra is more than three orders of magnitude (T 2min = 0.1 ms and T 2max = 10000 ms).As time increases, the signal will attenuate to none and the T 2 component will reduce too.Thus we choose with logarithms equally spaced from 0.1 ms to 10000 ms [11][12][13].

Simulated Ideal Model
Simulated data were created for two components T 2 distribution with value of 5 ms and 100 ms.The algorithm is tested with simulated data, the figure below shows the inversion results.
Comparing with the configured spectrum from Figure 2, we can see that inversion spectrum is in conformance with the configured one.So the result authenticates that our method is feasible.
Figure 3 shows the procedure of regularization parameter selection.Here the final regularization parameter value is 0.06 while chi-square value is 1017 approximately equal to the number of echo strings.

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 inversion results 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 than 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.

Experimental Validation
In order to test the feasibility of the procedure, we collected 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 yoghurt and another soybean oil.Meanwhile we make the T 2 spectrums derived from the Niumag inversion software as referenced spectrums.Through running the new algorithm in MATLAB with the measured three echo trains gathered from the PQ-001 NMR analyzer, 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 results 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.

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 inputting them into the algorithm.Here, we use the following two methods for data processing: Firstly, the front data of echo train contains short T 2 relaxation information, and its frequency is high while the followed data include long T 2 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     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 measured signal data should minus the baseline value.Although a lot of efforts are being spent on improving those methods, the efficient and effective inversion algorithm has yet to be developed.

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.Example analyses are presented using both synthetic data and experimental results of NMR studies on the liquid samples 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 distributions and help us reduce the uncertainty in the estimated distribution.Meanwhile, this study also provides guidance to get a reliable NMR T 2 inversion distribution in the field of MRI research and test.

Figure 1 .
Figure 1.The flow diagram of the regularization inversion.

Figure 2 .
Figure 2. T 2 inversion spectrum by the algorithm (circle) and the configured spectrum (cross).

Figure 3 .
Figure 3.The most reasonable way to choose the regularization parameter.

Figure 5 .
Figure 5. T 2 inversion spectrum (point) and the referenced spectrum (circle) obtained by pure vegetable oil added a few drops of water.

Figure 7 .
Figure 7. T 2 inversion spectrum (point) and the referenced spectrum (circle) obtained by soybean oil.