CN113009561A - DEM model-based seismic wave velocity parameter determination method, device and equipment - Google Patents
DEM model-based seismic wave velocity parameter determination method, device and equipment Download PDFInfo
- Publication number
- CN113009561A CN113009561A CN202110313216.3A CN202110313216A CN113009561A CN 113009561 A CN113009561 A CN 113009561A CN 202110313216 A CN202110313216 A CN 202110313216A CN 113009561 A CN113009561 A CN 113009561A
- Authority
- CN
- China
- Prior art keywords
- modulus
- pressure
- rock sample
- fracture density
- shear modulus
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 39
- 239000011435 rock Substances 0.000 claims abstract description 135
- 239000011148 porous material Substances 0.000 claims abstract description 111
- 238000009826 distribution Methods 0.000 claims abstract description 52
- 238000001228 spectrum Methods 0.000 claims abstract description 45
- 239000006185 dispersion Substances 0.000 claims abstract description 44
- 230000001186 cumulative effect Effects 0.000 claims abstract description 33
- 238000004364 calculation method Methods 0.000 claims abstract description 27
- 208000010392 Bone Fractures Diseases 0.000 claims description 115
- 206010017076 Fracture Diseases 0.000 claims description 115
- 208000013201 Stress fracture Diseases 0.000 claims description 33
- 229910052500 inorganic mineral Inorganic materials 0.000 claims description 29
- 239000011707 mineral Substances 0.000 claims description 29
- 230000006870 function Effects 0.000 claims description 22
- 239000011159 matrix material Substances 0.000 claims description 16
- 239000012530 fluid Substances 0.000 claims description 14
- 239000007787 solid Substances 0.000 claims description 11
- 238000005259 measurement Methods 0.000 claims description 8
- 230000008859 change Effects 0.000 claims description 6
- 150000001875 compounds Chemical class 0.000 claims description 6
- 238000004590 computer program Methods 0.000 claims description 6
- 238000005406 washing Methods 0.000 claims description 5
- 238000001035 drying Methods 0.000 claims description 4
- 150000003839 salts Chemical class 0.000 claims description 4
- 238000002441 X-ray diffraction Methods 0.000 claims description 3
- 239000000470 constituent Substances 0.000 claims description 2
- 229920000642 polymer Polymers 0.000 claims description 2
- 230000002194 synthesizing effect Effects 0.000 claims 1
- 238000011161 development Methods 0.000 abstract description 6
- 238000004519 manufacturing process Methods 0.000 abstract description 3
- 230000009286 beneficial effect Effects 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 12
- 238000002591 computed tomography Methods 0.000 description 9
- 230000008569 process Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 230000018109 developmental process Effects 0.000 description 5
- 235000019738 Limestone Nutrition 0.000 description 3
- 239000006028 limestone Substances 0.000 description 3
- 230000000704 physical effect Effects 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 238000009530 blood pressure measurement Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 229920006395 saturated elastomer Polymers 0.000 description 2
- 238000003860 storage Methods 0.000 description 2
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000005266 casting Methods 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000005055 memory storage Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
The embodiment of the specification provides a method, a device and equipment for determining seismic wave velocity parameters based on a DEM model. The method comprises the following steps: acquiring an actually measured bulk modulus and an actually measured shear modulus of a rock sample under at least two different pressures; calculating a high-pressure modulus by using the measured bulk modulus and the measured shear modulus; determining an accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on a modulus-fracture density relation under a DEM model; calculating microcrack porosity distribution spectra at the at least two different pressures from the cumulative fracture density; determining a seismic wave velocity parameter corresponding to a rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation. The method realizes accurate calculation of the seismic wave velocity dispersion and attenuation in geology dominated by the micro-pore structure, and is beneficial to production and development in practical application.
Description
Technical Field
The embodiment of the specification relates to the technical field of geological exploration and development, in particular to a method, a device and equipment for determining seismic wave velocity parameters based on a DEM model.
Background
In the field of oil and gas exploration and development, the physical properties of rocks and other seismic response parameters can be determined by obtaining elastic wave parameters of reservoir rocks, wherein the frequency dispersion characteristic and the attenuation characteristic of the velocity lay a key theoretical foundation in the aspects of reservoir detection, fluid identification and the like in a frequency domain, and the method can be used for solving the calibration problem when geophysical data of different frequency bands are jointly applied.
Velocity dispersion generally refers to the phenomenon that the velocity changes although the frequency changes when seismic waves propagate in an actual stratum, and the corresponding characteristic that the amplitude is weakened along with the increase of the distance is usually accompanied, namely the attenuation of the seismic waves. Fluid flow occurring on the pore scale is called jet flow when seismic waves induce velocity dispersion and attenuation over a wide frequency range and spatial scale. The jet flow is generally caused by a local pressure gradient formed by the difference in stiffness of adjacent pores, also referred to as local flow. Specifically, as the seismic wave passes through the rock, the "soft regions" in the rock (e.g., pore throats, microfractures, etc.) close, forcing the internal fluid into the "hard regions" (e.g., large aspect ratio pores), creating jets. The velocity dispersion and attenuation caused by the jet flow have a close relation with the rock pore structure parameters, especially the pore aspect ratio.
In the prior art, when velocity dispersion and attenuation parameters are obtained, generally, the velocity dispersion and attenuation of high-pressure saturated rocks in an ultrasonic frequency band can only be well described, but due to neglecting the influence of a micro-pore structure, the prediction under low effective pressure has a poor result. Generally, for interconnected pore systems, pores with smaller aspect ratios are closed earlier, and internal fluid is squeezed into adjacent pores with larger aspect ratios after closing, so that local fluid flow is generated among various pores, and the existing scheme cannot accurately and completely describe the velocity dispersion and attenuation dominated by the micro-pore structure. Therefore, a method for accurately and completely solving the seismic wave velocity dispersion and attenuation in geology dominated by a micro-pore structure is needed.
Disclosure of Invention
The embodiment of the specification aims to provide a method, a device and equipment for determining seismic wave velocity parameters based on a DEM (digital elevation model) model, so as to solve the problem of accurately solving the seismic wave velocity dispersion and attenuation.
In order to solve the technical problem, an embodiment of the present specification provides a method for determining a seismic wave velocity parameter based on a DEM model, including: acquiring an actually measured bulk modulus and an actually measured shear modulus of a rock sample under at least two different pressures; calculating a high-pressure modulus by using the measured bulk modulus and the measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores comprise pores that are not compressible under confining pressure; determining an accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on a modulus-fracture density relation under a DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample; calculating microcrack porosity distribution spectra at the at least two different pressures from the cumulative fracture density; determining a seismic wave velocity parameter corresponding to a rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
The embodiment of the present specification further provides a seismic wave velocity parameter determining device based on a DEM model, including: the modulus acquisition module is used for acquiring the actually measured bulk modulus and the actually measured shear modulus of the rock sample under at least two different pressures; the high-pressure modulus calculation module is used for calculating the high-pressure modulus by utilizing the actually measured bulk modulus and the actually measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores comprise pores that are not compressible under confining pressure; the accumulated fracture density determining module is used for determining the accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on the relation between the modulus and the fracture density under the DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample; the microcrack porosity distribution spectrum calculation module is used for calculating the microcrack porosity distribution spectrum under the at least two different pressures according to the accumulated fracture density; the seismic wave velocity parameter determining module is used for determining a seismic wave velocity parameter corresponding to the rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
The embodiment of the specification also provides a device for determining the seismic wave velocity parameters based on the DEM model, which comprises a memory and a processor; the memory to store computer program instructions; the processor to execute the computer program instructions to implement the steps of: acquiring an actually measured bulk modulus and an actually measured shear modulus of a rock sample under at least two different pressures; calculating a high-pressure modulus by using the measured bulk modulus and the measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores comprise pores that are not compressible under confining pressure; determining an accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on a modulus-fracture density relation under a DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample; calculating microcrack porosity distribution spectra at the at least two different pressures from the cumulative fracture density; determining a seismic wave velocity parameter corresponding to a rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
According to the technical scheme provided by the embodiment of the specification, firstly, the measured bulk modulus and the measured shear modulus of the rock sample under different pressures are obtained, the high-pressure modulus under high pressure is calculated according to the measured bulk modulus and the measured shear modulus, then the cumulative fracture density and the micro fracture porosity distribution spectrum are sequentially obtained, and the calculation of seismic wave velocity parameters including velocity dispersion and attenuation is realized according to the micro fracture porosity distribution spectrum. By the method, in the process of obtaining the seismic wave velocity parameter, the influence of the micro-cracks in the rock on the jet flow of the seismic wave is considered, the calculation of the seismic wave velocity parameter is further realized based on the aspect ratio and the porosity of the micro-pores, the accurate calculation of the seismic wave velocity dispersion and attenuation in geology dominated by the micro-pore structure is realized, and the production and development in practical application are facilitated.
Drawings
In order to more clearly illustrate the embodiments of the present specification or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments described in the specification, and other drawings can be obtained by those skilled in the art without creative efforts.
FIG. 1 is a flow chart of a method for determining seismic velocity parameters based on a DEM model according to an embodiment of the present disclosure;
figure 2 is a schematic illustration of a sandstone core sample in accordance with embodiments of the present disclosure;
FIG. 3A is a schematic view of a CT scan slice of a sandstone core in accordance with embodiments of the present disclosure;
fig. 3B is a schematic perspective view of a CT scan of a sandstone core according to an embodiment of the present disclosure;
FIG. 4 is a schematic diagram of a dry rock pressure change ultrasonic measurement of longitudinal and transverse wave velocity according to an embodiment of the present disclosure;
FIG. 5 is a graphical representation of modulus fit data in an embodiment of the present disclosure;
FIG. 6 is a schematic illustration of a cumulative fracture density fit, in accordance with an embodiment of the present disclosure;
FIG. 7 is a graph of intersection of microcrack aspect ratio and porosity at different pressures in accordance with an embodiment of the disclosure;
FIG. 8A is a schematic illustration of the bulk modulus of a "wet skeleton" according to an embodiment of the present disclosure;
FIG. 8B is a schematic illustration of the shear modulus of a "wet skeleton" according to an embodiment of the present disclosure;
FIG. 9A is a schematic diagram of longitudinal wave velocity dispersion according to an embodiment of the present disclosure;
FIG. 9B is a schematic diagram of a transverse wave velocity dispersion according to an embodiment of the present disclosure;
FIG. 9C is a schematic diagram of a longitudinal wave attenuation coefficient according to an embodiment of the present disclosure;
FIG. 9D is a schematic diagram of attenuation coefficient of transverse waves according to an embodiment of the present disclosure;
FIG. 10 is a schematic view of a limestone core sample according to an embodiment of the present disclosure;
FIG. 11 is a schematic illustration of a sheet of a limestone core casting according to an embodiment of the present disclosure;
FIG. 12 is a sectional view of a CT scan of a limestone core according to an embodiment of the present disclosure;
FIG. 13 is a schematic diagram of a dry rock pressure change ultrasonic measurement of longitudinal and transverse wave velocity in accordance with an embodiment of the present disclosure;
FIG. 14 is a graphical representation of modulus fit data in accordance with an embodiment of the present disclosure;
FIG. 15 is a schematic illustration of a cumulative fracture density fit, in accordance with an embodiment of the present disclosure;
FIG. 16 is a graph of intersection of microcrack aspect ratio and porosity at different pressures in accordance with an embodiment of the disclosure;
FIG. 17A is a schematic illustration of the bulk modulus of a "wet skeleton" according to an embodiment of the present disclosure;
FIG. 17B is a schematic illustration of the shear modulus of a "wet skeleton" according to an embodiment of the present disclosure;
FIG. 18A is a schematic diagram of longitudinal wave velocity dispersion according to an embodiment of the present disclosure;
FIG. 18B is a schematic diagram of a transverse wave velocity dispersion according to an embodiment of the present disclosure;
FIG. 18C is a schematic diagram of a longitudinal wave attenuation coefficient according to an embodiment of the present disclosure;
FIG. 18D is a schematic diagram of attenuation coefficient of transverse waves according to an embodiment of the present disclosure;
FIG. 19 is a block diagram of a seismic velocity parameter determining apparatus based on a DEM model according to an embodiment of the present disclosure;
fig. 20 is a block diagram of a seismic velocity parameter determination apparatus based on a DEM model according to an embodiment of the present disclosure.
Detailed Description
The technical solutions in the embodiments of the present disclosure will be clearly and completely described below with reference to the drawings in the embodiments of the present disclosure, and it is obvious that the described embodiments are only a part of the embodiments of the present disclosure, and not all of the embodiments. All other embodiments obtained by a person of ordinary skill in the art based on the embodiments in the present specification without any creative effort shall fall within the protection scope of the present specification.
In order to solve the technical problem, the embodiment of the specification provides a seismic wave velocity parameter determination method based on a DEM model. The execution subject of the DEM model-based seismic wave velocity parameter determination method is DEM model-based seismic wave velocity parameter determination equipment. As shown in fig. 1, the method for determining the seismic wave velocity parameter based on the DEM model may specifically include the following steps.
S110: and acquiring the measured bulk modulus and the measured shear modulus of the rock sample under at least two different pressures.
The rock sample may be a rock sample collected for a target reservoir or a laboratory prepared rock sample simulating the rock results of the corresponding formation. There may be no limitation on the lithology of the rock sample. In performing the experiment, in order to facilitate the handling of the rock sample, the rock sample may be processed into a specific shape, such as a cylindrical shape, so as to facilitate the subsequent operations.
The bulk modulus is a physical quantity reflecting the relationship between the bulk strain and the average stress of an object, and is used to reflect the macroscopic properties of a material. Shear modulus is the ratio of shear stress to strain and is used to characterize a material's ability to resist shear strain.
The measured bulk modulus and the measured shear modulus may be the moduli obtained after actual measurements are made on the rock sample. Specifically, the measured bulk modulus and the measured shear modulus corresponding to different pressures may be indirectly obtained through acquiring the longitudinal and transverse wave velocities corresponding to the rock sample under different pressures.
Specifically, the longitudinal wave velocity V can be obtained by first performing variable pressure measurement on the longitudinal wave velocity and the transverse wave velocity of the rock sample to which the variable pressure measurement is appliedP(P) and transverse wave velocity VS(P) thereafter, first using the formula
Calculating the measured bulk modulus, where KD_meas(P) is the measured bulk modulus at pressure P, ρ is the density of the rock sample, VP(P) is the longitudinal wave velocity under pressure P, VS(P) is the velocity of the transverse wave under pressure P, and the formula G is reusedD_meas(P)＝ρVS(P)2Calculating the measured shear modulus, wherein GD_meas(P) is the measured shear modulus at pressure P.
In some embodiments, the rock sample may be pre-processed before performing the respective experimental operation on the rock sample to enable the rock sample to better comply with subsequent experimental criteria. Specifically, the pretreatment may include at least one of a wash oil treatment, a salt washing treatment, and a drying treatment. The oil washing treatment is used for washing out oil in the rock sample, the salt washing treatment is used for removing salt in the rock sample, and the drying treatment is used for removing moisture in the rock sample. Through the pretreatment, a dry rock sample can be obtained, so that the subsequent experimental steps can be better executed.
In some embodiments, prior to obtaining the modulus, pore structure parameters corresponding to the rock sample may also be obtained. The pore structure parameter may include at least one of density, porosity, estimated bulk modulus, and estimated shear modulus.
When the estimated bulk modulus and the estimated shear modulus are obtained, XRD full-core analysis can be performed on the rock sample to obtain mineral components of the rock sample. Accordingly, each can be obtained during the analysisVolume fraction, bulk modulus, and shear modulus of the individual mineral components. Then, the formula can be utilized
Calculating the estimated bulk modulus, where KmFor estimating the bulk modulus, N is the amount of mineral constituent, fiIs the volume fraction of the ith mineral component, KiIs the bulk modulus of the ith mineral component, and using the formula
Calculation of the estimated shear modulus, wherein GmTo estimate the shear modulus, GiIs the shear modulus of the ith mineral component.
The estimated bulk modulus and the estimated shear modulus are results obtained after the bulk modulus and the shear modulus are estimated under the condition of combining mineral composition of a rock sample, and the calculation of the wet skeleton modulus can be realized in the subsequent steps by combining a microcrack porosity distribution spectrum.
S120: calculating a high-pressure modulus by using the measured bulk modulus and the measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores include pores that are not compressible under confining pressure.
The measured bulk modulus and the measured shear modulus may be used to calculate the high pressure modulus. The high pressure modulus may be equivalent to the modulus of a rock sample consisting of a solid mineral matrix and hard pores, i.e. the modulus corresponding to a rock sample at extreme pressures, wherein hard pores may refer to pores that have been incompressible by confining pressure with increasing pressure. The high pressure modulus may include a high pressure bulk modulus and a high pressure shear modulus.
Since the measured bulk modulus and shear modulus corresponding to the rock sample have been acquired in step S110, the high-pressure modulus can be derived based on a relation constructed under different pressure variations of modulus.
In some embodiments, the bulk modulus relationship and the shear modulus relationship may be constructed separately; the volume modulus relational expression comprises a relational expression constructed on the basis of theoretical volume modulus, pressure coefficient, high-pressure volume modulus and zero-pressure volume modulus; the shear modulus relational expression comprises a relational expression constructed on the basis of theoretical shear modulus, pressure coefficient, high-pressure shear modulus and zero-pressure shear modulus. And then combining the volume modulus relational expression and the shear modulus relational expression, and fitting by utilizing the actually measured volume modulus and the actually measured shear modulus to obtain the high-pressure modulus.
Specifically, the bulk modulus relationship may be constructed in general
In the formula, KD(P) is the theoretical bulk modulus of the rock sample under pressure P,
the bulk modulus at zero pressure is given,
in order to have a high bulk modulus at high pressure,
is the pressure coefficient. The shear modulus relationship constructed may be
In the formula, GD(P) is the theoretical shear modulus of the rock sample at pressure P,
is a shear modulus at zero pressure and a shear modulus,
high pressure shear modulus.
And under the condition that the measured bulk modulus and the measured shear modulus are obtained, combining the theoretical bulk modulus and the theoretical shear modulus in the relational expression, the formula OF can be constructed1＝∑[(KD_meas(P)-KD(P))2+(GD_meas(P)-GD(P))2]OF in the formula1Is a first objective function, KD_meas(P) is the measured bulk modulus under pressure P, GD_meas(P) is the measured shear modulus at pressure P. And performing nonlinear least square fitting on the above formula to minimize the first objective function, namely to make the measured bulk modulus and the measured shear modulus closest to each other, so as to obtain a theoretical bulk modulus and a theoretical shear modulus corresponding to the specific pressure. Then, the corresponding high-pressure bulk modulus and high-pressure shear modulus can be obtained according to the relational expression.
S130: determining an accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on a modulus-fracture density relation under a DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample.
The modulus and fracture density relations under the DEM model comprise the relations between Poisson's ratio and fracture density, the relations between bulk modulus and fracture density, and the relations between shear modulus and fracture density. Specifically, the relationship between the Poisson's ratio and the fracture density comprises
In the formula (I), the compound is shown in the specification,
equivalent Poisson's ratio, v, estimated for DEM modelb＝(3Kb-2Gb)/(6Kb+2Gb)，KbEquivalent modulus of elasticity, G, to the high-pressure bulk modulusbEquivalent elastic modulus of high-pressure shear modulus, gamma is fracture density; the bulk modulus to fracture density relationship comprises
In the formula (I), the compound is shown in the specification,
estimating the volume modulus of the equivalent medium for the DEM model; the shear modulus and fracture densityThe relationship includes
In the formula (I), the compound is shown in the specification,
and estimating the obtained equivalent medium shear modulus for the DEM model.
Because the DEM model already defines the relationship between modulus and fracture density. Under the condition that the measured bulk modulus and the measured shear modulus have been obtained in step S110, the corresponding fracture density can be obtained by fitting based on the equivalent medium bulk modulus and the equivalent medium shear modulus estimated based on the DEM model in the above relational expression. For example, a relational expression can be constructed
OF in the formula2In order to be the second objective function,
corresponding to pressure interval P estimated for DEM modeliEquivalent medium bulk modulus of (2), KD_meas(Pi) To correspond to the pressure interval PiThe measured bulk modulus of the polymer (a),
corresponding to pressure interval P estimated for DEM modeliEquivalent medium shear modulus of (1), GD_meas(Pi) To correspond to the pressure interval PiMeasured shear modulus of.
Under the condition that the measured bulk modulus and the measured shear modulus can change along with the change of the pressure, the fracture density under different pressures can be respectively determined, and then the accumulated fracture density is obtained through accumulation. In practical implementation, after the highest pressure corresponding to the high-pressure modulus is determined, at least two pressure intervals are determined based on the highest pressure and the zero pressure, and the density of the microcracks in each pressure interval is obtained respectively. Preferably, for ease of calculation, the pressure interval may be set based on the pressure to which the measured shear modulus and the measured bulk modulus are measured.
In particular, p at the highest voltageNIn the case of (3), p may be first introducedN→pN-1This pressure interval will correspond to the pressure pNAnd pN-1The elastic modulus is used as the elastic parameter of the rock background matrix and the rock equivalent medium, namely the high-pressure modulus
Substitute for [ Kb，Gb]By means of pressure pN-1Measured shear modulus and measured bulk modulus [ K ] measuredD_meas(P)，GD_meas(P)]Based on the relation, the volume modulus and the shear modulus of the equivalent medium in the pressure interval are obtained through fitting by controlling the second objective function to be the minimum value
And then obtaining a pressure interval p based on the relation between the modulus and the fracture density under the DEM modelN→pN-1Lower microcrack density f (P)N)。
Thereafter, other pressure intervals are possible, e.g. pi→pi-1Taking the density of the microcracks calculated in the previous step and the hard matrix pores as background media, and taking the pressure piTheoretical bulk modulus and theoretical shear modulus [ K ] ofD(P)，GD(P)]Substitution of [ Kb，Gb]In combination with a pressure pi-1Measured bulk modulus and measured shear modulus [ K ] atD_meas(P)，GD_meas(P)]Based on the relation, the volume modulus and the shear modulus of the equivalent medium in the pressure interval are obtained through fitting by controlling the second objective function to be the minimum value
Finally obtaining a pressure interval pi→pi-1Lower microcrack density gamma (P)i)。
Through the steps, the steps are carried out in sequenceThe micro-crack density corresponding to each pressure interval is obtained, and finally, a formula is utilized
The cumulative fracture density is calculated, where,
is a pressure PkThe sum of the fracture densities of all the lower open-hole microcracks, i.e. the cumulative fracture density, N is the total number of pressure intervals, Γ (P)i) Is a pressure interval pi→pi-1Lower microcrack density.
S140: and calculating the microcrack porosity distribution spectrum under the at least two different pressures according to the accumulated fracture density.
After the cumulative fracture density is obtained, a microcrack porosity profile at the different pressures can be calculated from the cumulative fracture density. The microfracture porosity distribution spectrum is used to represent porosity data of microfractures in the rock sample. Because the velocity dispersion and attenuation based on the jet flow have strong correlation with the microfracture and the porosity corresponding to the microfracture, after the microfracture porosity of the rock sample is obtained, the velocity dispersion and attenuation properties can be effectively obtained in the subsequent steps.
Before the microcrack porosity distribution spectrum is obtained, the microcrack aspect ratio distribution spectrum under different pressures can be obtained, and then the microcrack porosity distribution spectrum is determined according to the microcrack aspect ratio distribution spectrum.
When obtaining the microcrack porosity distribution spectrum, fracture density curves corresponding to the at least two different pressures may be obtained by first combining the fracture density fit. The fracture density curve can be a curve corresponding to a theoretical value of fracture density changing along with pressure. Specifically, the fracture density curve can be formulated
Is expressed, where ε (p) is the fracture density corresponding to pressure p, which is a theoretical formula fit0In order to obtain the initial fracture density,
is the undetermined coefficient.
Since the fracture density corresponding to each pressure interval has been calculated in the above-described step S130
Thus, a formula can be constructed
OF in the formula3Is a third objective function, ε (P)i) Corresponding to a pressure interval P fitted to a theoretical formulaiThe crack density of (a) is high,
corresponding to the pressure interval P obtained for calculationiThe fracture density of (a). By fitting the formula to minimize the third objective function, from the fitted epsilon (P)i) And determining undetermined coefficients in the previous relational expression so as to determine the fracture density curve.
Then, an initial pore aspect ratio corresponding to each microcrack can be calculated according to the fracture density curve; the initial pore aspect ratio comprises the microcrack aspect ratio at zero pressure at which a particular pressure point closes. In particular, a formula may be utilized
Calculating the initial pore aspect ratio, wherein0(Pi) Is a pressure PiInitial pore aspect ratio of lower microcracks, KD(P) is the theoretical bulk modulus at pressure P,
high pressure bulk modulus. Through the calculation process, the vector alpha of the initial aspect ratio of all the microcracks under zero pressure can be obtained0(P)＝[α0(P1)，α0(P2)，…，α0(PN)]。
Correspondingly, the microcracked porosity distribution spectrum is obtained by integrating the initial pore aspect ratios, specifically, the formula α (P) ═ α can be used0(P)-αp(P) calculating a microcracked porosity profile, wherein α (P) is the microcracked porosity profile and αp(P)＝[αp，αp，…，αp]1*N. Through the calculation process, the micro-fracture aspect ratio distribution spectrum under any pressure P is obtained.
After the micro fracture aspect ratio is obtained, calculating by utilizing the micro fracture pore aspect ratio distribution spectrum to obtain the micro fracture porosity distribution spectrum, and specifically determining the nth group of micro fracture porosity at the confining pressure of 0
Is calculated by the formula
Calculating the microcrack porosity at 0 confining pressure, wherein,
porosity of the nth set of microcracks at 0 confining pressure. Then, the formula can be utilized
Calculating the microcrack porosity distribution spectrum with an aspect ratio of alpha under an effective pressure P, wherein phi isnAnd (alpha, P) is the n < th > set of microcracked porosities with aspect ratio alpha under effective pressure P. Thereby completing the calculation of the microcrack porosity distribution spectrum.
S150: determining a seismic wave velocity parameter corresponding to a rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
After the micro-fracture porosity distribution spectrum is obtained, seismic wave velocity parameters, mainly velocity dispersion and attenuation, can be determined by combining the pore structure parameters of the rock sample. For the description and calculation method of the pore structure parameter, reference may be made to the description in step S110, and details are not repeated here.
In some embodiments, the wet skeleton modulus may be obtained based on the extended Gurevich jet model, and then the velocity dispersion and attenuation may be calculated based on Biot pore elasticity theory in combination with pore structure parameters.
Since the velocity dispersion and the magnitude of the attenuation caused by the jet flow are closely related to the rock pore structure parameters, especially the pore aspect ratio, significant attenuation is caused only by cracks or microcracks with a small aspect ratio, while pores with a large aspect ratio (such as isodiametric pores) contribute little to the attenuation. In terms of simulation of the microfracture jet effect, a common approach is to assume that the rock pore space consists of two parts: hard pores dominated by volume content and microcracks or soft pores that are very sensitive to pressure changes. Based on this assumption, a "wet skeleton" model is proposed, i.e. a rock with soft pores saturated with fluid and hard pores empty, to quantify the elastic response of the jet under very high frequency conditions.
The wet skeleton modulus is the modulus corresponding to the "wet skeleton" model, and specifically, the wet skeleton modulus may further include a wet skeleton bulk modulus and a wet skeleton shear modulus. The extended Gurevich jet flow model contains a calculation formula corresponding to the wet skeleton bulk modulus
In the formula, Kwf(P, w) is the wet skeleton bulk modulus,
is the high-pressure bulk modulus, N is the number of pressure intervals, phin(P) is a microcracked porosity profile corresponding to pressure P,
wherein, KD(P) is the theoretical bulk modulus, α, corresponding to the pressure Pn(P) is the microcrack aspect ratio profile corresponding to pressure P,
to correspond to the cumulative fracture density of the pressure P,
J1(xi) is a first order Bessel function,
where ω is the phase, a is the aperture parameter, ρf1Is fluid density, η is pore fluid viscosity, J0(xi) is a zero order Bessel function, KfIs the bulk modulus of fluid, KmIs the mineral bulk modulus; also includes a formula corresponding to the shear modulus of the wet skeleton
In the formula, Gwf(P, w) is the wet skeleton shear modulus, GD(P) is the theoretical shear modulus corresponding to the pressure P. Through the formula and by combining various parameters obtained by calculation, the calculation of the wet framework volume modulus and the wet framework shear modulus can be realized.
After the wet skeleton modulus is obtained, calculations for velocity dispersion and attenuation can be done based on Biot pore elasticity theory. Specifically, the formula V ═ 1/Re (1/V) can be usedc) Calculating velocity dispersion, where V is phase velocity and VcFor complex velocities calculated by the Biot model, including shear complex velocities
And longitudinal wave complex velocity
Complex velocity of transverse wave
Where Δ ═ P ρ22+Rρ11-2Qρ12，
M＝[(τ-φ)/Km+φ/Kf]-1，τ＝1-Kwf/Km，KwfIs the wet skeleton bulk modulus, phi is the porosity, GwfIs the wet skeleton shear modulus, ρ22＝aφρfWhere a is the tortuosity and the pores are randomly distributed, a may be 3, ρfFor pore fluid density, R-M φ2，ρ11＝(1-φ)ρm-ρ12，ρmThe modulus of the rock matrix, i.e. the modulus, rho, taking into account only the mineral content, without taking into account the porosity and spatial structure12＝(1-a)φρf，Q＝M(τ-φ)φ，τ＝1-Kwf/KmVelocity of complex longitudinal wave
Where ρ is ρD+φρf，ρDDensity of the dry rock sample. After the phase velocity is obtained through calculation, the velocity dispersion can be obtained according to the change condition of the phase velocity. The formula Q can also be utilized-1＝Im(Vc2)/Re(Vc2) Calculating the reciprocal of the quality factor, wherein Q-1Is the reciprocal of the quality factor. The inverse quality factor is used to reflect the speed attenuation situation. The porosity and density obtained by the above formula can be obtained by referring to the description in step S110, and will not be described herein again.
The above method is explained below using a specific scenario example. In this scenario example, a piece of fine grained feldspathic rock debris sandstone was utilized as the rock sample, the appearance of which is shown in fig. 2. The internal pore structure distribution of the rock sample can be obtained by performing CT scanning on the rock sample, the slice result of the CT scanning is shown in figure 3A, and the stereogram result of the CT scanning is shown in figure 3B. The CT scanning result shows that the sandstone rock sample has very uniform pore structure distribution, is basically a round hole and almost has no crack development. And then, drying the sandstone rock sample, acquiring the porosity of the sandstone rock sample by using a pore-permeability tester, and acquiring the density of the sandstone rock sample by using an electronic balance.
Specific physical properties are shown in table 1 below.
| Porosity (%) | Density (g/cc) |
| 21.536 | 2.019615693 |
TABLE 1
Accordingly, other pore structure parameters can be obtained by XRD mineral analysis of the rock sample. Specific analysis results can be shown in table 2 below.
TABLE 2
Based on the pore structure parameters obtained from the above XRD analysis, the equivalent elastic modulus can be calculated by VRH (Voigt-reus-Hill) boundary averaging theory, in combination with the reference amounts for mineral modulus listed in the petrophysical handbook. The parameters relating specifically to the modulus of the rock matrix are shown in table 3 below.
TABLE 3
In the process of obtaining the actually measured bulk modulus and the actually measured shear modulus, the utilized longitudinal and transverse wave speeds can be completed under a dry condition by utilizing a high-temperature high-pressure ultrasonic measurement system. Specifically, the data corresponding to the velocity of the longitudinal and transverse waves may be as shown in FIG. 4. And fitting by using the volume modulus and the shear modulus calculated according to the actually measured longitudinal and transverse wave velocity verse to obtain modulus fitting data shown in figure 5. Accordingly, the high-pressure modulus obtained under extremely high pressure obtained from fitting the data can be shown in table 4 below.
TABLE 4
Substituting the high-pressure modulus and the elastic modulus obtained by fitting under each pressure into a formula corresponding to the single-hole DEM model, and inverting the accumulated fracture density under each pressure according to the target function to obtain the data of the accumulated fracture density obtained by inversion corresponding to the origin in the figure 6.
Fitting the cumulative fracture density epsilon under zero pressure according to the cumulative fracture density data obtained by inversion0And pressure parameter
The specific fitting function is shown in table 5 below, and the fitting results are shown in the fitted curve in fig. 6.
TABLE 5
The obtained accumulated fracture density can be used for obtaining the aspect ratio distribution of different microcracks under various pressures, and further obtaining the porosity distribution condition corresponding to the microcracks with various aspect ratios under different confining pressures. As shown in fig. 7, a graph of the microcrack aspect ratio versus the corresponding porosity at different pressures is shown.
After the aspect ratio and porosity of each type of soft pore were obtained, the "wet skeleton" modulus in the extended Gurevich-based jet model was calculated from the obtained parameters. As shown in fig. 8A and 8B, a wet skeleton bulk modulus and a wet skeleton shear modulus are respectively illustrated.
And after the wet skeleton modulus is obtained, calculating to obtain the velocity dispersion and the attenuation based on a formula corresponding to a Biot pore elasticity theory. As shown in fig. 9A and 9B, the calculated longitudinal wave velocity dispersion and transverse wave velocity dispersion; fig. 9C and 9D show the calculated longitudinal wave attenuation coefficient and the calculated transverse wave attenuation coefficient.
This is illustrated below using another scenario example. In this scenario example, tight carbonate rock was used as the rock sample, the appearance of which is shown in fig. 10. The internal pore structure is obtained by CT scanning, and the core cast body slice shown in figure 11 and the CT scanning segmentation perspective corresponding to figure 12 are obtained. According to the scanning result, the internal pore structure distribution is uneven, and the crack growth is obvious. Similarly, the porosity parameter is obtained from a porosimeter, the density is measured by an electronic balance, and the measurement is carried out under dry conditions. The measured physical property parameters are shown in Table 6 below.
| Porosity (%) | Density (g/cc) |
| 1.705083 | 2.66311 |
TABLE 6
Other pore structure parameters can be obtained by XRD mineral analysis of the rock sample. Specific analysis results can be shown in table 7 below.
TABLE 7
Based on the pore structure parameters obtained from the above XRD analysis, the equivalent elastic modulus can be calculated by VRH (Voigt-reus-Hill) boundary averaging theory, in combination with the reference amounts for mineral modulus listed in the petrophysical handbook. Parameters relating specifically to the modulus of the rock matrix are shown in table 8 below.
TABLE 8
In the process of obtaining the actually measured bulk modulus and the actually measured shear modulus, the utilized longitudinal and transverse wave speeds can be completed under a dry condition by utilizing a high-temperature high-pressure ultrasonic measurement system. Specifically, the data corresponding to the velocity of the longitudinal and transverse waves may be as shown in fig. 13. The bulk modulus and shear modulus calculated from the measured longitudinal and lateral wave velocity verses are used for fitting to obtain modulus fitting data as shown in fig. 14. Accordingly, the high-pressure modulus obtained under the extremely high-pressure condition obtained from the fitting data can be shown in table 9 below.
TABLE 9
Substituting the high-pressure modulus and the elastic modulus obtained by fitting under each pressure into a formula corresponding to the single-hole DEM model, and inverting the accumulated fracture density under each pressure according to the target function to obtain the data of the accumulated fracture density obtained by inversion corresponding to the origin in the figure 15.
Fitting the cumulative fracture density epsilon under zero pressure according to the cumulative fracture density data obtained by inversion0And pressure parameter
The specific fitting function is shown in table 10 below, and the fitting results are shown in the fitted curve in fig. 15.
The obtained accumulated fracture density can be used for obtaining the aspect ratio distribution of different microcracks under various pressures, and further obtaining the porosity distribution condition corresponding to the microcracks with various aspect ratios under different confining pressures. As shown in fig. 16, a graph of the microcrack aspect ratio versus the corresponding porosity at different pressures is shown.
After the aspect ratio and porosity of each type of soft pore were obtained, the "wet skeleton" modulus in the extended Gurevich-based jet model was calculated from the obtained parameters. As shown in fig. 17A and 17B, a wet skeleton bulk modulus and a wet skeleton shear modulus are respectively illustrated.
And after the wet skeleton modulus is obtained, calculating to obtain the velocity dispersion and the attenuation based on a formula corresponding to a Biot pore elasticity theory. As shown in fig. 18A and 18B, the calculated longitudinal wave velocity dispersion and transverse wave velocity dispersion are obtained; fig. 18C and 18D show the calculated longitudinal wave attenuation coefficient and the calculated transverse wave attenuation coefficient.
Through the introduction of the above embodiment and the scenario example, it can be seen that the method first obtains the measured bulk modulus and the measured shear modulus of the rock sample under different pressures, and after calculating the high-pressure modulus under high pressure according to the measured bulk modulus and the measured shear modulus, sequentially obtains the cumulative fracture density and the microcrack porosity distribution spectrum, thereby obtaining the seismic wave velocity parameters including velocity dispersion and attenuation according to the microcrack porosity distribution spectrum. By the method, in the process of obtaining the seismic wave velocity parameter, the influence of the micro-cracks in the rock on the jet flow of the seismic wave is considered, the calculation of the seismic wave velocity parameter is further realized based on the aspect ratio and the porosity of the micro-pores, the accurate calculation of the seismic wave velocity dispersion and attenuation in geology dominated by the micro-pore structure is realized, and the production and development in practical application are facilitated.
Based on the method for determining the seismic wave velocity parameter based on the DEM model, the specification also provides an embodiment of a device for determining the seismic wave velocity parameter based on the DEM model. As shown in fig. 19, the seismic wave velocity parameter determining apparatus based on the DEM model specifically includes the following modules.
A modulus obtaining module 1910 configured to obtain a measured bulk modulus and a measured shear modulus of the rock sample under at least two different pressures.
A high-pressure modulus calculation module 1920 configured to calculate a high-pressure modulus using the measured bulk modulus and the measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores include pores that are not compressible under confining pressure.
An accumulated fracture density determination module 1930, configured to determine, based on a relationship between modulus and fracture density under the DEM model, an accumulated fracture density corresponding to the rock sample through the high-pressure modulus; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample.
And a microfracture porosity distribution spectrum calculation module 1940, configured to calculate microfracture porosity distribution spectra at the at least two different pressures according to the cumulative fracture density.
A seismic velocity parameter determination module 1950, configured to determine a seismic velocity parameter corresponding to a rock sample in combination with a pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
Based on the seismic wave velocity parameter determination method based on the DEM model, the embodiment of the specification further provides seismic wave velocity parameter determination equipment based on the DEM model. As shown in fig. 20, the DEM model-based seismic wave velocity parameter determination apparatus includes a memory and a processor.
In this embodiment, the memory may be implemented in any suitable manner. For example, the memory may be a read-only memory, a mechanical hard disk, a solid state disk, a U disk, or the like. The memory may be used to store computer program instructions.
In this embodiment, the processor may be implemented in any suitable manner. For example, the processor may take the form of, for example, a microprocessor or processor and a computer-readable medium that stores computer-readable program code (e.g., software or firmware) executable by the (micro) processor, logic gates, switches, an Application Specific Integrated Circuit (ASIC), a programmable logic controller, an embedded microcontroller, and so forth. The processor may execute the computer program instructions to perform the steps of: acquiring an actually measured bulk modulus and an actually measured shear modulus of a rock sample under at least two different pressures; calculating a high-pressure modulus by using the measured bulk modulus and the measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores comprise pores that are not compressible under confining pressure; determining an accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on a modulus-fracture density relation under a DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample; calculating microcrack porosity distribution spectra at the at least two different pressures from the cumulative fracture density; determining a seismic wave velocity parameter corresponding to a rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
The systems, devices, modules or units illustrated in the above embodiments may be implemented by a computer chip or an entity, or by a product with certain functions. One typical implementation device is a computer. In particular, the computer may be, for example, a personal computer, a laptop computer, a cellular telephone, a camera phone, a smartphone, a personal digital assistant, a media player, a navigation device, an email device, a game console, a tablet computer, a wearable device, or a combination of any of these devices.
From the above description of the embodiments, it is clear to those skilled in the art that the present specification can be implemented by software plus a necessary general hardware platform. Based on such understanding, the technical solutions of the present specification may be essentially or partially implemented in the form of software products, which may be stored in a storage medium, such as ROM/RAM, magnetic disk, optical disk, etc., and include instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute the methods described in the embodiments or some parts of the embodiments of the present specification.
The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. In particular, for the system embodiment, since it is substantially similar to the method embodiment, the description is simple, and for the relevant points, reference may be made to the partial description of the method embodiment.
The description is operational with numerous general purpose or special purpose computing system environments or configurations. For example: personal computers, server computers, hand-held or portable devices, tablet-type devices, multiprocessor systems, microprocessor-based systems, set top boxes, programmable consumer electronics, network PCs, minicomputers, mainframe computers, distributed computing environments that include any of the above systems or devices, and the like.
This description may be described in the general context of computer-executable instructions, such as program modules, being executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. The specification may also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media including memory storage devices.
While the specification has been described with examples, those skilled in the art will appreciate that there are numerous variations and permutations of the specification that do not depart from the spirit of the specification, and it is intended that the appended claims include such variations and modifications that do not depart from the spirit of the specification.
Claims (10)
1. A seismic wave velocity parameter determination method based on a DEM model is characterized by comprising the following steps:
acquiring an actually measured bulk modulus and an actually measured shear modulus of a rock sample under at least two different pressures;
calculating a high-pressure modulus by using the measured bulk modulus and the measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores comprise pores that are not compressible under confining pressure;
determining an accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on a modulus-fracture density relation under a DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample;
calculating microcrack porosity distribution spectra at the at least two different pressures from the cumulative fracture density;
determining a seismic wave velocity parameter corresponding to a rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
2. The method of claim 1, wherein said obtaining a measured bulk modulus and a measured shear modulus of the rock sample at least two different pressures comprises:
acquiring longitudinal wave velocity and transverse wave velocity corresponding to the rock sample under at least two different pressures;
calculating the actual measurement volume modulus and the actual measurement shear modulus of the rock sample along with the pressure change by utilizing the longitudinal wave velocity and the transverse wave velocity; including, among other things, the use of formulas
Calculating the measured bulk modulus, where KD_meas(P) is the measured bulk modulus at pressure P, ρ is the density of the rock sample, VP(P) is the longitudinal wave velocity under pressure P, Vs(P) is the velocity of the transverse wave under pressure P; using formula GD_meas(P)＝ρVs(P)2Calculating the measured shear modulus, wherein GD_meas(P) is the measured shear modulus at pressure P.
3. The method of claim 1, wherein the high pressure modulus comprises a high pressure bulk modulus and a high pressure shear modulus; calculating the high-pressure modulus by using the measured bulk modulus and the measured shear modulus, comprising:
respectively constructing a volume modulus relational expression and a shear modulus relational expression; the bulk modulus relation comprising
In the formula, KD(P) is the theoretical bulk modulus of the rock sample under pressure P,
the bulk modulus at zero pressure is given,
in order to have a high bulk modulus at high pressure,
is the pressure coefficient; the shear modulus relationship comprising
In the formula, GD(P) is the theoretical shear modulus of the rock sample at pressure P,
is a shear modulus at zero pressure and a shear modulus,
high pressure shear modulus;
combining the volume modulus relational expression and the shear modulus relational expression, and fitting by utilizing the actually measured volume modulus and the actually measured shear modulus to obtain a high-pressure modulus; wherein, include: by the pair formula OF1＝∑[(KD_meas(P)-KD(P))2+(GD_meas(P)-GD(P))2]Performing nonlinear least square fitting to obtain high-pressure bulk modulus and high-pressure shear modulus when the first objective function is minimum, wherein OF1Is a first objective function, KD_meas(P) is the measured bulk modulus under pressure P, GD_meas(P) is the measured shear modulus at pressure P.
4. The method of claim 1, wherein determining the cumulative fracture density corresponding to the rock sample from the high pressure modulus based on the modulus versus fracture density relationship under the DEM model comprises:
determining the highest pressure corresponding to the high-pressure modulus;
determining at least two pressure intervals according to the highest pressure and the zero pressure;
by a pair of formulas
Fitting to obtain micro-fracture densities respectively corresponding to different pressure intervals; OF in the formula2In order to be the second objective function,
corresponding to pressure interval P estimated for DEM modeliEquivalent medium bulk modulus of (2), KD_meas(Pi) To correspond to the pressure interval PiThe measured bulk modulus of the polymer (a),
corresponding to pressure interval P estimated for DEM modeliEquivalent medium shear modulus of (1), GD_meas(Pi) To correspond to the pressure interval PiThe measured shear modulus of (a);
and accumulating the density of the micro fractures to obtain the accumulated fracture density.
5. The method of claim 4, wherein the pass pair formula
Fitting to obtain micro-fracture densities respectively corresponding to different pressure intervals, wherein the micro-fracture densities comprise:
minimizing the second objective function by fitting the formula;
obtaining fracture density corresponding to the second objective function as micro-fracture density corresponding to the pressure interval;
the modulus and fracture density relation under the DEM model comprises a Poisson ratio and fracture density relation, a volume modulus and fracture density relation and a shear modulus and fracture density relation;
the Poisson's ratio and fracture density relationship comprises
In the formula (I), the compound is shown in the specification,
equivalent Poisson's ratio, v, estimated for DEM modelb＝(3Kb-2Gb)/(6Kb+2Gb)，KbEquivalent modulus of elasticity, G, to the high-pressure bulk modulusbEquivalent elastic modulus of high-pressure shear modulus, gamma is fracture density;
the bulk modulus to fracture density relationship comprises
In the formula (I), the compound is shown in the specification,
estimating the volume modulus of the equivalent medium for the DEM model;
the shear modulus to fracture density relationship comprises
In the formula (I), the compound is shown in the specification,
and estimating the obtained equivalent medium shear modulus for the DEM model.
6. The method of claim 1, wherein calculating the microcrack porosity profile at the at least two different pressures from the cumulative fracture density comprises:
fitting in conjunction with the cumulative fracture density to obtain fracture density curves corresponding to the at least two different pressures; wherein, include: combination formula
And
fitting to obtain a fracture density curve, wherein the fracture density curve is constructed based on initial fracture density and undetermined coefficient, epsilon (p) is fracture density corresponding to pressure p and is fitted by a theoretical formula0In order to obtain the initial fracture density,
to be undetermined coefficients, OF3Is a third objective function, ε (P)i) Corresponding to a pressure interval P fitted to a theoretical formulaiThe crack density of (a) is high,
corresponding to the pressure interval P obtained for calculationiThe fracture density of (a);
calculating an initial pore aspect ratio corresponding to each microcrack according to the fracture density curve; the initial pore aspect ratio comprises a microfracture aspect ratio at zero pressure at which a particular pressure point closes; wherein, include: using formulas
Calculating the initial pore aspect ratio, wherein0(Pi) Is a pressure PiInitial pore aspect ratio of lower microcracks, KD(P) is the theoretical bulk modulus at pressure P,
high pressure bulk modulus;
synthesizing the aspect ratios of the initial pores to obtain a microcrack pore aspect ratio distribution spectrum; wherein, include: using the formula α (P) ═ α0(P)-αp(P) calculating a microcracked porosity profile, wherein α (P) is the microcracked porosity profile and αp(P)＝[αp，αp，…，αp]1*N；
Calculating by using the microcrack pore aspect ratio distribution spectrum to obtain a microcrack pore distribution spectrum; wherein, include: using formulas
Calculating the microcrack porosity at 0 confining pressure, wherein,
porosity of the nth set of microcracks at a confining pressure of 0; using formulas
Calculating the microcrack porosity distribution spectrum with an aspect ratio of alpha under an effective pressure P, wherein phi isnAnd (alpha, P) is the n < th > set of microcracked porosities with aspect ratio alpha under effective pressure P.
7. The method of claim 1, wherein determining a seismic velocity parameter corresponding to the rock sample in conjunction with the pore structure parameter of the rock sample and the microfracture porosity profile comprises:
acquiring a wet skeleton modulus based on the expanded Gurevich jet flow model; the wet skeleton modulus comprises a wet skeleton volume modulus and a wet skeleton shear modulus; wherein, include: using formulas
Calculating the wet skeleton bulk modulus, where Kwf(P, ω) is the wet skeleton bulk modulus,
is the high-pressure bulk modulus, N is the number of pressure intervals, phin(P) is a microcracked porosity profile corresponding to pressure P,
wherein, KD(P) is the theoretical bulk modulus, α, corresponding to the pressure Pn(P) is the microcrack aspect ratio profile corresponding to pressure P,
to correspond to the cumulative fracture density of the pressure P,
J1(xi) is a first order Bessel function,
where ω is the phase, a is the aperture parameter, ρflIs fluid density, η is pore fluid viscosity, J0(xi) is a zero order Bessel function, KfIs the bulk modulus of fluid, KmIs the mineral bulk modulus; using formulas
Calculating the wet skeleton shear modulus, wherein Gwf(P, ω) is the wet skeleton shear modulus, GD(P) is the theoretical shear modulus corresponding to pressure P;
calculating speed dispersion and attenuation by combining pore structure parameters and the wet skeleton modulus based on a Biot pore elasticity theory; wherein, include: using the formula V-1/Re (1/V)c) Calculating velocity dispersion, where V is phase velocity and VcFor complex velocities calculated by the Biot model, including shear complex velocities
And longitudinal wave complex velocity
Complex velocity of transverse wave
Where Δ ═ P ρ22+Rρ11-2Qρ12，
M＝[(τ-φ)/Km+φ/Kf]-1，τ＝1-Kwf/Km，KwfIs the wet skeleton bulk modulus, phi is the porosity, GwfIs the wet skeleton shear modulus, ρ22＝aφρfA is the tortuosity, ρfFor pore fluid density, R-M φ2，ρ11＝(1-φ)ρm-ρ12，ρmIs the rock matrix modulus, ρ12＝(1-a)φρf，Q＝M(τ-φ)φ，τ＝1-Kwf/KmVelocity of complex longitudinal wave
Where ρ is ρD+φρf，ρDDensity of the dry rock sample; using the formula Q-1＝Im(Vc2)/Re(Vc2) Calculating the reciprocal of the quality factor, wherein Q-1Is the reciprocal of the quality factor.
8. The method of claim 1, wherein the pore structure parameter comprises at least one of density, porosity, estimated bulk modulus, and estimated shear modulus;
the estimated bulk modulus and the estimated shear modulus are obtained by the following method:
performing core XRD analysis on the rock sample to obtain mineral components corresponding to the rock sample;
using formulas
Calculating the estimated bulk modulus, where KmFor estimating the bulk modulus, N is the amount of mineral constituent, fiIs the volume fraction of the ith mineral component, KiIs the bulk modulus of the ith mineral component;
using formulas
Calculation of the estimated shear modulus, wherein GmTo estimate the shear modulus, GiIs the shear modulus of the ith mineral component;
correspondingly, before obtaining the pore structure parameter, the method further includes:
pre-treating the rock sample; the pretreatment includes at least one of a wash oil treatment, a salt washing treatment and a drying treatment.
9. A seismic wave velocity parameter determination device based on a DEM model is characterized by comprising:
the modulus acquisition module is used for acquiring the actually measured bulk modulus and the actually measured shear modulus of the rock sample under at least two different pressures;
the high-pressure modulus calculation module is used for calculating the high-pressure modulus by utilizing the actually measured bulk modulus and the actually measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores comprise pores that are not compressible under confining pressure;
the accumulated fracture density determining module is used for determining the accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on the relation between the modulus and the fracture density under the DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample;
the microcrack porosity distribution spectrum calculation module is used for calculating the microcrack porosity distribution spectrum under the at least two different pressures according to the accumulated fracture density;
the seismic wave velocity parameter determining module is used for determining a seismic wave velocity parameter corresponding to the rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
10. A seismic wave velocity parameter determination device based on a DEM model comprises a memory and a processor;
the memory to store computer program instructions;
the processor to execute the computer program instructions to implement the steps of: acquiring an actually measured bulk modulus and an actually measured shear modulus of a rock sample under at least two different pressures; calculating a high-pressure modulus by using the measured bulk modulus and the measured shear modulus; the high-pressure modulus is used to represent the equivalent modulus of a rock sample consisting of a solid mineral matrix and hard pores; the hard pores comprise pores that are not compressible under confining pressure; determining an accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on a modulus-fracture density relation under a DEM model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample; calculating microcrack porosity distribution spectra at the at least two different pressures from the cumulative fracture density; determining a seismic wave velocity parameter corresponding to a rock sample by combining the pore structure parameter of the rock sample and the microcrack porosity distribution spectrum; the seismic wave velocity parameters include velocity dispersion and attenuation.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110313216.3A CN113009561A (en) | 2021-03-24 | 2021-03-24 | DEM model-based seismic wave velocity parameter determination method, device and equipment |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110313216.3A CN113009561A (en) | 2021-03-24 | 2021-03-24 | DEM model-based seismic wave velocity parameter determination method, device and equipment |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN113009561A true CN113009561A (en) | 2021-06-22 |
Family
ID=76406025
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110313216.3A Pending CN113009561A (en) | 2021-03-24 | 2021-03-24 | DEM model-based seismic wave velocity parameter determination method, device and equipment |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113009561A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113740145A (en) * | 2021-09-06 | 2021-12-03 | 中国工程物理研究院电子工程研究所 | Device and method for testing bulk modulus of elastomer material |
| CN114236609A (en) * | 2021-12-17 | 2022-03-25 | 河海大学 | Prediction method for longitudinal wave velocity and attenuation of partially saturated hole fracture medium |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104360383A (en) * | 2014-11-12 | 2015-02-18 | 中国石油大学(华东) | Method and system for predicting seismic wave attenuation |
| CN111208565A (en) * | 2020-03-04 | 2020-05-29 | 中国石油大学(北京) | KT model-based hole seam parameter inversion method and device and storage medium |
| CN112505772A (en) * | 2020-12-10 | 2021-03-16 | 中国石油大学(华东) | Method for inverting rock pore distribution characteristics by utilizing pore and fracture medium elastic wave theory |
-
2021
- 2021-03-24 CN CN202110313216.3A patent/CN113009561A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104360383A (en) * | 2014-11-12 | 2015-02-18 | 中国石油大学(华东) | Method and system for predicting seismic wave attenuation |
| CN111208565A (en) * | 2020-03-04 | 2020-05-29 | 中国石油大学(北京) | KT model-based hole seam parameter inversion method and device and storage medium |
| CN112505772A (en) * | 2020-12-10 | 2021-03-16 | 中国石油大学(华东) | Method for inverting rock pore distribution characteristics by utilizing pore and fracture medium elastic wave theory |
Non-Patent Citations (5)
| Title |
|---|
| E. C. DAVID AND R. W. ZIMMERMAN: "Pore structure model for elastic wave velocities in fluid-saturated sandstones", 《JOURNAL OF GEOPHYSICAL RESEARCH》 * |
| 李东庆: "致密砂岩储层岩石物理及AVO特征的实验研究", 《中国优秀博硕士学位论文全文数据库(博士)基础科学辑》 * |
| 欧阳芳等: "基于微观孔隙结构特征的速度频散和衰减模拟", 《地球物理学报》 * |
| 欧阳芳等: "基于等效介质理论的孔隙纵横比分布反演", 《地球物理学报》 * |
| 邓继新等: "基于储层砂岩微观孔隙结构特征的弹性波频散响应分析", 《地球物理学报》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113740145A (en) * | 2021-09-06 | 2021-12-03 | 中国工程物理研究院电子工程研究所 | Device and method for testing bulk modulus of elastomer material |
| CN113740145B (en) * | 2021-09-06 | 2023-05-05 | 中国工程物理研究院电子工程研究所 | Device and method for testing bulk modulus of elastomer material |
| CN114236609A (en) * | 2021-12-17 | 2022-03-25 | 河海大学 | Prediction method for longitudinal wave velocity and attenuation of partially saturated hole fracture medium |
| CN114236609B (en) * | 2021-12-17 | 2022-07-19 | 河海大学 | Prediction method for longitudinal wave velocity and attenuation of partially saturated hole fracture medium |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113009561A (en) | DEM model-based seismic wave velocity parameter determination method, device and equipment | |
| CN111425193B (en) | Reservoir compressibility evaluation method based on clustering analysis logging rock physical facies division | |
| CN113009565A (en) | Method, device and equipment for determining seismic wave velocity parameters based on SCA model | |
| CN109490965B (en) | Method and device for quantitatively evaluating formation heterogeneity | |
| CN112946783B (en) | Hydrate saturation determination method, device and equipment | |
| CN108181654A (en) | AVAF analogy methods and device based on multi-scale rock physical model | |
| CN107203005A (en) | A kind of method that quantification calculates crack characterising parameter | |
| CN105093331B (en) | The method for obtaining Rock Matrix bulk modulus | |
| CN111123354A (en) | Method and equipment for predicting dense gas layer based on frequency-dependent reflection amplitude attenuation | |
| CN110515126A (en) | A kind of velocity of sound calculation method of the transverse isotropy of crack containing random distribution rock | |
| CN113093276A (en) | Method, device, equipment and system for predicting seismic wave velocity dispersion and attenuation | |
| CN113009562A (en) | KT model-based seismic wave velocity parameter determination method, device and equipment | |
| CN114488302A (en) | In-situ anisotropic ground stress field prediction method and system | |
| CN111812709A (en) | Method, device and equipment for establishing multi-scale wave induced flow model | |
| Yang et al. | All-parameters Rayleigh wave inversion | |
| CN113009563A (en) | Seismic wave velocity parameter determination method, device and equipment based on MT model | |
| CN109709610B (en) | Rock crack detection method and system | |
| Jie et al. | Analysis and application of automatic deformation monitoring data for buildings and structures of mining area | |
| Telesca et al. | Investigating the dynamical features of the time distribution of the reservoir-induced seismicity in Enguri area (Georgia) | |
| CN109975189B (en) | Method and device for predicting productivity of pore type sandstone reservoir | |
| CN109490988A (en) | Establish the method for being suitable for the rock physics new model of hard rock | |
| CN117471540A (en) | Method, device and equipment for predicting physical parameters of reservoir based on post-stack seismic data | |
| CN112230278B (en) | Seepage field characteristic parameter determining method and device | |
| JP6778628B2 (en) | Underground structure estimator | |
| CN112782780A (en) | Reservoir evaluation method, device and equipment based on rock physical facies |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20210622 |