Gas hydrate three-phase coexistence region modeling method and system

Abstract

The invention relates to the field of environmental analysis, and provides a gas hydrate three-phase coexistence region modeling method and system. The method comprises the following steps: acquiring three-dimensional seismic data and original logging data, and preprocessing the original logging data to obtain corrected logging data; respectively performing logging interpretation and structural interpretation according to the corrected logging data and the three-dimensional seismic data to obtain interpretation results; converting a structure interpretation result in the interpretation results into a three-dimensional geological structure model, and performing grid division on the research area by adopting a corner grid system to obtain a structure grid model; and performing attribute modeling on the constructed grid model based on variation function analysis, and predicting porosity and saturation in a three-dimensional space through a Kriging interpolation algorithm to obtain a three-phase coexistence region model. According to the method, the influence of complex geological conditions on the accuracy of the model is reduced, and the precision of the three-phase coexistence region model is remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of environmental analysis technology, and in particular to a method and system for modeling the three-phase coexistence region of gas hydrates. Background Technology

[0002] With the continued growth of global energy demand and the gradual depletion of traditional fossil fuel resources, natural gas hydrates, as an abundant and clean alternative energy source, have attracted widespread attention. Natural gas hydrates are mainly found in deep-sea sediments or permafrost regions, and their unique physicochemical properties make extraction complex and risky. Under actual geological conditions, natural gas hydrates, free gas, and water often coexist in the same reservoir, forming a three-phase coexistence zone. This three-phase coexistence phenomenon is prevalent in many sea areas around the world, becoming an important direction for international energy development research.

[0003] Current modeling techniques primarily focus on single-phase regions of hydrates or free gas, with relatively little research on modeling three-phase coexisting regions. Existing methods, such as vector modeling, sedimentary facies-constrained modeling, and layered modeling, while achieving some success in single-phase modeling, have significant limitations in handling the complex spatial distribution of three-phase coexisting regions. Particularly in three-phase coexisting regions, the interactions and spatial distribution patterns of hydrates, free gas, and water are complex, making it difficult for existing modeling methods to accurately characterize the dynamic boundaries and saturation distribution characteristics between the three phases. This results in insufficient model accuracy and an inability to provide a reliable basis for resource assessment and development in three-phase coexisting regions. Summary of the Invention

[0004] This invention provides a method and system for modeling the three-phase coexistence region of gas hydrates, in order to overcome the shortcomings of the prior art.

[0005] This invention provides a method for modeling the three-phase coexistence region of gas hydrates, comprising: S1: Acquire 3D seismic data and raw well logging data, preprocess the raw well logging data to obtain calibrated well logging data; S2: Perform well logging interpretation and structural interpretation based on the corrected well logging data and the three-dimensional seismic data respectively to obtain the interpretation results; S3: Transform the structural interpretation results in the above interpretation results into a three-dimensional geological structural model, and use a corner grid system to divide the study area into grids to obtain a structural grid model; S4: Based on the variogram analysis, the constructed mesh model is modeled for properties, and the porosity and saturation are predicted in three-dimensional space by the Kriging interpolation algorithm to obtain a three-phase coexistence region model.

[0006] According to the method for modeling the three-phase coexistence region of gas hydrate provided by the present invention, step S1 further includes: S11: Acquire 3D seismic data using seismic acquisition equipment and raw well logging data using well logging tools; S12: For the original logging data, obtain the original resistivity measurement value, the real-time temperature of the downhole instrument, and the laboratory calibration temperature; S13: Perform temperature compensation correction on the original resistivity measurement value to obtain corrected resistivity data; S14: The density curve is corrected based on the phase control density prediction method, and the measurement error is reduced by bandpass filtering to obtain the corrected density data; S15: For non-equidistant logging data, generate equidistant curves through interpolation to ensure consistency with the time sequence of the three-dimensional seismic data, and obtain corrected logging data.

[0007] According to the method for modeling the three-phase coexistence region of gas hydrate provided by the present invention, the interpretation results in step S2 include: Well logging interpretation results, including formation dip data, fracture development data, hydrate saturation data, and permeability data; The structural interpretation results include fault locations and lithological interfaces.

[0008] According to the method for modeling a three-phase coexistence zone of gas hydrate provided by the present invention, step S2, the step of interpreting the logging data based on the corrected logging data, further includes: S211: Extract the dip angle of the formation by resistivity imaging to obtain the dip angle data of the formation; S212: Identify the development of fractures and micro-faults through acoustic longitudinal wave time difference analysis and obtain the fracture development data; S213: Based on the formation dip data and the fracture development data, perform conventional calculations on the hydrate saturation and permeability to obtain preliminary hydrate saturation data and preliminary permeability data; S214: Extract acoustic transverse wave time difference and elemental energy spectrum data, analyze nuclear magnetic resonance data, and obtain porosity data and fine permeability data; S215: Establish a lithological model by using elemental capture energy spectrum, calculate hydrate saturation by combining density data and nuclear magnetic resonance data, obtain the stratigraphic mineral type and dry weight, and obtain the hydrate saturation data and the permeability data.

[0009] According to the method for modeling a three-phase coexistence region of gas hydrate provided by the present invention, step S2, the step of performing structural interpretation based on the three-dimensional seismic data, further includes: S221: Identify the location where the reflected wave is interrupted, trace the interruption location to determine the existence and orientation of the fault, and obtain fault orientation data; S222: Analyze the degree of misalignment of the reflected wave, estimate the dip angle and displacement of the fault, infer the nature of the fault, and obtain the location of the fault; S223: Identify lithological changes through high-reflection-wave contrast interfaces, trace continuous reflection-wave sequences, determine the distribution and thickness variations of different lithologies, and obtain the lithological interfaces; S224: Identify the bending or truncation changes of the reflected wave, compare the reflection characteristics of different sequences, confirm the location of the unconformity, and obtain the structural interpretation results.

[0010] According to the method for modeling the three-phase coexistence region of gas hydrate provided by the present invention, step S3 further includes: S31: Design the grid precision based on the geological complexity of the study area, and select the grid precision in the vertical direction according to different geological layers to obtain the grid division parameters; S32: Based on the grid division parameters, establish a three-dimensional spatial grid based on the corner grid system, incorporate the fault location in the structural interpretation results into the three-dimensional spatial grid, establish the cutting and dislocation relationship between the fault and the rock strata, and obtain the fault grid model; S33: Map the lithological interfaces in the structural interpretation results to the fault grid model, and calibrate and verify the stratigraphic division, lithological distribution and fluid filling through seismic impedance data to obtain the structural grid model.

[0011] According to the modeling method for a three-phase coexistence region of gas hydrate provided by the present invention, in the mesh division parameters in step S31, a 10m×10m mesh is used in the plane of the three-phase coexistence region, and a 0.1m vertical mesh is used in the hydrate-bearing sand layer.

[0012] According to the method for modeling a three-phase coexistence region of gas hydrate provided by the present invention, in step S4, when performing property modeling on the constructed mesh model based on variogram analysis, the expression of the variogram function is as follows:

[0013] in, The lag distance, It is a variation function. Distance Data pairs of quantities The sample point numbers selected in constructing the mesh model. To construct the sample point locations in the mesh model, For sample point locations The attribute value at that location, For sample point locations The attribute value at that location.

[0014] According to the method for modeling the three-phase coexistence region of gas hydrate provided by the present invention, in step S4, the prediction expression for porosity and saturation in three-dimensional space using the Kriging interpolation algorithm is as follows:

[0015] in, To predict the obtained attribute values, The number of sample points selected in constructing the mesh model, For the first The weighting coefficients for each sample point.

[0016] The present invention also provides a system for modeling a three-phase coexistence region of gas hydrates, for performing a method for modeling a three-phase coexistence region of gas hydrates as described in any of the preceding claims, comprising: Acquisition module: used to acquire 3D seismic data and raw well logging data, and to preprocess the raw well logging data to obtain calibrated well logging data; Interpretation module: used to perform well logging interpretation and structural interpretation based on the corrected well logging data and the three-dimensional seismic data, respectively, to obtain interpretation results; The partitioning module is used to transform the structural interpretation results in the interpretation results into a three-dimensional geological structural model. It uses a corner grid system to partition the study area into a structural grid model. Modeling module: used to perform attribute modeling on the constructed mesh model based on variogram analysis, and to predict porosity and saturation in three-dimensional space using the Kriging interpolation algorithm to obtain a three-phase coexistence region model.

[0017] The method for modeling the three-phase coexistence region of gas hydrates provided by this invention can effectively solve the technical problems of insufficient modeling accuracy and inaccurate spatial distribution characterization in the prior art through a systematic data processing flow and refined modeling technology. First, the method for modeling three-phase coexistence zones of gas hydrates provided by this invention effectively reduces measurement errors by performing environmental correction and time conversion on well logging data, providing a high-quality data foundation for subsequent modeling. Subsequently, by combining well logging interpretation and structural interpretation techniques, it can comprehensively acquire key geological information such as formation dip angle, fracture development, fault location, and lithological interfaces, accurately depicting the geological structural characteristics of the three-phase coexistence zone. Second, by employing high-resolution grid generation and a corner grid system, the model's ability to represent complex geological structures is significantly improved, enabling accurate reproduction of fault cutting relationships and strata displacement characteristics. Furthermore, this invention also uses variogram analysis and Kriging interpolation algorithms for attribute modeling, which can reasonably characterize the distribution patterns of porosity and saturation in three-dimensional space. Compared with traditional modeling methods, it has higher prediction accuracy and geological rationality. The final three-phase coexistence zone model can accurately reflect the spatial distribution characteristics of the hydrate phase, free gas phase, and water phase, providing reliable technical support for the assessment and development of natural gas hydrate resources. Attached Figure Description

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

[0019] Figure 1 A schematic diagram of a method for modeling a three-phase coexistence region of gas hydrates provided by the present invention; Figure 2 A schematic diagram of a modeling system for a three-phase coexistence region of gas hydrates provided by the present invention; Figure 3 This diagram illustrates a comparison of verification results for various modeling methods provided by this invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, embodiments of this invention, and should not be construed as limiting the invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention. In the description of this invention, it should be understood that the terminology used is for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0021] The embodiments of the present invention are described below with reference to the figures.

[0022] like Figure 1 As shown, this invention provides a method for modeling a three-phase coexistence region of gas hydrates, comprising: S1: Acquire 3D seismic data and raw well logging data, preprocess the raw well logging data to obtain calibrated well logging data.

[0023] Step S1 further includes: S11: Acquire 3D seismic data through seismic acquisition equipment and obtain raw well logging data through well logging tools.

[0024] In step S11, the present invention first acquires seismic and well logging data. The main seismic acquisition equipment of the present invention includes an 8x20 cubic inch air gun array, a Seal digital recording system, a 24-bit Seal digital cable, and a SYS3 depth control device. The frequency range of the raw seismic data is 6-20Hz. The seismic data is processed using Focus software, preserving the true amplitude during processing to clearly interpret seismic characteristics associated with natural gas hydrates, especially BSR. Its dominant frequency is 40Hz, with online and cross-line spacing of 50 and 6.25m, respectively, including sonic logging, high-resolution resistivity imaging, and density data. For well logging data, the present invention uses logging-while-drilling tools MicroScopeHD, NeoScope, Sonicscope, and ProVISION to acquire gamma ray, wellbore, resistivity, high-resolution resistivity images, neutrons, density, P-wave transit time, S-wave transit time, T2 spectrum, and pore structure data.

[0025] S12: For the original logging data, obtain the original resistivity measurement value, the real-time temperature of the downhole instrument, and the laboratory calibration temperature.

[0026] In step S12, the present invention extracts the original resistivity measurement value from the original logging data and records the real-time temperature of the downhole instrument at the measurement time. This temperature is obtained in real time through a downhole temperature sensor. The laboratory calibration temperature is the ambient temperature when the instrument is calibrated under surface laboratory conditions, typically set to 25 degrees Celsius. Since there is a significant difference between the downhole ambient temperature and the laboratory calibration temperature, temperature changes can cause systematic deviations in the resistivity measurement value. Therefore, in step S12, the present invention needs to obtain these three temperature-related parameters for subsequent correction.

[0027] S13: Perform temperature compensation correction on the original resistivity measurement value to obtain corrected resistivity data.

[0028] In step S13, this invention uses a temperature compensation formula to correct the original resistivity measurement value. The temperature compensation coefficient is selected based on the formation fluid characteristics and downhole temperature conditions. Specifically, firstly, the temperature difference between the real-time temperature of the downhole instrument and the laboratory calibration temperature is calculated. Then, the temperature difference is multiplied by the temperature compensation coefficient to obtain the temperature compensation factor. Next, the temperature compensation factor is calculated with the original resistivity measurement value to obtain the corrected resistivity data. This invention eliminates the influence of temperature changes on resistivity measurement through correction, ensuring that the corrected resistivity data accurately reflects the resistivity characteristics of the formation.

[0029] S14: The density curve is corrected based on the phase control density prediction method, and the measurement error is reduced by bandpass filtering to obtain the corrected density data.

[0030] In step S14, the present invention corrects the density curve based on the phase-controlled density prediction method. Specifically, firstly, the lithofacies types in the well logging data are identified, including sandstone, mudstone, and transitional facies. Then, a corresponding density prediction relationship is established for each lithofacies type. Subsequently, by comparing the measured density data with the predicted density data, outliers and deviations generated during the measurement process are identified. Next, bandpass filtering is used to filter the density curve. The passband frequency range set by the bandpass filter can retain the true change signal of formation density while filtering out high-frequency noise and low-frequency drift introduced during the measurement process. Finally, after filtering, the corrected density data is obtained.

[0031] S15: For non-equidistant logging data, generate equidistant curves through interpolation to ensure consistency with the time sequence of the three-dimensional seismic data, and obtain corrected logging data.

[0032] In step S15, this invention performs interpolation processing on non-equidistant logging data. The original logging data has uneven sampling intervals along the depth direction, exhibiting varying density. The interpolation process employs linear interpolation or cubic spline interpolation methods. Based on known non-equidistant sampling point data, the data values ​​at equidistant sampling points are calculated. The sampling interval of the equidistant curve is determined according to the time sampling rate of the 3D seismic data, ensuring that the sampling density of the logging data in the time domain remains consistent with the seismic data. Subsequently, this invention performs time-depth conversion on the equidistant logging data. A time-depth relationship is established based on the sonic logging data, converting the depth-domain logging data to the time domain. This ensures that the converted logging data and the 3D seismic data accurately correspond on the time axis, ultimately yielding the corrected logging data.

[0033] S2: Perform well logging interpretation and structural interpretation based on the corrected well logging data and the three-dimensional seismic data respectively to obtain the interpretation results.

[0034] The interpretation results in step S2 include: well logging interpretation results, which include formation dip data, fracture development data, hydrate saturation data, and permeability data; and structural interpretation results, which include fault locations and lithological interfaces.

[0035] Specifically, step S2, the step of interpreting the logging data based on the corrected logging data, further includes: S211: Extract the dip angle of the formation by resistivity imaging to obtain the dip angle data of the formation.

[0036] In step S211, this invention acquires a high-resolution resistivity image of the wellbore using MicroScopeHD resistivity imaging technology. This image presents the resistivity distribution characteristics of the wellbore in the form of a two-dimensional unfolded diagram. In the resistivity image, geological structures such as formation bedding planes, fractures, and faults exhibit linear or curvilinear resistivity characteristics. This invention extracts the apparent dip angle and apparent dip direction data of the formation bedding planes by identifying the tilt angle and azimuth of these linear features in the image. Subsequently, based on the wellbore trajectory azimuth and inclination angle, the apparent dip angle and apparent dip direction are converted into the true dip angle and true dip direction of the formation, obtaining formation dip angle data, which includes the dip azimuth and dip angle values ​​of the formation.

[0037] S212: Identify the development of fractures and micro-faults through acoustic longitudinal wave time difference analysis, and obtain the fracture development data.

[0038] In step S212, this invention acquires P-wave transit time data using the Sonicscope acoustic logging tool. P-wave transit time reflects the velocity characteristics of acoustic waves propagating in the formation. In intact formations, P-wave transit time exhibits a relatively stable numerical range. When fractures or micro-faults exist in the formation, the acoustic wave propagation path is disrupted, leading to sudden increases or abnormal fluctuations in P-wave transit time. This invention identifies points of sudden increases and abnormal fluctuation segments by analyzing the changing characteristics of the P-wave transit time curve. For points of sudden increases, the location of fracture development is determined; for abnormal fluctuation segments, resistivity images and wellbore data are combined to determine whether they are micro-fault development segments, and the distribution depth, fracture aperture, and development density of fractures and micro-faults in the wellbore are statistically analyzed to obtain fracture development data.

[0039] S213: Based on the formation dip data and the fracture development data, perform conventional calculations on the hydrate saturation and permeability to obtain preliminary hydrate saturation data and preliminary permeability data.

[0040] In step S213, this invention performs conventional calculations of hydrate saturation and permeability based on formation dip data and fracture development data. The formation dip data reflects the structural morphology of the reservoir, and the fracture development data reflects the heterogeneity of the reservoir. The calculation method used in this invention is the Archie formula for hydrate saturation. Based on corrected resistivity data, porosity data, and formation water resistivity, preliminary hydrate saturation data is calculated using the saturation index and cementation index parameters of the Archie formula. Permeability calculation employs the Karman-Kozani equation. Based on porosity data and grain size analysis data, combined with fracture development data, the permeability is corrected to obtain preliminary permeability data.

[0041] S214: Extract acoustic transverse wave time difference and elemental energy spectrum data, analyze nuclear magnetic resonance data, and obtain porosity data and fine permeability data.

[0042] In step S214, this invention extracts shear wave transit time data using the Sonicscope tool. Shear wave transit time is related to the shear modulus and consolidation degree of the formation. Simultaneously, elemental energy spectrum data is obtained using the NeoScope tool. This data records the content distribution of various elements in the formation, including the weight percentage of elements such as silicon, calcium, iron, and sulfur. Specifically, this invention obtains T2 spectrum data using the ProVISION nuclear magnetic resonance logging tool. The T2 spectrum reflects the relaxation time distribution of fluids in the formation pores. Subsequently, based on the distribution characteristics of the T2 spectrum, the pores are divided into three categories: micropores, mesopores, and macropores. The volume percentage of each type of pore is calculated to obtain porosity data. Then, based on the geometric mean of the T2 spectrum and the porosity data, the formation permeability is calculated using an SDR model or a Coates model to obtain refined permeability data.

[0043] S215: Establish a lithological model by using elemental capture energy spectrum, calculate hydrate saturation by combining density data and nuclear magnetic resonance data, obtain the stratigraphic mineral type and dry weight, and obtain the hydrate saturation data and the permeability data.

[0044] Further, in step S215, this invention analyzes the content of various mineral components in the formation using elemental capture energy spectroscopy (EVS). EDS measures the gamma-ray energy spectrum generated by the reaction of neutrons with the nuclear elements in the formation to quantitatively analyze the elemental composition. In the analysis, based on the content of major elements such as silicon, calcium, and iron, this invention employs the ELAN multi-mineral interpretation model to convert the elemental content into the volume fractions of minerals such as quartz, feldspar, calcite, dolomite, and clay, establishing a lithological model. The obtained lithological model provides the dry weight and volume percentage of each mineral component in the formation. Subsequently, hydrate saturation is calculated by combining density data and NMR data. Density data reflects the bulk density of the formation, while NMR data reflects the type of pore fluid. Since the density of hydrates is between that of water and free gas, this invention identifies the presence of hydrates through abnormal changes in density values. In the NMR T2 spectrum, hydrates occupy pore space, leading to a decrease in effective porosity. Therefore, this invention calculates hydrate saturation data by comparing the difference between NMR porosity and conventional porosity. Finally, by combining the mineral composition information from the lithology model, the matrix density and porosity parameters in the hydrate saturation calculation were corrected to obtain more accurate hydrate saturation and permeability data.

[0045] In step S2, the step of constructing and interpreting the three-dimensional seismic data further includes: S221: Identify the location where the reflected wave is interrupted, trace the location of the interruption to determine the existence and orientation of the fault, and obtain fault orientation data.

[0046] In step S221, the present invention identifies the interruption location of reflected waves in the profile of three-dimensional seismic data. Reflected waves are the reflected signals generated by seismic waves at different lithological interfaces. On the seismic profile, continuous reflected wave layers represent continuous stratigraphic interfaces. When the strata are cut by a fault, the reflected waves are interrupted at the fault location, manifested as a sudden termination or displacement of the reflected waves. By tracing the continuity of reflected waves one by one, the coordinates of the location where the reflected waves are interrupted are marked. Subsequently, the interruption location is traced in three-dimensional space, and the interruption points on different profiles are connected to form the spatial distribution pattern of the fault, determine the strike and azimuth of the fault, and obtain the fault strike data.

[0047] S222: Analyze the degree of misalignment of the reflected wave, estimate the dip angle and displacement of the fault, infer the nature of the fault, and obtain the location of the fault.

[0048] In step S222, this invention analyzes the degree of displacement of reflected waves on both sides of the fault. Reflected wave displacement refers to the difference in vertical position of reflected waves from the same stratigraphic interface on both sides of the fault. This invention obtains the fault's displacement value by measuring the time difference or depth difference of the same reflecting layer on the hanging wall and footwall. Based on the dip angle and displacement direction of the reflected waves, the dip angle of the fault is determined; the dip angle is the angle between the fault plane and the horizontal plane. By analyzing the relative motion relationship between the strata on the hanging wall and footwall, the nature of the fault is determined, including normal faults, reverse faults, or strike-slip faults. Normal faults are characterized by a relative subsidence of the hanging wall, reverse faults by a relative uplift of the hanging wall, and strike-slip faults by horizontal displacement. Finally, by comprehensively considering the fault's strike, dip angle, displacement, and nature, the three-dimensional spatial coordinates and geometric parameters of the fault location are determined.

[0049] S223: Identify lithological changes through high-reflection-wave contrast interfaces, trace continuous reflection wave sequences, determine the distribution and thickness variations of different lithologies, and obtain the lithological interfaces.

[0050] In step S223, this invention identifies lithological changes by recognizing high-contrast interfaces in seismic profiles. The amplitude intensity of the reflected waves is determined by the difference in wave impedance between the rock layers on both sides of the interface. Significant differences in wave impedance exist between sandstone and mudstone, and between hydrated and non-hydrated layers, which manifest as strong-amplitude reflected waves in seismic profiles. This invention identifies the distribution of lithological interfaces by calibrating the spatial location of strong-amplitude reflected waves; subsequently, it traces continuous reflected wave sequences, which are a series of parallel or subparallel combinations of reflected waves representing stratigraphic assemblages of a specific depositional period; finally, by comparing the thickness variations of reflected wave sequences at different locations, it determines the lateral distribution characteristics and vertical stacking relationships of the strata, obtaining data on the distribution range and thickness variations of different lithologies, and ultimately determining the three-dimensional spatial morphology of the lithological interfaces.

[0051] S224: Identify the bending or truncation changes of the reflected wave, compare the reflection characteristics of different sequences, confirm the location of the unconformity, and obtain the structural interpretation results.

[0052] Furthermore, this invention identifies the location of unconformities by recognizing the bending or truncation changes of reflected waves. On seismic profiles, reflected waves from the underlying strata are truncated by the unconformity, manifesting as abrupt termination of the reflected waves; reflected waves from the overlying strata exhibit erosion or overlapping geometry, showing bending towards the unconformity. This invention marks the location of the unconformity by identifying these truncation and bending features; subsequently, it compares the differences in reflected wave characteristics between the strata above and below the unconformity, including the amplitude, frequency, and continuity of the reflected waves. Significant differences in the reflection characteristics between the upper and lower strata confirm the existence of the unconformity; finally, the identified unconformity is comprehensively analyzed in conjunction with fault locations and lithological interfaces to form a complete structural interpretation result. This result includes spatial distribution information of all important geological interfaces and structural elements within the study area.

[0053] S3: The structural interpretation results in the above interpretation results are transformed into a three-dimensional geological structural model. The study area is divided into grids using a corner grid system to obtain a structural grid model.

[0054] Step S3 further includes: S31: Design the grid precision according to the geological complexity of the study area, and select the grid precision in the vertical direction according to different geological layers to obtain the grid division parameters; In the grid division parameters in step S31, a 10m×10m grid is used in the plane of the three-phase coexistence zone, and a 0.1m vertical grid is used in the hydrate-bearing sand layer.

[0055] Furthermore, in step S31, the present invention designs the grid precision based on the geological complexity of the study area. The geological complexity is comprehensively evaluated by the number of faults, the frequency of lithological changes, and the amplitude of stratigraphic dip changes. In the three-phase coexistence zone, due to the drastic spatial distribution changes of hydrates, free gas, and water, the planar grid precision is designed to be 10m × 10m, which can capture the boundary positions and saturation gradient changes between the three phases. Vertically, the present invention selects the grid precision based on the thickness and internal heterogeneity of different geological layers. In hydrate-bearing sand layers, since the distribution of hydrates is controlled by bedding and pore structure, the vertical variation scale is small, and a vertical grid precision of 0.1m is adopted. For mudstone caprock and underlying strata, the vertical variation is relatively gentle, and a vertical grid precision of 0.5m to 1m is adopted. Finally, the obtained grid division parameters include planar grid size, vertical grid size, total number of grids, and grid distribution scheme, which are recorded in the grid parameter table.

[0056] S32: Based on the grid division parameters, establish a three-dimensional spatial grid based on the corner grid system, incorporate the fault location in the structural interpretation results into the three-dimensional spatial grid, establish the cutting and dislocation relationship of the fault on the rock strata, and obtain the fault grid model.

[0057] Furthermore, in step S32, this invention establishes a three-dimensional spatial grid of a corner grid system based on grid division parameters. The corner grid system is a method that describes the grid geometry by defining the coordinates of eight corner points of a grid cell. The corner grid can accurately represent complex structures such as fault cutting and stratum deformation. Specifically, firstly, planar grid nodes are generated within the study area according to a planar grid size of 10m × 10m, and the node coordinates adopt the geodetic coordinate system of the study area. Then, vertical grid layers are generated in the depth direction according to the vertical grid accuracy, and the depth value of each layer is determined according to the bedding depth data in the structural interpretation results. Next, the fault location data is imported into the three-dimensional spatial grid, and the fault location includes the set of spatial coordinate points of the fault plane. At the fault plane location, this invention cuts continuous grid cells into two independent grid blocks along the fault plane, with the hanging wall grid block and the footwall grid block being discontinuous at the fault plane. Subsequently, this invention adjusts the relative positions of the hanging wall and footwall grid blocks according to the fault drop value to realize the misalignment relationship of the strata on both sides of the fault. The grid cells on both sides of the fault plane are connected by the fault transfer coefficient, which reflects the fault's ability to block or conduct fluid flow. After fault cutting, the final fault mesh model is obtained, which includes the corner coordinates of all mesh elements, mesh connection relationships, and fault transmission coefficients.

[0058] S33: Map the lithological interfaces in the structural interpretation results to the fault grid model, and calibrate and verify the stratigraphic division, lithological distribution and fluid filling through seismic impedance data to obtain the structural grid model.

[0059] In step S33, the present invention maps the lithological interfaces in the structural interpretation results to a fault grid model. The lithological interface data includes a set of spatial coordinate points for each lithological interface. The mapping process involves matching the spatial coordinates of the lithological interfaces with the depth range of the grid cells to determine the lithological type of each grid cell. For grid cells that cross lithological interfaces, the grid cell is divided into upper and lower sub-cells based on the proportion of the interface within the grid cell, and each sub-cell is assigned different lithological properties.

[0060] Subsequently, the stratigraphic division was calibrated and verified using seismic impedance data. Seismic impedance is equal to the product of stratigraphic density and seismic wave velocity, reflecting the lithology and fluid properties of the strata. This invention extracts seismic impedance volume data from three-dimensional seismic data, converting the position coordinates of each grid cell in the grid model into the time coordinates of the seismic data volume, and extracting the seismic impedance value at the corresponding location.

[0061] After extraction, this invention compares the extracted seismic impedance values ​​with the lithological types of the grid cells. Sandstone typically has a higher seismic impedance than mudstone; the presence of hydrates increases seismic impedance, while the presence of free gas decreases it. Through comparative analysis, regions in the grid model where the lithological distribution is inconsistent with the seismic impedance data are identified, and the lithological interface positions in these regions are adjusted and corrected. Simultaneously, based on the numerical range of the seismic impedance, the fluid filling type within the grid cells is determined, distinguishing between hydrate layers, free gas layers, and aquifers. Finally, after calibration and verification, a structural grid model is obtained, whose stratigraphic division, lithological distribution, and fluid filling are consistent with the seismic data.

[0062] S4: Based on the variogram analysis, the constructed mesh model is modeled for properties, and the porosity and saturation are predicted in three-dimensional space by the Kriging interpolation algorithm to obtain a three-phase coexistence region model.

[0063] In step S4, when performing attribute modeling on the constructed mesh model based on variogram analysis, the expression of the variogram function used is as follows:

[0064] in, The lag distance, It is a variation function. Distance Data pairs of quantities The sample point numbers selected in constructing the mesh model. To construct the sample point locations in the mesh model, For sample point locations The attribute value at that location, For sample point locations The attribute value at that location.

[0065] In step S4, this invention performs attribute modeling on the constructed mesh model based on variogram analysis, aiming to extend discrete wellpoint measurement data to the entire three-dimensional mesh space. Variation function analysis is used to quantify the spatial variation characteristics of attribute values. In the variogram expression, The sample point number selected in the construction grid model is the sample point sourced from the well logging interpretation results at the well location, including porosity data and hydrate saturation data. Indicates the first The spatial coordinates of each sample point, including planar coordinates and depth coordinates; Indicates the location of the sample point The attribute value at that location is either porosity or saturation. The lag distance represents the spatial distance between two sample points. The lag distance starts from 0 and increases with a certain step size, which is set to an integer multiple of the planar grid size.

[0066] In variogram analysis, for each lag distance This invention searches for the distance between all sample points. The sample point pairs, which include location and location Two sample points, Indicates the location of the sample point The attribute value at the location; count all distances The number of sample point pairs, denoted as Then, the attribute value difference for each pair of sample points is calculated; the difference equals... minus Square the differences and sum them up. Divide the sum by 2. , to obtain the lag distance Corresponding variation function value .

[0067] Subsequently, the lag distance was changed. By taking the value of , repeating the above calculation process, a series of lag distances and their corresponding variogram values ​​are obtained. Plot the variogram curve. The obtained variogram curve shows that it rises with the increase of lag distance and tends to stabilize. The rate of rise of the curve reflects the spatial continuity of the attribute value, and the stable value of the curve reflects the overall degree of variation of the attribute value.

[0068] In step S4, the prediction expression for porosity and saturation in three-dimensional space using the Kriging interpolation algorithm is as follows:

[0069] in, To predict the obtained attribute values, The number of sample points selected in constructing the mesh model, For the first The weighting coefficients for each sample point.

[0070] In step S4, the present invention further predicts porosity and saturation in three-dimensional space using the Kriging interpolation algorithm. Kriging interpolation is an optimal unbiased linear estimation method based on the variogram function. In the prediction expression, This indicates the location of the point to be predicted in the constructed mesh model, which corresponds to the center coordinates of the mesh cell; This indicates the points to be predicted obtained from the prediction. The attribute value at the location; This indicates the number of sample points selected in constructing the mesh model. The selection rule is based on the number of points to be predicted. Search all sample points within a certain range around the variogram. The search range is determined based on the range parameter of the variogram, which is usually set to 1.5 to 2 times the range. Indicates the first The weighting coefficients for each sample point reflect their contribution to the prediction result. These weighting coefficients are calculated by solving the Kriging equations, a system of linear equations established based on the unbiased estimation condition and the minimum estimation variance condition. The coefficient matrix in the equations consists of the variogram values ​​between the sample points, and the right-hand side consists of the variogram values ​​between the sample points and the point to be predicted. Solving the equations yields... Weight coefficients And a Lagrange multiplier. Finally, the weighting coefficients... Attribute values ​​of corresponding sample points Multiply the results by summing them over all sample points to obtain the predicted attribute value at the point to be predicted. .

[0071] Subsequently, the above prediction process is repeated for each grid cell in the constructed grid model, extending the porosity data to three-dimensional space to form a porosity model, which records the porosity value of each grid cell. Similarly, hydrate saturation data is extended to three-dimensional space to form a hydrate saturation model, which records the hydrate saturation value of each grid cell. Simultaneously, this invention establishes a free gas saturation model based on free gas saturation data from well logging interpretation results using a Kriging interpolation algorithm. Within each grid cell, the hydrate saturation... Free gas saturation and water saturation Satisfying the constraint condition that the sum equals 1, calculate the water saturation based on the predicted hydrate saturation and free gas saturation. .

[0072] Ultimately, this invention integrates the porosity model, hydrate saturation model, free gas saturation model, and water saturation model to obtain a three-phase coexistence region model, which fully describes the distribution characteristics of the hydrate phase, free gas phase, and water phase in three-dimensional space.

[0073] like Figure 2 As shown, the present invention also provides a gas hydrate three-phase coexistence region modeling system for performing a gas hydrate three-phase coexistence region modeling method as described in any of the above claims, comprising: Acquisition module 100: used to acquire three-dimensional seismic data and raw well logging data, and to preprocess the raw well logging data to obtain calibrated well logging data; Interpretation module 200: used to perform well logging interpretation and structural interpretation based on the corrected well logging data and the three-dimensional seismic data respectively, and obtain interpretation results; Module 300: Used to convert the structural interpretation results in the interpretation results into a three-dimensional geological structural model, and to divide the study area into grids using a corner grid system to obtain a structural grid model; Modeling module 400: Used to perform attribute modeling on the constructed mesh model based on variogram analysis, and to predict porosity and saturation in three-dimensional space using the Kriging interpolation algorithm to obtain a three-phase coexistence region model.

[0074] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0075] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the gas hydrate three-phase coexistence region modeling method described in various embodiments or some parts of embodiments.

[0076] In addition, the present invention also conducted error analysis and verification on the proposed model.

[0077] In the verification process, this invention utilizes seismic data from well locations in the work area to conduct a preliminary assessment of the subsurface structure. Due to the limited number of control wells intersecting with the triphase zone within the hydrate study area, some differences were observed in the distribution of model data, coarsened data, and raw data. However, the model effectively captures high porosity and high saturation regions, as demonstrated by the seismic profiles, indicating a relatively high overall reliability of the newly constructed model. When comparing models generated using three different modeling methods, values ​​within the same model were selected for analysis. The comparison shows that the saturation data obtained by the deterministic complex morphology geological modeling method matches the variation range of the oil well data better. In contrast, the porosity values ​​obtained by Kriging interpolation and sequential Gaussian models exhibit greater fluctuations.

[0078] Specific analysis results are as follows Figure 3 As shown, Figure 3 (a) shows a comparison of saturation data established using three modeling methods. Figure 3 (b) shows a comparison of porosity data established using three modeling methods. Figure 3 (c) shows the saturation probability distribution histograms established by the three modeling methods. Figure 3 In the middle (d), the porosity probability distribution histograms are established using the three modeling methods. Figure 3 Analysis shows that the three-dimensional multi-scale model established through deterministic complex morphological geological modeling is basically consistent with general geological understanding and actual well logging data. In addition, the results confirm the superiority of layered modeling over models obtained from a single modeling method.

[0079] In the well logging interpretation stage, this invention extracts formation dip data using resistivity imaging technology and identifies fracture and micro-fault development through acoustic P-wave time-of-flight analysis. This multi-parameter comprehensive analysis method comprehensively reveals the internal structural characteristics of the reservoir, providing detailed geological evidence for understanding the formation mechanism and fluid migration patterns of three-phase coexistence zones. Compared to single-parameter analysis methods, this method offers greater reliability and more comprehensive interpretation capabilities. Furthermore, this invention further extracts acoustic S-wave time-of-flight and elemental energy spectrum data, and conducts in-depth analysis of nuclear magnetic resonance (NMR) data to obtain more refined porosity and permeability parameters. These refined reservoir properties are crucial for accurately assessing the reservoir capacity and development potential of three-phase coexistence zones. Moreover, this invention establishes a lithological model using elemental capture energy spectrum and calculates hydrate saturation by combining density and NMR data, obtaining detailed petrological information such as formation mineral types and dry weight. This comprehensive analysis method based on multiple well logging parameters offers higher accuracy compared to traditional single-parameter calculation methods, effectively reducing the inherent ambiguity of single well logging methods and making the calculation results of hydrate saturation and permeability more accurate and reliable.

[0080] In the structural interpretation stage, this method accurately characterizes the fracture system within the study area through systematic fault identification and characterization techniques. Faults, as crucial channels for fluid migration, are essential for understanding the formation process of three-phase coexistence zones and predicting fluid distribution patterns. This method identifies lithological variations through high-contrast interfaces and traces continuous reflective wave sequences to determine the distribution and thickness variations of different lithologies. This enables precise delineation of the vertical layering structure of reservoirs, and the detailed sequence delineation provides an accurate geological framework for subsequent grid and attribute modeling. In the gridding stage, this method employs a high-resolution grid design, significantly improving grid accuracy compared to traditional methods. It accurately captures the spatial variation details within the three-phase coexistence zone, allowing the model to more realistically reflect the boundary positions and saturation gradient changes between the three phases.

[0081] This invention also identifies the spatial variation characteristics of porosity and saturation through variogram analysis, and then performs three-dimensional spatial prediction by combining it with the Kriging interpolation algorithm. This can reasonably characterize the continuity and variability of attribute parameters in space, so that the prediction results not only meet the constraints of known data, but also reflect the spatial distribution characteristics of geological laws. Compared with simple interpolation methods, it has stronger geological rationality and higher prediction accuracy. The three-phase coexistence zone model constructed in the end can accurately reflect the spatial distribution law of hydrate phase, free gas phase and water phase, and provide reliable geological model support for the development and utilization of natural gas hydrate resources.

[0082] Overall, the high-resolution 3D geological modeling technology employed in this invention, targeting the three-phase coexistence region of natural gas hydrates (NGHs), demonstrates the enormous potential and advantages of this technology in well logging. It provides powerful assistance for geological steering of exploration and appraisal wells and development horizontal wells, optimizing and adjusting drilling trajectories, thereby improving the success rate and economic benefits of oil and gas development. Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for modeling a three-phase coexistence region of gas hydrates, characterized in that, include: S1: Acquire 3D seismic data and raw well logging data, preprocess the raw well logging data to obtain calibrated well logging data; S2: Perform well logging interpretation and structural interpretation based on the corrected well logging data and the three-dimensional seismic data respectively to obtain the interpretation results; S3: Transform the structural interpretation results in the above interpretation results into a three-dimensional geological structural model, and use a corner grid system to divide the study area into grids to obtain a structural grid model; S4: Based on the variogram analysis, the constructed mesh model is modeled for properties, and the porosity and saturation are predicted in three-dimensional space by the Kriging interpolation algorithm to obtain a three-phase coexistence region model.

2. The method for modeling a three-phase coexistence region of gas hydrates according to claim 1, characterized in that, Step S1 further includes: S11: Acquire 3D seismic data using seismic acquisition equipment and raw well logging data using well logging tools; S12: For the original logging data, obtain the original resistivity measurement value, the real-time temperature of the downhole instrument, and the laboratory calibration temperature; S13: Perform temperature compensation correction on the original resistivity measurement value to obtain corrected resistivity data; S14: The density curve is corrected based on the phase control density prediction method, and the measurement error is reduced by bandpass filtering to obtain the corrected density data; S15: For non-equidistant logging data, generate equidistant curves through interpolation to ensure consistency with the time sequence of the three-dimensional seismic data, and obtain corrected logging data.

3. The method for modeling a three-phase coexistence region of gas hydrates according to claim 1, characterized in that, The interpretation results in step S2 include: Well logging interpretation results, including formation dip data, fracture development data, hydrate saturation data, and permeability data; The structural interpretation results include fault locations and lithological interfaces.

4. The method for modeling a three-phase coexistence region of gas hydrate according to claim 3, characterized in that, Step S2, the step of interpreting the logging data based on the corrected logging data, further includes: S211: Extract the dip angle of the formation by resistivity imaging to obtain the dip angle data of the formation; S212: Identify the development of fractures and micro-faults through acoustic longitudinal wave time difference analysis and obtain the fracture development data; S213: Based on the formation dip data and the fracture development data, perform conventional calculations on the hydrate saturation and permeability to obtain preliminary hydrate saturation data and preliminary permeability data; S214: Extract acoustic transverse wave time difference and elemental energy spectrum data, analyze nuclear magnetic resonance data, and obtain porosity data and fine permeability data; S215: Establish a lithological model by using elemental capture energy spectrum, calculate hydrate saturation by combining density data and nuclear magnetic resonance data, obtain the stratigraphic mineral type and dry weight, and obtain the hydrate saturation data and the permeability data.

5. The method for modeling a three-phase coexistence region of gas hydrates according to claim 3, characterized in that, Step S2, the step of constructing and interpreting the three-dimensional seismic data, further includes: S221: Identify the location where the reflected wave is interrupted, trace the interruption location to determine the existence and orientation of the fault, and obtain fault orientation data; S222: Analyze the degree of misalignment of the reflected wave, estimate the dip angle and displacement of the fault, infer the nature of the fault, and obtain the location of the fault; S223: Identify lithological changes through high-reflection-wave contrast interfaces, trace continuous reflection-wave sequences, determine the distribution and thickness variations of different lithologies, and obtain the lithological interfaces; S224: Identify the bending or truncation changes of the reflected wave, compare the reflection characteristics of different sequences, confirm the location of the unconformity, and obtain the structural interpretation results.

6. The method for modeling a three-phase coexistence region of gas hydrates according to claim 1, characterized in that, Step S3 further includes: S31: Design the grid precision based on the geological complexity of the study area, and select the grid precision in the vertical direction according to different geological layers to obtain the grid division parameters; S32: Based on the grid division parameters, establish a three-dimensional spatial grid based on the corner grid system, incorporate the fault location in the structural interpretation results into the three-dimensional spatial grid, establish the cutting and dislocation relationship between the fault and the rock strata, and obtain the fault grid model; S33: Map the lithological interfaces in the structural interpretation results to the fault grid model, and calibrate and verify the stratigraphic division, lithological distribution and fluid filling through seismic impedance data to obtain the structural grid model.

7. The method for modeling a three-phase coexistence region of gas hydrates according to claim 6, characterized in that, In the grid division parameters mentioned in step S31, a 10m×10m grid is used in the plane of the three-phase coexistence zone, and a 0.1m vertical grid is used in the hydrate-bearing sand layer.

8. The method for modeling a three-phase coexistence region of gas hydrates according to claim 1, characterized in that, In step S4, when performing attribute modeling on the constructed mesh model based on variogram analysis, the expression of the variogram function used is as follows: in, The lag distance, It is a variation function. Distance Data pairs of quantities The sample point numbers selected in constructing the mesh model. To construct the sample point locations in the mesh model, For sample point locations The attribute value at that location, For sample point locations The attribute value at that location.

9. The method for modeling a three-phase coexistence region of gas hydrates according to claim 8, characterized in that, In step S4, the prediction expressions for porosity and saturation in three-dimensional space using the Kriging interpolation algorithm are as follows: in, To predict the obtained attribute values, The number of sample points selected in constructing the mesh model, For the first The weighting coefficients for each sample point.

10. A system for modeling a three-phase coexistence region of gas hydrates, used to execute a method for modeling a three-phase coexistence region of gas hydrates as described in any one of claims 1 to 9, characterized in that, include: Acquisition module: used to acquire 3D seismic data and raw well logging data, and to preprocess the raw well logging data to obtain calibrated well logging data; Interpretation module: used to perform well logging interpretation and structural interpretation based on the corrected well logging data and the three-dimensional seismic data, respectively, to obtain interpretation results; The partitioning module is used to transform the structural interpretation results in the interpretation results into a three-dimensional geological structural model. It uses a corner grid system to partition the study area into a structural grid model. Modeling module: used to perform attribute modeling on the constructed mesh model based on variogram analysis, and to predict porosity and saturation in three-dimensional space using the Kriging interpolation algorithm to obtain a three-phase coexistence region model.

Images