A method for joint inversion of reservoir parameters using resistivity and nuclear logging under mud invasion

Abstract

The present invention is a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion, which can be used to simulate logging curves and reservoir physical parameter inversion at any time under mud invasion conditions. First, a wellbore numerical model is constructed, and then the wellbore numerical model is assigned values ​​according to the formation reservoir characteristics, and then the assigned numerical model is numerically simulated to obtain the resistivity response, neutron response and density response curves of the reservoir, and the measured logging curve is compared with the numerical simulation response curve, and the assignment of the numerical model parameters is continuously adjusted until the numerical simulation curve is basically consistent with the measured logging curve, and finally the real reservoir parameters under formation conditions are obtained. This method is based on rock physical parameters combined with wellbore information and phase permeability characteristic curves, and integrates seepage characteristics, electromagnetic fields and nuclear logging principles based on instrument principles. The actual reservoir parameters of the reservoir are inverted through multi-time, multi-physical field and multi-dimensional logging response information, thereby improving the accuracy of reservoir fluid property identification and parameter evaluation.

Description

Technical Field

[0001] The invention relates to the technical field of reservoir logging evaluation in oil and gas exploration and development, and in particular to a method for jointly inverting reservoir parameters using resistivity and nuclear logging under mud invasion. Background Art

[0002] At present, the reserves of low-permeability oil and gas fields in China are large, but such oil and gas reservoirs often have poor reservoir properties due to heavy local mud, heavy gray or complex pore structure, so the natural production capacity is low. During artificial mining, after the drilling is invaded by mud, the logging curve response is affected by many factors, resulting in a small difference in resistivity between oil and gas and water layers, which brings great difficulties to reservoir evaluation.

[0003] At present, conventional fluid identification and reservoir parameter measurement are mainly used to simulate different logging projects for different targets. For example, an existing joint inversion method for resistivity logging while drilling based on parallel computing technology includes the following steps: A. Loading resistivity logging while drilling data; B. Smooth filtering and automatic stratification of logging curves; C. Selecting the invasion depth of the formation model, the resistivity of the invasion zone, and the initial value of the formation resistivity according to the degree of curve separation; D. Using the finite element method, based on multi-core parallel or GPU parallel computing technology, calculate the resistivity logging response of the corresponding depth point of the measured curve; E. Calculating the relative error between the forward simulation response value and the measured data. If the relative error is less than the set threshold, the model value is output as the inversion result. If not, the model change is calculated according to the Marquette iterative algorithm, the formation model parameters are reset, and the step D is returned until the inversion result is output.

[0004] However, the fluid properties and pore structure parameters obtained by simulating different logging projects for different targets using conventional fluid identification and reservoir parameter measurement have low accuracy. Summary of the invention

[0005] The present invention provides a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion to solve the problems of the above technical solutions. This solution can carry out joint simulation of resistivity and nuclear logging based on rock physical parameters combined with wellbore information and phase permeability characteristic curves, and obtain reliable reservoir parameters through joint inversion of multi-dimensional logging response information such as multi-time and multi-physical fields, thereby improving the accuracy of reservoir fluid property identification and parameter evaluation, and providing technical support for oil and gas exploration and development.

[0006] The technical solution adopted by the present invention is: a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion, comprising the following steps:

[0007] Step 1: construct a wellbore numerical model; layer the underground reservoir measured logging curve model according to the reservoir gamma, resistivity, neutron, and density curve characteristics, and then square the underground reservoir measured logging curve by layer, each layer represents a rock physics model;

[0008] Step 2: Setting the initial value of the wellbore numerical model; inputting the rock physics model parameters, seepage information and wellbore information parameters into the rock physics model in step 1, wherein the rock physics model parameters and seepage information parameters are assigned by layer, and each layer inputs the parameters into the rock physics model in step 1 according to the formation reservoir characteristics, and the wellbore information parameters are processed according to the whole well section, that is, a unified value is input;

[0009] Step 3: Dynamic numerical simulation: According to the relevant porous media seepage theory and diffusion theory, the initial numerical model of the wellbore obtained in step 2 is subjected to a dynamic invasion simulation of drilling fluid, and the dynamic fluid distribution profile of the reservoir at any time is obtained, including the dynamic saturation and salinity profile of the reservoir. The obtained dynamic fluid distribution profile of the reservoir is converted by a relevant formula to obtain the corresponding radial distribution profiles of resistivity, migration length, and density at different times;

[0010] Step 4: Numerical simulation of resistivity, neutron and density measurement responses; obtain the numerically simulated resistivity logging response based on the electromagnetic field finite element simulation method combined with the instrument structure parameters and working principle, and obtain the numerically simulated neutron and density measurement curves based on the instrument structure parameters and nuclear logging fast simulation method;

[0011] Step 5: Invert reservoir parameters; compare the measured logging curve of the underground reservoir in step 1 with the numerical simulation curve obtained in step 4, and make the numerical simulation curve consistent with the measured logging curve through model updating and iteration, and finally obtain the reservoir parameters of the inverted rock physical structure.

[0012] In the above technical scheme, the wellbore numerical model is first constructed, and then the wellbore numerical model is assigned and numerically simulated according to the formation reservoir characteristics, and then the measured logging curve is compared with the numerical simulation curve, and the assignment of the numerical model parameters is continuously adjusted until the numerical simulation curve is basically consistent with the measured logging curve, thereby obtaining the real reservoir parameters under formation conditions. The processing method of this scheme can carry out resistivity and nuclear logging joint dynamic simulation and inversion to obtain accurate fluid properties and pore structure parameters of the reservoir.

[0013] Preferably, in step 1, the logging curve is a record of the change of the physical properties of the formation with the depth of the well, and the top and bottom interfaces of each sublayer are divided by the natural gamma half-amplitude point, the reference resistivity, neutron and density curves. Squaring means averaging the numerical segments of the logging curve in a graphical display mode to form a sawtooth shape, and the average value of the logging curve of each sublayer is used to represent the sublayer.

[0014] Preferably, in step 2, the input model parameters include rock physics model parameters, seepage information and wellbore information parameters. The rock physics model parameters and seepage parameters are processed by layer, and the wellbore information parameters are processed according to the whole well section. The specific parameters of the rock physics model include rock skeleton parameters, formation physical property parameters, formation pore fluid parameters, and mud content. Rock skeleton parameters include mineral components, content and density; formation physical property parameters include saturation, permeability, porosity and pore structure; formation pore fluid parameters include fluid components, content, density, mineralization and viscosity; seepage information includes phase permeability curve and capillary pressure curve; wellbore information parameters include wellbore diameter, mud information, drill string pressure and temperature. Rock skeleton parameters are obtained through X-ray diffraction data, pore fluid information including fluid composition, density, mineralization, and viscosity is obtained through water analysis experimental materials, seepage information is obtained through regional displacement experiments, wellbore information is obtained through drilling logs and geological daily reports, formation pore structure parameters are obtained through resistivity core experiments, and shale content, porosity, permeability, and saturation are obtained through the average value of the logging curve of each small layer using the following logging interpretation formula. The formula for obtaining the initial value of shale content is:

[0015]

[0016]

[0017] In the formula, I sh Indicates the mud content index, decimal; V sh Indicates the formation mud content; decimal; GR indicates the formation natural gamma measurement value, gAPI; GR max Indicates the maximum value of natural gamma of pure mudstone in this well, gAPI; GR min The minimum natural gamma value of pure sandstone in this well, gAPI.

[0018] The formula for obtaining the initial porosity value is:

[0019]

[0020]

[0021]

[0022] Where: Indicates the calculated density porosity, v / v; DEN indicates the density curve logging reading, g / cm3; DEN ma Indicates the density response parameter of the skeleton, g / cm3; DEN f Indicates the density response parameter of the fluid, g / cm3; V sh Indicates the mud content of the formation, v / v; DEN sh Represents the density response parameter of mudstone, g / cm3; Indicates the calculated neutron porosity, v / v; CNC indicates the neutron curve logging reading, %; CNC ma Indicates the neutron response parameter of the skeleton, %; CNC f Indicates the neutron response parameter of the fluid, %; CNC sh represents the neutron response parameter of mudstone, %; Represents formation porosity, v / v.

[0023] The permeability calculation can be done using the regional core porosity-permeability regression model. For example, the model formula for the study area is:

[0024]

[0025] Where: K represents the formation permeability, mD; Represents the formation porosity, v / v.

[0026] The formula for obtaining the initial value of water saturation is:

[0027]

[0028] Where: S w Indicates the water saturation of the formation; R t Indicates true resistivity of formation, Ω.m; V sh Indicates the mud content of the formation, decimal; R sh Represents the resistivity of pure mudstone, Ω.m; R w Indicates the resistivity of formation water, Ω.m; It represents the formation porosity, a decimal; a represents the lithology coefficient; m represents the cementation index; n represents the saturation index. a, m, and n are obtained using block empirical values ​​or other logging interpretation methods.

[0029] Preferably, in step three, the phase permeability curve and the capillary pressure curve are obtained through core experiments. According to the relevant permeability theory and ion diffusion theory including the phase permeability curve and the capillary pressure curve, the dynamic distribution profile of the reservoir fluid including the formation saturation and salinity at any time is simulated, and the corresponding radial distribution profiles of resistivity, migration length, and density at different times are obtained through conversion of resistivity and nuclear logging related formulas. The resistivity radial profile conversion formula is:

[0030] Where: R t Represents formation resistivity; R w Represents the resistivity of formation water; S w Indicates the water saturation of the formation; Represents formation porosity; R sh Represents the resistivity of mud; V sh represents the shale content; a represents the lithology coefficient; m represents the cementation index; n represents the saturation index; c represents the shale content index coefficient. The general calculation formula is: a, m, and n are obtained using block experience or other well logging interpretation methods.

[0031] The density radial profile conversion formula is:

[0032]

[0033] Where: b represents the formation density; ρ w represents the formation water density; ρ h Indicates the density of oil and gas fluid; represents the formation porosity; ρ sh Indicates the density of mud; C sh Indicates the shale porosity; S w Indicates the formation water saturation; ρ i The density of the i-th component in the skeleton; C i represents the content of the i-th component in the skeleton; n c Represents the total number of component types in the skeleton.

[0034] The calculation formula of the radial profile of neutron migration length is:

[0035]

[0036] Where: ξ represents the formation migration length; S w Indicates the water saturation of the formation; ξ w represents the migration length of formation water; ξ h It represents the migration length of oil and gas fluid; represents the formation porosity; ξ m represents the skeleton migration length; η represents the exponential coefficient.

[0037] Preferably, in step 4, the resistivity logging response is obtained according to the radial distribution profile of resistivity converted in step 3, and the corresponding resistivity logging response is obtained mainly by using the electromagnetic field finite element simulation method combined with the instrument structure parameters and the numerical simulation of the working principle. The three-dimensional vector finite element method is used to spatially discretize the logging physical model, establish the unit electromagnetic field discrete equation, install, eliminate and solve all units, and then obtain the electrical logging instrument response. Resistivity logging simulation is to solve the Maxwell equations under given boundary conditions. Starting from Maxwell's equations, the electromagnetic field satisfies the following equation:

[0038]

[0039]

[0040] Where: E represents the electric field intensity; H represents the magnetic field intensity; J is the source current density; ω is the angular frequency; σ is the electrical conductivity; μ is the magnetic permeability.

[0041] From the above two equations, we can deduce that the vector wave equation satisfied by the electric field is:

[0042]

[0043] is the complex dielectric constant; ε = ε r ε0, where ε0 is the dielectric constant of vacuum; ε r is the relative dielectric constant.

[0044] make

[0045] E=E p +E s

[0046] Where: Background field E p is the electric field when the entire space is filled with a medium with conductivity σ0, which satisfies the equation:

[0047]

[0048] in, By changing the above three equations, we can get the vector wave equation:

[0049]

[0050] Background Field E p By analytical calculation, the secondary field E s The finite element method is used for calculation. The solution of the above equation changes smoothly, and a sparse grid can be used to solve it, reducing the computational workload. Select a sufficiently large area so that the electric field on the boundary decays to approximately 0, then the above equation only needs to satisfy the boundary condition equation:

[0051]

[0052] Where: To solve the boundary of the region ω, n is the direction of its normal.

[0053] Considering the boundary condition equation, the vector wave equation is transformed into its weak product form:

[0054]

[0055] Where: N is the vector basis function, V is a single solution element, Ω is the entire solution space, E p represents the background electric field, E s represents the secondary electric field, ε c0 represents the dielectric constant of the background field, ε represents the dielectric constant of the secondary field, ω is the angular frequency, and μ is the magnetic permeability.

[0056] Preferably, in step 4, according to the radial distribution profile of formation density and the radial distribution profile of migration length converted in step 3, the neutron and density logging responses are obtained by numerical simulation in combination with instrument structural parameters and nuclear logging fast simulation method.

[0057] The basic principle is that the rate of change of neutrons and gamma photons over time in a certain volume is equal to the generation rate minus the leakage rate and absorption rate. The following formula can be used to describe this process. The flux at the detector can be expressed as the sum of the contributions of each particle emitted from the source to the detector.

[0058] For density logging, the gamma source emits gamma photons, which undergo Compton scattering and photoelectric scattering with atomic nuclei and produce scattered gamma rays. The detector mainly detects the count rate of scattered gamma rays. For neutron porosity logging, the pulsed neutron source emits fast neutrons, which undergo inelastic scattering and elastic scattering reactions with atomic nuclei. The neutron energy is reduced and becomes thermal neutrons. The detector mainly detects the flux of thermal neutrons. According to the above theory, the background area flux can be expressed as:

[0059] N B (r R )=∫dr∫dE∫dΩψ B (r s ,r,E,Ω)S(r R ,r,E,Ω)

[0060] Where N B (r R ) represents the flux of the background area detector, n; r s Represents the position of the source in space, dimensionless; r Rrepresents the position of the detector in space, dimensionless; r represents the distance between a point in space and the source, cm; Ω represents the collision angle of the particle, °; E represents the energy of the particle, MeV; ψ B Indicates that it is located at r s The flux of particles emitted by the radiation source at r, n / (cm 2 ·s). S represents r R The detector response function at can be expressed as:

[0061] S(r R ,r,E,Ω)=∫dE′∫dΩ′∑ s (r,E,Ω→E′Ω′)Ψ + (r R ,r,E′,Ω′)

[0062] Where ∑ S represents the macroscopic scattering cross section of the particle, 1 / cm; Ω′ represents the angle after the particle collision, °; ​​E′ represents the energy after the particle collision, MeV. Ψ + is the importance equation. The background flux sensitivity equation is defined as:

[0063]

[0064] The contribution of the changing area to the detector flux can be expressed as:

[0065]

[0066] Where Δ∑ is the particle cross section in the changing area, 1 / cm; ∑ B is the particle cross section in the background area, 1 / cm. The detector flux is the superposition of the background flux and the flux in the changing area, expressed as:

[0067]

[0068] For density logging, the variables are set to the density values ​​as shown below, C p is the density coefficient, dimensionless; ρ is the density, g / cm 3 ρ B Indicates background density, g / cm 3 .

[0069]

[0070] For neutron logging, the variable is set to the migration length derivative as shown below, C N is the neutron coefficient; Lm is the migration length, dimensionless; LmB is the background migration length, dimensionless.

[0071]

[0072] Preferably, in step 5, the actual logging curve of the underground reservoir is compared with the numerical simulation curve, and through nonlinear iterative inversion and model updating, the forward simulation process of step 3 and step 4 is continuously called to make the numerical simulation curve consistent with the measured logging curve, and finally the reservoir parameters including porosity, permeability and saturation under formation conditions are obtained. The nonlinear iterative inversion is a set of model parameters that is obtained by finding the minimum value of the target equation, and the target equation is:

[0073]

[0074] The objective equation is constrained by the following formula

[0075]

[0076] Where: W d is the data weight; d i is a well logging response obtained by simulation; d m is the measured logging response; α is the damping factor; n c is the number of pre-set mineral components; x is the rock physics model parameter vector; it can be expressed as:

[0077]

[0078] S w is the formation water saturation, represents the formation porosity, C i Content of the i-th mineral component, C sh represents the shale porosity, and the simulated logging response d can be expressed as:

[0079] d=[CNC,DEN,1 / R] T

[0080] CNC, DEN and R are the values ​​of neutron logging, density logging and resistivity logging respectively.

[0081] Compared with the prior art, the beneficial effects of the present invention are as follows: this method is based on rock physical parameters combined with wellbore information and phase permeability characteristic curves, integrates seepage characteristics, electromagnetic fields and nuclear logging principles based on instrument principles, and inverts actual reservoir parameters through multi-time, multi-physical field and multi-dimensional logging response information, thereby improving the accuracy of reservoir fluid property identification and parameter evaluation, and providing technical support for oil and gas exploration and development. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1 It is a flow chart of a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to the present invention.

[0083] Figure 2It is a schematic diagram of the working process of a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to the present invention.

[0084] Figure 3 It is a dynamic cross-sectional diagram of mud invasion according to the present invention.

[0085] Figure 4 It is the resistivity and nuclear logging response diagram of the present invention.

[0086] Figure 5 It is a schematic diagram of the processing results of the actual well based on numerical simulation of the present invention. DETAILED DESCRIPTION

[0087] The drawings are only for illustrative purposes and cannot be construed as limiting the present invention. To better illustrate the present embodiment, some parts in the drawings may be omitted, enlarged, or reduced, and do not represent the size of the actual product. For those skilled in the art, it is understandable that some well-known structures and their descriptions in the drawings may be omitted. The positional relationships described in the drawings are only for illustrative purposes and cannot be construed as limiting the present invention.

[0088] The same or similar numbers in the drawings of the embodiments of the present invention correspond to the same or similar parts; in the description of the present invention, it should be understood that if the terms "upper", "lower", "left", "right", "long", "short" and the like indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings, they are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, the terms describing the positional relationship in the drawings are only used for illustrative purposes and cannot be understood as limitations on this patent. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances.

[0089] The technical solution of the present invention is further described in detail below through specific embodiments and in conjunction with the accompanying drawings:

[0090] Example 1

[0091] like Figure 1-Figure 5 The figure shows a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion, Example 1.

[0092] Step 1: construct a wellbore numerical model; layer the underground reservoir measured logging curve model according to the reservoir gamma, resistivity, neutron, and density curve characteristics, and then square the underground reservoir measured logging curve by layer, each layer represents a rock physics model;

[0093] Step 2: Setting the initial value of the wellbore numerical model; inputting the rock physics model parameters, seepage information and wellbore information parameters into the rock physics model in step 1, wherein the rock physics model parameters and seepage information parameters are assigned by layer, and each layer inputs the parameters into the rock physics model in step 1 according to the formation reservoir characteristics, and the wellbore information parameters are processed according to the whole well section, that is, a unified value is input;

[0094] Step 3: Dynamic numerical simulation: According to the relevant porous media seepage theory and diffusion theory, the initial numerical model of the wellbore obtained in step 2 is subjected to a dynamic invasion simulation of drilling fluid, and the dynamic fluid distribution profile of the reservoir at any time is obtained, including the dynamic saturation and salinity profile of the reservoir. The obtained dynamic fluid distribution profile of the reservoir is converted by a relevant formula to obtain the corresponding radial distribution profiles of resistivity, migration length, and density at different times;

[0095] Step 4: Numerical simulation of resistivity, neutron and density measurement responses; obtain the numerically simulated resistivity logging response based on the electromagnetic field finite element simulation method combined with the instrument structure parameters and working principle, and obtain the numerically simulated neutron and density measurement curves based on the instrument structure parameters and nuclear logging fast simulation method;

[0096] Step 5: Invert reservoir parameters; compare the measured logging curve of the underground reservoir in step 1 with the numerical simulation curve obtained in step 4, and make the numerical simulation curve consistent with the measured logging curve through model updating and iteration, and finally obtain the reservoir parameters of the inverted rock physical structure.

[0097] Beneficial effects of this embodiment: The processing method of this solution can carry out resistivity and nuclear logging joint dynamic simulation and inversion to obtain accurate fluid properties and pore structure parameters of the reservoir.

[0098] Example 2

[0099] Embodiment 2 of a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion, based on embodiment 1, as Figure 1-Figure 5 As shown, steps one to four are further limited.

[0100] Specifically, in step 1, the well logging curve is a record of the change of the physical properties of the formation with the depth of the well. The top and bottom interfaces of each sublayer are divided by the natural gamma half-amplitude point and the reference resistivity, neutron and density curves. Squaring means averaging the numerical segments of the well logging curve in the graphical display mode to form a sawtooth shape, and the average value of the well logging curve of each sublayer is used to represent the sublayer.

[0101] Specifically, in step 2, the input model parameters include rock physics model parameters, seepage information and wellbore information parameters. The rock physics model parameters and seepage parameters are processed by layer, and the wellbore information parameters are processed according to the whole well section. The specific parameters of the rock physics model include rock skeleton parameters, formation physical property parameters, formation pore fluid parameters, and mud content. Rock skeleton parameters include mineral components, content and density; formation physical property parameters include saturation, permeability, porosity and pore structure; formation pore fluid parameters include fluid components, content, density, mineralization and viscosity; seepage information includes phase permeability curve and capillary pressure curve; wellbore information parameters include wellbore diameter, mud information, drill string pressure and temperature. Rock skeleton parameters are obtained through X-ray diffraction data, pore fluid information including fluid composition, density, mineralization, and viscosity is obtained through water analysis experimental materials, seepage information is obtained through regional displacement experiments, wellbore information is obtained through drilling logs and geological daily reports, formation pore structure parameters are obtained through resistivity core experiments, and shale content, porosity, permeability, and saturation are obtained through the average value of the logging curve of each small layer using the following logging interpretation formula. The formula for obtaining the initial value of shale content is:

[0102]

[0103]

[0104] In the formula, I sh Indicates the mud content index, V sh Indicates the formation mud content; decimal; GR indicates the formation natural gamma measurement value, gAPI; GR max Indicates the maximum value of natural gamma in the pure mudstone formation of this well, gAPI; GR min It indicates the minimum natural gamma value of the pure sandstone formation in this well, gAPI.

[0105] The formula for obtaining the initial porosity value is:

[0106]

[0107]

[0108]

[0109] Where: Indicates the calculated density porosity, v / v; DEN indicates the density curve logging reading, g / cm3; DEN ma Indicates the density response parameter of the skeleton, g / cm3; DEN f Indicates the density response parameter of the fluid, g / cm3; V sh Indicates the mud content of the formation, v / v; DEN sh Represents the density response parameter of mudstone, g / cm3; Indicates the calculated neutron porosity, v / v; CNC indicates the neutron curve logging reading, %; CNC ma Indicates the neutron response parameter of the skeleton, %; CNC f Indicates the neutron response parameter of the fluid, %; CNC sh represents the neutron response parameter of mudstone, %; Represents the formation porosity, v / v.

[0110] The permeability calculation can be done using the regional core porosity-permeability regression model. For example, the model formula for the study area is:

[0111]

[0112] Where: K represents the formation permeability, mD; Represents the formation porosity, v / v.

[0113] The formula for obtaining the initial value of water saturation is:

[0114]

[0115] Where: S w Indicates the water saturation of the formation; R t Indicates true resistivity of formation, Ω.m; V sh Indicates the mud content, decimal; R sh Represents the resistivity of pure mudstone, Ω.m; R w Indicates the resistivity of formation water, Ω.m; represents porosity, a decimal; a represents lithology coefficient; m represents cementation index; n represents saturation index, a, m, n are obtained by block experience or other logging interpretation methods.

[0116] Specifically, in step three, the phase permeability curve and capillary pressure curve are obtained through core experiments. According to the relevant permeability theory and ion diffusion theory including the phase permeability curve and capillary pressure curve, the dynamic distribution profile of the reservoir fluid including the formation saturation and salinity at any time is simulated, and the corresponding radial distribution profiles of resistivity, migration length and density at different times are obtained through resistivity and nuclear logging related formula conversion. The resistivity radial profile conversion formula is:

[0117] Where: R t Represents formation resistivity; R w Represents the resistivity of formation water; S w Indicates the water saturation of the formation; Represents formation porosity; R sh Represents the resistivity of mud; V shrepresents the shale content; a represents the lithology coefficient; m represents the cementation index; n represents the saturation index; c represents the shale content index coefficient. The general calculation formula is: a, m, and n are obtained using block experience or other well logging interpretation methods.

[0118] The density radial profile conversion formula is:

[0119]

[0120] Where: b represents the formation density; ρ w represents the formation water density; ρ h Indicates the density of oil and gas fluid; represents the formation porosity; ρ sh Indicates the density of mud; C sh Indicates the shale porosity; S w Indicates the formation water saturation; ρ i The density of the i-th component in the skeleton; C i represents the content of the i-th component in the skeleton; n c Represents the total number of component types in the skeleton.

[0121] The calculation formula of the radial profile of neutron migration length is:

[0122]

[0123] Where: ξ represents the formation migration length; S w Indicates the water saturation of the formation; ξ w represents the migration length of formation water; ξ h It represents the migration length of oil and gas fluid; represents the formation porosity; ξ m represents the skeleton migration length; η represents the exponential coefficient.

[0124] Specifically, in step 4, the resistivity logging response is obtained according to the radial distribution profile of resistivity converted in step 3, and the corresponding resistivity logging response is obtained mainly by using the electromagnetic field finite element simulation method combined with the instrument structure parameters and the numerical simulation of the working principle. The three-dimensional vector finite element method is used to spatially discretize the logging physical model, establish the unit electromagnetic field discrete equation, install, eliminate and solve all units, and then obtain the electrical logging instrument response. Resistivity logging simulation is to solve the Maxwell equations under given boundary conditions. Starting from Maxwell's equations, the electromagnetic field satisfies the following equation:

[0125]

[0126]

[0127] Where: E represents the electric field intensity; H represents the magnetic field intensity; J is the source current density; ω is the angular frequency; σ is the electrical conductivity; μ is the magnetic permeability.

[0128] From the above two equations, we can deduce that the vector wave equation satisfied by the electric field is:

[0129]

[0130] is the complex dielectric constant, ε=ε r ε0 where ε0 is the dielectric constant of vacuum, ε r is the relative dielectric constant.

[0131] make

[0132] E=E p +E s

[0133] Where, the background field E p is the electric field when the entire space is filled with a medium with conductivity σ0, which satisfies the equation:

[0134]

[0135] in, By changing the above three equations, we can get the vector wave equation:

[0136]

[0137] The background field is calculated by analytical methods, and the secondary field is calculated by the finite element method. The solution of the above equation changes slowly, and a sparser grid can be used to solve it, reducing the computational workload. Select a large enough area so that the electric field on the boundary decays to approximately 0, then the above equation only needs to satisfy the boundary condition equation:

[0138]

[0139] Where: To solve the boundary of the region ω, n is the direction of its normal.

[0140] Considering the boundary condition equation, the vector wave equation is transformed into its weak product form:

[0141]

[0142] Where: N is the vector basis function, V is a single solution element, Ω is the entire solution space, E p represents the background electric field, E s represents the secondary electric field, ε c0 represents the dielectric constant of the background field, ε represents the dielectric constant of the secondary field, ω is the angular frequency, and μ is the magnetic permeability.

[0143] Specifically, in step 4, according to the radial distribution profile of formation density and the radial distribution profile of migration length converted in step 3, the neutron and density logging responses are obtained by numerical simulation in combination with instrument structural parameters and nuclear logging fast simulation method.

[0144] The basic principle is that the rate of change of neutrons and gamma photons over time in a certain volume is equal to the generation rate minus the leakage rate and absorption rate. The following formula can be used to describe this process. The flux at the detector can be expressed as the sum of the contributions of each particle emitted from the source to the detector.

[0145] For density logging, the gamma source emits gamma photons, which undergo Compton scattering and photoelectric scattering with atomic nuclei and produce scattered gamma rays. The detector mainly detects the count rate of scattered gamma rays. For neutron porosity logging, the pulsed neutron source emits fast neutrons, which undergo inelastic scattering and elastic scattering reactions with atomic nuclei. The neutron energy is reduced and becomes thermal neutrons. The detector mainly detects the flux of thermal neutrons. According to the above theory, the background area flux can be expressed as:

[0146] N B (r R )=∫dr∫dE∫dΩψ B (r s ,r,E,Ω)S(r R ,r,E,Ω)

[0147] Where N B (r R ) represents the flux of the background area detector, n; r s Represents the position of the source in space, dimensionless; r R represents the position of the detector in space, dimensionless; r represents the distance between a point in space and the source, cm; Ω represents the collision angle of the particle, °; E represents the energy of the particle, MeV; ψ B Indicates that it is located at r s The flux of particles emitted by the radiation source at r, n / (cm 2 ·s). S represents r R The detector response function at can be expressed as:

[0148] S(r R ,r,E,Ω)=∫dE′∫dΩ′∑ s (r,E,Ω→E′Ω′)Ψ + (r R ,r,E′,Ω′)

[0149] Where ∑ S represents the macroscopic scattering cross section of the particle, 1 / cm; Ω′ represents the angle after the particle collision, °; ​​E′ represents the energy after the particle collision, MeV. Ψ+ is the importance equation. The background flux sensitivity equation is defined as:

[0150]

[0151] The contribution of the changing area to the detector flux can be expressed as:

[0152]

[0153] Where Δ∑ is the particle cross section in the changing area, 1 / cm; ∑ B is the particle cross section in the background area, 1 / cm. The detector flux is the superposition of the background flux and the flux in the changing area, expressed as:

[0154]

[0155] For density logging, the variables are set to the density values ​​as shown below, C p is the density coefficient, dimensionless; ρ is the density, g / cm 3 ρ B Indicates background density, g / cm 3 .

[0156]

[0157] For neutron logging, the variable is set to the migration length derivative as shown below, C N is the neutron coefficient; Lm is the migration length, dimensionless; LmB is the background migration length, dimensionless.

[0158]

[0159] Beneficial effects of this embodiment: In step 1, the well logging curve is squared and the numerical segmentation of the well logging curve is averaged and layered, so that the well logging curve of each small layer can be more representative. In step 4, the three-dimensional vector finite element method is used to spatially discretize the well logging physical model, establish the unit electromagnetic field discrete equation, install, eliminate and solve all units, and accurately obtain the response of the electric logging instrument.

[0160] Example 3

[0161] Embodiment 3 of a method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion, based on embodiment 1 or 2, as Figure 1-Figure 5 As shown, step five is further limited.

[0162] Specifically, in step 5, the underground reservoir logging curve is compared with the numerical simulation curve. Through nonlinear iterative inversion and model update, the forward simulation process of step 3 and step 4 is continuously called to make the numerical simulation curve consistent with the measured logging curve, and finally the reservoir parameters including porosity, permeability and saturation under formation conditions are obtained. Among them, nonlinear iterative inversion is to find a set of model parameters that minimize the target equation. The target equation is:

[0163]

[0164] The objective equation is constrained by the following formula

[0165]

[0166] Where: W d is the data weight; d i is a well logging response obtained by simulation; d m is the measured logging response; α is the damping factor; n c is the number of pre-set mineral components; x is the rock physics model parameter vector; it can be expressed as:

[0167]

[0168] S w is the formation water saturation, represents the formation porosity, C i Content of the i-th mineral component, C sh represents the shale porosity, and the simulated logging response d can be expressed as:

[0169] d=[CNC,DEN,1 / R] T

[0170] CNC, DEN and R are the values ​​of neutron logging, density logging and resistivity logging respectively.

[0171] The beneficial effects of this embodiment are as follows: in step five, through nonlinear iterative inversion and model updating, the forward simulation process of steps three and four can be quickly called to make the numerical simulation curve consistent with the measured logging curve, and finally obtain the reservoir parameters including porosity, permeability and saturation under formation conditions.

[0172] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is not necessary and impossible to list all the embodiments here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. A method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion, characterized in that: The steps include: Step 1: construct a wellbore numerical model; layer the underground reservoir measured logging curve model according to the reservoir gamma, resistivity, neutron, and density curve characteristics, and then square the underground reservoir measured logging curve by layer, each layer represents a rock physics model; Step 2: Setting the initial value of the wellbore numerical model; inputting the rock physics model parameters, seepage information and wellbore information parameters into the rock physics model in step 1, wherein the rock physics model parameters and seepage information parameters are assigned by layer, and each layer inputs the parameters into the rock physics model in step 1 according to the formation reservoir characteristics, and the wellbore information parameters are processed according to the whole well section, that is, a unified value is input; Step 3: Dynamic numerical simulation: According to the relevant porous media seepage theory and diffusion theory, the drilling fluid dynamic invasion simulation is performed on the wellbore numerical model obtained in the step 2 to obtain the reservoir dynamic fluid distribution profile at any time, including the reservoir dynamic saturation and salinity profile, and the obtained reservoir dynamic fluid distribution profile is converted by relevant formulas to obtain the corresponding radial distribution profiles of resistivity, migration length, and density at different times; Step 4: Numerical simulation of resistivity, neutron and density measurement responses; obtain the numerically simulated resistivity logging response based on the electromagnetic field finite element simulation method combined with the instrument structure parameters and working principle, and obtain the numerically simulated neutron and density measurement curves based on the instrument structure parameters and nuclear logging fast simulation method; Step 5: Invert reservoir parameters; compare the measured logging curve of the underground reservoir in step 1 with the numerical simulation curve obtained in step 4, and make the numerical simulation curve consistent with the measured logging curve through model updating and iteration, and finally obtain the reservoir parameters of the inverted rock physical structure.

2. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 1, characterized in that: In step 2, the rock physics model parameters include rock skeleton parameters, formation physical property parameters, formation pore fluid parameters, and mud content; the seepage information includes phase permeability curves and capillary pressure curves; and the wellbore information parameters include wellbore diameter, mud information, drill string pressure, and temperature.

3. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 2, characterized in that: The rock skeleton parameters include mineral composition, content and density; the formation physical property parameters include saturation, permeability, porosity and pore structure; the formation pore fluid parameters include fluid composition, content, density, mineralization and viscosity; the mud information includes drilling fluid mineralization, mud solid content and drilling fluid viscosity.

4. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 3, characterized in that: In the step 2, the initial value is set, the rock skeleton parameters are obtained through X-diffraction data, including fluid composition, density, mineralization, viscosity, pore fluid information is obtained through water analysis experimental materials, the seepage information is obtained through core displacement experiments, the wellbore information is obtained through drilling logs and geological daily reports, the formation pore structure parameters are obtained through resistivity core experiments, the mud content, the porosity, the permeability, and the saturation are obtained through the average value of the logging curve of each small layer using conventional logging interpretation formulas.

5. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 1, characterized in that: In step three, the resistivity radial distribution profile is obtained by converting the saturation and mineralization dynamic profiles using the Indonesian formula, and the density radial distribution profile is obtained by calculating the rock physics volume model based on the saturation dynamic profile, porosity, pore fluid components and rock skeleton component density values.

6. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 1, characterized in that: In step three, the migration length radial distribution profile is calculated through a rock physics volume model based on the saturation dynamic profile, porosity, pore fluid and rock skeleton component migration length values.

7. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 1, characterized in that: In step 4, the instrument is a resistivity logging instrument, a neutron logging instrument, or a density logging instrument.

8. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 1, characterized in that: In the step 4, based on the radial distribution profile of the formation resistivity converted in the step 3, the three-dimensional vector finite element method is used to spatially discretize the logging physical model, establish the unit electromagnetic field discrete equation, install, eliminate and solve all units, and then obtain the resistivity logging response.

9. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 1, characterized in that: In the step 4, the density and neutron response curves are obtained based on the radial distribution profile of the formation density and the radial distribution profile of the migration length converted in the step 3, combined with the instrument structure parameters and the nuclear logging fast simulation method.

10. The method for joint inversion of reservoir parameters by resistivity and nuclear logging under mud invasion according to claim 1, characterized in that: In step five, the underground reservoir logging curve is compared with the numerical simulation curve. Through nonlinear iterative inversion and model updating, the forward simulation process of step two and step four is continuously called to make the numerical simulation curve consistent with the measured logging curve, and finally the reservoir parameters under formation conditions are obtained.

Images