MUPen2DTool: a Matlab Tool for 2D Nuclear Magnetic Resonance relaxation data inversion

Accurate and efficient analysis of materials properties from Nuclear Magnetic Resonance (NMR) relaxation data requires robust and efficient inversion procedures. Despite the great variety of applications requiring to process two-dimensional NMR data (2DNMR), a few software tools are freely available. The aim of this paper is to present MUPen2DTool an open-source MATLAB based software tool for 2DNMR data inversion. The user can choose among several types of NMR experiments, and the software provides codes that can be used and extended easily. Furthermore, a MATLAB interface makes it easier to include users own data. The practical use is demonstrated in the reported examples of both synthetic and real NMR data. Introduction Nuclear Magnetic Resonance relaxometry of H nuclei ( H NMR) can give crucial information about the properties of many materials, ranging from cement [1] to biological tissues [2]. So, for example in case of porous media, H NMR permits the accurate estimate of important petrophysical parameters, such as porosity, saturation and permeability. For example, borehole H NMR is extensively used in oil and gas reservoir characterisation, and recent technological advances have led to tools suitable for environmental applications (see details in [3]). Usually the NMR parameters investigated are relaxation times (longitudinal T1 and/or transverse T2) and self-diffusion coefficient (D) as they are sensitive to the local physical environment and can also provide some chemical information. As there may be a range of the NMR parameter values that characterize a given system, in order to correctly interpret NMR experiment results and compute the distributions of these parameters, it is necessary to processes the experimental data with a robust and accurate inversion procedure. Recently, two-dimensional NMR (2DNMR) techniques are gaining increasing importance in analysing different porous media [4]. Moreover, new computationally intensive applications, such as multidimensional logging and general January 19, 2022 1/15 3DNMR data inversion [5], require efficient methods for the inversion of 2DNMR data. For all these reasons, there is an increasing request of software that can be easily applied to process 2DNMR data to compute 2D parameter distribution (NMR maps). Even if a considerable amount of inversion methods has been proposed in literature, the software tools implementing such methods are seldom freely available for testing. We strongly believe that open source software gives a determinant contribution to the progress of knowledge by making it possible for scholars to compare and improve their achievements. Therefore, starting from 2009, we released Upenwin [6] and recently Upen2dTool [7], an open-source software tool for 2DNMR data inversion. From a mathematical point of view, the problem of computing the two-dimensional relaxation time distributions from NMR data is a linear ill-posed problem modelled by a Fredholm integral equation with separable kernel. The strong ill-conditioning and the presence of data noise make the inverse problem very challenging. A regularisation technique is applied to reformulate the inversion problem as a non-negatively constrained optimisation problem, whose objective function contains a data fitting term and a regularisation term. In [8, 9] the Uniform Penalty (UPEN) principle, a multiple-parameters locally adapted Tikhonov-like regularization method, has been stated for one-dimensional NMR data and has been implemented in Upenwin software. Successively, in 2016, such principle has been extended for two-dimensions data [10], and further analysed and improved in [11, 12]. In 2019, Upen2dTool has been made available [7]. Both Upenwin and Upen2dTool implement a L2−norm locally adapted regularisation (in one and two dimensions, respectively), where the automatic computation of the regularisation parameters follows the UPEN principle. Although such methods compute very accurate distributions, their computational cost may be high since they require the solution of several non-negatively constrained least-squares problems. For this reason, a new method has been studied and proposed in [13] which represents a substantial change in the inversion strategy, consisting of adding an L1 penalty term to the locally adapted L2 term and removing the non-negativity constraint. The new method allows us to substantially improve the computation efficiency of the inversion process and for some kind of 2D data, to obtain even more robust and accurate NMR maps. We now release the Multiple Uniform Penalty 2D Tool (MUPen2DTool) open source software implementing the method proposed in [13]. MUPen2DTool consists of the source code, software documentation and a user guide which contains an installation guide, a technical description of synthetic NMR tests and input data format.MUPen2DTool comes also with a user friendly GUI (Graphical User Interface) that guides the user in the different steps of the inversion process. Moreover NMR data of several representative examples are available to help the interested user to assess the toolbox efficiency and effectiveness. The GUI makes it easy to handle the inversion parameters and inspect the loaded data, therefore our software can be flexibly used in the analysis of different types of samples. Being MUPen2DTool open source and user friendly, we believe that it has a large number of potential users from different application fields. In this paper we describe the software through a brief overview of the implemented algorithms and a detailed analysis of the 2D distributions computed from the data set enclosed in the software package. The paper has the following structure. We first introduce the general structure of MUPen2DTool describing its key features. Then, we present the problem of NMR data inversion and the characteristics of the implemented regularization algorithm. Finally, we report the software validation on several representative NMR relaxometry data relative to different types of samples. January 19, 2022 2/15 The problem of NMR data inversion In this section, we describe the mathematical model for NMR data inversion and the numerical scheme used by MUPen2DTool for its solution. The continuous model In MUPen2DTool, we consider 2DNMR maps corresponding to T1-T2, T2-T2 and T2-D relaxation data; in all these cases, the measured NMR signal is supposed to be related to an underlying distribution function by a Fredholm integral equation of the first kind. T1-T2 case. In a conventional Inversion-Recovery (IR) or Saturation Recovery (SR) experiment detected by a Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence [14], the relaxation data S(t1, t2) depending on t1, t2 evolution times can be expressed as: S(t1, t2) = ∫∫ ∞ 0 k1(t1, T1)k2(t2, T2)F (T1, T2)dT1dT2 + e(t1, t2) (1) where F (T1, T2) is the unknown distribution of T1 and T2 relaxation times and the kernels k1 and k2 have the expression k1(t1, T1) = { 1− 2 exp(−t1/T1), for IR sequence 1− exp(−t1/T1), for SR sequence and k2(t2, T2) = exp(−t2/T2). (2) Here and henceforth, the function e(·, ·) represents Gaussian additive noise. T2-T2 case. In a CPMG-CPMG experiment, the measured data S(t1, t2) is related to the underlying distribution F (T21, T22) by the integral equation S(t1, t2) = ∫∫ ∞ 0 k1(t1, T21)k2(t2, T22)F (T21, T22)dT21dT22 + e(t1, t2) (3) where both kernels k1 and k2 refer to transversal relaxation times T21, T22 and are defined as k1(t1, T21) = exp(−t1/T21), k2(t2, T22) = exp(−t2/T22). (4) Diffusion-T2 case. In a Stimulated Echo-CPMG experiment, the acquired echo amplitude S(t1, t2) can be expressed as S(t1, t2) = ∫∫ ∞ 0 k1(t1, T2)k2(t2, D)F (T2, D)dT2dD + e(t1, t2) (5) where the kernels k1 and k2 are k1(t1, T2) = exp(−t1/T2), k1(t2, D) = exp(−t2 ·D)). (6) The objective is to estimate the T1-T2, T2-T2 or D − T2 map from the measured data; this inversion is an ill-posed problem, which means that small noise in the data can cause significant changes in the computed 2D distribution. The discrete model The discretization of the linear integral equations (1), (3) and (5) leads to the liner system Kf + e = s (7) January 19, 2022 3/15 where K = K2 ⊗K1 is the Kronecker product of the discretized kernels K1 ∈ R11 and K2 ∈ R22 , the vector s ∈ R , M = M1 ·M2, represents the measured noisy signal, f ∈ R , N = N1 ·N2, is the vector reordering of the 2D distribution to be computed and e ∈ R represents the additive Gaussian noise. The linear problem (7) is typically ill-conditioned and regularization strategies are necessary to obtain stable discrete distributions. A recent review of regularization techniques used for 2DNMR data inversion can be found in [13]. The minimization problem MUPen2DTool uses a multipenalty approach based on both L1 and L2 regularization with locally adapted regularization parameters. In this regularization framework, the NMR data inversion problem is reformulated as the unconstrained minimization problem min f {