Pore structure parameter calculation method based on pore type identification

Abstract

The present application relates to the technical field of oil and gas exploration, and particularly relates to a pore structure parameter calculation method based on pore type identification. First, pore types and face porosities are identified through well logging data and core experiments, and a nuclear magnetic T2 spectrum based on the face porosity is constructed. Then, based on the well logging data, pore types are identified and porosities of different pore types are calculated. Then, based on the calibrated nuclear magnetic T2 spectrum, matrix pore, dissolution pore and fracture pore nuclear magnetic T2 spectra are divided. Finally, based on the spherical tube pore structure model, the nuclear magnetic T2 spectra of different pore types are split into corresponding spherical pore spectra and tubular throat spectra, and corresponding pore structure parameters are calculated. Through identification of pore types and calculation of pore structure parameters, the technical scheme solves the problems that the laboratory image method is mostly qualitative analysis, the mercury injection experiment cannot distinguish pores and throats, the experimental measurement cost is high, data points are discrete, and the current nuclear magnetic resonance logging cannot distinguish different pore type nuclear magnetic resonance T2 spectra.

Description

Technical Field

[0001] This invention relates to the field of oil and gas exploration technology, and in particular to a method for calculating pore structure parameters based on pore type identification. Background Technology

[0002] In recent years, with the increasing demand for oil and natural gas, the difficulty of exploration and development has also increased. Oil and gas exploration and development has shifted from conventional oil and gas to unconventional oil and gas such as shale oil and tight oil. Since the pore structure of reservoir rocks has a significant impact on reservoir properties, oil and gas migration and storage capacity, and oilfield development, the pore structure of reservoirs is currently the focus and challenge of research on unconventional oil and gas exploration and development.

[0003] Currently, methods that can intuitively reflect and characterize the microscopic pore structure of reservoirs include laboratory measurement methods and nuclear magnetic resonance logging methods. Laboratory measurement methods include laboratory imaging methods and mercury intrusion porosimetry (MIP) analysis methods. However, laboratory imaging analysis methods are qualitative analyses and lack quantitative characterization of three-dimensional pore structures. The mercury intrusion curves obtained by mercury intrusion porosimetry cannot distinguish between pores and throats, nor can they reflect the heterogeneity of pore-throat distribution. Furthermore, laboratory measurement methods also suffer from high costs, data dispersion, and the inability to provide continuous evaluation. On the other hand, nuclear magnetic resonance logging methods cannot achieve nuclear magnetic T2 spectrum segmentation for different pore types and have relatively few studies on reservoir pore connectivity.

[0004] In summary, among the traditional methods for studying microscopic pore structures, laboratory imaging methods are mostly qualitative analyses, mercury intrusion porosimetry cannot distinguish between pores and throats, experimental measurement is costly and data points are discrete, and current nuclear magnetic resonance logging cannot distinguish between different pore types in its T2 spectra. Summary of the Invention

[0005] The purpose of this invention is to provide a method for calculating pore structure parameters based on pore type identification, which solves the problems of laboratory image methods being mostly qualitative analysis, mercury intrusion porosimetry being unable to distinguish between pores and throats, high experimental measurement costs and discrete data points, as well as the current inability of nuclear magnetic resonance logging to distinguish nuclear magnetic resonance T2 spectra of different pore types.

[0006] To achieve the above objectives, the present invention provides a method for calculating pore structure parameters based on pore type identification. The method includes the following steps:

[0007] S1: Basic processing of core data and well logging data to obtain NMR T2 spectrum, total porosity, matrix porosity, fracture NMR T2 spectrum, dissolution cavity NMR T2 spectrum and core surface relaxation rate.

[0008] S2: Based on well logging data, identify pore types and calculate the porosity of different pore types to obtain fracture porosity and dissolution cavity porosity;

[0009] S3: Based on the initial NMR T2 spectrum calibrated with total porosity, obtain the NMR T2 spectrum NT2 after total porosity calibration. Use the NMR T2 spectrum of dissolution pores calibrated with porosity to establish the NMR T2 spectrum of dissolution pores calibrated with porosity VT2. Use the NMR T2 spectrum of fractures calibrated with porosity to establish the NMR T2 spectrum of fractures calibrated with porosity FT2. Based on the NMR T2 spectrum calibrated with total porosity, the NMR T2 spectrum of fractures calibrated with porosity, and the NMR T2 spectrum of dissolution pores calibrated with porosity, obtain the NMR T2 spectrum BT2 reflecting the porosity of the matrix.

[0010] S4: Using the spherical tube pore structure model, NT2, FT2, VT2 and BT2 are decomposed into corresponding spherical pore spectra and tubular throat spectra, respectively. Based on the spherical pore spectra and tubular throat spectra of each pore type, the pore structure parameters of spherical pores and tubular throats for each pore type are calculated.

[0011] The specific content of step S1 is as follows:

[0012] Based on one-dimensional nuclear magnetic resonance logging data, the singular value decomposition method was used to obtain the nuclear magnetic resonance T2 spectrum;

[0013] Total porosity, secondary porosity, and matrix porosity are obtained based on conventional logging data including acoustic waves, density, and neutrons.

[0014] Based on electrical imaging logging data, pore types are identified, and fracture NMR T2 spectra based on fracture porosity are constructed by calculating fracture porosity, fracture opening, fracture porosity, and the actual size of pixels in electrical imaging images. Dissolution pore porosity NMR T2 spectra based on dissolution pore porosity are constructed by calculating the porosity, radius, and pore size distribution of dissolution pores.

[0015] The relaxation rate of the core surface was obtained based on the nuclear magnetic resonance experiment of the core.

[0016] Before the step of "identifying pore types based on electrical imaging logging data", it is necessary to fill the blank zone of electrical imaging. Specifically, for the static data of electrical imaging logging, the blank zone of electrical imaging is filled by the multi-point geostatistical FilterSim algorithm to fully reflect the distribution of wellbore fractures and dissolution pores.

[0017] In the step of "identifying pore types based on electrical imaging logging data", the specific identification principle is as follows: in electrical imaging logging data, the high-resistivity rock skeleton is displayed in a bright color, and the formation water exists in the dissolution pores and fracture pores, which have low resistivity and are displayed in a dark color in the electrical imaging data. The fracture image is a dark sine curve, and the dissolution pores appear as irregular clumps, spots, or near-circular shapes in the electrical imaging logging data.

[0018] The specific content of step S2 is as follows: based on the total porosity, matrix porosity and secondary porosity calculated from conventional logging data, the pore type is identified by combining electrical imaging logging data, and the porosity of fracture porosity and dissolution pores is calculated. The secondary porosity is divided according to the ratio of fracture porosity to dissolution pore porosity, and the porosity of fracture porosity and dissolution pores is calculated.

[0019] The specific steps of step S3 are as follows:

[0020] Based on the nuclear magnetic resonance T2 spectrum and total porosity, the amplitude of each relaxation time in the nuclear magnetic resonance T2 spectrum is calibrated to establish the nuclear magnetic resonance T2 spectrum NT2 after total porosity calibration.

[0021] Based on the NMR T2 spectrum of dissolution pores and the porosity of dissolution pores, the porosity of dissolution pores is calibrated by the NMR T2 spectrum of dissolution pores, thereby establishing the NMR T2 spectrum VT2 of dissolution pores after porosity calibration.

[0022] Based on the crack NMR T2 spectrum and crack porosity, the crack porosity is calibrated by the crack NMR T2 spectrum, thereby establishing the crack porosity calibrated crack porosity NMR T2 spectrum FT2.

[0023] The NMR T2 spectrum BT2, which reflects the matrix porosity, is obtained by subtracting the fracture porosity-calibrated fracture porosity NMR T2 spectrum FT2 and the dissolution cavity porosity-calibrated dissolution cavity NMR T2 spectrum VT2 from the total porosity-calibrated NMR T2 spectrum NT2.

[0024] The specific content of step S4 is as follows:

[0025] Based on the NMR T2 spectrum VT2 of the dissolved pores after porosity calibration, a spherical tube pore structure model is adopted. The pore space of the dissolved pores is decomposed into spherical pore bodies and tubular throats using the principle of equal division of surface area. The spherical pore body spectrum VTP2 and the tubular throat spectrum VTT2 of the NMR T2 spectrum of the dissolved pores are obtained. Then, based on the spherical pore body spectrum VTP2 and the tubular throat spectrum VTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats of the dissolved pores are calculated.

[0026] Based on the NMR T2 spectrum BT2 of the matrix pores, a spherical tube pore structure model is adopted. The pore space of the matrix pores is decomposed into spherical pore bodies and tubular throats by using the principle of equal division of surface area. Thus, the spherical pore body spectrum BTP2 and the tubular throat spectrum BTT2 of the matrix pore NMR T2 spectrum are obtained. Then, based on the spherical pore body spectrum BTP2 and the tubular throat spectrum BTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats in the matrix pores are calculated.

[0027] Based on the fracture porosity-calibrated fracture pore NMR T2 spectrum FT2, a spherical tube pore structure model is used. The pore space of the fracture pore is decomposed into spherical pore bodies and tubular throats using the principle of equal surface area division. This yields the spherical pore body spectrum FTP2 and the tubular throat spectrum FTT2 of the fracture pore NMR T2 spectrum. Then, based on the spherical pore body spectrum FTP2 and the tubular throat spectrum FTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats in the fracture pore are calculated.

[0028] This invention discloses a method for calculating pore structure parameters based on pore type identification. First, pore types and porosity are identified using well logging data and core experiments, and a NMR T2 spectrum based on porosity is constructed. Then, based on well logging data, pore types are identified and porosity is calculated for different pore types. Next, based on the porosity-calibrated NMR T2 spectrum, matrix pores, dissolution pores, and fracture pores are segmented into NMR T2 spectra. Finally, based on a pneumatic tube pore structure model, the NMR T2 spectra of different pore types are decomposed into corresponding spherical pore spectra and tubular throat spectra, and the corresponding pore structure parameters are calculated. This technical solution, by identifying pore types and calculating pore structure parameters, solves the problems of laboratory image methods being mostly qualitative analysis, mercury intrusion porosimetry being unable to distinguish between pores and throats, high experimental measurement costs, discrete data points, and the inability of current NMR logging to distinguish NMR T2 spectra of different pore types. Attached Figure Description

[0029] To more clearly illustrate the technical solutions 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.

[0030] Figure 1 This is a flowchart of the steps of the pore structure parameter calculation method based on pore type identification provided by the present invention.

[0031] Figure 2 The present invention provides a method for picking up dissolution cavities based on electrical imaging logging data and calculating the pore size distribution (first channel electrical imaging logging data, second channel picking up dissolution cavities, and third channel pore size distribution).

[0032] Figure 3 This invention provides the total NMR T2 spectrum at a certain depth, the fracture NMR T2 spectrum and the dissolution cavity T2 spectrum after porosity calibration based on the pore size distribution of fractures and dissolution pores, and the matrix porosity NMR T2 spectrum (a is the total porosity calibrated NMR T2 spectrum NT2 after NMR logging data despectroscopy; b is the dissolution cavity NMR T2 spectrum VT2 after dissolution cavity porosity calibration; c is the fracture porosity calibrated fracture porosity NMR T2 spectrum FT2; d is the matrix porosity NMR T2 spectrum BT2).

[0033] Figure 4 This is a schematic diagram of the pore structure model of the X-ray tube provided by the present invention.

[0034] Figure 5 The present invention provides the NMR T2 spectrum VT2 of the dissolution pores, and obtains the NMR T2 spectrum VTP2 of the spherical pore body and the NMR T2 spectrum VTT2 of the tubular throat through the optimal configuration relationship (a is the NMR T2 spectrum NT2 of the dissolution pores; b is the optimal configuration relationship of each group of spherical tube pore structure models; c is the NMR T2 spectrum VTP2 of the spherical pore body; d is the NMR T2 spectrum VTT2 of the tubular throat).

[0035] Figure 6 The matrix pore NMR T2 spectrum BT2 provided by this invention obtains the spherical pore body NMR T2 spectrum BTP2 and the tubular throat NMR T2 spectrum BTT2 through the optimal configuration relationship (a is the matrix pore NMR T2 spectrum BT2; b is the optimal configuration relationship of each group of spherical tube pore structure models; c is the spherical pore body NMR T2 spectrum BTP2; d is the tubular throat NMR T2 spectrum BTT2).

[0036] Figure 7 The present invention provides NMR T2 spectrum FT2 of crack pores, and obtains NMR T2 spectrum FTP2 of spherical pore bodies and NMR T2 spectrum FTT2 of tubular throats through optimal configuration relationship (a is NMR T2 spectrum FT2 of crack pores; b is the optimal configuration relationship of each group of spherical tube pore structure models; c is NMR T2 spectrum FTP2 of spherical pore bodies; d is NMR T2 spectrum FTT2 of tubular throats).

[0037] Figure 8 This invention provides a method for identifying pore type and pore structure parameters in the 6435-6485 section of a carbonate reservoir in an oilfield example well based on electrical imaging and nuclear magnetic resonance T2 spectroscopy. Detailed Implementation

[0038] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0039] Please see Figures 1 to 8 Taking the identification of pore type and calculation of pore structure parameters in a carbonate reservoir at a depth of 6435-6485 m in a typical well of an oilfield in China as an example, this invention provides a method for calculating pore structure parameters based on pore type identification. The method includes the following steps (e.g.) Figure 1 As shown):

[0040] S1: Basic processing of core data and well logging data to obtain NMR T2 spectrum, total porosity, matrix porosity, fracture NMR T2 spectrum, dissolution cavity NMR T2 spectrum and core surface relaxation rate.

[0041] S2: Based on well logging data, identify pore types and calculate the porosity of different pore types to obtain fracture porosity and dissolution cavity porosity;

[0042] S3: Based on the initial NMR T2 spectrum calibrated with total porosity, obtain the NMR T2 spectrum NT2 after total porosity calibration. Use the NMR T2 spectrum of dissolution pores calibrated with porosity to establish the NMR T2 spectrum of dissolution pores calibrated with porosity VT2. Use the NMR T2 spectrum of fractures calibrated with porosity to establish the NMR T2 spectrum of fractures calibrated with porosity FT2. Based on the NMR T2 spectrum calibrated with total porosity, the NMR T2 spectrum of fractures calibrated with porosity, and the NMR T2 spectrum of dissolution pores calibrated with porosity, obtain the NMR T2 spectrum BT2 reflecting the porosity of the matrix.

[0043] S4: Using the spherical tube pore structure model, NT2, FT2, VT2 and BT2 are decomposed into corresponding spherical pore spectra and tubular throat spectra, respectively. Based on the spherical pore spectra and tubular throat spectra of each pore type, the pore structure parameters of spherical pores and tubular throats for each pore type are calculated.

[0044] In this embodiment, pore types and porosity are first identified using well logging data and core experiments, and a nuclear magnetic resonance (NMR) T2 spectrum based on porosity is constructed. Then, based on well logging data, pore types are identified and porosity for different pore types is calculated. Next, based on the porosity-calibrated NMR T2 spectrum, matrix porosity, dissolution pores, and fracture pores are segmented into NMR T2 spectra. Finally, based on a pneumatic tube pore structure model, the NMR T2 spectra for different pore types are decomposed into corresponding spherical pore spectra and tubular throat spectra, and the corresponding pore structure parameters are calculated. This technical solution, by identifying pore types and calculating pore structure parameters, solves the problems of laboratory image methods being mostly qualitative analysis, mercury intrusion porosimetry being unable to distinguish between pores and throats, high experimental measurement costs and data point dispersion, as well as the inability of current NMR logging to distinguish NMR T2 spectra for different pore types. It plays an important role in reservoir evaluation, including the identification of physical properties and pore types in unconventional oil and gas reservoirs, and the characterization of pore microstructure, thus laying the foundation for evaluating oil and gas migration and storage capacity and oilfield development.

[0045] Furthermore, the specific content of step S1 is as follows:

[0046] S1.1. Based on one-dimensional nuclear magnetic resonance logging data, the nuclear magnetic resonance T2 spectrum is obtained using the singular value decomposition method;

[0047] S1.2. Obtain total porosity, secondary porosity, and matrix porosity based on conventional logging data including acoustic waves, density, and neutrons.

[0048] S1.3. Based on electrical imaging logging data, identify pore types, and construct a fracture NMR T2 spectrum based on fracture porosity by calculating fracture porosity, fracture opening, fracture porosity, and the actual size of pixels in the electrical imaging image. Construct a dissolution cavity NMR T2 spectrum based on dissolution cavity porosity by calculating the porosity, radius, and pore size distribution of the dissolution cavity.

[0049] S1.4. Based on core nuclear magnetic resonance experiments, obtain the relaxation rate of the core surface.

[0050] In this embodiment, the specific formula for obtaining the nuclear magnetic resonance T2 spectrum using singular value decomposition (SVD) in step S1.1 is as follows:

[0051] T2=(T 2i Am i ),(i=1,2,L,n)(0)

[0052] In the formula, T 2i The time interval (ms) or the value at the i-th placement point in the T2 spectrum is given by the time interval (m). i Let be the spectral amplitude corresponding to the i-th relaxation time.

[0053] In step S1.2, the three-porosity logging data in the conventional logging data provides formation porosity information. Among them, neutron logging porosity φN and density logging porosity φD reflect the total formation porosity, while sonic velocity logging mainly reflects the primary intergranular porosity φs.

[0054]

[0055] In the formula, φ D Porosity measured by density logging, decimal; ρ ma , ρ f , ρ b These are density values ​​measured from rock skeleton, pore fluid, and density logging, respectively, in g / cm³. 3 ;

[0056]

[0057] In the formula, φ N Porosity measured by neutron logging, decimal; Φ ma , Φ f , Φ N The figures represent the hydrogen content indices measured by rock skeleton, pore fluid, and neutron logging, respectively, in %;

[0058]

[0059] In the formula, φ s Porosity measured by neutron logging, decimal; Δt ma , Δt f Δ t These represent the time differences for rock skeleton, pore fluid, and logging time difference measurements, respectively, in μs / m.

[0060]

[0061] In the formula, φ t For total porosity, φ D φ N For density, porosity measured by neutron logging, φ fv For secondary (fractures and caves) porosity, φ b This refers to the matrix porosity.

[0062] Furthermore, before the step of "identifying porosity types based on electrical imaging logging data", it is necessary to fill the blank areas of electrical imaging. Specifically, for the static data of electrical imaging logging, the blank areas of electrical imaging are filled by the multi-point geostatistical FilterSim algorithm, so as to fully reflect the distribution of wellbore fractures and dissolution cavities.

[0063] In this embodiment, the FilterSim algorithm for multi-point geostatistics is a filter-based multi-point geostatistics method. The filter is a data template where each pixel location has a weight. The algorithm works as follows: static data from electrical imaging logging is used as image patches with known pixels, called training images. After filtering these image patches, k-means clustering is used to form multiple image patch template categories. During infilling, the nearest module category is selected by comparing the metric distance between the current missing image patch and each category of image patch templates. If there are no known pixels in the current missing image patch, a random pixel is selected from the image patch templates for infilling.

[0064] Furthermore, in the step of "identifying pore types based on electrical imaging logging data", the specific identification principle is as follows: in electrical imaging logging data, the high-resistivity rock skeleton is displayed in bright color, and the formation water exists in the dissolution pores and fracture pores, which have low resistivity and are displayed in dark color on the electrical imaging data. The fracture image is a dark sine curve, and the dissolution pores appear as irregular clumps, spots, or near-circular shapes in the electrical imaging logging data.

[0065] In this embodiment, the method for picking up cracks and dissolution cavities in electrical imaging data is as follows:

[0066] Median filtering is used to suppress image noise.

[0067] The Otsu method in image processing was used to identify conductivity anomalies corresponding to annotations and cavities in conductivity images, and to extract crack and dissolution cavity regions.

[0068] The cracks and dissolution cavities were picked up, and Hough and sine function fitting were used to identify the cracks, laying the foundation for the extraction of crack parameters.

[0069] After removing cracks, the developed areas of dissolution cavities were effectively divided on the image using the 8-connected domain labeling method, and each cavity was assigned a unique identifier, laying the foundation for the subsequent extraction of cavity parameters.

[0070] Furthermore, in step S1.3, by calculating the crack porosity, crack opening, crack porosity, and the actual size of the pixels in the electro-imaging image, the specific content of constructing the crack NMR T2 spectrum based on the crack porosity is as follows:

[0071] (1) Based on the electrical imaging data, the black sine curve is regarded as an open fracture. The fracture is picked in the electrical imaging logging software ciflog. The dip angle qng (°) and dip direction azi_ang (°) of the fracture are directly obtained by fitting the Hough transform and the sine function.

[0072] (2) Based on the resistivity anomaly area A of the crack in the electrical imaging data, the crack opening is calculated according to the following formula:

[0073]

[0074] In the formula, W represents the crack width in mm; R represents the crack width in mm. xo - Resistivity of the intrusion zone, Ω·m; A - Area of ​​conductivity anomaly caused by the crack, mm / Ω·m; R m - Mud resistivity, Ω.m; C, b - Parameters related to instrument structure. In this embodiment, C = 0.004801, b = 0.863, which are related to the instrument structure; Z is the displacement perpendicular to the crack trajectory direction; Z0 is the basic half-width, which is generally greater than half of the crack width shown in the figure; V e I represents the plate potential value. a I is the electrode current (mA), a function of the instrument's vertical position Z as it crosses the crack. b It is the electrode current in the original strata or framework.

[0075] (3) Based on the resistivity anomaly range and crack opening of the crack in the electro-imaging data, calculate the actual size Pixs of each pixel in the electro-imaging image:

[0076]

[0077] In the formula, W is the fracture width in mm; pixel is the number of pixels of the abnormal current area A in the Z direction of the electrical imaging logging data; and pixel is the length of each voxel in the electrical imaging image in mm / pixel.

[0078] (4) Calculate the porosity or porosity of the crack based on the electrical imaging data. The porosity of the i-th crack is expressed as:

[0079]

[0080] In the formula, Sφ fi For i crack porosity or porosity, W i L is the average width of the i-th crack; i Let L be the length of the i-th crack within the statistical window; L is generally chosen as 1m or 0.6096m, and D is the well diameter.

[0081] The total porosity or void fraction of the cracks is expressed as:

[0082]

[0083] In the formula, Sφ f denoted as crack porosity; M represents the total number of cracks.

[0084] (5) Calculate the crack radius and tortuosity based on the crack opening and the arc length of the crack sine curve;

[0085] The crack tortuosity τwfp and crack radius Radwfp are respectively:

[0086]

[0087] In the formula, τ wfp qng represents the crack tortuosity, qng represents the crack dip angle (°), and W represents the crack opening; Rad wfp Where is the crack radius.

[0088] (6) Based on the crack radius Rad wfp And the relationship between nuclear magnetic resonance T2 relaxation time and pore size (Equation 14), calculate the nuclear magnetic relaxation time corresponding to the crack radius, and the crack aperture Sφ f (Formula 10) is used as the amplitude of the corresponding relaxation time, thereby establishing the NMR T2 spectrum corresponding to the crack (Formula 15);

[0089]

[0090] In the formula, T2 wfpi To calculate the T2 relaxation time (ms) of the NMR for the i-th slit within a statistical window length L; Rad wfpi Let be the radius of the i-th crack within the statistical window length L, in μm; ρ be the surface relaxation rate, in μm / ms; F s F is a geometric shape factor, which is related to the pore shape. For spherical pores, F s =3; for columnar pores, F s =2;

[0091] T2f0=(T2 wfpi ,Sφ fi ), i = 1, 2, L, M (15)

[0092] Where M is the total number of cracks within the statistical window length L.

[0093] Furthermore, in step S1.3, the specific content of constructing the NMR T2 spectrum of the dissolution pores based on the porosity, radius, and pore size distribution of the dissolution pores is as follows:

[0094] (1) Calculation of porosity (pore ratio) and equivalent radius of dissolution cavities:

[0095] For a specific identified dissolution cavity, within the range of abnormal conductivity distribution of the cavity, the total number of pixels within the distribution range and the number of pixels at the edges are calculated using point statistics. This allows for the calculation of the porosity (pore density) and the equivalent radius of the dissolution cavity. The specific formulas are as follows:

[0096]

[0097] In the formula, Sφ vi To calculate the porosity or porosity of the i-th dissolution cavity within a window length L, VP i TP represents the number of pixels within the i-th dissolution hole; TP represents the total number of pixels in the electrical imaging within the statistical window length L of the dissolution holes; L is generally chosen as 1m or 0.6096m.

[0098]

[0099] In the formula, Rad wvpi To calculate the equivalent radius of the i-th dissolution cavity within the statistical window length L; CP i is the number of pixels at the edge of the i-th dissolution hole, reflecting the perimeter of the dissolution hole; pixels is the actual size of each pixel in the electro-imaging, mm / pixel.

[0100] (2) Characterization of the total porosity of dissolution pores within the statistical window length: The porosity of dissolution pores within the statistical window length is equal to the sum of the porosities of all dissolution pores within that window length, i.e.:

[0101]

[0102] In the formula, Sφ v To statistically analyze the porosity or void ratio of dissolution cavities within a window length L.

[0103] (3) Characterization of the pore size distribution of dissolution cavities within the statistical window length, based on the porosity Sφ of all dissolution cavities within the statistical window length. vi And the equivalent radius Rad of the dissolution cavity wvpi Radv: Pore size distribution of dissolution cavities within the statistical window length

[0104] Radv = (Rad wvpi ,Sφ vi ), i = 1, 2, L, N (19)

[0105] Where N is the total number of dissolution cavities within the statistical window length L.

[0106] (4) The NMR T2 distribution corresponding to the pore size distribution of dissolution cavities within the statistical window length. Based on the pore size distribution of dissolution cavities within the statistical window length, Radv = (Rad wvpi ,Sφ vi ), i=1,2,L,N, combining the rock surface relaxation rate ρ and formula (20), establish the radius Rad of the i-th dissolution cavity. wvpi With NMR T2 wvpiThe corresponding relationship of relaxation time is used to establish the T2 spectrum distribution of dissolution pores within the statistical window length (Equation 21):

[0107]

[0108] Where T2 is the nuclear magnetic resonance relaxation time, ms; ρ is the surface relaxation rate, μm / ms; and S is the pore surface area, μm. 2 V is the pore volume, in μm. 3 Fs is the geometric shape factor, which is related to the pore shape. For spherical pores, Fs = 3; for columnar pores, Fs = 2; r is the equivalent radius of the pore, in μm.

[0109] T2v0=(T2 wvpi ,Sφ vi ), i = 1, 2, L, N (21)

[0110] Among them, T2 wvpi The NMR T2 relaxation time is used to calculate the i-th dissolution cavity within the statistical window length L.

[0111] Furthermore, in step S1.4, based on the core nuclear magnetic resonance experiment, the specific details of the relaxation rate of the core surface are as follows:

[0112] (1) ln(I / Io) and -γ were obtained through PGSE experiments. 2 g 2 δ 2 The (Δ-δ / 3) relationship is used to obtain the diffusion coefficient D under different conditions by applying formula (22):

[0113]

[0114] In the formula, I is the echo amplitude; I0 is the initial echo amplitude; D is the diffusion coefficient; γ is the gyromagnetic ratio; g is the magnetic field gradient; δ is the gradient pulse width; Δ is the propagation time; in the PGSE sequence experiment, a series of echo signals can be obtained by changing g, δ or Δ.

[0115] (2) The specific surface area (S / V) is obtained by fitting the following formula (23):

[0116]

[0117] In the formula, D(t) is the apparent diffusion coefficient; D0 is the free diffusion coefficient; It is a higher-order quantity of diffusion time, which is related to diffusion time t, surface relaxation rate ρ, and tortuosity τ.

[0118] (3) Calculate the harmonic mean (1 / T2) of T2 using the T2 spectrum of the CPMG sequence.

[0119] (4) The surface relaxation rate ρ is obtained by substituting (S / V) into formula (24);

[0120]

[0121] In the formula, ρ is the surface relaxation rate, S / V is the specific surface area, and T2 is the NMR relaxation time.

[0122] Furthermore, the specific content of step S2 is as follows: based on the total porosity, matrix porosity and secondary porosity calculated from conventional logging data, the pore type is identified by combining electrical imaging logging data, and the porosity of fracture porosity and dissolution pores is calculated. The secondary porosity is divided according to the ratio of fracture porosity to dissolution pore porosity, and the porosity of fracture porosity and dissolution pores is calculated.

[0123] In this embodiment, based on the secondary porosity φ at the same depth fv Crack surface ratio Sφ f and the porosity Sφ of dissolution pores v Secondary porosity is classified as fracture porosity φ based on the ratio of the porosity of cracks to that of dissolution cavities. f and the porosity φ of the dissolution pores v (Formula 25):

[0124]

[0125] Furthermore, the specific steps of step S3 are as follows:

[0126] S3.1 Based on the nuclear magnetic resonance T2 spectrum and total porosity, the amplitude of each relaxation time in the nuclear magnetic resonance T2 spectrum is calibrated to establish the nuclear magnetic resonance T2 spectrum NT2 after total porosity calibration.

[0127] S3.2. Based on the NMR T2 spectrum of the dissolution pores and the porosity of the dissolution pores, the porosity of the dissolution pores is calibrated by the NMR T2 spectrum of the dissolution pores, thereby establishing the NMR T2 spectrum VT2 of the dissolution pores after the porosity calibration.

[0128] S3.3. Based on the crack NMR T2 spectrum and crack porosity, the crack porosity is calibrated on the crack NMR T2 spectrum, thereby establishing the crack porosity-calibrated crack porosity NMR T2 spectrum FT2.

[0129] S3.4 Subtract the fracture porosity-calibrated fracture porosity NMR T2 spectrum FT2 and the dissolution cavity porosity-calibrated dissolution cavity NMR T2 spectrum VT2 from the total porosity-calibrated NMR T2 spectrum NT2 to obtain the NMR T2 spectrum BT2 reflecting the matrix porosity.

[0130] In this embodiment, step S3.1 specifically involves: based on the nuclear magnetic resonance T2 spectrum, T2 = (T 2iAm i ), (i = 1, 2, ..., n), T 2i The time interval (ms) or the value at the i-th placement point in the T2 spectrum is given by the time interval (m). i The spectral amplitude corresponding to the i-th relaxation time; the T2 spectrum after porosity calibration, NT2 = (T 2i ,NAm i ), (i = 1, 2, ..., n), T 2i The i-th relaxation time or the i-th placement point value in the T2 spectrum, in milliseconds; NAM i denoted as the porosity-calibrated spectral amplitude corresponding to the i-th relaxation time, and n is the number of NMR relaxation times in the T2 spectrum.

[0131] Step S3.2 specifically includes:

[0132] (1) Based on the T2 NMR spectrum of the dissolved pores, T2v0, T2v0=(T2 wvpi ,Sφ vi ), i = 1, 2, ..., N, (where T2 wvpi To calculate the NMR T2 relaxation time of the i-th dissolution cavity within the statistical window length L; Sφ vi To determine the porosity or porosity of the i-th dissolution pore within the statistical window length L, where N is the total number of dissolution pores within the statistical window length L, cubic spline interpolation is performed to establish the NMR T2 spectrum T2v1 (Formula 26) for dissolution pores with the same NMR relaxation time distribution range (T2i, i = 1, 2, ..., n) as the well logging NMR T2 spectrum relaxation time distribution range.

[0133] T2v1=(T 2i ,Sφ vi ),i=1,2,L,n (26);

[0134] In the formula, T 2i Sφ is the nuclear magnetic relaxation time after cubic spline interpolation. vi T2i represents the amplitude corresponding to the i-th relaxation time after cubic spline interpolation; n is the number of relaxation times after cubic spline interpolation.

[0135] (2) Recalibrate the amplitude Sφ of nuclear magnetic relaxation time in the T2 NMR spectrum. vi (i = 1, 2, ..., n) is Vm i (Formula 27) After establishing the porosity calibration of the dissolved pores, the NMR T2 spectrum is established, VT2=(T 2i Vm i ), (i = 1, 2, ..., n), T 2i Vm represents the i-th relaxation time or the i-th placement point value in the T2 spectrum, in milliseconds; i The spectral amplitude after porosity calibration corresponds to the i-th relaxation time;

[0136] Vm i =fvam×Sφ vi (i = 1, 2, ..., n)(27)

[0137] In the formula, Vm i Sφ represents the spectral amplitude corresponding to the i-th relaxation time in the NMR T2 spectrum after porosity calibration. vi denoted as fvam, where fvam is the conversion coefficient and n is the number of NMR relaxation times in the T2 spectrum of the dissolution pores.

[0138] Step S3.3 specifically includes:

[0139] (1) Based on the crack NMR T2 spectrum T2f0, T2f0=(T2 wfpi ,Sφ fi ), i = 1, 2, L, M, (where T2 wfpi To calculate the T2 relaxation time of the NMR spectrum for the i-th crack within the statistical window length L; Sφ fi To determine the porosity or void fraction of the i-th fracture within a statistical window length L, where M is the total number of fractures within L, cubic spline interpolation is performed to establish that the NMR relaxation time distribution range is the same as that of the well logging NMR T2 spectrum relaxation time distribution range (T2...). i The NMR T2 spectrum T2f1 of the dissolution cracks (i = 1, 2, ..., n) (Formula 28);

[0140] T2f1=(T 2i ,Sφ fi ), i = 1, 2, ..., n (28)

[0141] In the formula, T 2i Sφ is the nuclear magnetic relaxation time after cubic spline interpolation. fi The relaxation time T after cubic spline interpolation is the i-th time. 2i The corresponding amplitude; n is the number of relaxation times after cubic spline interpolation.

[0142] (2) Recalibrate the amplitude Sφ of nuclear magnetic relaxation time in the T2 NMR spectrum. fi (i = 1, 2, ..., n) is Fm i (Formula 29) Establish the NMR T2 spectrum after crack porosity calibration, FT2=(T 2i ,Fm i ), (i = 1, 2, ..., n), T 2i Fm represents the i-th relaxation time or the i-th placement point value in the T2 spectrum, in milliseconds (ms); i The spectral amplitude after porosity calibration corresponds to the i-th relaxation time;

[0143] Fm i =ffam×Sφfi (i = 1, 2, ..., n) (29)

[0144] In the formula, Fm i Sφ represents the spectral amplitude corresponding to the i-th relaxation time in the NMR T2 spectrum after porosity calibration. fi denoted as , where is the spectral amplitude corresponding to the i-th relaxation time in the T2 spectrum of the crack pore nuclear magnetic resonance (NMR) T2f1, ffam is the conversion coefficient, and n is the number of NMR relaxation times in the T2 spectrum.

[0145] Furthermore, the specific content of step S3.4 is as follows: based on the porosity-calibrated nuclear magnetic resonance T2 spectrum, NT2, NT2=(T 2i ,NAm i ), (i=1,2,L,n), fracture T2 spectrum after fracture porosity calibration, FT2, FT2=(T 2i ,Fm i ), (i=1,2,L,n) and the T2 spectrum of the dissolved pores after the porosity calibration, VT2, VT2=(T 2i Vm i ), (i=1,2,L,n), obtain the matrix porosity T2 spectrum by subtracting the fracture T2 spectrum and the dissolution pore T2 spectrum from the total NMR T2 spectrum, BT2, BT2=(T 2i ,Bm i ), (i = 1, 2, ..., n);

[0146] Bm i =NAM i -Fm i -Vm i (i = 1, 2, ..., n) (30)

[0147] In the formula, Bm i The spectral amplitude corresponding to the i-th NMR relaxation time in the matrix pore T2 spectrum; NAM i The i-th relaxation time T2 in the NMR T2 spectrum i Corresponding spectral amplitude; Fm i The i-th relaxation time T2 in the crack porosity T2 spectrum i Corresponding spectral amplitude; Vm i The i-th relaxation time T2 in the T2 spectrum of the dissolution cavity i The corresponding spectral amplitude.

[0148] Furthermore, the specific content of step S4 is as follows:

[0149] S4.1: Based on the NMR T2 spectrum VT2 of the dissolved pores after the porosity calibration, a spherical tube pore structure model is adopted. The pore space of the dissolved pores is decomposed into spherical pore bodies and tubular throats using the principle of equal division of surface area. The spherical pore body spectrum VTP2 and the tubular throat spectrum VTT2 of the NMR T2 spectrum of the dissolved pores are obtained. Then, based on the spherical pore body spectrum VTP2 and the tubular throat spectrum VTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats of the dissolved pores are calculated.

[0150] S4.2: Based on the NMR T2 spectrum BT2 of the matrix pores, a spherical tube pore structure model is adopted. The pore space of the matrix pores is decomposed into spherical pore bodies and tubular throats using the principle of equal division of surface area. Thus, the spherical pore body spectrum BTP2 and the tubular throat spectrum BTT2 of the matrix pore NMR T2 spectrum are obtained. Then, based on the spherical pore body spectrum BTP2 and the tubular throat spectrum BTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats in the matrix pores are calculated.

[0151] S4.3: Based on the fracture porosity-calibrated fracture porosity NMR T2 spectrum FT2, a spherical tube pore structure model is adopted. The pore space of the fracture pore is decomposed into spherical pore bodies and tubular throats using the principle of equal surface area division. Thus, the spherical pore body spectrum FTP2 and the tubular throat spectrum FTT2 of the fracture pore NMR T2 spectrum are obtained. Then, based on the spherical pore body spectrum FTP2 and the tubular throat spectrum FTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats in the fracture pore are calculated.

[0152] The basic unit of pore space is the sphere-tube pore structure model, whose length is twice the equivalent pore radius Rad. e The pore structure model of the cyclone tube consists of two parts: a spherical pore body and a cylindrical throat. The radius of the spherical pore is Rad. s The radius of the tubular larynx is Rad c (like Figure 4 Formula (31) gives the pore radius Rad of the sphere in the sphere pore structure model. s With equivalent pore radius Rad e The relationship, formula (32) gives the equivalent pore radius Rad of the i-th group of sphere pore structure model. ei With the i-th relaxation time T 2i The quantitative relationship; Cd represents the ratio of the tubular throat radius to the sphere pore radius in the sphere pore structure model (Formula 33), reflecting the spatial configuration relationship of the sphere pore structure model; the surface area of ​​the sphere pore is SFs, the surface area of ​​the tubular throat is SFc, and the total pore space surface area of ​​the sphere pore structure model is SFt;

[0153]

[0154] In the formula, Cs represents the relationship between the average radius of the grouped pores and the radius of the spherical pores, and is generally taken as 3;

[0155] Rad ei =3ρ×T 2i (32)

[0156] In the formula, T 2i Rad represents a relaxation time in the T2 spectrum, in milliseconds (ms); ei ρ is the equivalent spherical aperture radius, in μm; ρ is the surface relaxation rate, in μm / ms.

[0157] Cd = Rad c / Rad s , Cd∈[0,1] (33)

[0158] In the formula, Cd represents the ratio of the tubular throat radius to the pore radius of the sphere in the PV tube pore structure model, which is called the configuration relationship of the PV tube pore structure model. It is usually taken as six discrete values: 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0; Rad c Rad is the radius of the tubular larynx. s Let be the radius of the pore in the sphere.

[0159] Furthermore, the specific content of step S4.1 includes:

[0160] Step 4.1.1: Based on the NMR T2 spectrum VT2 of the dissolution pores, VT2 = (T 2i Vm i ), (i=1,2,L,n), and calculate the equivalent pore radius representing the pore space in the n groups of sphere tube pore structure models according to formula (32), and calculate the spherical pore radius and tubular throat radius under all spatial configuration relationships Cd in each group of sphere tube pore structure models according to formula (31) and formula (33).

[0161] Step 4.1.2: Based on the equivalent pore radius, spherical pore radius and tubular throat radius under all spatial configuration relationships Cd in the n-group spherical tube pore structure model, calculate the surface area of ​​the spherical pore and the surface area of ​​the tubular throat, and divide the spectral amplitude of the NMR T2 spectrum VT2 of the dissolution pore into the spectral amplitudes corresponding to the spherical pore and the tubular throat according to the ratio of the surface area of ​​the spherical pore to the surface area of ​​the tubular throat (Formulas 34 and 35);

[0162]

[0163] In the formula, Vm i The relaxation time T2 in the NMR T2 spectrum of the dissolution pores is indicated. i The corresponding spectral amplitude, SF vcik SF vsik SF vtikThese represent the surface areas of the tubular throats, spherical pores, and total pore space in the i-th group of dissolution pore sphere pore structure models under the current k-th configuration relationship, respectively (i = 1, 2, ..., n; k = 1, 2, ..., 6n); i is the number of groups of sphere pore structures; k is the number of configuration relationships between spherical pores and tubular throats under each group of sphere pore structure models, which is usually 6.

[0164] Step 4.1.3: Based on the radii of the spherical pores and tubular throats of the n-group X-ray tube models, calculate the corresponding NMR T2 time for the spherical pores and the NMR T2 time for the tubular throats according to formula (14 or 20).

[0165] Step 4.1.4: Based on the T2 time of NMR of the spherical pore and the T2 time of NMR of the tubular throat of the n-group X-ray tube models, as well as the corresponding spectral amplitude of the spherical pore and the tubular throat, the spherical pore volume spectrum and the tubular throat spectrum are established.

[0166] VT 2sk =(T 2vsik Vm sik (i = 1, 2, ..., n; k = 1, 2, ..., 6) n (36)

[0167] VT 2ck =(T 2vcik Vm cik (i = 1, 2, ..., n; k = 1, 2, ..., 6) n (37)

[0168] In the formula, VT 2sk The spherical pore volume spectrum of the PV tube pore structure model, VT 2ck Tubular throat spectrum of the pore structure model of a PV tube; T2 vsik Vm represents the NMR relaxation time of the spherical pore under the configuration relationship Cd of the k-th pore structure in the i-th PV tube pore structure model. sik T2 represents the spectral amplitude corresponding to the NMR relaxation time of the spherical pores under the configuration relationship Cd in the i-th group of X-ray tube pore structure models; vcik Vm represents the NMR relaxation time of the tubular throat under Cd in the configuration relationship of the k-th pore structure of the i-th PV tube in the PV tube pore structure model. cik The spectral amplitude corresponding to the NMR relaxation time of the tubular throat under the configuration relationship Cd in the i-th group of X-ray tube pore structure models is represented by: i represents the number of X-ray tube models; k represents the number of configuration relationships in each group of X-ray tube pore structure models.

[0169] Step 4.1.5: Based on the NMR T2 spectrum of the spherical pores according to the configuration relationship between the spherical pores and the tubular throats in all the spherical tube pore structure models, (V T2s(spectrum) and the nuclear magnetic resonance T2 spectrum characterizing the tubular larynx, (V T2c (spectrum), calculate the echo train characterizing the pore space of dissolution cavities, VM sck (Formula 38);

[0170]

[0171] In the formula, t is the echo train time; VM sck This is the echo train calculated based on the sphere pore spectrum and tubular throat spectrum under the current k-th configuration relationship, using the X-ray tube pore structure model.

[0172] Step 4.1.6: Based on the T2 NMR spectrum and VT2 of the dissolution pores, calculate the echo train VM that characterizes the pore space of the dissolution pores. ori (Formula 39), and calculate the difference between the two characterizations of the echo train in the pore space of the dissolution cavity (Formula 39).

[0173] Formula 40);

[0174]

[0175] In the formula, t is the echo train time; VM ori Echo train calculated from the inverse NMR T2 spectrum of the dissolution pores;

[0176] VM errork =|VM sck -VM ori |,(k=1,2,L,6 n (40)

[0177] In the formula, k represents the k-th configuration relationship, VM sck This represents the echo train calculated from the pore spectrum of the dissolved pore sphere and the pore spectrum of the tubular throat obtained based on the pore structure model of the sphere tube under the current k-th configuration relationship.

[0178] Step 4.1.7: Select the spherical pore volume spectrum and tubular throat channel spectrum corresponding to the minimum echo train error as the optimal spherical pore spectrum VTP2 and the optimal tubular throat channel spectrum VTT2 of the NMR T2 spectrum of the dissolution pores.

[0179] VTP2=(T 2vsi * Vm si * ),(i=1,2,L,n)(41)

[0180] VTT2=(T 2vci * Vm ci * ),(i=1,2,L,n) (42)

[0181] In the formula, T2vsi * represents the pore size V of the sphere. T2s The i-th relaxation time or the i-th placement point value in the spectrum, in milliseconds; V msi * represents the spectral amplitude corresponding to the i-th relaxation time; T 2vci * indicates the tubular larynx (V) T2s The i-th relaxation time or the i-th placement point value in the spectrum, in milliseconds; Vm ci * represents the spectral amplitude corresponding to the i-th relaxation time.

[0182] Step 4.1.8: Based on the optimal configuration relationship of the spherical pore spectrum VTP2 and the optimal tubular throat spectrum VTT2, calculate the pore radius and tortuosity of the spherical pores in the dissolution cavity, and the pore radius and tortuosity of the tubular throat.

[0183] (1) Based on the spherical pore spectrum VTP2 of the optimal configuration relationship, calculate the arithmetic mean VT of the T2 relaxation time. 2SAM and the geometric mean of T2 relaxation time VT 2SGM ;

[0184] (2) Based on the tubular laryngeal spectrum VTT2 of the optimal configuration relationship, calculate the arithmetic mean VT of the T2 relaxation time. 2CAM and the geometric mean of T2 relaxation time VT 2CGM ;

[0185] (3) Arithmetic mean VT of T2 relaxation time based on spherical pores 2SAM and the geometric mean of T2 relaxation time VT 2SGM The radius of the spherical pores in the dissolution cavity is calculated using formula (14 or 20), and the tortuosity of the spherical pores in the dissolution cavity is calculated based on the ratio of the arithmetic mean to the geometric mean.

[0186]

[0187] In the formula, VT 2SAM The T2 relaxation time of the pore spectrum of the dissolved cavity sphere is expressed in ms; VT 2SGM Vτ represents the geometric mean of the T2 relaxation time of the pore spectrum of the dissolved cavity sphere, in milliseconds; wp The tortuosity of the sphere's pores is characterized by the x-tube model; vca and vcb are the parameters of the model.

[0188] (4) Arithmetic mean VT of T2 relaxation time based on tubular larynx 2CAM and the geometric mean of T2 relaxation time VT 2CGM The pore radius of the tubular throat in the dissolution cavity is calculated using formula (14 or 20), and the tortuosity of the tubular throat in the dissolution cavity is calculated based on the ratio of the arithmetic mean to the geometric mean.

[0189]

[0190] In the formula, VT 2CAM The arithmetic mean of the T2 relaxation times of the tubular laryngeal tract spectrum, in milliseconds (ms); VT 2CGM Vτ represents the geometric mean of the T2 relaxation time of the tubular laryngeal tract spectrum, in milliseconds (ms); wt VCA represents the tortuosity of the tubular throat in the X-ray tube model; VCA and VCB are the parameters of the model (as shown in Table 1, where VCA = 2.798 and VCB = 0.3072).

[0191] Furthermore, the specific content of step S4.2 includes:

[0192] Step 4.2.1: Based on the matrix porosity NMR T2 spectrum BT2, BT2 = (T 2i ,Bm i ), (i=1,2,L,n), and calculate the equivalent pore radius representing the pore space in the n groups of sphere tube pore structure models according to formula (32), and calculate the spherical pore radius and tubular throat radius of all configuration relationships Cd in each group of sphere tube pore structure models according to formula (31) and formula (33).

[0193] Step 4.2.2: Based on the equivalent pore radius, spherical pore radius, and tubular throat radius under all spatial configuration relationships Cd in the n sets of X-ray tube pore structure models, calculate the surface area of ​​the spherical pores and the surface area of ​​the tubular throats, and then calculate the spectral amplitude Bm of the matrix pore NMR T2 spectrum BT2 according to the ratio of the surface area of ​​the spherical pores to the surface area of ​​the tubular throats. i The spectral amplitudes are divided into those corresponding to spherical pores and tubular throats (Formulas 45 and 46);

[0194]

[0195] In the formula, Bm i SF represents the spectral amplitude corresponding to the relaxation time T2i in the matrix pore NMR T2 spectrum. bcik SF bsik SF btik These represent the surface areas of the tubular throats, spherical pores, and total pore space in the i-th group of matrix pore sphere pore structure models under the current k-th configuration relationship, respectively (i = 1, 2, ..., n; k = 1, 2, ..., 6n); i is the number of groups of sphere pore structures; k is the number of configuration relationships between spherical pores and tubular throats under each group of sphere pore structure models, which is usually 6.

[0196] Step 4.2.3: Based on the radii of the spherical pores and tubular throats of the n sets of X-ray tube models, calculate the corresponding NMR T2 time for the spherical pores and the NMR T2 time for the tubular throats according to formula (14 or 20).

[0197] Step 4.2.4: Based on the T2 time of NMR for the spherical pores and the T2 time of NMR for the tubular throats of the n sets of X-ray tube models, as well as the corresponding spectral amplitudes of the spherical pores and tubular throats, the spherical pore volume spectrum and the tubular throat spectrum are established.

[0198] BT 2sk =(T 2bsik ,Bm sik (i = 1, 2, ..., n; k = 1, 2, ..., 6) n (47)

[0199] BT 2ck =(T 2bcik ,Bm cik (i = 1, 2, ..., n; k = 1, 2, ..., 6) n (48)

[0200] In the formula, BT 2sk The spherical pore volume spectrum of the PV tube pore structure model, BT 2ck Tubular throat spectrum for the pore structure model of a PV tube; T 2bsik Bm represents the NMR relaxation time of the spherical pore under the configuration relationship Cd of the k-th sphere pore structure in the i-th sphere pore structure model. sik T represents the spectral amplitude corresponding to the NMR relaxation time of the spherical pores under the configuration relationship Cd in the i-th group of X-ray tube pore structure models; 2bcik Bm represents the NMR relaxation time of the tubular throat under Cd in the configuration relationship of the k-th pore structure of the i-th PV tube in the PV tube pore structure model. cik The spectral amplitude corresponding to the NMR relaxation time of the tubular throat under the configuration relationship Cd in the i-th group of X-ray tube pore structure models is represented by: i represents the number of X-ray tube models; k represents the number of configuration relationships in each group of X-ray tube pore structure models.

[0201] Step 4.2.5: Based on the NMR T2 spectrum of the spherical pores according to the configuration relationship between the spherical pores and the tubular throats in all the spherical tube pore structure models, (B T2s (spectrum) and the nuclear magnetic resonance T2 spectrum characterizing the tubular larynx, (B T2c (spectrum), calculate the echo train characterizing the matrix pore space, BM sck (Formula 49);

[0202]

[0203] In the formula, t is the echo train time; BM sck This is the echo train calculated based on the sphere pore spectrum and tubular throat spectrum obtained from the X-ray tube model under the current k-th configuration relationship.

[0204] Step 4.2.6: Based on the matrix pore NMR T2 spectrum and BT2, calculate the echo train BM characterizing the matrix pore space. ori (Formula 50), and calculate the difference in pore space echo trains between the two characterizing matrix pores (Formula 51);

[0205]

[0206] In the formula, t is the echo train time; BM ori The echo train calculated from the inverse T2 NMR spectrum of the matrix pores;

[0207] BM errork =|BM sck -BM ori |,(k=1,2,L,6 n (51)

[0208] In the formula, k represents the k-th configuration relationship, BM sck This represents the echo train calculated from the matrix pore sphere pore spectrum and the tubular throat pore spectrum obtained based on the sphere pore structure model under the current k-th configuration relationship.

[0209] Step 4.2.7: Select the spherical pore volume spectrum and the tubular throat channel spectrum corresponding to the minimum echo train error as the optimal spherical pore spectrum BTP2 and the optimal tubular throat channel spectrum BTT2 of the matrix pore NMR T2 spectrum;

[0210] BTP2=(T 2bsi * ,Bm si * ),(i=1,2,L,n) (52)

[0211] BTT2=(T 2bci * ,Bm ci * ),(i=1,2,L,n) (53)

[0212] In the formula, T 2bsi * represents the i-th relaxation time or the i-th point value in the BTP2 spectrum of the spherical pores, in milliseconds; B msi * represents the spectral amplitude corresponding to the i-th relaxation time; T 2bci * represents the i-th relaxation time or the i-th location value in the BTT2 spectrum of the tubular larynx, in milliseconds; Bm ci * represents the spectral amplitude corresponding to the i-th relaxation time.

[0213] Step 4.2.8: Based on the optimal configuration relationship of the spherical pore spectrum BTP2 and the optimal tubular throat spectrum BTT2, calculate the pore radius and tortuosity of the spheres and the pore radius and tortuosity of the tubular throat in the matrix pore space.

[0214] (1) Calculate the arithmetic mean BT of the T2 relaxation time based on the spherical pore spectrum BTP2 of the optimal configuration relationship. 2SAM Geometric mean of T2 relaxation time BT 2SGM ;

[0215] (2) Based on the tubular laryngeal tract spectrum BTT2 of the optimal configuration relationship, calculate the arithmetic mean BT of the T2 relaxation time. 2CAM Geometric mean of T2 relaxation time BT 2CGM ;

[0216] (3) Arithmetic mean BT of T2 relaxation time based on spherical pores 2SAM Geometric mean of T2 relaxation time BT 2SGM The radius of the spherical pores in the matrix pores is calculated using formula (14 or 20), and the tortuosity of the spherical pores in the matrix pores is calculated based on the ratio of the arithmetic mean to the geometric mean.

[0217]

[0218] In the formula, BT 2SAM The BT value represents the arithmetic mean of the T2 relaxation time of the pore spectrum of the matrix pore spheres, in ms; 2SGM Bτ represents the geometric mean of the T2 relaxation time of the pore spectrum of matrix-porous spheres, in ms; wp The tortuosity of the sphere's pores is characterized by the cycloidal tube model; bca and bcb are the parameters of the model.

[0219] (4) Arithmetic mean BT of T2 relaxation time based on tubular larynx 2CAM Geometric mean of T2 relaxation time BT 2CGM The pore radius of the tubular throat in the matrix pores is calculated using formula (14 or 20), and the tortuosity of the tubular throat in the matrix pores is calculated based on the ratio of the arithmetic mean to the geometric mean.

[0220]

[0221] In the formula, BT 2CAM The arithmetic mean of the T2 relaxation time of the tubular laryngeal tract spectrum, in milliseconds (ms); BT 2CGM Bτ represents the geometric mean of the T2 relaxation time of the tubular laryngeal tract spectrum, in milliseconds; wt denoted as the tortuosity of the tubular throat in the X-ray tube model; bca and bcb are the parameters of the model (as shown in Table 1, where bca = 2.798; bcb = 0.3072).

[0222] Furthermore, the specific content of step S4.3 includes:

[0223] Step 4.3.1: Based on the fracture pore NMR T2 spectrum FT2, FT2 = (T 2i ,Fm i ), (i=1,2,L,n), and calculate the equivalent pore radius representing the pore space in the n groups of sphere tube pore structure models according to formula (32), and calculate the spherical pore radius and tubular throat radius of all configuration relationships Cd in each group of sphere tube pore structure models according to formula (31) and formula (33).

[0224] Step 4.3.2: Based on the equivalent pore radius, spherical pore radius, and tubular throat radius calculated from all configuration relationships Cd in the n sets of X-ray tube pore structure models, calculate the surface area of ​​the spherical pore and the surface area of ​​the tubular throat, and then calculate the spectral amplitude Fm of the crack pore NMR T2 spectrum FT2 according to the ratio of the surface area of ​​the spherical pore to the surface area of ​​the tubular throat. i The spectral amplitudes are divided into those corresponding to spherical pores and tubular throats (Formulas 56 and 57);

[0225]

[0226] In the formula, Fm i SF represents the spectral amplitude corresponding to T2i during relaxation in the T2 NMR spectrum of the crack pores. fcik SF fsik SF ftik These represent the surface areas of the tubular throat, the spherical pores, and the total pore space in the i-th group of the fracture-pore sphere pore structure model under the current k-th configuration relationship, respectively (i = 1, 2, ..., n; k = 1, 2, ..., 6n); i is the number of groups of sphere pore structures; k is the number of configuration relationships between spherical pores and tubular throats under each group of sphere pore structure models, which is usually 6.

[0227] Step 4.3.3: Based on the radii of the spherical pores and tubular throats of the n sets of X-ray tube models, calculate the corresponding NMR T2 time for the spherical pores and the NMR T2 time for the tubular throats according to formula (14 or 20).

[0228] Step 4.3.4: Based on the T2 time of NMR of the spherical pore and the T2 time of NMR of the tubular throat of the n-group X-ray tube models, as well as the corresponding spectral amplitude of the spherical pore and the tubular throat, the spherical pore volume spectrum and the tubular throat spectrum are established.

[0229] FT 2sk =(T 2fsik ,Fm sik (i = 1, 2, ..., n; k = 1, 2, ..., 6) n (58)

[0230] FT 2ck =(T 2fcik ,Fmcik (i = 1, 2, ..., n; k = 1, 2, ..., 6) n (59)

[0231] In the formula, FT 2sk The spherical pore volume spectrum of the PV tube pore structure model, FT 2ck Tubular throat spectrum of the pore structure model of a PV tube; T2 fsik Fm represents the NMR relaxation time of the spherical pore under the configuration relationship Cd of the k-th pore structure in the i-th PV tube pore structure model. sik T2 represents the spectral amplitude corresponding to the NMR relaxation time of the spherical pores under the configuration relationship Cd in the i-th group of X-ray tube pore structure models; fcik Fm represents the NMR relaxation time of the tubular throat under Cd in the configuration relationship of the k-th pore structure of the i-th PV tube in the PV tube pore structure model. cik The spectral amplitude corresponding to the NMR relaxation time of the tubular throat under the configuration relationship Cd in the i-th group of X-ray tube pore structure models is represented by: i represents the number of X-ray tube models; k represents the number of configuration relationships in each group of X-ray tube pore structure models.

[0232] Step 4.3.5: Based on the NMR T2 spectrum of the spherical pores according to the configuration relationship between the spherical pores and the tubular throats in all the spherical tube pore structure models, (F T2s (F) and the nuclear magnetic resonance T2 spectrum characterizing the tubular larynx, (F) T2c FM spectrum, used to calculate the echo train characterizing the fracture pore space. sck (Formula 60);

[0233]

[0234] In the formula, t is the echo train time; FM sck This is the echo train calculated from the spherical pore spectrum and tubular throat spectrum obtained based on the X-ray tube model under the current k-th configuration relationship.

[0235] Step 4.3.6: Based on the fracture pore NMR T2 spectrum and FT2, calculate the echo train FM characterizing the fracture pore space. ori (Formula 61), and calculate the difference in pore space echo trains between the two characterizing crack pores (Formula 62);

[0236]

[0237] In the formula, t is the echo train time; FM ori The echo train calculated from the T2 NMR spectrum of the crack pores;

[0238] FM errork =|FM sck -FM ori|,(k=1,2,L,6 n (62)

[0239] In the formula, k represents the k-th configuration relationship, FM sck This represents the echo train calculated from the crack pore sphere pore spectrum and the tubular throat pore spectrum obtained based on the sphere pore structure model under the current k-th configuration relationship.

[0240] Step 4.3.7: Select the spherical pore volume spectrum and the tubular throat channel spectrum corresponding to the minimum echo train error as the optimal spherical pore volume spectrum FTP2 and the optimal tubular throat channel spectrum FTT2 of the crack NMR T2 spectrum;

[0241] FTP2 = (T 2fsi * ,Fm si * ),(i=1,2,L,n)(63)

[0242] FTT2=(T 2fci * ,Fm ci * ),(i=1,2,L,n)(64)

[0243] In the formula, T2 fsi * represents the i-th relaxation time or the i-th point value in the FTP2 spectrum of the spherical pores, in milliseconds; Fm si * represents the spectral amplitude corresponding to the i-th relaxation time; T 2fci * represents the i-th relaxation time or the i-th position value in the FTT2 spectrum of the tubular larynx, in milliseconds; Fm ci * represents the spectral amplitude corresponding to the i-th relaxation time.

[0244] Step 4.3.8: Based on the optimal configuration relationship of the spherical pore spectrum FTP2 and the optimal tubular throat spectrum FTT2, calculate the pore radius and tortuosity of the spherical pores and the pore radius and tortuosity of the tubular throat in the crack pores.

[0245] (1) Calculate the arithmetic mean FT of the T2 relaxation time based on the spherical pore spectrum FTP2 of the optimal configuration relationship. 2SAM and the geometric mean of T2 relaxation time (FT) 2SGM ;

[0246] (2) Based on the tubular throat spectrum FTT2 of the optimal configuration relationship, calculate the arithmetic mean of the T2 relaxation time FT. 2CAM and the geometric mean of T2 relaxation time (FT) 2CGM ;

[0247] (3) Arithmetic mean of T2 relaxation time based on spherical pores (FT) 2SAMand the geometric mean of T2 relaxation time (FT) 2SGM The radius of the spherical pore in the crack is calculated using formula (14 or 20), and the tortuosity of the spherical pore in the crack is calculated based on the ratio of the arithmetic mean to the geometric mean.

[0248]

[0249] In the formula, FT 2SAM The arithmetic mean of the T2 relaxation time of the pore spectrum of the fractured sphere, in milliseconds; FT 2SGM Fτ represents the geometric mean of the T2 relaxation time of the pore spectrum of the cracked sphere, in milliseconds. wp The tortuosity of the sphere's pores is characterized by the x-tube model; fca and fcb are the parameters of the model.

[0250] (4) Arithmetic mean of T2 relaxation time based on tubular larynx (FT) 2CAM and the geometric mean of T2 relaxation time (FT) 2CGM The pore radius of the tubular throat in the fracture pore is calculated using formula (14 or 20), and the tortuosity of the tubular throat in the fracture pore is calculated based on the ratio of the arithmetic mean to the geometric mean.

[0251]

[0252] In the formula, FT 2CAM The arithmetic mean of the T2 relaxation times of the tubular laryngeal tract spectrum, in milliseconds; FT 2CGM Fτ represents the geometric mean of the T2 relaxation time of the tubular laryngeal tract spectrum, in milliseconds (ms); wp fca represents the tortuosity of the tubular throat in the X-ray tube model; fca and fcb are the parameters of the model (as shown in Table 1, where fca = 2.798; fcb = 0.3072).

[0253] like Figure 8 As shown, channel 3 contains conventional acoustic, density, and neutron logging data; channel 6 contains electrical imaging logging data; channels 7 and 8 contain the distribution of dissolution cavities identified by electrical imaging and their calculated pore size distribution; channel 9 contains the NMR T2 spectrum; channel 10 contains the porosity calculated based on electrical imaging logging data; channel 11 contains the porosity curve; and channel 11 contains the calculated pore structure parameters: the radius of spherical pores in dissolution cavities and the radius of spherical pores in matrix pores.

[0254] Table 1. Pore structure parameters calculated based on the NMR T2 spectra of dissolution pores (VT2), matrix pores (BT2), and fracture pores (FT2).

[0255]

[0256] In summary:

[0257] (1) To address the problem that laboratory image methods are qualitative analyses and lack quantitative characterization of three-dimensional pore structures, this technical solution is adopted: combining nuclear magnetic resonance T2 spectroscopy with electrical imaging logging to achieve quantitative calculation of pore structure parameters;

[0258] (2) To address the problem that the mercury intrusion curves obtained by the mercury intrusion method cannot distinguish between pores and throats and cannot reflect the heterogeneity of pore and throat distribution, this technical solution is adopted: nuclear magnetic resonance data and the pore structure model of the tube are used to calculate the spectrum of spherical pore bodies and tubular throats of the same pore type, and further calculate the pore structure parameters of spherical pores and tubular throats of the same pore type.

[0259] (3) To address the problems of high cost, discrete data, and inability to conduct continuous evaluation using laboratory measurement methods, this technical solution is proposed: a well logging evaluation method based on nuclear magnetic resonance T2 spectrum combined with electrical imaging data is proposed to realize pore type identification and pore structure parameter calculation.

[0260] (4) In view of the current problem that nuclear magnetic resonance logging cannot achieve the nuclear magnetic T2 spectrum segmentation of different pore types, this technical solution is adopted: combining nuclear magnetic resonance and electrical imaging data to realize the construction of nuclear magnetic T2 spectrum of dissolution pores and fracture pores, and to segment the matrix pores, dissolution pores and fracture pores from the total nuclear magnetic T2 spectrum.

[0261] (5) In view of the relatively limited research on reservoir pore connectivity in current nuclear magnetic resonance logging, this technical solution is adopted: combining nuclear magnetic resonance and electrical imaging data to obtain spherical pore spectrum and tubular throat spectrum of different pore types, and calculating the tortuosity of spherical pores and throat tortuosity to reflect pore connectivity characteristics.

[0262] Therefore, this technical solution, based on electrical imaging and nuclear magnetic resonance (NMR) T2 spectroscopy, achieves pore type classification and pore structure parameter calculation, solving the problems of laboratory imaging methods being mostly qualitative analysis, mercury intrusion porosimetry being unable to distinguish between pores and throats, high experimental measurement costs and data point dispersion, and the inability of current NMR logging to distinguish NMR T2 spectra of different pore types. With the continuous development of oil and gas exploration, complex reservoirs have gradually become the focus of oil and gas exploration and development. Pore structure is an important factor affecting reservoir fluid flow, restricting the accumulation and flow of oil and gas in the reservoir. Therefore, the study of reservoir pore structure is currently a key focus of oil and gas development research in complex reservoirs. This invention identifies pore types and further classifies the same pore type into spherical pores and tubular throats, calculating the corresponding pore structure parameters, thereby achieving pore type classification and fine characterization of pore structure. It plays an important role in reservoir evaluation, including physical properties, pore type identification, and pore microstructure characterization, both now and in the future, thus laying the foundation for evaluating oil and gas migration and accumulation capacity and oilfield development.

[0263] The above description discloses only one preferred embodiment of the present invention, and should not be construed as limiting the scope of the present invention. Those skilled in the art will understand that all or part of the processes of the above embodiments can be implemented, and equivalent changes made in accordance with the claims of the present invention are still within the scope of the invention.

Claims

1. A method for calculating pore structure parameters based on pore type identification, characterized in that, Includes the following steps: S1: Basic processing of core data and well logging data to obtain NMR T2 spectrum, total porosity, matrix porosity, fracture NMR T2 spectrum, dissolution cavity NMR T2 spectrum and core surface relaxation rate. S2: Based on well logging data, identify pore types and calculate the porosity of different pore types to obtain fracture porosity and dissolution cavity porosity; S3: Based on the initial NMR T2 spectrum calibrated with total porosity, obtain the NMR T2 spectrum NT2 after total porosity calibration. Use the NMR T2 spectrum of dissolution pores calibrated with porosity to establish the NMR T2 spectrum of dissolution pores calibrated with porosity VT2. Use the NMR T2 spectrum of fractures calibrated with porosity to establish the NMR T2 spectrum of fractures calibrated with porosity FT2. Based on the NMR T2 spectrum calibrated with total porosity, the NMR T2 spectrum of fractures calibrated with porosity, and the NMR T2 spectrum of dissolution pores calibrated with porosity, obtain the NMR T2 spectrum BT2 reflecting the porosity of the matrix. S4: Using the spherical tube pore structure model, NT2, FT2, VT2 and BT2 are decomposed into corresponding spherical pore spectra and tubular throat spectra, respectively. Based on the spherical pore spectra and tubular throat spectra of each pore type, the pore structure parameters of spherical pores and tubular throats for each pore type are calculated.

2. The pore structure parameter calculation method based on pore type identification as described in claim 1, characterized in that, The specific content of step S1 is as follows: Based on one-dimensional nuclear magnetic resonance logging data, the singular value decomposition method was used to obtain the nuclear magnetic resonance T2 spectrum; Total porosity, secondary porosity, and matrix porosity are obtained based on conventional logging data including acoustic waves, density, and neutrons. Based on electrical imaging logging data, pore types are identified, and fracture NMR T2 spectra based on fracture porosity are constructed by calculating fracture porosity, fracture opening, fracture porosity, and the actual size of pixels in electrical imaging images. Dissolution pore porosity NMR T2 spectra based on dissolution pore porosity are constructed by calculating the porosity, radius, and pore size distribution of dissolution pores. The relaxation rate of the core surface was obtained based on the nuclear magnetic resonance experiment of the core.

3. The pore structure parameter calculation method based on pore type identification as described in claim 2, characterized in that, Before the step "Identify porosity type based on electrical imaging logging data", it is necessary to fill the blank zone of electrical imaging. Specifically, for static data of electrical imaging logging, the blank zone of electrical imaging is filled by the multi-point geostatistical FilterSim algorithm to fully reflect the distribution of wellbore fractures and dissolution cavities.

4. The pore structure parameter calculation method based on pore type identification as described in claim 3, characterized in that, In the step "Identifying Pore Types Based on Electro-Imaging Logging Data", the specific identification principle is as follows: In the electro-imaging logging data, the high-resistivity rock skeleton is displayed in a bright color. Formation water exists in the dissolution pores and fracture pores, which have low resistivity and are displayed in a dark color in the electro-imaging data. Fractures are imaged as dark sinusoidal curves. Dissolution pores appear as irregular clumps, spots, or near-circular shapes in the electro-imaging logging data.

5. The pore structure parameter calculation method based on pore type identification as described in claim 4, characterized in that, The specific content of step S2 is as follows: Based on the total porosity, matrix porosity and secondary porosity calculated from conventional logging data, the porosity type is identified by combining electrical imaging logging data, and the porosity ratios of fracture porosity and dissolution porosity are calculated. The secondary porosity is divided according to the ratio of fracture porosity to dissolution porosity, and the porosity of fracture porosity and dissolution porosity are calculated.

6. The pore structure parameter calculation method based on pore type identification as described in claim 5, characterized in that, The specific steps of step S3 are as follows: Based on the nuclear magnetic resonance T2 spectrum and total porosity, the amplitude of each relaxation time in the nuclear magnetic resonance T2 spectrum is calibrated to establish the nuclear magnetic resonance T2 spectrum NT2 after total porosity calibration. Based on the NMR T2 spectrum of dissolution pores and the porosity of dissolution pores, the porosity of dissolution pores is calibrated by the NMR T2 spectrum of dissolution pores, thereby establishing the NMR T2 spectrum VT2 of dissolution pores after porosity calibration. Based on the crack NMR T2 spectrum and crack porosity, the crack porosity is calibrated by the crack NMR T2 spectrum, thereby establishing the crack porosity calibrated crack porosity NMR T2 spectrum FT2. The NMR T2 spectrum BT2, which reflects the matrix porosity, is obtained by subtracting the fracture porosity-calibrated fracture porosity NMR T2 spectrum FT2 and the dissolution cavity porosity-calibrated dissolution cavity NMR T2 spectrum VT2 from the total porosity-calibrated NMR T2 spectrum NT2.

7. The pore structure parameter calculation method based on pore type identification as described in claim 6, characterized in that, The specific content of step S4 is as follows: Based on the NMR T2 spectrum VT2 of the dissolved pores after porosity calibration, a spherical tube pore structure model is adopted. The pore space of the dissolved pores is decomposed into spherical pore bodies and tubular throats using the principle of equal division of surface area. The spherical pore body spectrum VTP2 and the tubular throat spectrum VTT2 of the NMR T2 spectrum of the dissolved pores are obtained. Then, based on the spherical pore body spectrum VTP2 and the tubular throat spectrum VTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats of the dissolved pores are calculated. Based on the NMR T2 spectrum BT2 of the matrix pores, a spherical tube pore structure model is adopted. The pore space of the matrix pores is decomposed into spherical pore bodies and tubular throats by using the principle of equal division of surface area. Thus, the spherical pore body spectrum BTP2 and the tubular throat spectrum BTT2 of the matrix pore NMR T2 spectrum are obtained. Then, based on the spherical pore body spectrum BTP2 and the tubular throat spectrum BTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats in the matrix pores are calculated. Based on the fracture porosity-calibrated fracture pore NMR T2 spectrum FT2, a spherical tube pore structure model is used. The pore space of the fracture pore is decomposed into spherical pore bodies and tubular throats using the principle of equal surface area division. This yields the spherical pore body spectrum FTP2 and the tubular throat spectrum FTT2 of the fracture pore NMR T2 spectrum. Then, based on the spherical pore body spectrum FTP2 and the tubular throat spectrum FTT2, the pore radius and tortuosity of the spherical pore bodies and tubular throats in the fracture pore are calculated.

Images