A Novel Method to Enhance the Inversion Speed and Precision of the NMR T₂ Spectrum by the TSVD Based Linearized Bregman Iteration

1  Introduction

The time domain signal determined by the low field nuclear magnetic resonance (NMR) can be used to detect many features of the porous rock, such as the pore size distribution, pore connectivity, wettability, viscosity, as well as fluid saturation [1–4]. The decaying signals are inverted to the time domain spectrums, which are frequently related to the formation evaluation and petrophysical study. However, the inversion of the NMR data is usually ill-posed, resulting in significant difficulty in the computation precision and efficiency. To deal with this problem, many attempts and developments have been tried in recent decades, which can be summarized into three categories: iterative algorithms, regularization algorithms and intelligent algorithms [5–7]. However, there are many drawbacks for these existing algorithms. Firstly, for the iterative algorithms, such as the truncated singular value decomposition (TSVD), the nonnegative least squares (NNLS) and the simultaneous iterative reconstruction technique (SIRT), the initial guess and the iteration criterion are hard to be determined [8–10]. Moreover, some unwanted signals are introduced during the iteration process, leading to the misinterpretation of the results [11]. Secondly, for the regularized algorithms, such as the Butler-Reeds-Dawson (BRD) algorithm and the maximum entropy algorithm, the penalty term and the regularization parameters are difficult to estimate [12,13]. In addition, it has the risk of trapping in local minimal values. Thirdly, although some intelligent algorithms, such as the simulated annealing and the differential evolution can reach favorable inversion result at certain conditions, they are immature and not used massively. In recent years, some researchers utilized the L1 regularization to invert the raw decaying data to spectrums and achieved favorable results [14–17].

The Bregman iteration is considered as an effective method for solving constrained optimization problems. It has been widely applied in many fields such as seismic prospecting, image reconstruction, as well as noise reduction [18–21]. The best advantage of the Bregman iteration is the easy convergence [20]. Therefore, it is introduced into the NMR inversion, aiming to get the stable solution at a favorable iteration time. In this paper, the novel method composed by the combination of improved TSVD and the linearized Bregman iteration is investigated. And then, it is validated by comparing its results with those of other previous algorithms.

2  NMR Relaxation Theory

Based on the basic principle NMR theory, the longitudinal relaxation time (T1) and the transversal relaxation time (T2) are used to depict the polarizing and decaying behaviors during the nuclear precession process, which are given by [1,2]:

1T1=1T1b+1T1s=1T1b+ρ1sSV(1)

1T2=1T2b+1T2s+1T2d=1T2b+ρ2sSV+1T2d(2)

where T1b and T2b are the bulk relaxation times of T1 and T2, determined by the intrinsic properties of the fluids such as viscosity and chemical compositions; T1s and T2s are the surface relaxation times of T1 and T2, controlled by the properties of fluid-rock interface such as pore size and surface relaxivity; T2d is the diffusional relaxation time of T2, affected by the diffusion coefficient, temperature and pressure; S/V is the pore surface to volume ratio; ρ1s and ρ2s are the surface relaxivities of T1 and T2, respectively [2].

T2 is the most frequently measured parameter in conventional low-field NMR experiments and logging, which are obtained by Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence. Therefore, only the inversion problem of T2 is examined in this study. The discrete relationship between the magnetization intensity and T2 satisfies the following equation:

f(t)=f0∑i=1naiexp⁡(−t/T2i)+ε(t).(3)

where ai is the volume fraction of the i-th relaxation component decaying with relaxation time of T2; f0 is the initial magnetization intensity; t is the evolution time; n is the number of T2, which is always predefined as the exponent of 2, such as 32, 64, and 128; ε(t) is the recorded noise.

In the conventional inversion algorithms, Eq. (3) is expressed as the optimization form with nonnegative constraints.

min||F−AX||2s.t.x(j)≥0forj=1,2,…,n(4)

where F is recorded magnetic intensity with the size of m×1, A is the kernel matrix with the size of m×n, and X is the T2 with the size of n×1.

Since A is ill-conditioned, the optimization problem is an ill-posed problem. The solution with the minimal fitting residual is easily contaminated by noise.

3  Methodology

It is often the case that most elements in the spectrum are zeros, so L1 penalty term is added into the function to enforce the sparsity of the spectrum. Therefore, the cost function is expressed as [18–21]:

minx(j)≥012||F−AX||22+λ||X||1(5)

where ||X||1 denotes the one-norm of X; λ is the regularization factor, which provides a tradeoff between fidelity-to-data and the noise sensitivity.

3.1 The Linearized Bregman Iteration

The Bregman iteration is very popular in many optimization problems due to its very nice convergence properties, including monotonic decrease in the residual term, convergence to the original signal, and convergence in terms of Bregman distance to the original signal with noisy data [22].

Assuming J(x) is a differentiable convex function (Bregman function), the Bregman distance between two points is given by [18–22]:

DJ(x,y)=J(x)−J(y)−<p,x−y>.(6)

where p∈∂J(x) is the subgradient of J(x).

The Bregman distance has several nice properties that make it an efficient tool for solving L1 regularization problems. Using the Bregman function, the optimization problem is transformed as:

minx(j)≥0{J(x)+H(x,f)}(7)

where J(x)=λ||X||1 and H(x,f)=12||F−AX||22.

Consequently, the key steps of the Bregman iteration are expressed as follows:

xk+1=argmin{DJpk(x,xk)+H(x,f)}(8)

pk+1=pk−∇H(xk+1,f)∈∂J(xk+1)(9)

where k denotes the iteration time.

The main function of Eq. (8) is to update the iterative solution, and the main function of Eq. (9) is to update the search direction for a new iteration. Due to the complexity of the computation of the subgradient, the simplified iteration operators are given by [23]:

fk+1=f−(Axk−fk)(10)

xk+1=minJ(x)+12||Ax−fk+1||22(11)

Then the problem is transformed to find the optimal values to satisfy Eqs. (8) or (11), which can be computationally expensive. Here we adopt the linearized version of the Bregman iteration to accelerate the computation since the minimization of Eq. (8) can be replaced by a minimization step that can be solved exactly.

The first order Taylor expansion of H(x,f) is expressed as:

H(x,f)=H(xk,f)+<∇H(xk,f),xk−x>(12)

By adding the penalty term of 12δ||x−xk||22 , the iteration in Eq. (8) is then transformed as:

xk+1=argminDJpk(x,xk)+H(xk,f)+<∇H(xk,f),xk−x>+12δ||x−xk||22(13)

Substituting Eq. (4) into Eq. (13), the iteration form is derived as:

xk+1=argminDJpk(x,xk)+12δ||x−(xk−1δAT(Axk−f))||22(14)

Eq. (9) is equivalent to:

pk+1=pk−δ(xk+1−(xk−1δAT(Axk−f)))(15)

Combining vk+1=∑j=0kAT(f−Axj) and Eq. (14) gives:

xik+1={1δ(vik−λ)vik∈(λ,+∞)0vik∈[−λ,λ]1δ(vik+λ)vik∈(−∞,−λ)(16)

Hence, the implementation of the linearized Bregman algorithm is as follows: (1) Initialize the values of x and v, iteration time of k and the maximal iteration time T, and the deviation of the stop criterion ξ; (2) Do the iterations of Eqs. (15) and (16) when ||Axk−f||>ξ or k<T; (3) Break the iteration by the stop criterion, and output the solution.

3.2 The Empirical TSVD

The above section gives the full part of the linearized Bregman iteration. However, direct application of this algorithm is impossible since the ill-conditioned kernel matrix A. So some pretreatments are necessary to ensure suitability. SVD is the effective method to decompose the ill-conditioned matrix, which is expressed as [23–25]:

A=USVT(17)

where U=[u1,u2,…,um] and V=[v1,v2,…,vn] are orthogonal matrices with the size of m×m and n×n; S=diag(s1,s2,…,sr) is the diagonal matrix with the size of m×n, and r=min(m,n).

Singular values of S are ordered as:

s1≥s2≥⋅⋅⋅≥sr=0(18)

It should be noted that singular values decay quickly to zero, and small singular values are often considered harmful components since they will result in the large error of the solutions. Conventional methods are to set a threshold and drop values smaller than the threshold such as the generalized cross validation (GCV), L curve, and signal to noise ratio (SNR) based methods [26].

Through many times of simulations trials, an empirical equation is established to predict the proper truncated position and can be expressed by:

q=a×SNRb(19)

where a and b are fitting parameters. In this study, they are 2.869 and 0.438, respectively.

Therefore, the kernel matrix is expressed as [27,28]:

Aq=UqSqVqT(20)

where Uq=(u1,u2,…,uq), Sq=(s1,s2,…,sq) and Vq=(v1,v1,…,vq).

Therefore, the iteration function for Eq. (15) is expressed as:

pk+1=pk−δ(xk+1−(xk−1δ(xk−VqSq−1UqTf))(21)

The space complexity, time complexity, and computational complexity for Aq are 3q2, 4q3, and O(q3), respectively. Hence, the computational complexity is greatly decreased by O(m2n)/O(q3) since q is much smaller than m and n. The flowchart of the proposed inversion method is shown in Fig. 1.

Figure 1: The flowchart of the proposed inversion algorithm

4  Numerical Examples

The simulation tests are conducted to verify the proposed novel inversion method using two typical distributions of the T2 spectrum, including the unimodal spectrum, the bimodal spectrum. Moreover, numerical data with different SNRs will be used to investigate the noise tolerance of this method. In our simulation, the noise is added using the in-built ‘awgn’ function of the Matlab software. It is typical additive Gaussian white noise. The SNR is specified with dB. All the simulations are conducted on a computer labeled the ‘OptiPlex 7050’ with the Intel(R) Core(TM) i7-7700 CPU @ 3.60 GHz. In the simulation, the echo spacing (TE) and the number of echoes are set as 0.2 ms and 8000, respectively. After many runs of simulations, λ and δ are found to be 0.01 and 0.5, respectively. The maximal iterations are set as 100. The total duration for each inversion is only about 1.2 s, which is suitable for both laboratory investigations and field applications. For comparison, the conventionally used method named ‘BRD’ is also used.

4.1 The Unimodal Spectrum

The Gaussian’s function is used to establish the forward unimodal distributed T2 spectrum with a peak of 10 ms and the decaying signals with SNR of 100, 50, 30, 20, 10, and 5, respectively. It should be noted that for simplicity, the normalized amplitude less than 0.003 is dropped and the spectrum is then renormalized. The simulated echo trains are shown in Fig. 2. The inverted spectrums are presented in Fig. 3. It can be seen that for the echo with SNR higher than 20, the inverted spectrum is strongly consistent with the forward model. However, if the SNR is lower than 10, the inverted results deviate from the model, but the peak position is only slightly moved. The inversion inaccuracy can be attributed to the contamination of the noise. However, the inversion results obtained by the BRD method deviate from the true value obviously at the SNR lower than 20.

Figure 2: Simulated NMR trains of the unimodal spectrum with different SNRs

Figure 3: The inverted unimodal spectrums of the signals with different SNRs

4.2 The Bimodal Spectrum

Considering the most common case, the bimodally distributed spectrum is constructed. The positions of the two peaks are 10 ms and 60. Similar to the previous case, Fig. 4 shows the simulated Gaussian’s noises with SNR of 100, 50, 30, 20, 10, and 5. The normalized amplitude of less than 0.001 is dropped and the data is renormalized. The inversion results are shown in Fig. 5. The results show that the inverted spectrum is coherent with the model for echo trains with SNR higher than 30. For echo trains with SNR lower than 20, the deviation between the inverted spectrum and the model increases with the decrease of the SNR. The inversion results are unacceptable when the SNR is less than 10, since the position and the shape are totally changed. The inversion results obtained by the BRD algorithm are similar to the unimodal spectrum.

Figure 4: Simulated NMR trains of the bimodal spectrum with different SNRs

Figure 5: The inverted bimodal spectrums of the signals with different SNRs

Based on the above simulation results, it is concluded that the proposed novel algorithm works very well for the inversion of the T2 spectrum. The proposed algorithm can achieve fast convergence and its inversion results are only slightly influenced by the noise when the SNR is higher than 20.

4.3 Comparison of the Computation Efficiency

Figs. 3 and 5 show the comparisons of the forward model, the inverted spectrums by the new method, and the inverted spectrums by the BRD algorithm. It is obvious that at high SNR (approximately larger than 50), both methods achieve favorable inversion results. However, the inversion results ofthe new method show better performance than the BRD algorithm, when the echo trains contain noises (particularly for the SNR lower than 20). The MSE is defined by:

MSE=1N∑i=1N(yi−inverted−yi−forward)2×10000(22)

where yi−inverted is the amplitude for i-th component of the inverted spectrum; yi−forward is the amplitude for the i-th component of the forward spectrum; N is the number of the components of the spectrum.

Fig. 6 gives the mean square error (MSE) between the forward model and the inversion results for different cases. It is observed that the new method achieves lower MSE than the BRD algorithm, especially for the unimodal distributed spectrum. Fig. 7 compares the time duration between the new method and the BRD inversion for our numerical cases. It is seen that the computation efficiency is greatly enhanced using the new method. Both the results of MSE and the time duration reveal the superiority of the new method.

Figure 6: The relationship between the MSE and the SNR at different numerical cases

Figure 7: The relationship between the time duration and the SNR at different numerical cases

5  Conclusions

A novel, efficient and accurate algorithm is developed for the inversion of the NMR T2 spectrum. In this method, the empirical TSVD and linearized Bregman iteration are used to enhance the speed and accuracy of the numerical inversion of the NMR T2 spectrum for the first time. The results of the numerical inversion study show that the linearized Bregman iteration can obtain quick and effective performance in the solution of ill-posed and over-determined problems. Moreover, this novel method can work well for data with the SNR higher than 20.

However, much research should be conducted to generalize this method to field applications and experiments. The effects of the iteration parameters and the acquisition parameters on inversion instability are also very important and can be investigated in the future.

Funding Statement: The authors are grateful for the financial support by the National Nature Science Foundation of China (42174142), CNPC Innovation Found (2021DQ02-0402), and National Key Foundation for Exploring Scientific Instrument of China (2013YQ170463).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.