Joint inversion method for multi-field coupling parameters of submarine muddy silt hydrate reservoir and based on multi-source monitoring data

Abstract

The present invention belongs to the technical field of natural gas hydrate development and monitoring. Disclosed is a joint inversion method for multi-field coupling parameters of a submarine muddy silt hydrate reservoir and based on multi-source monitoring data. The method comprises: acquiring logging data and laboratory test data, and determining the mud content and pore water salinity of a reservoir therefrom; determining a natural gas phase state of the reservoir by means of the data and a natural gas hydrate phase equilibrium model; constructing a resistivity model and an acoustic wave velocity model which take the influence of the mud content into consideration, adding a physical constraint, using a differential optimization algorithm to invert the porosity and hydrate saturation of the reservoir, and then verifying the reliability of the models by means of forward modeling; and on the basis of an inversion result, using a thermodynamic model, a permeability model and a mechanical strength model, which take the mud content into consideration, to perform forward modeling, so as to acquire multi-field coupling parameters of the reservoir. The method is mainly used for parameter inversion of submarine muddy silt hydrate reservoirs, and can improve the parameter inversion accuracy.

Description

Multi-field coupling parameter joint inversion method for hydrate reservoir of submarine argillaceous silt type based on multi-source monitoring data TECHNICAL FIELD

[0001] The application belongs to the technical field of natural gas hydrate development monitoring, and particularly relates to the technical field of multi-field coupling parameter joint inversion of hydrate reservoir of submarine argillaceous silt type. BACKGROUND

[0002] Natural gas hydrate is an ice-like crystalline compound formed under low temperature and high pressure conditions, and is mostly present in submarine shallow argillaceous silt or sandy loose sediment layers. It is considered as an important clean alternative energy source due to its large reserves, high energy density, and non-polluting combustion. In the process of exploitation, the hydrate in the reservoir pore will undergo phase transition, resulting in complex changes in the reservoir pore structure, permeability, mechanical properties and heat transfer characteristics, and these changes interact with each other, which has an important influence on the safety and efficiency of exploitation. Therefore, it is necessary to effectively monitor the multi-field coupling process of heat flow and force of the reservoir.

[0003] At present, the main ways to obtain the characteristic data of natural gas hydrate reservoir and monitor the multi-field coupling parameters include exploration, well logging, sampling analysis, etc. The methods used cover multiple fields such as geophysics and geochemistry. Although these single monitoring technologies can characterize part of the characteristics and dynamic changes of the reservoir from different angles, the information provided by each technology is limited, and it is difficult to fully reflect the complex situation of multi-field coupling of the reservoir. At the same time, there are differences in the space-time conditions of different test processes, so the reservoir characteristic parameters inverted from different sources of monitoring data often have large deviations, which affects the judgment of the true state of the reservoir.

[0004] Especially crucial is that the submarine hydrate reservoir is mostly of argillaceous silt type, and the argillaceous content will have a significant impact on the electrical, acoustic, thermal, mechanical and other characteristics of the reservoir. Most of the existing reservoir parameter inversion models do not fully consider the influence of argillaceous content, or the correction of its influence is not perfect, resulting in insufficient accuracy of the inversion results, and it is difficult to accurately obtain key parameters such as reservoir porosity, hydrate saturation, permeability, and mechanical strength. How to effectively integrate multi-source monitoring data, build a multi-field coupling inversion model that fully considers the influence of argillaceous content, and comprehensively and systematically invert key parameters of the reservoir to improve the inversion accuracy has become a technical problem to be solved in the field of dynamic monitoring of natural gas hydrate exploitation reservoir. SUMMARY

[0005] The purpose of the present application is to provide a multi-field coupling parameter joint inversion method for hydrate reservoir of submarine argillaceous silt type based on multi-source monitoring data, which solves the problem of low accuracy of reservoir parameter inversion caused by insufficient consideration of the influence of argillaceous content and insufficient integration of multi-source monitoring data in the prior art, and realizes integrated inversion of multi-field coupling parameters.

[0006] To achieve the above object, the present application provides the following technical solutions.

[0007] A multi-field coupling parameter joint inversion method for a seabed argillaceous silt type hydrate reservoir, comprising the steps of:

[0008] Obtaining logging data and laboratory test data;

[0009] Determining the argillaceous content and the pore water salinity of the reservoir from the logging data and the laboratory test data;

[0010] Judging the natural gas phase state of the reservoir according to a natural gas hydrate phase equilibrium model in combination with the argillaceous content and the pore water salinity;

[0011] Constructing a resistivity model and a sonic velocity model considering the influence of the argillaceous content, constructing an objective function based on the resistivity data and the sonic velocity data in the logging data, and performing global optimization inversion by using a differential optimization algorithm after adding physical constraints to obtain the reservoir porosity and the hydrate saturation;

[0012] Performing forward calculation on the resistivity model and the sonic velocity model by using the reservoir porosity and the hydrate saturation to complete model reliability verification;

[0013] Constructing a thermodynamic model, a permeability model and a mechanical strength model considering the influence of the argillaceous content, and performing forward calculation by using the reservoir porosity, the hydrate saturation and the argillaceous content to obtain reservoir parameters.

[0014] In a possible implementation manner, the reservoir parameters include reservoir thermal conductivity, thermal capacity, absolute permeability, gas-water relative permeability, dynamic-static elastic parameters, compressive strength and internal friction angle;

[0015] When the logging data is obtained, the temperature, the pore pressure, the resistivity, the sonic velocity and the natural gamma data of the reservoir are collected;

[0016] When the laboratory test data is obtained, the mineral composition, the mechanical properties and the thermal conduction properties related data of the reservoir rock sample are collected.

[0017] In a possible implementation manner, when the argillaceous content is determined, the natural gamma data in the logging data is used to obtain the reservoir argillaceous content through natural gamma normalization calculation and a related empirical formula;

[0018] When the pore water salinity is determined, it is directly obtained through the laboratory test data or indirectly calculated by using related response parameters in the logging data.

[0019] In one possible implementation, when determining the phase state of the natural gas in the reservoir, a hydrate thermodynamic model based on the fugacity equilibrium theory is used, which is modified by combining the equation of state, the salinity enhancement effect of the clay reservoir, and the pore size effect to determine the stable region of the reservoir hydrate, and then the phase state of the natural gas is determined.

[0020] In one possible implementation, when constructing the resistivity model, the Waxman-Smits model is improved by incorporating the influence of clay content on reservoir conductivity, and the relationship between reservoir resistivity and water saturation and hydrate saturation is correlated through the model.

[0021] When constructing the acoustic velocity model, the Thomas-Stieber rock physics hybrid model was improved, and the acoustic properties of sandstone, mudstone, hydrate phase and fluid phase were combined to establish the correlation between mud content and longitudinal wave velocity and transverse wave velocity.

[0022] In one possible implementation, during the global optimization inversion, the deviation between the calculated values ​​of the resistivity model and the sonic velocity model and the measured values ​​from the well logging is used as the objective function. Physical range constraints on the reservoir parameters are added, and the optimal solution is iteratively searched using a differential optimization algorithm to obtain the reservoir porosity and hydrate saturation.

[0023] In one possible implementation, when constructing the thermodynamic model, the calculation relationship between reservoir thermal conductivity and heat capacity is established by incorporating parameters such as clay content, porosity, and hydrate saturation based on the volume-weighted mixing law.

[0024] When performing forward modeling calculations using the aforementioned thermodynamic model, calibration is performed by combining laboratory measurement data, well logging statistics, and dynamic correction parameters to improve the accuracy of thermal conductivity and heat capacity calculations.

[0025] In one possible implementation, when constructing the permeability model, the absolute permeability of the reservoir is calculated by using the improved Kozeny-Carman equation and incorporating a clay content correction term; the relationship between water saturation, hydrate saturation and gas-water relative permeability is correlated by combining the Brooks-Corey model.

[0026] When performing forward modeling calculations using the aforementioned permeability model, the parameters of reservoir fluid flow characteristics are determined.

[0027] In one possible implementation, when constructing the mechanical strength model, the correlation between reservoir compressive strength, internal friction angle and the above parameters is established by combining porosity, hydrate saturation and clay content through an empirical strength model.

[0028] In the process of constructing the dynamic-static elastic parameter calculation relationship, the dynamic elastic parameters are calculated through the longitudinal wave velocity and the transverse wave velocity, and the static elastic parameters are obtained by using the dynamic-static elastic parameter conversion model.

[0029] In a possible implementation, in the process of calculating the dynamic elastic parameters, the dynamic Young's modulus and the dynamic Poisson's ratio are calculated by using the related formula of elastic mechanics based on the longitudinal wave velocity and the transverse wave velocity obtained by the acoustic wave velocity model and in combination with the reservoir saturation density.

[0030] In the process of dynamic-static elastic parameter conversion, the empirical conversion model containing a conversion coefficient and a conversion intercept is used to obtain the static Young's modulus and the static Poisson's ratio.

[0031] Compared with the prior art, the present application has the beneficial effects that, compared with the prior single monitoring technology with limited information and insufficient data integration, the present application simultaneously obtains logging data and laboratory test data, covers multiple types of information such as reservoir temperature, pore pressure, resistivity, acoustic wave velocity, and rock sample mineral composition, effectively solves the parameter inversion deviation problem caused by the dispersion of data from different sources, and makes the inversion result more consistent with the real state of the reservoir.

[0032] In view of the limitation that the existing inversion model does not fully consider the influence of shale content, the present application integrates the shale content factor in all core models. The improved Waxman-Smits resistivity model and the Thomas-Stieber acoustic wave velocity model can accurately associate the shale content with the electrical and acoustic properties of the reservoir, and improve the inversion accuracy of porosity and hydrate saturation; the corrected Kozeny-Carman permeability model, the thermodynamic model based on the volume-weighted mixing law, and the empirical mechanical strength model respectively adapt to the characteristics of the shale silt type reservoir from the dimensions of fluid flow, heat conduction, and mechanical stability, avoid the parameter calculation deviation caused by ignoring the influence of shale, and make each reservoir parameter more practically valuable.

[0033] In the inversion calculation process, the present application constructs an objective function and adds physical constraints, and performs global optimization by combining a differential optimization algorithm, compared with the existing inversion method lacking constraints or local optimization, can effectively avoid local optimal solution, and ensure the rationality of the inversion of core parameters. At the same time, the newly added forward verification step uses the porosity and hydrate saturation obtained by inversion to perform forward calculation on the resistivity model and the acoustic wave velocity model, corrects the model deviation through the result verification, and further improves the reliability of the inversion result.

[0034] Existing technologies can only obtain partial reservoir parameters, while this invention achieves integrated inversion of multi-field coupled parameters. Through the inversion of basic physical property parameters and subsequent forward modeling calculations of thermodynamics, permeability, and mechanical strength models, key parameters such as reservoir thermal conductivity, heat capacity, absolute permeability, dynamic-static elastic parameters, and compressive strength are comprehensively obtained, fully covering the thermo-fluidomechanical multi-field coupled characteristics. This allows for more refined and comprehensive reservoir characterization, solving the problem of insufficient targeted development planning caused by incomplete parameter supply in existing technologies.

[0035] In practical applications, through the integration of multi-source data and correction of mud content, the consistency between the inversion results of various parameters and the measured data has been significantly improved. As can be seen from the comparison between the observed values ​​and the predicted values ​​in the attached figure, the inversion results of key parameters such as resistivity and longitudinal wave velocity can accurately match the actual monitoring situation. Attached Figure Description

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

[0037] Figure 1 is a schematic diagram of the technical process of the joint inversion method of multi-field coupled parameters of seabed muddy silt-type hydrate reservoir based on multi-source monitoring data according to an embodiment of the present invention.

[0038] Figure 2 is a comparative verification analysis of the resistivity inversion calculation of the SH-XX well in an embodiment of the present invention;

[0039] Figure 3 is a comparative verification analysis diagram of the longitudinal wave velocity inversion calculation of the SH-XX well according to an embodiment of the present invention;

[0040] Figure 4 is a diagram showing the inversion calculation and analysis of shear wave velocity in the SH-XX well according to an embodiment of the present invention;

[0041] Figure 5 is a calculation and analysis diagram of the thermal conductivity inversion of the SH-XX well according to an embodiment of the present invention;

[0042] Figure 6 is a thermal capacity inversion calculation and analysis diagram of the SH-XX well according to an embodiment of the present invention;

[0043] Figure 7 is a diagram showing the absolute permeability inversion calculation and analysis of the SH-XX well according to an embodiment of the present invention;

[0044] Figure 8 is a diagram showing the inversion calculation and analysis of the elastic modulus of the SH-XX well according to an embodiment of the present invention.

[0045] Figure 9 is a diagram showing the Poisson's ratio inversion calculation and analysis of the SH-XX well according to an embodiment of the present invention;

[0046] Fig. 10 is a SH-XX well compressive strength inversion calculation analysis diagram of the embodiment of the present application. DETAILED DESCRIPTION

[0047] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all the other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present application.

[0048] EMBODIMENT

[0049] It should be noted that the terms "comprising" and "having" and any variations thereof in the embodiments of the present application are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device including a series of steps or units does not have to be limited to only those steps or units clearly listed, but can include other steps or units not clearly listed or inherent to the process, method, product or device.

[0050] The multi-field coupling parameter joint inversion method for a submarine muddy silt type hydrate reservoir in the embodiments of the present application comprises the following steps:

[0051] Step 101: Obtain logging data and laboratory test data.

[0052] Specifically, the logging data can be temperature, pore pressure, resistivity, acoustic velocity and natural gamma data of the reservoir; and the laboratory test data can be mineral composition, mechanical property and thermal conductivity data of the reservoir rock sample.

[0053] The reservoir parameters include reservoir thermal conductivity, heat capacity, absolute permeability, gas-water relative permeability, dynamic-static elastic parameters, compressive strength and internal friction angle; when the logging data is obtained, the temperature, pore pressure, resistivity, acoustic velocity and natural gamma data of the reservoir are collected; and when the laboratory test data is obtained, the mineral composition, mechanical property and thermal conductivity of the reservoir rock sample are collected.

[0054] Specifically, the reservoir thermal conductivity can be a parameter characterizing the heat transfer capacity of the reservoir; the heat capacity can be the heat required for the temperature of the reservoir per unit volume to rise by 1K; the absolute permeability can be the inherent ability of the reservoir to allow fluid to pass through; the gas-water relative permeability can be the ratio of the gas phase or water phase to the absolute permeability; the dynamic-static elastic parameters can be the dynamic Young's modulus, the static Young's modulus, the dynamic Poisson's ratio, the static Poisson's ratio; the compressive strength can be the ability of the reservoir to resist pressure damage; the internal friction angle can be a parameter characterizing the shear resistance of the reservoir; the mineral composition data can be the proportion of quartz, clay minerals in the rock sample; the mechanical property data can be the compressive strength, the elastic modulus of the rock sample; the thermal conduction characteristic data can be the thermal conductivity of the rock sample.

[0055] When acquiring the logging data, the temperature data, the pore pressure data, the resistivity data, the longitudinal wave velocity data, the transverse wave velocity data, and the natural gamma data of the target reservoir at different depths are collected by measuring downhole with the logging instrument; when acquiring the laboratory test data, the rock sample obtained from the target reservoir is analyzed, the mineral composition is tested by the X-ray diffractometer to obtain the percentage data of quartz, clay, and other minerals, the mechanical properties are tested by the pressure testing machine to obtain the compressive strength and the elastic modulus data of the rock sample, and the thermal conduction characteristics are tested by the thermal conductivity instrument to obtain the thermal conductivity and the heat capacity data of the rock sample. All the collected data are classified and stored in the corresponding database.

[0056] Step 102: determining the shale content and the pore water salinity of the reservoir from the logging data and the laboratory test data.

[0057] Specifically, the shale content can be the volume proportion of shale components in the reservoir; the pore water salinity can be the weight percentage of salt in the pore water.

[0058] Wherein, when determining the shale content, the natural gamma data in the logging data is used to obtain the shale content of the reservoir through natural gamma normalization calculation and a relevant empirical formula; when determining the pore water salinity, the laboratory test data is directly acquired or indirectly calculated by using the relevant response parameters in the logging data.

[0059] Specifically, the natural gamma normalization calculation can be to convert the measured natural gamma value to a normalized value between 0 and 1, for example, I GR =(GR-GRminmax min ; the empirical formula can be a shale content calculation formula based on the HILCHI coefficient V sh =(2 CGUR×IGR -1) / (2 CGUR -1); the relevant response parameters can be resistivity, acoustic wave velocity, and other logging parameters related to the pore water salinity.

[0060] When determining the pore water salinity, if there is laboratory test data, the test result of the pore water of the rock sample is directly used; if there is no laboratory data, the pore water salinity is indirectly calculated by using the resistivity response parameter in the logging data, combining the Arps salinity-temperature formula, and the correlation between the resistivity and the salinity.

[0061] Step 103: judging the reservoir natural gas phase state according to the natural gas hydrate phase equilibrium model in combination with the shale content and the pore water salinity.

[0062] Specifically, the natural gas hydrate phase equilibrium model can be a CSMHYD thermodynamic equilibrium model.

[0063] In the process of judging the reservoir natural gas phase state, a hydrate thermodynamic model based on the fugacity equilibrium theory is used, and the state equation, the salinity enhancement effect of the shale reservoir and the pore size effect are combined for correction to determine the hydrate stable region of the reservoir and further judge the natural gas phase state.

[0064] Specifically, the fugacity equilibrium theory can be the equilibrium conditions of the hydrate phase, the gas phase, the water phase and the like; the state equation can be the Peng-Robinson state equation; the salinity enhancement effect can be the phenomenon that the adsorption of ions on the surface of clay minerals leads to the increase of the equivalent salinity; the pore size effect can be the phenomenon that the decrease of the effective pore radius caused by the shale leads to the change of the capillary pressure; and the hydrate stable region can be the region that meets the hydrate formation in the temperature-pressure coordinate system.

[0065] In the process of judging the reservoir natural gas phase state, the CSMHYD hydrate thermodynamic model based on the fugacity equilibrium theory is used, and the equilibrium conditions are f hydrate =f gas =f water , wherein f hydrate , f gas and f water are the fugacities of the hydrate phase, the gas phase and the water phase respectively; the fugacities of the gas phase and the water phase are calculated by combining the Peng-Robinson state equation;

[0066] The calculation formula of the hydrate phase fugacity is as follows:

[0067] wherein v i is the proportion of the i-th gas molecule occupying the hydrate cage; φ i is the fugacity coefficient of the cage structure; y i is the mole fraction of component i in the gas phase; R is the gas constant; T is the temperature (℃); Δμ empty is the chemical potential difference of the empty hydrate lattice, which is a physical parameter of the hydrate phase equilibrium model and determines the relative stability of the empty lattice;

[0068] The calculation formula of the gas phase fugacity is as follows:

[0069] where P is the system pressure; is the fugacity coefficient of component i in the gas phase;

[0070] The water phase fugacity calculation formula is as follows:

[0071] x i : the mole fraction of component i dissolved in water; γ i : the activity coefficient; Psat: the saturated vapor pressure of pure water; v i : the partial molar volume of component i gas;

[0072] Considering the salinity enhancement effect of the argillaceous reservoir, the ion adsorbed on the surface of clay minerals makes the equivalent salinity increase, and the water phase fugacity is modified; considering the pore size effect, the argillaceous content reduces the effective pore radius, and the capillary pressure term is modified; for example, the reservoir temperature, pore pressure, argillaceous content, and pore water salinity are substituted into the modified model to determine the hydrate stable region of the reservoir; the positional relationship between the actual temperature and pressure of the reservoir and the stable region is compared to determine the phase state of the natural gas as solid hydrate.

[0073] Step 104: constructing a resistivity model and a sonic velocity model considering the influence of argillaceous content, constructing a target function based on resistivity data and sonic velocity data in the logging data, adding physical constraints and using a differential optimization algorithm for global optimization inversion to obtain reservoir porosity and hydrate saturation.

[0074] Specifically, the resistivity model can be an improved Waxman-Smits model; the sonic velocity model can be an improved Thomas-Stieber rock physics hybrid model; the target function can be the sum of squares of deviations between model calculation values and logging measured values; and the physical constraint can be a value range of 10% to 40% of porosity.

[0075] In the global optimization inversion, the deviation between the calculation values of the resistivity model and the sonic velocity model and the logging measured values is taken as the target function, the physical value range constraint of the reservoir parameters is added, and the optimal solution is iteratively searched by a differential optimization algorithm to obtain the reservoir porosity and hydrate saturation.

[0076] Specifically, the improved Waxman-Smits model can be a resistivity calculation model incorporating the contribution term of argillaceous conductivity; the influence of argillaceous on the conductivity of the reservoir can be the cation exchange of argillaceous to enhance the conductivity; and the water saturation can be the volume ratio of water in the reservoir pores.

[0077] When constructing the resistivity model, the Waxman-Smits model is improved to incorporate the influence of argillaceous content on the conductivity of the reservoir, and the model expression is where R t is the formation resistivity (Ω·m) ; φ is the porosity; S w is the water saturation, which is the inversion target parameter; a is the cementation factor; m is the porosity exponent; n is the saturation exponent; R w is the formation water resistivity (Ω·m) ; Q v is the cation exchange capacity (meq / mL) ; V sh is the shale content; C sh is the shale conductivity contribution (S / m) ; R sh is the shale resistivity (Ω·m) ; the model is used to correlate the reservoir resistivity with the water saturation and the hydrate saturation.

[0078] When constructing the acoustic velocity model, the Thomas-Stieber petrophysical hybrid model is improved, and the acoustic characteristics of the sandstone skeleton P-wave velocity, the shale skeleton P-wave velocity, the hydrate P-wave velocity and the water P-wave velocity are combined to establish a P-wave velocity expression V p,water =1357.0+4.9T-0.05T 2 +1.3(S salt -35)

[0079] where V p is the P-wave velocity (m / s) ; V s is the P-wave velocity (m / s) ; φ is the porosity; V sh is the shale content; S h is the hydrate saturation; S w is the water saturation; V p,sand is the sandstone skeleton P-wave velocity (m / s) ; V p,hydrate is the hydrate P-wave velocity (m / s) ; V p,water is the water P-wave velocity (m / s) ; V p,shale is the shale skeleton P-wave velocity (m / s) ; T is the temperature (℃) ; S salt is the salinity percentage; the S-wave velocity expression is established by combining the shale content, the porosity, the hydrate saturation and the water saturation, the correlation between the shale content and the P-wave velocity and the S-wave velocity is established by the model, and the S-wave velocity equation adopts the BGTL theoretical equation:

[0080] where k: bulk modulus (unit: Pa) ; μ eff : effective shear modulus (unit: Pa) ; ρ sat : saturated rock density (unit: kg / m 3 ).

[0081] Bulk modulus: k=k ma(1-β)+β 2 M

[0082] where k ma is the bulk modulus of the rock matrix; k fl is the bulk modulus of the fluid in the pore; β is the Biot coefficient, representing the ratio of the fluid volume change to the rock volume change, which is related to the porosity of the sedimentary medium; and M is a modulus, representing the increment of the hydrostatic pressure required to press a certain amount of water into the sedimentary medium under the condition that the volume of the sedimentary medium is constant.

[0083] Effective shear modulus:

[0084] where μ ma is the shear modulus of the rock matrix (GPa);

[0085] k ma and μ ma are calculated by the Hill average equation:

[0086] where m is the number of minerals in the solid phase of the rock; f i is the volume fraction of the i-th mineral in the solid phase; k i and μ i are the bulk modulus and shear modulus of the i-th mineral, respectively.

[0087] Saturation density: ρ sat = (1-φ-V sh )·ρ sand +φ·(S h ·ρ hydrate +S w ·ρ water )+V sh ·ρ shale

[0088] where ρ sand is the density of the sandstone (kg / m 3 ); ρ hydrate is the density of the hydrate (kg / m 3 ); and ρ water is the density of water (kg / m 3 ).

[0089] Step 105: Perform forward calculation on the resistivity model and acoustic velocity model using the reservoir porosity and hydrate saturation to complete the model reliability verification.

[0090] Step 106: build a thermodynamic model, a permeability model and a mechanical strength model considering the influence of shale content, perform forward calculation by using the reservoir porosity, hydrate saturation and shale content to obtain reservoir parameters.

[0091] Wherein, when building the thermodynamic model, based on the volume-weighted mixing law, the shale content, porosity and hydrate saturation parameters are integrated to establish the calculation relationship of reservoir thermal conductivity and heat capacity; when performing forward calculation by the thermodynamic model, the calculation accuracy of thermal conductivity and heat capacity is improved by combining laboratory measurement data, logging statistical results and dynamic correction parameters for calibration. Specifically, the volume-weighted mixing law can be a law of calculating the overall characteristics by volume ratio of each component; the thermal conductivity calculation relationship can be λ = (1-φ-V sh )·λ sand +φ·(S h ·λ hydrate +S w ·λ water )+V sh ·λ shale

[0092] Wherein, λ is the reservoir thermal conductivity W / (m·K); φ is the porosity; V sh is the shale content; λ sand is the sandstone skeleton thermal conductivity W / (m·K); S h is the hydrate saturation; λ hydrate is the hydrate thermal conductivity W / (m·K); S w is the water saturation; λ water is the water thermal conductivity W / (m·K); λ shale is the mudstone thermal conductivity W / (m·K);

[0093] The heat capacity calculation relationship can be the heat capacity mixing law: C p =(1-φ-V sh )·C p,sand +φ·(S h ·C p,hydrate +S w ·C p,water )+V sh ·C p,shale

[0094] Wherein, C p is the reservoir heat capacity J / (m 3 ·K); φ is the porosity; V sh is the shale content; C p,sand is the sandstone skeleton heat capacity J / (m 3 ·K); S h is the hydrate saturation; C p,hydrate is the hydrate heat capacity J / (m 3·K);S w C represents water saturation. p,water Hydrothermal capacity J / (m 3 ·K); C p,shale Heat capacity of mudstone (J / (m³)) 3 ·K).

[0095] Specifically, when constructing the permeability model, the absolute permeability of the reservoir is calculated by incorporating a modified Kozeny-Carman equation with a clay content correction term; the relationship between water saturation, hydrate saturation, and gas-water relative permeability is correlated using the Brooks-Corey model; and the reservoir fluid flow characteristic parameters are clarified through forward modeling of the permeability model. Specifically, the modified Kozeny-Carman equation can be: φ eff =φ·(1-S h -S w,irr )

[0096] Where k0 is the baseline permeability (corrected for clay content) (mD); φ eff c is the effective porosity; d is the porosity index; d is the hydrate blockage coefficient; S h S represents the hydrate saturation level. h,crit φ represents the critical saturation degree of the hydrate; φ represents the porosity; S w,irr To bind water saturation;

[0097] The Brooks-Corey model can be:

[0098] Where, k rw The relative permeability of the aqueous phase; S represents the relative permeability of the water phase at the endpoint. w Water saturation; S w,irr S represents the bound water saturation. g,res n represents the residual gas saturation. w k is the water phase permeability index. rg This refers to the relative permeability of the gas phase. S represents the relative permeability of the gas phase at the endpoint. g n represents the gas saturation level. g It is the gas phase permeability index.

[0099] In constructing the mechanical strength model, the correlation between reservoir compressive strength, internal friction angle and the above parameters is established by combining porosity, hydrate saturation and clay content through an empirical strength model; in constructing the dynamic-static elastic parameter calculation relationship, the dynamic elastic parameters are first calculated by P-wave velocity and S-wave velocity, and then the static elastic parameters are obtained by using the dynamic-static elastic parameter conversion model.

[0100] Specifically, the empirical strength model can be:

[0101] Compressive strength: σ c =α·φ -β ·(1-S h ) γ ·(1-δV sh )

[0102] Where, φ f φ is the internal friction angle. f0 The reference internal friction angle is (°); η is the hydrate strengthening coefficient (° / saturation); S h θ represents hydrate saturation; θ represents the clay weakening coefficient (° / clay content); V sh The content of clay;

[0103] Angle of internal friction: φ f =φ f0 +ηS h -θV sh

[0104] Where φ f φ is the internal friction angle. f0 The reference internal friction angle is (°); η is the hydrate strengthening coefficient (° / saturation); S h θ represents hydrate saturation; θ represents the clay weakening coefficient (° / clay content); V sh This refers to the mud content.

[0105] Furthermore, when calculating the dynamic elastic parameters, the dynamic Young's modulus and dynamic Poisson's ratio are calculated using relevant formulas of elasticity based on the longitudinal wave velocity and transverse wave velocity obtained from the acoustic velocity model and combined with the reservoir saturation density; when performing the dynamic-static elastic parameter conversion, an empirical conversion model containing conversion coefficients and conversion intercepts is used to obtain the static Young's modulus and static Poisson's ratio.

[0106] Specifically, the relevant formulas in elasticity mechanics can be:

[0107] Among them, E dyn V is the dynamic Young's modulus; ρ is the reservoir density; V p V is the longitudinal wave velocity (m / s); s The longitudinal wave velocity (m / s); ν dyn This is the dynamic Poisson's ratio.

[0108] Dynamic-to-static elastic parameter conversion: E static =a·E dyn +b,ν static =c·ν dyn +d

[0109] Among them, E static ν is the static Young's modulus; a is the dynamic-to-static Young's modulus conversion factor; b is the conversion intercept; static is the static Poisson ratio; c is the dynamic-static Poisson ratio conversion coefficient; d is the conversion intercept.

[0110] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.

[0111] The above embodiments are merely illustrative of the technical concept and features of the present invention, and are intended to enable those skilled in the art to understand the content of the present invention and implement it accordingly. They should not be construed as limiting the scope of protection of the present invention. All equivalent changes or modifications made based on the essence of the content of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A multi-field coupling parameter joint inversion method for a seabed muddy silt type hydrate reservoir, characterized in that, The method comprises the steps of: obtaining well logging data and laboratory test data; determining the shale content and the pore water salinity of the reservoir from the well logging data and the laboratory test data; judging the gas phase state of the reservoir according to a natural gas hydrate phase equilibrium model in combination with the shale content and the pore water salinity; constructing a resistivity model and a sonic velocity model considering the influence of the shale content, constructing an objective function based on the resistivity data and the sonic velocity data in the well logging data, and performing global optimization inversion by using a differential optimization algorithm after adding physical constraints to obtain the reservoir porosity and the hydrate saturation; performing forward calculation on the resistivity model and the sonic velocity model by using the reservoir porosity and the hydrate saturation to complete model reliability verification; constructing a thermodynamic model, a permeability model and a mechanical strength model considering the influence of the shale content, and performing forward calculation by using the reservoir porosity, the hydrate saturation and the shale content to obtain reservoir parameters.

2. The method according to claim 1, wherein the reservoir parameters include reservoir thermal conductivity, heat capacity, absolute permeability, gas-water relative permeability, dynamic-static elastic parameters, compressive strength and internal friction angle; the temperature, pore pressure, resistivity, sonic velocity and natural gamma data of the reservoir are collected when the well logging data is obtained; and the mineral composition, mechanical properties and thermal conduction characteristics related data of the reservoir rock sample are collected when the laboratory test data is obtained.

3. The method according to claim 1, wherein the shale content is determined by using the natural gamma data in the well logging data, and the shale content of the reservoir is obtained through natural gamma normalization calculation and related empirical formula; and the pore water salinity is determined by directly obtaining the laboratory test data or indirectly calculating the pore water salinity by using related response parameters in the well logging data.

4. The method according to claim 1, wherein the gas phase state of the reservoir is judged by using a hydrate thermodynamic model based on the fugacity equilibrium theory, and the hydrate stable region of the reservoir is determined by combining the state equation, the shale reservoir salinity enhancement effect and the pore size effect for correction, so as to further judge the gas phase state.

5. The method according to claim 1, wherein the resistivity model is constructed by improving the Waxman-Smits model and integrating the influence of the shale content on the reservoir conductivity, and the relationship between the reservoir resistivity and the water saturation and the hydrate saturation is associated by using the model; and the sonic velocity model is constructed by improving the Thomas-Stieber rock physics hybrid model, and the correlation between the shale content and the longitudinal wave velocity and the transverse wave velocity is established by combining the acoustic characteristics of the sandstone, the shale, the hydrate phase and the fluid phase.

6. The method according to claim 1, wherein the deviation between the calculated values of the resistivity model and the sonic velocity model and the measured values of the well logging is taken as the objective function when the global optimization inversion is performed, the physical value range constraint of the reservoir parameters is added, the optimal solution is iteratively searched by using the differential optimization algorithm, and the reservoir porosity and the hydrate saturation are obtained.

7. The method according to claim 1, wherein ​ ​ ​ ​ ​ ​ ​ ​ ​ In constructing the thermodynamic model, the volume-weighted mixing law is used to incorporate the parameters of shale content, porosity and hydrate saturation to establish the calculation relationship of reservoir thermal conductivity and heat capacity; In the forward calculation of the thermodynamic model, the calculation accuracy of thermal conductivity and heat capacity is improved by combining with the laboratory measurement data, logging statistical results and dynamic correction parameters.

8. The method of claim 1, wherein, In constructing the permeability model, the improved Kozeny-Carman equation is used to calculate the absolute permeability of the reservoir by incorporating the shale content correction term; and the Brooks-Corey model is used to correlate the relationship between water saturation, hydrate saturation and gas-water relative permeability; In the forward calculation of the permeability model, the flow characteristics parameters of the reservoir fluid are determined.

9. The method of claim 1, wherein, In constructing the mechanical strength model, the empirical strength model is used to establish the correlation between the reservoir compressive strength, internal friction angle and the parameters of porosity, hydrate saturation and shale content; In constructing the dynamic-static elastic parameter calculation relationship, the dynamic elastic parameters are calculated through the longitudinal wave velocity and the transverse wave velocity, and then the static elastic parameters are obtained by using the dynamic-static elastic parameter conversion model.

10. The method of claim 9, wherein, In calculating the dynamic elastic parameters, the dynamic Young's modulus and the dynamic Poisson's ratio are calculated by using the relevant formulas of elastic mechanics based on the longitudinal wave velocity and the transverse wave velocity obtained from the acoustic velocity model and in combination with the reservoir saturated density; In the dynamic-static elastic parameter conversion, the empirical conversion model with conversion coefficient and conversion intercept is used to obtain the static Young's modulus and the static Poisson's ratio.

Images