Academic Journal of Science and Technology ISSN: 2771-3032 | Vol. 3, No. 2, 2022 125 Study on Inversion Algorithm of NMR T2 Spectrum Yong Deng*, Lin Xu School of Mechanical and Electrical Engineering, Southwest Petroleum University, Chengdu 610500, China * Corresponding author Abstract: Nuclear magnetic resonance (NMR) logging technology is a kind of technology developed gradually in recent decades, which has great significance for the development of logging technology. Nuclear magnetic resonance logging technology compared to traditional logging technology, not only has the ideal signal - hydrogen nuclei, and the exploration of the original signal is not affected by the influence of the rock skeleton, is the only provides information that is unrelated to formation rock porosity, can better provide the formation permeability, movable fluid, oil and gas content, such as information, It is a kind of logging method suitable for low porosity and low permeability and complex special lithology, and overcomes the defects of traditional logging methods .With the development of NMR theory, NMR has been applied to the analysis of well logging data. The original NMR logging data are the spin-echo string observed by the spin echo string, and its amplitude decays with time. The T2 spectrum can be inversely calculated from the measured logging data, and then the T2 spectrum can be further solved to obtain key data such as permeability and porosity. In this paper, a simulated T2 spectrum is constructed by forward modeling, and the inversion is carried out under different SNR. On this basis, the specific relationship between the optimal number of singular values retained and the change of SNR is summarized, and an improved singular value decomposition algorithm is proposed. Inverse problem can be converted to nonlinear fitting optimization problems with nonnegative constraints, linear truncation algorithm on calculation precision has increased significantly, compared with the traditional singular value decomposition inversion algorithm, a new algorithm of signal-to-noise ratio on the solution of smaller and has better stability, and through the improvement of nonnegative constraint processing structure was more rapid inversion iteration algorithm. Keywords: Nuclear magnetic resonance, T2 spectrum, The inversion, Singular value. 1. Introduction NMR was discovered in 1946 by two professors in the United States. The inversion of NMR T2 spectrum is a muruhe lti-exponential fitting of echo string. NMR has been active in the field of logging technology for more than 50 years, and has accumulated rich experience in basic theories such as logging principles. Noise is inevitably generated by logging instruments. In general, the signal-to-noise ratio (SNR) is calculated from the echo signal collected in the pore. The signal is extracted from the noise interference and the expected T2 spectrum is inverted. For the distribution of T2 transverse relaxation time, since the entire T2 spectrum distribution is equal to the initial echo size, too many distribution points will increase the uncertainty of spectral solution results, thus increasing the error of inversion results. If the number of distribution points is too small, most useful data will be lost, and the authenticity of the inversion image cannot be guaranteed. In the current research of NMR logging inversion, singular value decomposition (SVD) is the most widely used method. The new algorithm has strong adaptability to SNR and can get better inversion results under different distribution methods and different distribution points. The inversion method is based on the research of the traditional singular value decomposition inversion algorithm, analyzes the main factors affecting the inversion results, on this basis, the algorithm is optimized and improved, and the feasibility of the algorithm is proved by simulation experiment. In addition to the application value, the inversion algorithm research and experimental research can further improve the understanding of the nature of the inversion algorithm, and lay the foundation for the research of a more robust and practical new algorithm. 2. Theoretical Basis 2.1. Nuclear magnetic resonance theory As for the nucleus, it consists mainly of protons and neutrons. And there are big differences and differences between the two, one of the main aspects is whether charged, and protons and neutrons are called nucleon. The electric charge and the mass of the nucleus fully show its basic characteristics. The number of protons in the nucleus determines the charge of the nucleus, and the sum of the number of protons and neutrons is the mass of the nucleus. When the charge and mass of the nucleus are taken as the starting point in the analysis process, it is obvious that there is a strong interaction between the nucleus and the surrounding particles, but it is difficult to explain some weak interactions, such as the now well-known nuclear magnetic resonance phenomenon. Understanding the spin in the nucleus, which is produced by the momentum of another internal Angle, will explain the phenomena of NMR. At present, when we divide nuclei, we can divide nuclei with spin and nuclei without spin on the basis of spin status. Existing research has found that nuclei with odd atomic numbers have "spin". When we set up the corresponding magnetic field around the core in the current state, its motion state is similar to that of a gyro. In fact, one of the things that we're going to focus on in our analysis of NMR is the spin nucleus. 2.2. The relaxation process When no external magnetic field is applied, the nucleus population is in equilibrium, and its magnetization vector is in the same direction as the static magnetic field. When an external RF pulse is applied, the equilibrium of the nuclear 126 system is destroyed, the nucleus precession occurs, the magnetization vector direction of the RF pulse is shifted, and the nuclear system is in an unbalanced state. When the RF pulse is removed, the nuclear system is only affected by the static magnetic field, and the magnetization vector will gradually recover in the direction of the static magnetic field. This process of nuclear system gradually recovering from non-equilibrium to equilibrium state is called relaxation. Figure 1 shows the magnetization in equilibrium state, and Figure 2 shows the magnetization in imbalance state. The relaxation can be divided into transverse relaxation and longitudinal relaxation. When an external RF pulse is applied, the direction of the macroscopic magnetization vector M will shift, which can be decomposed into transverse component along the horizontal plane and longitudinal component along the longitudinal axis. Figure 1. Magnetization at equilibrium The recovery process along the horizontal plane of the transverse relaxation component M to the initial state value of zero is called the transverse relaxation process, and the time required for complete recovery is called the transverse relaxation time, denoted by T2. The recovery process of the longitudinal relaxation component M along the longitudinal axis to the initial state value M is called the longitudinal relaxation process, and the time required for complete recovery is called the longitudinal relaxation time, denoted by T1. Figure 2. Magnetization in an unbalanced state 3. NMR Relaxation Signal Model From the perspective of the composition of reservoir pores, they contain many fluid components, and for the spin echo series obtained with the support of CPMG pulse train, it mainly refers to the multi-exponential decay relaxation signals formed by the joint action of many single exponential decay signals. It can be obviously found that the component of the ith echo signal amplitude b mainly includes the response signals of various relaxation components. Based on this, the corresponding discrete model can be obtained as follows: 𝑏 βˆ‘ π‘₯ 𝑒 / πœ€ , 𝑖 1,2, … 𝑛 (1) In Equation 1, B 𝑏 , 𝑏 , … 𝑏 refers to the echo signal collected by the measurement, n is the number of echoes collected in the actual logging. X π‘₯ , π‘₯ , … π‘₯ is the contribution of different fluid components to the echo signal, that is, the relaxation spectrum amplitude of different fluid components to be calculated by inversion, and m is the number of relaxation components. 𝑑 𝑖. 𝜏 , Ο„ is the echo interval; T_2j is the relaxation time for the JTH relaxation component; πœ€ is the noise signal amplitude of the ith echo string. Thus can analyze the vector form of the discrete model is A B = AX + epsilon, in which A 𝑒 , 𝑖 1,2, … , 𝑛; 𝑗 1,2, … , π‘š. Equation 2 is the Fredholm integral equation of the first kind when the echo signal and the fluid component are regarded as continuous rather than finite: 𝑏 𝐾 𝑑, 𝑇 π‘₯ 𝑇 𝑑𝑇 (2) Where 𝐾 𝑑, 𝑇 𝑒 is called the integral kernel. T is the transverse relaxation time. π‘₯ 𝑇 refers to the curve formed by the transverse relaxation times of all fluid components, that is, the T2 spectrum. You solve integral equation 2 to get π‘₯ 𝑇 , and that's actually the T2 spectrum inversion problem. Discrete solution of continuous integral equation is usually the most common method to solve the above inversion problems. The discrete acquisition of n echo signals is the discretization of the echo signals. Specifically, the discretization of 𝑋 𝑇 requires the determination of m discrete points on the T2 time axis, which are the scatter points. Since most of the echo signals are in exponential decay state, the distribution of the distribution points is generally in the form of logarithmic equal spacing. In order to transform the complex integral equation into a simple system of linear equations 𝐡 𝐴𝑋 πœ€, we only need to discretize the complex integral equation. Based on the discretization of the first kind of Fredholm integral equation, it is found that most of the linear equations generated are ill-conditioned. Specifically, it means that the solution of the equation is more sensitive to the input, and a slight adjustment of the input may result in a big difference. In the vector form 𝐡 𝐴𝑋 πœ€ of the discrete model of the first kind of Fredholm integral equation, there is A noisy signal Ξ΅, so that the solutionX 𝐴 𝐡 of the system of equations will fluctuate to some extent. In the practical application process, because of the low SNR obtained, it is not obvious to solve 𝐡 𝐴𝑋 directly. The most widely used methods are singular value decomposition method and non-negative least square method. However, most of the above methods have some limitations, including high SNR requirements, slow computing speed, and spectral line discontinuity. 4. Inversion of T2 Spectrum Based on Singular Value Decomposition 4.1. Basic principle of T2 spectrum inversion The inversion process of T2 spectrum actually refers to the 127 process of solving the system 𝐡 𝐴𝑋 πœ€ . Since π‘›π‘˜ 𝐴, 𝐡 π‘Ÿπ‘Žπ‘›π‘˜ 𝐴 , the system 𝐡 𝐴𝑋 πœ€ is incompatible. For the incompatible system of equations, we want to obtain the approximate solution of the system of equations -least squares solution. Generally speaking, the solution of the system of equations is not unique. In engineering, we want to obtain the best least squares solution, that is, the least modular solution of the general solution. 4.2. Principle of singular value decomposition Gaussian iteration and LU decomposition methods are often used to deal with ill-conditioned or nearly ill- conditioned values, but in many cases, these methods have many shortcomings and are difficult to give satisfactory solutions. The singular value decomposition (SVD) method can effectively make up for the shortcomings of LU decomposition method, and is a more scientific and effective solution. Singular value decomposition theory is based on the following linear algebra theory, A M N represents M row N column matrix, equivalent to matrix UW𝑉 multiplication form, U is M*N column orthogonal matrix, W is diagonal element positive or zero diagonal matrix, 𝑉 is N *N column orthogonal matrix V transpose. Calculate the generalized inverse matrix from the M*N matrix, U and V are column orthogonal, the inverse matrix is equal to their transpose; W is a diagonal matrix, and the inverse matrix is the inverse of diagonal elements, which can be written as Equation 3: 𝐴 𝑉. π‘‘π‘–π‘Žπ‘” 1/𝑀 . π‘ˆ (3) If one or more Ο‰ zero or numerical value is too small, formula 3 mistakes will happen, in the computer rounding errors or numerical recognition difficulty formula 3 will be an error occurs, when the formula of zero increase matrix is singular, the matrix is called a singular matrix, with the singular value of the maximum and the minimum ratio increases, matrix is singular, If the value is too large for the computer to compute, it is called ill-conditioned matrix. The value domain of the null space of singular matrix has great influence on the matrix. The biggest advantage of singular value decomposition method is that it can accurately determine the value domain and null space of matrix. For example, the matrix A is A linear mapping equation Ax b from the vector space x to the vector space b, the null space dimension is called nullspace, is A subspace of x mapped to zero Ax 0 ), b is the range of A, calculated by linear combinations of the columns of A subspace, and the dimension of b is the rank of A. The range of A is the combination of columns in U corresponding to Ο‰ 0 , and the null space of A is the column in V corresponding to Ο‰ 0. When b is equal to 0, linear combinations of the null space are solutions. If bβ‰ 0, the equation has A solution if b and a have the same rank. And since there's an infinite number of linear combinations of the null space column vectors and solutions, there's an infinite number of solutions to this equation. We usually obtain the least squares solution formula by singular value decomposition method, and obtain certain values, as shown in Equation 4: 𝑓 𝑉. π‘‘π‘–π‘Žπ‘” . π‘ˆ . 𝑏 (4) If Ο‰ 0, then 0. If you compare the rank of matrix b to the rank of matrix A, if the rank of b is greater than the rank of A, then this equation has no solution. However, at this time, the singular value decomposition will still get a solution to the equation, which is not the exact solution of Ax b, but it is the closest solution in the sense of least squares. This approach doesn't take into account some of the linear components involved in the problem, but from a practical point of view, the overlooked variable is the exact solution including rounding errors, which is close to or even close to the null space vector. In the inversion of T2 spectrum, the number of rows of the system of equations is much larger than the number of columns, W will not have singular condition, and Ο‰ will not be zero. However, matrix A may have matrix degradation, so it is necessary to zero some small Ο‰ 5. Improved Singular Value Decomposition Algorithm 5.1. Optimal number of retained singular values For singular value decomposition (SVD), it is the key to select the appropriate number of retained singular values. When the SNR is too low, the solution is closely related to the number of retained singular values. In this paper, the minimum residual sum of squares is used as the measurement standard, and the optimal number of retained singular values is determined by calculating the optimal number of retained singular values under different SNR. The optimal number of retained singular values l satisfies: min 𝜎 βˆ‘ π‘₯ π‘₯ (5) π‘₯ : Construct the ith relaxation component of the spectrum; π‘₯ : The ith relaxation component of T2 spectrum obtained by inversion with l singular values retained; 𝜎 : Residual of the inversion results when the singular value of l is preserved. Taking the constructed echo string as the experimental data, the optimal number of preserved singular values is obtained by calculation, and compared with the number of singular values obtained by the traditional SVD algorithm using SNR as the condition number. The results are shown in Table 1. From the table, we can see that the number of singular values retained by the traditional SVD algorithm is generally less than the best number of retained values calculated. 128 Table 1. Number of retained singular values under different SNR SNR 5 10 20 30 50 70 90 100 SVD 保留ε₯‡ εΌ‚ε€ΌδΈͺζ•° 2 2 3 4 5 5 5 5 ζœ€δ½³δΏη•™ε₯‡ εΌ‚ε€ΌδΈͺζ•° 4 4 4 5 5 6 6 6 It can be seen from Figure 3 and Figure 4 that it is easy to invert the multi-peak T2 spectrum into a single-peak curve by using the traditional SVD algorithm, which is more serious when the SNR is low. Figure 3. Comparison results when SNR = 5. Figure 4. Comparison results when SNR = 30 Therefore, when measuring the number of retained singular values, only referring to the maximum singular value ratio (SVR) signal-to-noise ratio (SNR) cannot meet the information demand, and how to choose the appropriate condition number is the key to solve the problem. Because the condition number itself is continuous values, but to choose the singular value decomposition method in singular value reservation number when it is discrete, general best reserved to the singular values the number of l and the best condition number is a step function, both correspondence and are continuous, so the singular value decomposition method is the best condition number of a range. Therefore, when confirming the value interval of the best condition number, it can be obtained by the number of the best singular values retained,Ο‰ /Ο‰ ,Ο‰ /Ο‰ . Based on this, the singular value truncation algorithm is optimized: X Vdiag , , … , . , 0, … ,0 π‘ˆ π‘Œ (5) It can be called linear truncation algorithm. The traditional SVD algorithm using SNR as condition number is a special case of linear truncation algorithm a=1,b=0. The coefficients a and b used in the actual logging are determined by the above method, based on the different number of spots and echoes. Assumption conditions: Number of points m=32, number of echoes n=500. The experimental results show that when a=1.45 and b=16, the SVD algorithm can achieve the best effect under any SNR condition. The number of scatter points and echo numbers can be obtained through experiments by taking the values of linear truncation coefficients corresponding to other values, as shown in Table 2. There is an inverse relationship between the difference of the results of different truncation methods and the SNR. Compared with the traditional SVD algorithm, the linear truncation algorithm is more suitable for low SNR. When the SNR is higher than 100, the inversion results of both the traditional SVD and the linear truncation algorithm are basically consistent with the real spectrum, while the linear truncation algorithm is not only suitable for high SNR, When the SNR is not less than 10, the T2 spectrum can still be accurately reflected. 129 εΈƒη‚ΉδΈͺζ•° m ε›žζ³’δΈͺζ•° n η³»ζ•° a η³»ζ•° b 13 125 1 31 13 500 1 48 30 125 1.1 65 30 500 1.1 74 50 125 1.25 78 50 500 1.25 80 Combined with Table 2, it can be seen that the linear truncation coefficient a is not much affected by changes in the number of distribution points. When the number of spots increases by 4 times, the coefficient a only shows an increase of 1/4. The linear truncation coefficient b changes greatly with the increase of linear truncation coefficient. Meanwhile, the linear truncation coefficient a is not significantly affected by the variation of the number of echoes, while the coefficient b is on the contrary, which is greatly affected. It can be seen that the coefficient matrix A will change greatly when the number of scatter points and echo numbers rise, resulting in the continuous rise of the optimal condition number. Based on the experimental data, the relationship between the coefficient and the number of scatter points and the number of echoes was clarified. a and b represented the coefficient, m represented the number of scatter points, and n represented the number of echo: a 1 0.005m b 35 ln π‘š 58 0.16 𝑛 125 /βˆšπ‘š 5.2. Handling of nonnegative constraints In the process of non-negative constraint inversion of singular value, the minimum solution X needs to be found. If the solution is less than zero, the corresponding column of the solution needs to be deleted from matrix A and solved again until the minimum X is no longer less than zero, or the upper limit of forced decomposition times 𝐾 is met. The specific process is shown in Figure 5. Figure 5. Flow chart of constraints of conventional SVD algorithm This nonnegative constraint method is easy to control and simple to implement. However, in the practical application process, it is also found that this situation sometimes occurs. When n(nβ‰₯1) negative solutions are obtained after the KTH solution, the minimum value in matrix A is found, and then its corresponding column in matrix A is forcibly removed. In the K 1 solution process, the negative solutions in X will exceed n. In this case, it can be called "nonnegative constraint excess". If the columns in matrix A are deleted improperly, the continuity band of T2 spectrum distribution will be damaged greatly, and even the T2 spectrum will be deformed and distorted. It is also difficult to control when determining the upper limit of the number of forced decomposition 𝐾 . Once the number of decomposition is insufficient, there will be too much gap between the inversion spectrum and the real spectrum. On the contrary, too many decomposition times will affect the inversion speed greatly. In order to solve the above problems, Lin Feng et al. proposed an improved nonnegative constraint method through research. Only for matrix A on A singular value decomposition, the number of iterations assumed 𝐾 , for solving X, at the same time to get the minimum X min X , once X not greater than zero, is forced to make X_k zero and gain new X, called X ', After X' is obtained, Y AX is continued to be solved until all solutions are greater than zero or reach the upper limit of iteration 𝐾 . The flow diagram of specific methods is shown in FIG. 6. Figure 6. The improved inversion algorithm is a non- negative constraint scheme The improved algorithm does not need to do singular value decomposition repeatedly, so the inversion speed is accelerated. However, the term less than zero in X is forced to bereset to zero and then Y is recalculated, which will increase the amplitude of echo signal.For many iterations, the inversion result may be far from the original echo signal. In this way, the calculated echo will be above the spectral graph of the original echo signal, resulting in a big deviation of the result. In view of the above phenomena, this paper starts from the non-negative constraint mechanism to optimize. After the K 1 solution, the following process is performed: in the first step, the number of negative solutions in X is taken as the object to determine whether the number of negative solutions in X is less than the number of negative 130 solutions in X for the KTH solution. If it is less than the number of negative solutions in X for the KTH solution, the cycle of solving is stopped, and the negative solutions in X for the KTH solution are set to zero and the result is output. Second, take the minimum negative solution as the object, judge whether the minimum negative solution in X is greater than the minimum negative solution in X, if so, stop the cycle of solving, set the negative solution in the KTH solution to zero and output the result. The optimized non- negative constraint scheme can greatly improve the stability of inversion and greatly reduce the occurrence of malformations in the spectrogram. At the same time, because the possibility of excessive nonnegative constraints is reduced, the number of singular value decomposition also decreases, so the inversion speed is improved. In the inversion process of SVD, SVD of matrix A takes A lot of time. The complexity of decompressing matrix A into matrices U, D and V can be expressed as 4m n 8mn 9n flops. Therefore, in addition to the size of matrix A, the inversion speed is also affected by the number of singular value decomposition. If the number of echoes and the distribution scheme are not changed, it is necessary to reduce the number of singular value decomposition (SVD) or to speed up the efficiency of SVD in order to improve the inversion speed. In the process of processing the actual log data, the results of the first singular value decomposition (SVD) were consistent for all depth data after the placement scheme was determined. Therefore, the first result of all points can be represented by U,D and V obtained by the first decomposition. In this way, the number of singular value decomposition will decrease and the inversion speed will be improved. If H is used to represent the well depth and L is used to represent the sampling interval, the number of decompositions can be reduced to H/L. Through numerical simulation experiments, it is found that for most data points, the number of singular value decomposition is between 0 and 5, and the number of singular value decomposition can be reduced by at least 20%. At the same time, the upper limit of the number of mandatory iterations is set as K M method can also achieve the purpose of reducing the number of singular value decomposition, when the number of matrix factorization is higher than M times, the obtained solution will be forced to zero. Fast inversion is one of the best methods for most operations, even though it may eventually result in a reduction in the resolution of the results. Figure 7. Flow chart of improved algorithm 6. Conclusion By introducing the basic principle of T2 spectrum, the singular value decomposition method which is widely used is chosen to retrieve NMR logging signal. In essence, singular value decomposition (SVD) is to determine an appropriate linear cut-off value to help determine the condition number of linear equations, and then change the stability of the solution. How to choose the best condition number is the key to solve the inversion problem and obtain good inversion results.When the SNR is greater than 5, the linear truncation algorithm can be used, and when the SNR is very low, the authenticity of the relaxation spectrum distribution can be well maintained. When the SNR is greater than 30, the linear truncation algorithm should be used for better results. When the SNR is very low, the improved algorithm in this paper can still keep the authenticity of the relaxation spectrum distribution and accurately reflect the true information of the stratum. When the SNR is lower than 5, the fixed number of retained singular values is selected, and the inversion effect is good. By improving the treatment of nonnegative constraints, the stability of the adjusted inversion is significantly improved, which effectively improves the problem of malformation in the spectrum and controls it in the minimum range. In addition, the number of singular value decomposition decreases, so the inversion speed will be effectively improved. References [1] Zou Y L. Inversion method of NMR logging data and uncertainty of T_2 spectrum [D]. China University of Petroleum (Beijing), 2016. [2] Zhang Yan. Application of wavelet denoising in NMR logging [D]. Huazhong University of Science and Technology, 2015. [3] Guo Zhong-hua, LI Shu-Qing, WANG Lei, TANG Yan-wei. Improved Wavelet denoising algorithm based on Adaptive threshold [J]. Journal of Chongqing University of Posts and Telecommunications (Natural Science Edition), 2015, 27(06):740-744+750. 131 [4] Xiao Lizhi, ZHANG Hengrong, Liao Guangzhi, Fu Shaoqing, Li Kui. NMR relaxation inversion of porous media based on Backus-Gilbert theory [J]. Chinese Journal of Geophysics, 2012, 55(11):3821-3828. [5] Gao Yang, XIAO Lizhi. Inversion of NMR relaxation time using improved truncated singular value decomposition [J]. Oil Geophysical Prospecting, 2015, 50(02):376-381+8. [6] Chen Yuan. Inversion algorithm of T2 spectrum distribution in NMR logging [D]. Huazhong University of Science and Technology, 2009. [7] Li Jing. Research on T2 spectrum inversion algorithm based on singular value decomposition [D]. Northeast Petroleum University, 2018. [8] Zhu B. Overview of development, application and physical basis of nuclear magnetic resonance (NMR) [J]. Science and Technology Innovation and Application, 2013(5):11-11. [9] Anand C S,Sahambi J S. Wavelet domain non-linear filtering for MRI denoising.Magnetic Resonance Imaging, 2010,28(6):842-861. [10] Zwickl, Derrick Joel, PhD. Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criterion[J]. ProQuest Dissertations&Theses,2006. [11] Ma S, Kong L, Chen J. An improved NMR signal de-noising algorithm based on wavelet transform[J].Journal of Computational Information Systems, 2011, 7(13):4651- 4659. [12] Michael Prange, Yi-Qiao Song. Quantifying uncertainity in NMR T2 spectra using Monte Carlo inversion. Journal of Magnetic Resonance, 2009, 196:54-60. [13] Supriyo Ghosh, Kevin M. Keener, Yong Pan. A simulation based method to assess inversion algorithm for transverse relaxation data. Journal of Magnetic Resonance, 2008, 191:226-230.