Shale in-situ macropore fracture discrimination method based on nuclear magnetic-mercury injection combination

Abstract

The present application relates to the technical field of shale pore and fracture analysis, and particularly relates to a shale in-situ macropore and fracture discrimination method based on nuclear magnetic resonance-mercury injection. The shale sample in an in-situ state is first obtained; then, two-dimensional nuclear magnetic resonance testing is performed on the shale sample, and two-dimensional relaxation spectrum is inverted based on a truncated singular value decomposition method; in the inversion process, the value range of the singular value cutoff value is determined based on the two-dimensional data matrix under different signal-to-noise ratios, and the optimal cutoff value is determined through global optimization; then, high-pressure mercury injection testing is performed on the shale sample, and finally, the macropore and fracture region in the shale sample is determined by comprehensively considering the high-pressure mercury injection testing result and the two-dimensional relaxation spectrum. The present application determines the value range of the singular value cutoff value based on the nuclear magnetic resonance testing result under different signal-to-noise ratios, and determines the optimal cutoff value through global optimization, so as to improve the inversion effect, obtain more reliable two-dimensional relaxation spectrum, and then comprehensively discriminate the macropore and fracture by combining the testing results of nuclear magnetic resonance-mercury injection, thereby improving the discrimination accuracy.

Description

Technical Field

[0001] This invention relates to the field of shale pore and fracture analysis technology, specifically to a method for identifying macropores and fractures in shale based on nuclear magnetic resonance-mercury intrusion porosimetry. Background Technology

[0002] As an important carrier of unconventional oil and gas resources, the accurate characterization of the pore structure of shale reservoirs is crucial for reservoir evaluation and oil and gas development. Common pore structure analysis techniques include mercury intrusion porosimetry (MIP) and nuclear magnetic resonance (NMR). MIP establishes the relationship between injection pressure and pore size using the Washburn equation, enabling the measurement of pore distribution from nanometer to hundreds of micrometer scales. NMR assesses pore size by analyzing the lateral relaxation time distribution, offering advantages such as non-destructive testing and in-situ fluid sensitivity.

[0003] However, mercury intrusion porosimetry carries the risk of structural damage. To detect large pores, high pressure is usually applied during the test, which can easily cause shale samples to fracture, generating experimentally induced cracks and interfering with the true pore signal (such as producing false large pores). In the process of nuclear magnetic resonance (NMR) detection, the conversion of transverse relaxation time into pore size depends on the accurate calibration of surface relaxation intensity. However, the complex composition of shale minerals results in significant heterogeneity in surface relaxation intensity, which can easily lead to deviations in pore analysis. At the same time, pore volume is proportional to NMR signal intensity, and large pore signals can easily be masked by small pore signals, resulting in missed detections. Mercury intrusion porosimetry and NMR technology have low accuracy in detecting large pores in shale in situ. Summary of the Invention

[0004] To address the technical problem of low accuracy in detecting in-situ macropore fractures in shale, the present invention aims to provide a method for identifying in-situ macropore fractures in shale based on nuclear magnetic resonance-mercury intrusion porosimetry (NMR-MILI). The specific technical solution adopted is as follows:

[0005] A method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry, the method comprising:

[0006] Obtain shale samples in situ;

[0007] Two-dimensional nuclear magnetic resonance (NMR) tests were performed on the shale sample based on preset NMR test parameters. A two-dimensional data matrix was constructed using the NMR test results, and a two-dimensional relaxation spectrum was inverted based on the truncated singular value decomposition method. In the inversion process, the range of singular value cutoff values ​​was determined based on different signal-to-noise ratios, and the optimal cutoff value within the range was obtained based on a global optimization algorithm. The inversion result under the optimal cutoff value was used as the two-dimensional relaxation spectrum.

[0008] High-pressure mercury intrusion testing was performed on the shale sample based on preset mercury intrusion testing parameters. The macropore and fracture regions in the shale sample were determined by combining the high-pressure mercury intrusion testing results with the two-dimensional relaxation spectrum.

[0009] Furthermore, the preset two-dimensional NMR scan parameters include at least the inversion recovery sequence parameters, the spin echo sequence parameters, and the number of scans; wherein, the inversion recovery sequence parameters include the inversion time and the waiting time, and the spin echo sequence parameters include the echo interval and the number of echoes.

[0010] Furthermore, the method for obtaining the two-dimensional data matrix includes:

[0011] In the two-dimensional nuclear magnetic resonance test, the shale sample is scanned a preset number of times to obtain the CPMG echo string at each reversal time; the CPMG echo string at each reversal time is used as a row vector to construct a two-dimensional data matrix.

[0012] Furthermore, the method for obtaining the two-dimensional data matrix under different signal-to-noise ratios includes:

[0013] By changing the preset number of scans, the two-dimensional data matrix constructed under each preset number of scans is used as a two-dimensional data matrix with a signal-to-noise ratio.

[0014] Furthermore, the method for obtaining the range of singular value cutoff values ​​includes:

[0015] The two-dimensional data matrix at each signal-to-noise ratio is used as the test matrix; for each test matrix, a linear model is constructed based on the attenuation characteristics of the CPMG echo train with waiting time and acquisition time, and the coefficient matrix in the linear model is determined; wherein, the dependent variable of the linear model is the test matrix, and the independent variable of the linear model is the two-dimensional relaxation spectrum corresponding matrix to be solved.

[0016] Based on a preset algorithm, the singular value cutoff value of the coefficient matrix corresponding to each test matrix is ​​determined, and the initial inversion of the corresponding linear model is performed based on the singular value cutoff value to obtain the proportion of negative components in the initial inversion result; based on the proportion of negative components corresponding to all test matrices, the range of values ​​for the singular value cutoff value is determined.

[0017] Furthermore, based on the proportion of the negative components corresponding to all the test matrices, the range of singular value cutoff values ​​is determined, including:

[0018] All test matrices corresponding to the proportion of negative components are sorted and fitted in ascending order of signal-to-noise ratio to obtain the curve of change of negative component proportion; the extreme value among all abrupt change points in the curve of change of negative component proportion is determined, and the singular value cutoff value of the coefficient matrix of the test matrix corresponding to the extreme value is successively used as the upper and lower limits to obtain the range of singular value cutoff value.

[0019] Furthermore, the preset algorithm includes at least the L-curve method.

[0020] Furthermore, the method for obtaining the optimal cutoff value includes:

[0021] A linear model of the two-dimensional data matrix is ​​constructed, and the linear model is inverted based on the truncated singular value decomposition method. Based on the difference between the inversion results before and after truncating the singular values ​​within the range of values, and the negative component in the inversion results after truncation, the target optimization function is determined.

[0022] The optimal cutoff value within the range of values ​​is determined using the particle swarm optimization algorithm. The optimal cutoff value is the cutoff value corresponding to the minimum value of the objective optimization function.

[0023] Furthermore, the preset mercury intrusion porosimetry test parameters include at least the termination condition during the high-pressure mercury intrusion porosimetry test, the pressure change step size of the applied pressure, and the pressure holding time at each pressure; wherein, the termination condition is that the applied pressure reaches the preset safe pressure threshold.

[0024] Furthermore, by combining the results of high-pressure mercury intrusion porosimetry and the two-dimensional relaxation spectrum, the macroporous fracture regions in the shale sample were determined, including:

[0025] A pore size-pore volume reference map is fitted based on the high-pressure mercury intrusion porosimetry test results, and a pore size-pore volume reference map is fitted based on the two-dimensional relaxation spectrum. When both the mercury intrusion porosimetry reference map and the NMR reference map have pore sizes larger than a preset size, it is determined that there are large pores in the shale sample.

[0026] The present invention has the following beneficial effects:

[0027] This invention first obtains shale samples in situ to preserve the original reservoir state of the shale samples to the greatest extent possible and avoid interference with subsequent test results. Then, based on preset NMR test parameters, two-dimensional NMR tests are performed on the shale samples to assess porosity distribution by analyzing relaxation time. A two-dimensional data matrix is ​​then constructed using the two-dimensional NMR test results, and a two-dimensional relaxation spectrum is inverted using the truncated singular value decomposition method to suppress noise while preserving important signals. During the inversion process, the influence of the signal-to-noise ratio on the inversion results under the traditional truncation value is analyzed based on the two-dimensional data matrix under different signal-to-noise ratios to determine the range of singular value cutoff values. An optimal cutoff value within this range is obtained using a global optimization algorithm to avoid local convergence of cutoff values. The inversion result under the optimal cutoff value is then used as the two-dimensional relaxation spectrum. High-pressure mercury intrusion porosimetry (HIP) is performed on the shale samples based on preset mercury intrusion porosimetry parameters for combined NMR-HIP testing. Finally, the macropore and fracture regions in the shale samples are determined by combining the HIP results and the two-dimensional relaxation spectrum. This invention determines the range of singular value cutoff values ​​based on NMR test results under different signal-to-noise ratios, and globally optimizes to determine the optimal cutoff value, thereby improving the inversion effect and obtaining a more reliable two-dimensional relaxation spectrum. Then, it combines the NMR-mercury intrusion porosimetry test results to comprehensively identify large pores, thus improving the identification accuracy. Attached Figure Description

[0028] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0029] Figure 1 A flowchart illustrating a method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry, provided as an embodiment of the present invention;

[0030] Figure 2 This is a flowchart illustrating a method for obtaining the range of singular value cutoff values ​​according to an embodiment of the present invention. Detailed Implementation

[0031] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the following, in conjunction with the accompanying drawings and preferred embodiments, details the specific implementation, structure, features, and effects of a method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry (NMR-MILI) proposed according to the present invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0032] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0033] The following description, in conjunction with the accompanying drawings, details a specific scheme for an in-situ macropore identification method for shale based on nuclear magnetic resonance-mercury intrusion porosimetry provided by this invention.

[0034] Please see Figure 1 The document illustrates a flowchart of a method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry (NMR-MJI) according to an embodiment of the present invention, specifically including:

[0035] Step S1: Obtain shale samples in situ.

[0036] To identify large pores and fractures in shale, one embodiment of the present invention first drills a core from the shale reservoir to obtain a shale sample in its in-situ state. The shale sample is a core sample that can represent the original underground reservoir state and has an intact physical structure. During the drilling process, areas with natural fracture development should be avoided to prevent interference from natural fractures and reduce the risk of damage during core processing or mercury intrusion porosimetry.

[0037] It should be noted that the specific procedures for obtaining shale samples in situ are well-known to those skilled in the art, and will be briefly described here:

[0038] Shale samples were obtained by cutting the core using a diamond wire cutter. During the cutting process, the cutting speed was controlled at 0.5-1.0 mm / min to avoid thermal damage. Deionized water was used as the coolant to avoid ion contamination of the sample, which could affect the relaxation signal in subsequent nuclear magnetic resonance (NMR) tests.

[0039] The diameter of the shale sample was controlled to be 25.0±0.1mm, the height to be 30.0±0.2mm, and the end face parallelism error to be ≤0.05mm, in order to meet the testing requirements of subsequent nuclear magnetic resonance testing and mercury intrusion porosimetry testing. The implementer may also adjust it according to the actual situation.

[0040] Immediately after cutting, the shale samples are placed in a constant temperature chamber at the formation temperature to prevent fluid evaporation or oxidation, so as to preserve the original reservoir state of the shale samples to the greatest extent and avoid interfering with subsequent test results.

[0041] Step S2: Two-dimensional nuclear magnetic resonance (NMR) testing is performed on the shale sample based on preset NMR test parameters. A two-dimensional data matrix is ​​constructed using the two-dimensional NMR test results, and a two-dimensional relaxation spectrum is inverted based on the truncated singular value decomposition method. During the inversion process, the range of singular value cutoff values ​​is determined based on the two-dimensional data matrix under different signal-to-noise ratios, and the optimal cutoff value within the range is obtained based on a global optimization algorithm. The inversion result under the optimal cutoff value is used as the two-dimensional relaxation spectrum.

[0042] Given the highly complex pore structure in shale, which includes nanoscale and microscale pores, and considering that two-dimensional nuclear magnetic resonance (NMR) is a non-destructive testing method that can distinguish different types of pores at multiple scales, it helps to comprehensively understand the distribution of pores.

[0043] Furthermore, considering that the test results of two-dimensional nuclear magnetic resonance are usually a time-domain signal, which is the superposition result of the hydrogen nucleus decay signals of the fluid in the pores of the shale sample, it is impossible to directly resolve the pore distribution information. It is necessary to convert the time-domain signal into a two-dimensional relaxation spectrum through mathematical inversion before the distribution characteristics of different pores can be analyzed.

[0044] Furthermore, it is considered that the truncated singular value decomposition (TSVD) method automatically filters out the signal components that contribute the most to the inversion results by truncating singular values ​​and ignoring noise-dominated components, which can effectively suppress the influence of noise on the inversion results and adaptively preserve the key features of pore components.

[0045] Based on this, the embodiments of the present invention first perform two-dimensional nuclear magnetic resonance (NMR) tests on shale samples based on preset NMR test parameters, construct a two-dimensional data matrix using the two-dimensional NMR test results, and invert the two-dimensional relaxation spectrum based on the truncated singular value decomposition method.

[0046] It should be noted that before performing two-dimensional nuclear magnetic resonance (NMR) testing on shale samples, the surface paraffin layer must be removed to avoid contaminating the NMR probe. This is a well-known technical method, and the specific processing procedure will not be described in detail here.

[0047] Among them, considering that there may be random noise in the test signal of nuclear magnetic resonance, and nuclear magnetic resonance inversion is very sensitive to noise, test signals with low signal-to-noise ratio are prone to generating a large number of negative components during inversion, which leads to unreliable inversion results; by repeatedly scanning nuclear magnetic resonance, the signal-to-noise ratio of the test signal can be improved to a certain extent, thereby reducing noise interference in the inversion process;

[0048] Furthermore, considering that relaxation time reflects the speed at which fluid molecules recover from the excited state to the ground state; for fluid molecules in pores, relaxation time is affected by their degrees of freedom of motion in the pores and their interaction with the pore walls; longitudinal relaxation time T1 reflects the interaction strength between the pore surface and the fluid, while transverse relaxation time T2 is related to the interaction or diffusion behavior between fluid molecules; therefore, analyzing longitudinal and transverse relaxation times can help assess pore distribution; and IR-CPMG can detect longitudinal relaxation time T1 by inverting pulses, while obtaining relevant information about transverse relaxation time T2 through CPMG sequences;

[0049] Based on this, in a preferred embodiment of the present invention, the preset two-dimensional nuclear magnetic resonance scanning parameters include at least inversion recovery sequence parameters, spin echo sequence parameters, and number of scans; wherein, the inversion recovery sequence parameters include inversion time and waiting time, and the spin echo sequence parameters include echo interval and number of echoes.

[0050] Then, based on the preset NMR test parameters, two-dimensional NMR tests were performed on the shale samples, and a two-dimensional data matrix was constructed using the two-dimensional NMR test results.

[0051] Preferably, in one embodiment of the present invention, the method for obtaining a two-dimensional data matrix includes: in a two-dimensional nuclear magnetic resonance test, scanning the shale sample a preset number of times to obtain the CPMG echo string at each reversal time; and using the CPMG echo string at each reversal time as a row vector to construct a two-dimensional data matrix.

[0052] Specifically, a nuclear magnetic resonance (NMR) spectrometer such as the Bruker Avance NEO is used, with the magnetic field strength set between 1T and 3T. In this embodiment, it is set to 0.47T. The NMR spectrometer is equipped with a gradient system and a temperature control unit. During the test, the formation temperature still needs to be simulated to maintain the in-situ state, and the temperature control accuracy is set to ±0.5 degrees Celsius to ensure stable test conditions and make the relaxation time measurement repeatable and comparable. The pulse sequence type of the NMR spectrometer is set to IR-CPMG sequence.

[0053] The IR sequence is used to measure the longitudinal relaxation time T1: the value range of T1 is set to 1ms-10s to cover the T1 values ​​corresponding to fluids in various pores in the shale sample; a preset number of reversal times TI are set, which is 16 in this embodiment, and is specifically determined by logarithmic intervals (a well-known technique, which will not be elaborated further); a waiting time TW (also known as polarization waiting time) is set, which satisfies TW≥5×T1max to ensure complete recovery of longitudinal magnetization, where T1max is the preset maximum value of longitudinal relaxation time, and TW is set to 60s in this embodiment; the CPMG sequence is used to measure the transverse relaxation time T2: the echo interval TE between two adjacent echoes in the CPMG echo train is set to 0.2ms to capture micropore signals, and the number of echoes NECH is set to 5000 to ensure T2 spectral resolution;

[0054] Two-dimensional nuclear magnetic resonance (NMR) testing is a well-known technique. Here is a brief description of the general process of obtaining test results and constructing a two-dimensional data matrix: First, the shale sample is placed in a static magnetic field for time TW. The magnetization vector M will reach its maximum value M0 along the magnetic field direction (+Z axis). For each reversal time TI, a 180° reversal pulse is applied to flip the magnetization vector M from the +Z axis to the -Z axis. The reversal pulse is then turned off. After waiting for another reversal time TI, a 90° pulse is applied to flip the current longitudinal magnetization to the XY plane. Then, a series of 180° pulses are immediately applied to acquire the CPMG echo train.

[0055] The number of scans N needs to be set to ensure that the signal-to-noise ratio is increased to more than 8 times the baseline noise. In this embodiment, the number of scans is set to 64. For each reversal time TI, the scan is repeated 64 times (a series of 180° pulses are repeatedly applied to acquire CPMG echo strings), resulting in 64 CPMG echo strings. The 64 CPMG echo strings are averaged and fused to obtain the final CPMG echo string at that reversal time. The same operation is performed at each reversal time to acquire 16 CPMG echo strings. The CPMG echo strings at each reversal time are used as row vectors to construct a two-dimensional data matrix. The dimension of the two-dimensional data matrix is ​​16×5000.

[0056] After constructing the two-dimensional data matrix, a linear model can be constructed by combining the attenuation characteristics of the CPMG echo train with the waiting time and acquisition time. This allows for the decoupling of the longitudinal relaxation time T1 and the transverse relaxation time T2, and the inversion of the two-dimensional relaxation spectrum (T1-T2 spectrum), which prepares for subsequent analysis of the porosity distribution under the nuclear magnetic resonance testing angle.

[0057] According to the Fredholm equation of the first kind, the attenuation of the CPMG echo train with waiting time and acquisition time can be expressed in integral form (a well-known technique, briefly described here); specifically:

[0058] ;

[0059] In the formula, The index of the reversal time in the IR sequence can also be called the index of the repolarization time of the magnetization vector. Indicates the sequence number of the echo time in the CPMG echo train; In a two-dimensional data matrix, the time inversion is... Echo time is The signal amplitude at that time, i.e. The first CPMG echo train acquired below Echo amplitude at each echo time; Indicates the initial magnetization intensity (signal amplitude after full polarization); It is a natural constant; For longitudinal relaxation time; This refers to the lateral relaxation time; The content of magnetic nuclei with transverse relaxation time T2 and longitudinal relaxation time T1 is represented, and it also represents the relaxation time distribution function; the integration ranges of the two integral symbols are [T2min, T2max] and [T1min, T1max], respectively; d is the differential symbol;

[0060] In the above formula, The kernel function is of dimension T1, representing the dynamic process of polarization recovery; Let be a kernel function of dimension T2, representing the dynamic process of signal attenuation; discretizing the above integral formula, it can be transformed into matrix form: Where M represents a two-dimensional data matrix; This represents the matrix corresponding to the kernel function in the T2 dimension; The matrix representing the content of magnetic cores; This is the matrix corresponding to the kernel function in dimension T1;

[0061] Based on the Kronecker product, the above matrix is... Rewritten as , and then simplified to ;in, ; ; ; For matrix , The Kronecker expression, that is, K is the kernel matrix composed of T1 and T2 dimensions; This indicates that the matrices within the parentheses are concatenated into a column vector in column order; where, If the linear model that needs to be inverted is a given m and K, then the inversion problem is transformed into the problem of finding s.

[0062] Then, according to the Lanczos singular value decomposition theory, the kernel matrix K can be transformed and decomposed into... Assuming K is an r-order m×n matrix, then Let m×r be a matrix that satisfies ( (where r is the identity matrix). Let n×r be a matrix that satisfies ; Represent a diagonal matrix, by or It consists of the positive square roots of the r non-zero eigenvalues, that is, it consists of the r singular values ​​of the kernel matrix K.

[0063] It should be noted that the above inversion process analysis is a well-known technique familiar to those skilled in the art.

[0064] Considering that small singular values ​​can be truncated during inversion based on the TSVD algorithm to effectively suppress noise interference; however, traditional TSVD inversion relies on a fixed singular value truncation value. If the truncation value is too small (i.e., too many singular values ​​are retained), small singular values ​​corresponding to noise will be retained, causing noise to be amplified in the inversion result, resulting in spurious oscillations and negative values, and unstable inversion results. If the truncation value is too large (i.e., too few singular values ​​are retained), the singular values ​​corresponding to the true signal may be lost, resulting in overly smoothed inversion results and loss of detailed information, especially weak but important signals such as signals from large apertures.

[0065] Furthermore, considering that the quality of the inversion results will change under different signal-to-noise ratios (SNRs), by systematically changing the SNR (e.g., changing the number of scans), two-dimensional data matrices under different SNRs are obtained. Then, an initial inversion is performed under each SNR. Based on the singular value cutoff value set in the initial inversion process where the quality of the inversion results is relatively stable, the range of values ​​for a singular value cutoff value can be determined. This allows the truncation strategy to be adaptively adjusted in subsequent inversion processes, enabling the inversion algorithm to automatically select the optimal cutoff value and obtain a more reliable two-dimensional relaxation spectrum.

[0066] Therefore, in the process of inverting the two-dimensional relaxation spectrum based on the truncated singular value decomposition method, the embodiments of the present invention determine the range of singular value cutoff values ​​based on the two-dimensional data matrix under different signal-to-noise ratios, and further obtain the optimal cutoff value within the range based on the global optimization algorithm. The inversion result under the optimal cutoff value is used as the two-dimensional relaxation spectrum to prepare for subsequent analysis of pore distribution.

[0067] Preferably, in one embodiment of the present invention, the method for obtaining two-dimensional data matrices under different signal-to-noise ratios includes:

[0068] By changing the number of scans, the two-dimensional data matrix constructed under each number of scans is used as a two-dimensional data matrix with a signal-to-noise ratio. Specifically, in this embodiment, the number of scans is set to 8, 16, 32, 64, and 128 in sequence. Based on the above two-dimensional data matrix construction method, the two-dimensional data matrix constructed under each number of scans is obtained in sequence, which will not be repeated here. Each two-dimensional data matrix corresponds to a signal-to-noise ratio; wherein, the signal-to-noise ratio increases with the increase of the number of scans.

[0069] Further, the range of singular value cutoff values ​​is determined based on two-dimensional data matrices under different signal-to-noise ratios.

[0070] Preferably, in one embodiment of the present invention, the method for obtaining the range of singular value cutoff values ​​includes:

[0071] Please see Figure 2The flowchart illustrates a method for obtaining the range of singular value cutoff values ​​according to an embodiment of the present invention, specifically including:

[0072] Step S201: Use the two-dimensional data matrix under each signal-to-noise ratio as the test matrix; for each test matrix, construct a linear model based on the attenuation characteristics of the CPMG echo train with waiting time and acquisition time, and determine the coefficient matrix in the linear model; where the dependent variable of the linear model is the test matrix, and the independent variable of the linear model is the two-dimensional relaxation spectrum corresponding matrix to be solved.

[0073] First, the two-dimensional data matrix at each signal-to-noise ratio is used as the test matrix for subsequent initial inversion to analyze the quality of the inversion results and to prepare for determining the range of singular value cutoffs; then, for each test matrix... Based on the above analysis and transformation of the integral form, a linear model is constructed. The specific construction process will not be repeated here; among them, For the dependent variable in a linear model; This is the coefficient matrix in the linear model; , where is the independent variable of the linear model, is the corresponding matrix of the two-dimensional relaxation spectrum to be solved, and is the inversion result to be solved in the end.

[0074] Step S202: Determine the singular value cutoff value of the coefficient matrix corresponding to each test matrix based on a preset algorithm, and perform initial inversion on the corresponding linear model based on the singular value cutoff value to obtain the proportion of negative components in the initial inversion result; determine the range of values ​​for the singular value cutoff value based on the proportion of negative components corresponding to all test matrices.

[0075] After determining the test matrices under different signal-to-noise ratios, the initial inversion can be performed based on the TSVD algorithm. First, the singular value cutoff value of the coefficient matrix corresponding to each test matrix is ​​determined based on the preset algorithm. The preset algorithm includes at least the L-curve method (a well-known technical means, which will not be described in detail). In other embodiments, the implementer can also use the BG compromise criterion to obtain the singular value cutoff value. Then, the initial inversion is performed based on the obtained singular value cutoff value to obtain the initial solution S(0). The proportion of negative components in S(0) is counted. The proportion of negative components reflects the quality of the inversion result to a certain extent. The lower the proportion of negative components, the higher the quality of the inversion result.

[0076] It should be noted that obtaining the proportion of negative components is a well-known technical method, and the specific acquisition process will not be described in detail. Implementers can also perform multiple initial inversions, such as 10, for each linear model, and average the proportion of negative components in all initial inversion results to obtain the final proportion of negative components, so as to avoid random errors in a single inversion process.

[0077] Considering that a higher proportion of negative components in the initial inversion results means a larger singular value cutoff value, retaining too many small singular values, which are usually highly correlated with noise; conversely, a lower proportion of negative components means fewer small singular values ​​are retained, resulting in better noise suppression; therefore, the range of singular value cutoff values ​​is further determined based on the proportion of negative components corresponding to all test matrices.

[0078] Preferably, in one embodiment of the present invention, by systematically changing the signal-to-noise ratio (SNR) and performing an initial inversion, a curve relating the SNR to the proportion of negative components is fitted. The abrupt change points on this curve mark the critical points where the inversion quality significantly changes, such as a shift from noise dominance to signal dominance, or from excessive smoothing to the point where effective signals are beginning to be retained. Furthermore, the range of values ​​can be determined based on the singular value cutoff value during the inversion process where the proportion of negative components is relatively stable between abrupt change points. Therefore, the range of values ​​for the singular value cutoff value is determined based on the proportion of negative components corresponding to all test matrices, including:

[0079] All test matrices are sorted and fitted in ascending order of signal-to-noise ratio to obtain the curve of negative component proportion change. The extreme values ​​among all abrupt change points in the curve of negative component proportion change are determined. The singular value cutoff values ​​of the coefficient matrix of the test matrix corresponding to the extreme values ​​are used as the upper and lower limits to obtain the range of singular value cutoff values.

[0080] Specifically, all test matrices are sorted in ascending order of signal-to-noise ratio for the proportion of negative components. Curve fitting is performed based on the least squares method to obtain the curve of the change in the proportion of negative components. Then, abrupt change point detection is performed on the curve of the change in the proportion of negative components, such as t-test and MK test, to determine the maximum and minimum values ​​among all abrupt change points in the curve of the change in the proportion of negative components.

[0081] The mutation point indicates a sudden and significant change in the proportion of negative components at the corresponding signal-to-noise ratio, until another significant change occurs at the next mutation point. The minimum value at the mutation point corresponds to the singular value truncation value of the coefficient matrix of the corresponding test matrix during the initial inversion process, which serves as the lower limit of the range of singular value cutoff values. The maximum value at the mutation point corresponds to the singular value truncation value of the coefficient matrix of the corresponding test matrix during the initial inversion process, which serves as the upper limit of the range of singular value cutoff values.

[0082] It should be noted that the least squares method for curve fitting and mutation point detection is a well-known technique and will not be elaborated further.

[0083] Considering that the core role of TSVD in inversion is to reduce the sensitivity and interdependence of model parameters by truncation, it completely discards the subspace information corresponding to the truncated singular values, which often leads to local optima and failure to converge to the global minimum, resulting in local convergence. However, the population cooperation mechanism or random perturbation mechanism in the global optimization mode can solve the problem of local convergence caused by truncating singular values ​​during the inversion process. Therefore, this embodiment of the invention further obtains the optimal cutoff value within the range of values ​​based on the global optimization algorithm, and uses the inversion result under the optimal cutoff value as a two-dimensional relaxation spectrum.

[0084] Preferably, in one embodiment of the present invention, the method for obtaining the optimal cutoff value includes:

[0085] A linear model of a two-dimensional data matrix is ​​constructed. The linear model is inverted based on the truncated singular value decomposition method. The objective optimization function is determined based on the difference between the inversion results before and after truncating the singular values ​​within the range of values, as well as the negative components in the inversion results after truncation. The optimal cutoff value within the range of values ​​is determined using the particle swarm optimization algorithm. The optimal cutoff value is the cutoff value corresponding to the minimum value of the objective optimization function.

[0086] Specifically, to determine the optimal cutoff value for inverting the two-dimensional data matrix, a linear model of the two-dimensional data matrix can first be constructed based on the analytical transformation of the above integral form. The details will not be elaborated further; then, the linear model is inverted based on the truncated singular value decomposition method, and the optimal cutoff value within the range is determined by the particle swarm optimization algorithm during the inversion process.

[0087] The objective optimization function is set as follows: , For any singular value within the range of values, truncate the value and change It can obtain the inversion results after truncating several singular values ​​to perform target optimization; Singular value cutoff value The truncated inversion results; The inversion result without singular value truncation; the numerator in the objective optimization function. The difference between the inversion results before and after truncating the singular values ​​within the range of values ​​is used to measure the amount of information loss; the denominator As a normalization factor, it avoids the magnitude of the solution affecting the value of the objective optimization function;

[0088] It should be noted that particle swarm optimization is a well-known technique familiar to those skilled in the art; the general optimization process is briefly described here:

[0089] First, the particle swarm is initialized: the number of particles is set to 1000, and the positions of 70% of the total number of particles are set as the initial solution S(0), that is, the element values ​​of the initial solution are assigned to the parameter vector of the particles. Within the range of the singular value cutoff value determined above, the initial positions of the remaining 30% of particles are randomly generated; the position of each particle represents a cutoff value to avoid the search from stalling prematurely due to all particle positions being completely identical. The particle swarm optimization parameters are set as follows: the inertia weight w ranges from [0.1, 1.5], and in this embodiment, it is linearly decreased from 0.9 to 0.4 (large inertia at the beginning is beneficial for global exploration; small inertia in the later stage is beneficial for local search); the learning factors c1 and c2 both range from [1.1, 2], and in this embodiment, they are both set to 1.5; the number of iterations is set to 2000.

[0090] Then, within the range of singular value cutoff values, iterative optimization is performed using particle swarm optimization. The fitness of each particle is evaluated through the objective optimization function, and the cutoff value corresponding to the minimum value of the objective optimization function is found and taken as the optimal cutoff value.

[0091] Obtaining the inversion result at the optimal cutoff value yields a two-dimensional relaxation spectrum, a well-known technique that will not be elaborated further.

[0092] Thus, two-dimensional relaxation spectra of shale samples were obtained based on nuclear magnetic resonance (NMR) testing, preparing for further analysis of porosity by combining mercury intrusion porosimetry (MIP) results.

[0093] Step S3: Perform high-pressure mercury intrusion testing on the shale sample based on preset mercury intrusion testing parameters, and determine the macropore fracture region in the shale sample by combining the high-pressure mercury intrusion test results and two-dimensional relaxation spectrum.

[0094] It should be noted that the shale sample needs to be dried before the mercury intrusion porosimetry test to avoid splitting caused by subsequent pressure.

[0095] Specifically, the shale samples were dried in a 60°C vacuum oven for 48 hours until constant weight was achieved (mass change <0.1%). Then, a fully automatic mercury porosimeter (Micromeritics AutoPore V series) was used for mercury porosimeter testing. The fully automatic mercury porosimeter is also equipped with a high-precision pressure sensor (range 0.1-50 MPa) and a mercury intrusion detection unit.

[0096] In a preferred embodiment of the present invention, the preset mercury intrusion porosimetry test parameters include at least the termination condition during the high-pressure mercury intrusion porosimetry test, the pressure change step size of the applied pressure, and the pressure holding time at each pressure; wherein, the termination condition is that the applied pressure reaches a preset safe pressure threshold.

[0097] As an example, the applied pressure during the high-pressure mercury intrusion porosimetry test is controlled in segments, divided into low-pressure and high-pressure segments. The low-pressure segment has a pressure range of 0.1-10 MPa, a pressure change step of 0.5 MPa, and a pressure holding time of 30 seconds at each pressure level. That is, after each pressure change, the pressure is held for 30 seconds before the next change. The high-pressure segment has a pressure range of 10-50 MPa, a pressure change step of 5 MPa, and a pressure holding time of 60 seconds at each pressure level. The termination condition is that the applied pressure reaches the preset safe pressure threshold of 50 MPa. In other examples, the implementer can also customize the termination condition, such as stopping when the mercury injection rate increases by less than 0.1% for three consecutive levels, or adjust the pressure control scheme according to their own preferences.

[0098] During the mercury intrusion porosimetry test, the volume of mercury injected at each pressure level is recorded, the cumulative amount of mercury injected is calculated, and the pore size at each pressure level is calculated to generate a pore size-pore volume distribution curve. The above operations are well-known technical methods and will not be described in detail here.

[0099] Further analysis of high-pressure mercury intrusion testing results and two-dimensional relaxation spectra determined the macropore fracture regions in the shale samples.

[0100] Preferably, in one embodiment of the present invention, determining the macroporous fracture region in a shale sample by combining high-pressure mercury intrusion porosimetry results and two-dimensional relaxation spectroscopy includes:

[0101] A pore size-pore volume reference map is fitted based on the high-pressure mercury intrusion porosimetry test results, and a pore size-pore volume reference map is fitted based on the two-dimensional relaxation spectrum. When the pore size in both the mercury intrusion porosimetry reference map and the NMR reference map is larger than the preset size, it is determined that there are large pores in the shale sample.

[0102] Considering that the pore structure in some complex rock materials may be highly irregular, specific models such as power functions may not accurately describe the relationship between T2 and pore size. However, the frequency alignment operation in Origin fitting is not limited by a model and can find a mapping relationship between T2 measured by NMR and pore size obtained by mercury intrusion porosimetry. This mapping relationship is not based on a specific, assumed model (e.g., a power function model), but on the common physical essence (statistical distribution of pore volume) reflected by the two measurement methods, resulting in better physical consistency. Therefore, in one embodiment of this invention, NMR-mercury intrusion porosimetry pore size mapping is further performed based on Origin fitting, thereby fitting a pore size-pore volume NMR reference image based on a two-dimensional relaxation spectrum. This is a well-known technique, and the general process is briefly described here:

[0103] First, select the data: In Origin, select data AX (cumulative distribution frequency of mercury intrusion porosimetry), BY (pore size of mercury intrusion porosimetry), and CY (cumulative distribution frequency data of T2). Then, use the "Analyze → Mathematics → Interpolate / Extrapolate Y from X" tool to interpolate or extrapolate the Y data from the given set of X values. The X value used for interpolation is the cumulative distribution frequency data of T2 in the CY column. Then, input [Bookl]Sheetl!(A,B), where A is the cumulative distribution frequency data of mercury intrusion porosimetry in column AX, and B is the pore size data of mercury intrusion porosimetry in column BY. In this embodiment, linear extrapolation is selected, and the generated new column is the equivalent pore size distribution after T2 relaxation time conversion.

[0104] Then, the obtained pore size-pore volume distribution curve is directly used as the pore size-pore volume mercury intrusion porosimetry reference map; the equivalent pore size after T2 relaxation time conversion is used as the X-axis parameter, and the sum of signal amplitudes in the pore size region corresponding to T2 relaxation time is used as the relative content of fluid in the pores as the Y-axis parameter, thereby fitting the pore size-pore volume NMR reference map; the mercury intrusion porosimetry reference map and the NMR reference map are compared, and pores with a pore size greater than >1μm in each map are considered large pores; when the distribution of large pores in the two maps is consistent, that is, the pore size and corresponding pore volume of large pores are consistent, it is determined that large pores have high confidence; otherwise, it is determined that the NMR-mercury intrusion porosimetry data conflict, and the experiment needs to be repeated.

[0105] The implementer can also analyze the pore volume and fluid content within a specific pore size based on mercury intrusion porosimetry reference maps and NMR reference maps. This is a well-known technique and will not be elaborated further.

[0106] In another embodiment of the present invention, the relaxation time distribution in the two-dimensional relaxation spectrum can be directly analyzed, and partitioning rules can be set to identify large pores; further, a pore size threshold can be set, and the pore size-pore volume distribution curve in the mercury intrusion porosimetry test results can be analyzed to identify large pores; finally, the dual discrimination results of NMR and mercury intrusion porosimetry are combined to obtain the final discrimination result.

[0107] For the two-dimensional relaxation spectrum obtained by inversion, a partitioning rule can be set: Free oil phase: distributed in the region of T1>100ms and T2>10ms, corresponding to free-flowing oil in microcracks or macropores; Bound water phase: distributed in the region of T1<10ms and T2<1ms, corresponding to bound water in clay or organic nanopores; Transitional fluid: the region of T1 in 10-100ms and T2 in 1-10ms, corresponding to light oil or mixed-phase fluid in organic pores; The sum of signal amplitudes corresponding to the free oil region is calculated by integration, and the sum of signal amplitudes corresponding to the free oil region is divided by the sum of all signal amplitudes to characterize the relative content of fluid in macropores;

[0108] For the mercury intrusion porosimetry results, the pore size threshold is set to 1 μm, and pores with a pore size r > 1 μm are considered as macropores. The cumulative mercury intrusion volume in the region with a pore size r > 1 μm is calculated by integration, and the cumulative mercury intrusion volume in the region with a pore size r > 1 μm is divided by the total mercury intrusion volume to characterize the structural volume of the macropores.

[0109] When both the two-dimensional relaxation spectrum obtained from nuclear magnetic resonance inversion and the mercury porosimetry results indicate the presence of macropores, it is determined that macropores exist in the shale sample.

[0110] In summary, this invention first obtains shale samples in situ; then, it performs two-dimensional nuclear magnetic resonance (NMR) testing on the shale samples based on preset NMR testing parameters, constructs a two-dimensional data matrix using the NMR test results, and inverts the two-dimensional relaxation spectrum using the truncated singular value decomposition (SVD) method. During the inversion process, the range of singular value cutoffs is determined based on the two-dimensional data matrix under different signal-to-noise ratios, and the optimal cutoff value within this range is obtained using a global optimization algorithm. The inversion result at the optimal cutoff value is used as the two-dimensional relaxation spectrum. High-pressure mercury intrusion porosimetry (HSI) is then performed on the shale samples based on preset mercury intrusion porosimetry (HIP) parameters. The macropore and fracture regions in the shale samples are determined by combining the HIP results and the two-dimensional relaxation spectrum. This invention improves the inversion effect by determining the range of singular value cutoffs based on NMR test results under different signal-to-noise ratios and using global optimization to determine the optimal cutoff value, thus obtaining a more reliable two-dimensional relaxation spectrum. Finally, the combined NMR-HIP test results are used to comprehensively identify macropores and fractures, improving the accuracy of identification.

[0111] It should be noted that the order of the above embodiments of the present invention is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. The processes depicted in the accompanying drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0112] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

Claims

1. A method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry, characterized in that, The method includes: Obtain shale samples in situ; Two-dimensional nuclear magnetic resonance (NMR) tests were performed on the shale sample based on preset NMR test parameters. A two-dimensional data matrix was constructed using the NMR test results, and a two-dimensional relaxation spectrum was inverted based on the truncated singular value decomposition method. In the inversion process, the range of singular value cutoff values ​​was determined based on different signal-to-noise ratios, and the optimal cutoff value within the range was obtained based on a global optimization algorithm. The inversion result under the optimal cutoff value was used as the two-dimensional relaxation spectrum. High-pressure mercury intrusion testing was performed on the shale sample based on preset mercury intrusion testing parameters. The macroporous fracture region in the shale sample was determined by combining the high-pressure mercury intrusion test results and the two-dimensional relaxation spectrum. The preset two-dimensional NMR scan parameters include at least the inversion recovery sequence parameters, the spin echo sequence parameters, and the number of scans; among which, the inversion recovery sequence parameters include the inversion time and the waiting time, and the spin echo sequence parameters include the echo interval and the number of echoes; Based on the combined results of high-pressure mercury intrusion porosimetry and the aforementioned two-dimensional relaxation spectrum, the macroporous fracture regions in the shale sample were determined, including: A pore size-pore volume reference map is fitted based on the high-pressure mercury intrusion porosimetry test results, and a pore size-pore volume reference map is fitted based on the two-dimensional relaxation spectrum. When both the mercury intrusion porosimetry reference map and the NMR reference map have pore sizes larger than a preset size, it is determined that there are large pores in the shale sample.

2. The method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry according to claim 1, characterized in that, The method for obtaining the two-dimensional data matrix includes: In the two-dimensional nuclear magnetic resonance test, the shale sample is scanned a preset number of times to obtain the CPMG echo string at each reversal time; the CPMG echo string at each reversal time is used as a row vector to construct a two-dimensional data matrix.

3. The method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry according to claim 2, characterized in that, The methods for obtaining the two-dimensional data matrix under different signal-to-noise ratios include: By changing the preset number of scans, the two-dimensional data matrix constructed under each preset number of scans is used as a two-dimensional data matrix with a signal-to-noise ratio.

4. The method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry according to claim 3, characterized in that, The method for obtaining the range of singular value cutoff values ​​includes: The two-dimensional data matrix at each signal-to-noise ratio is used as the test matrix; for each test matrix, a linear model is constructed based on the attenuation characteristics of the CPMG echo train with waiting time and acquisition time, and the coefficient matrix in the linear model is determined; wherein, the dependent variable of the linear model is the test matrix, and the independent variable of the linear model is the two-dimensional relaxation spectrum corresponding matrix to be solved. Based on a preset algorithm, the singular value cutoff value of the coefficient matrix corresponding to each test matrix is ​​determined, and the initial inversion of the corresponding linear model is performed based on the singular value cutoff value to obtain the proportion of negative components in the initial inversion result; based on the proportion of negative components corresponding to all test matrices, the range of values ​​for the singular value cutoff value is determined.

5. The method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry according to claim 4, characterized in that, Based on the proportion of the negative components corresponding to all the test matrices, the range of singular value cutoff values ​​is determined, including: All test matrices corresponding to the proportion of negative components are sorted and fitted in ascending order of signal-to-noise ratio to obtain the curve of change of negative component proportion; the extreme value among all abrupt change points in the curve of change of negative component proportion is determined, and the singular value cutoff value of the coefficient matrix of the test matrix corresponding to the extreme value is successively used as the upper and lower limits to obtain the range of singular value cutoff value.

6. The method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry according to claim 4, characterized in that, The preset algorithm includes at least the L-curve method.

7. The method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry according to claim 1, characterized in that, The method for obtaining the optimal cutoff value includes: A linear model of the two-dimensional data matrix is ​​constructed, and the linear model is inverted based on the truncated singular value decomposition method. Based on the difference between the inversion results before and after truncating the singular values ​​within the range of values, and the negative component in the inversion results after truncation, the target optimization function is determined. The optimal cutoff value within the range of values ​​is determined using the particle swarm optimization algorithm. The optimal cutoff value is the cutoff value corresponding to the minimum value of the objective optimization function.

8. The method for identifying macropores in shale based on nuclear magnetic resonance-mercury intrusion porosimetry according to claim 1, characterized in that, The preset mercury intrusion porosimetry test parameters include at least the termination condition during the high-pressure mercury intrusion porosimetry test, the pressure change step size of the applied pressure, and the pressure holding time at each pressure; wherein, the termination condition is that the applied pressure reaches the preset safe pressure threshold.

Images