CN113009562A - KT model-based seismic wave velocity parameter determination method, device and equipment - Google Patents

KT model-based seismic wave velocity parameter determination method, device and equipment Download PDF

Info

Publication number
CN113009562A
CN113009562A CN202110313218.2A CN202110313218A CN113009562A CN 113009562 A CN113009562 A CN 113009562A CN 202110313218 A CN202110313218 A CN 202110313218A CN 113009562 A CN113009562 A CN 113009562A
Authority
CN
China
Prior art keywords
modulus
pressure
rock sample
fracture density
calculating
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
Application number
CN202110313218.2A
Other languages
Chinese (zh)
Inventor
赵建国
闫秀懿
欧阳芳
赵皓
李智
肖增佳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum Beijing
Original Assignee
China University of Petroleum Beijing
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Petroleum Beijing filed Critical China University of Petroleum Beijing
Priority to CN202110313218.2A priority Critical patent/CN113009562A/en
Publication of CN113009562A publication Critical patent/CN113009562A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282Application of seismic models, synthetic seismograms
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/303Analysis for determining velocity profiles or travel times
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/307Analysis 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 KT 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 relationship under a KT 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

KT model-based seismic wave velocity parameter determination method, device and equipment
Technical Field
The embodiment of the specification relates to the technical field of geological exploration and development, in particular to a seismic wave velocity parameter determination method, device and equipment based on a KT 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
An object of an embodiment of the specification is to provide a method, a device and equipment for determining seismic wave velocity parameters based on a KT model, so as to solve the problem of how to accurately solve 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 KT 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 relationship under a KT 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 device for determining seismic wave velocity parameters based on a KT 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 determination module is used for determining accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on the relation between the modulus and the fracture density in a KT 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 further provides a device for determining the seismic wave velocity parameters based on the KT 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 relationship under a KT 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 flowchart of a method for determining seismic wave velocity parameters based on a KT 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 determination apparatus based on a KT model according to an embodiment of the present disclosure;
fig. 20 is a block diagram of a seismic velocity parameter determination device based on a KT 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, an embodiment of the specification provides a method for determining seismic wave velocity parameters based on a KT model. An execution body of the seismic wave velocity parameter determination method based on the KT model is seismic wave velocity parameter determination equipment based on the KT model. As shown in FIG. 1, the KT model-based seismic wave velocity parameter determination method specifically comprises 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
Figure BDA0002990116950000041
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, the individual mineral constituents can be obtained during the analysisVolume fraction, bulk modulus, and shear modulus. Then, the formula can be utilized
Figure BDA0002990116950000051
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
Figure BDA0002990116950000052
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
Figure BDA0002990116950000061
In the formula, KD(P) is the theoretical bulk modulus of the rock sample under pressure P,
Figure BDA0002990116950000062
the bulk modulus at zero pressure is given,
Figure BDA0002990116950000063
in order to have a high bulk modulus at high pressure,
Figure BDA0002990116950000064
is the pressure coefficient. The shear modulus relationship constructed may be
Figure BDA0002990116950000065
In the formula, GD(P) is the theoretical shear modulus of the rock sample at pressure P,
Figure BDA0002990116950000066
is a shear modulus at zero pressure and a shear modulus,
Figure BDA0002990116950000067
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 relationship under a KT model; the cumulative fracture density represents a sum of fracture densities of open-hole microfractures in the rock sample.
The relation between the modulus and the fracture density in the KT model comprises the relation between the bulk modulus and the fracture density and the relation between the shear modulus and the fracture density. In particular, the bulk modulus versus fracture density relationship can include
Figure BDA0002990116950000068
In the formula (I), the compound is shown in the specification,
Figure BDA0002990116950000069
equivalent medium bulk modulus, K, estimated for KT modelbEquivalent modulus of elasticity, G, to the high-pressure bulk modulusbEquivalent modulus of elasticity, upsilon, which is the high-pressure shear modulusb=(3Kb-2Gb)/(6Kb+2Gb) Gamma-ray is fracture density; the shear modulus to fracture density relationship comprises
Figure BDA00029901169500000610
In the formula (I), the compound is shown in the specification,
Figure BDA00029901169500000611
the equivalent medium shear modulus estimated for the KT model,
Figure BDA00029901169500000612
since the KT model already defines the relation between modulus and flaw density. In the case that the actually measured bulk modulus and the actually measured shear modulus have been obtained in step S110, the corresponding fracture density may be obtained by a fitting manner based on the equivalent medium bulk modulus and the equivalent medium shear modulus estimated based on the KT model in the above relational expression. For example, a relational expression can be constructed
Figure BDA00029901169500000613
OF in the formula2In order to be the second objective function,
Figure BDA00029901169500000614
corresponding to a pressure interval P estimated for the KT modeliEquivalent medium bulk modulus of (2), KD_meas(Pi) To correspond to the pressure interval PiThe measured bulk modulus of the polymer (a),
Figure BDA0002990116950000071
corresponding to a pressure interval P estimated for the KT 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-1Elasticity of lower partModulus as elastic parameter of rock background matrix and rock equivalent medium, namely high-pressure modulus
Figure BDA0002990116950000072
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
Figure BDA0002990116950000073
And then obtaining a pressure interval p based on the relation between the modulus and the fracture density under the KT modelN→pN-1Lower microcrack density gamma (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
Figure BDA0002990116950000074
Finally obtaining a pressure interval pi→pi-1Lower microcrack density gamma (P)i)。
Sequentially solving the micro-crack density corresponding to each pressure interval through the steps, and finally utilizing a formula
Figure BDA0002990116950000075
The cumulative fracture density is calculated, where,
Figure BDA0002990116950000076
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
Figure BDA0002990116950000081
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,
Figure BDA0002990116950000082
is the undetermined coefficient.
Since the values corresponding to the respective ones have been calculated in the above-described step S130Fracture density in pressure interval
Figure BDA0002990116950000083
Thus, a formula can be constructed
Figure BDA0002990116950000084
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,
Figure BDA0002990116950000085
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
Figure BDA0002990116950000086
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,
Figure BDA0002990116950000087
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 distribution spectrum, α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
Figure BDA0002990116950000088
Is calculated by the formula
Figure BDA0002990116950000089
Calculating the microcrack porosity at 0 confining pressure, wherein,
Figure BDA00029901169500000810
porosity of the nth set of microcracks at 0 confining pressure. Then, the formula can be utilized
Figure BDA00029901169500000811
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
Figure BDA0002990116950000091
In the formula, Kwf(P, ω) is the wet skeleton bulk modulus,
Figure BDA0002990116950000092
is the high-pressure bulk modulus, N is the number of pressure intervals, phin(P) is a microcracked porosity profile corresponding to pressure P,
Figure BDA0002990116950000093
wherein, KD(P) is the theoretical bulk modulus, α, corresponding to the pressure Pn(P) is the microcrack aspect ratio profile corresponding to pressure P,
Figure BDA0002990116950000094
to correspond to the cumulative fracture density of the pressure P,
Figure BDA0002990116950000095
J1(xi) is a first order Bessel function,
Figure BDA0002990116950000096
where, ω is the phase, a is the microcrack pore aspect ratio distribution spectrum, ρ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
Figure BDA0002990116950000097
In the formula, Gwf(P, ω) 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
Figure BDA0002990116950000098
And longitudinal wave complex velocity
Figure BDA0002990116950000099
Complex velocity of transverse wave
Figure BDA0002990116950000101
Where Δ ═ P ρ22+Rρ11-2Qρ12
Figure BDA0002990116950000102
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-φ)ρm12,ρ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
Figure BDA0002990116950000103
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.
Figure BDA0002990116950000104
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.
Figure BDA0002990116950000111
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.
Figure BDA0002990116950000112
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 KT model, and inverting the accumulated fracture density under each pressure according to the objective 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
Figure BDA0002990116950000113
The specific fitting function is shown in table 5 below, and the fitting results are shown in the fitted curve in fig. 6.
Figure BDA0002990116950000114
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.
Figure BDA0002990116950000121
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.
Figure BDA0002990116950000122
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.
Figure BDA0002990116950000123
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 KT model, and inverting the accumulated fracture density under each pressure according to the objective 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
Figure BDA0002990116950000125
The specific fitting function is shown in table 10 below, and the fitting results are shown in the fitted curve in fig. 15.
Figure BDA0002990116950000124
Watch 10
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 parameters based on the KT model, the specification further provides an embodiment of a device for determining the seismic wave velocity parameters based on the KT model. As shown in fig. 19, the apparatus for determining seismic wave velocity parameters based on KT 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.
A cumulative fracture density determination module 1930 configured to determine a cumulative fracture density corresponding to the rock sample from the high-pressure modulus based on a modulus-to-fracture density relationship in the KT model; 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 method for determining the seismic wave velocity parameters based on the KT model, the embodiment of the specification further provides a device for determining the seismic wave velocity parameters based on the KT model. As shown in fig. 20, the KT 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 relationship under a KT 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 KT 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 relationship under a KT 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
Figure FDA0002990116940000011
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
Figure FDA0002990116940000012
In the formula, KD(P) is the theoretical bulk modulus of the rock sample under pressure P,
Figure FDA0002990116940000013
the bulk modulus at zero pressure is given,
Figure FDA0002990116940000014
in order to have a high bulk modulus at high pressure,
Figure FDA0002990116940000015
is the pressure coefficient; the shear modulus relationship comprising
Figure FDA0002990116940000016
In the formula, GD(P) is the theoretical shear modulus of the rock sample at pressure P,
Figure FDA0002990116940000017
is a shear modulus at zero pressure and a shear modulus,
Figure FDA0002990116940000018
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 applying formula PF1=∑[(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 under pressure PMeasured shear modulus of.
4. The method of claim 1, wherein determining a cumulative fracture density corresponding to the rock sample from the high-pressure modulus based on a modulus versus fracture density relationship under the KT 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
Figure FDA0002990116940000021
Fitting to obtain micro-fracture densities respectively corresponding to different pressure intervals; OF in the formula2In order to be the second objective function,
Figure FDA0002990116940000022
corresponding to a pressure interval P estimated for the KT modeliEquivalent medium bulk modulus of (2), KD_meas(Pi) To correspond to the pressure interval PiThe measured bulk modulus of the polymer (a),
Figure FDA0002990116940000023
corresponding to a pressure interval P estimated for the KT 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
Figure FDA0002990116940000024
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 relation between the modulus and the fracture density under the KT model comprises a relation between the bulk modulus and the fracture density and a relation between the shear modulus and the fracture density;
the bulk modulus to fracture density relationship comprises
Figure FDA0002990116940000025
In the formula (I), the compound is shown in the specification,
Figure FDA0002990116940000026
equivalent medium bulk modulus, K, estimated for KT modelbEquivalent modulus of elasticity, G, to the high-pressure bulk modulusbEquivalent modulus of elasticity, v, to high shear modulusb=(3Kb-2Gb)/(6Kb+2Gb) Gamma-ray is fracture density;
the shear modulus to fracture density relationship comprises
Figure FDA0002990116940000027
In the formula (I), the compound is shown in the specification,
Figure FDA0002990116940000028
the equivalent medium shear modulus estimated for the KT model,
Figure FDA0002990116940000029
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
Figure FDA0002990116940000031
And
Figure FDA0002990116940000032
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,
Figure FDA0002990116940000033
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,
Figure FDA0002990116940000034
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
Figure FDA0002990116940000035
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,
Figure FDA0002990116940000036
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
Figure FDA0002990116940000037
Calculating the microcrack porosity at 0 confining pressure, wherein,
Figure FDA0002990116940000038
porosity of the nth set of microcracks at a confining pressure of 0; using formulas
Figure FDA0002990116940000039
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
Figure FDA00029901169400000310
Calculating the wet skeleton bulk modulus, where Kwf(P, w) is the wet skeleton bulk modulus,
Figure FDA00029901169400000311
is the high-pressure bulk modulus, N is the number of pressure intervals, phin(P) is a microcracked porosity profile corresponding to pressure P,
Figure FDA00029901169400000312
it is composed of
In, KD(P) is a radical ofTheoretical bulk modulus of pressure P, alphan(P) is the microcrack aspect ratio profile corresponding to pressure P,
Figure FDA0002990116940000041
to correspond to the cumulative fracture density of the pressure P,
Figure FDA0002990116940000042
J1(xi) is a first order Bessel function,
Figure FDA0002990116940000043
where, ω is the phase, a is the microcrack pore aspect ratio distribution spectrum, ρ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
Figure FDA0002990116940000044
Calculating the wet skeleton shear modulus, wherein Gwf(P, w) 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
Figure FDA0002990116940000045
And longitudinal wave complex velocity
Figure FDA0002990116940000046
Complex velocity of transverse wave
Figure FDA0002990116940000047
Where Δ ═ P ρ22+Rρ11-2Qρ12
Figure FDA0002990116940000048
KwfWet skeleton bulk modulus, porosity, GwfIs the wet skeleton shear modulus, ρ22=aφρfA is the tortuosity, ρfFor pore fluid density, R-M φ2,ρ11=(1-φ)ρm12,ρmIs the rock matrix modulus, ρ12=(1-a)φρf,Q=M(τ-φ)φ,τ=1-Kwf/KmVelocity of complex longitudinal wave
Figure FDA0002990116940000049
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
Figure FDA00029901169400000410
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
Figure FDA00029901169400000411
Calculation of the estimated shear modulus, wherein GmTo estimateShear 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 KT model-based seismic wave velocity parameter determination apparatus, 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 determination module is used for determining accumulated fracture density corresponding to the rock sample through the high-pressure modulus based on the relation between the modulus and the fracture density in a KT 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 KT model-based seismic wave velocity parameter determination device 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 relationship under a KT 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.
CN202110313218.2A 2021-03-24 2021-03-24 KT model-based seismic wave velocity parameter determination method, device and equipment Pending CN113009562A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110313218.2A CN113009562A (en) 2021-03-24 2021-03-24 KT model-based seismic wave velocity parameter determination method, device and equipment

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110313218.2A CN113009562A (en) 2021-03-24 2021-03-24 KT model-based seismic wave velocity parameter determination method, device and equipment

Publications (1)

Publication Number Publication Date
CN113009562A true CN113009562A (en) 2021-06-22

Family

ID=76406026

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110313218.2A Pending CN113009562A (en) 2021-03-24 2021-03-24 KT model-based seismic wave velocity parameter determination method, device and equipment

Country Status (1)

Country Link
CN (1) CN113009562A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114924318A (en) * 2022-05-06 2022-08-19 中国石油大学(华东) Seismic rock physical modeling method for stable prediction of mineral modulus

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105425280A (en) * 2015-11-21 2016-03-23 西南石油大学 Prediction method for mineral modulus and pore structure
CN109471168A (en) * 2018-11-06 2019-03-15 河海大学 The prediction technique of velocity of longitudinal wave and decaying in a kind of hole fissuted medium
CN110276091A (en) * 2019-04-26 2019-09-24 中国石油化工股份有限公司 Elastic wave response model modelling approach based on rock multi-modal pore system structure
CN111208565A (en) * 2020-03-04 2020-05-29 中国石油大学(北京) KT model-based hole seam parameter inversion method and device and storage medium
CN111208566A (en) * 2020-03-04 2020-05-29 中国石油大学(北京) Hole seam parameter inversion method and device based on SCA model and storage medium

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105425280A (en) * 2015-11-21 2016-03-23 西南石油大学 Prediction method for mineral modulus and pore structure
CN109471168A (en) * 2018-11-06 2019-03-15 河海大学 The prediction technique of velocity of longitudinal wave and decaying in a kind of hole fissuted medium
CN110276091A (en) * 2019-04-26 2019-09-24 中国石油化工股份有限公司 Elastic wave response model modelling approach based on rock multi-modal pore system structure
CN111208565A (en) * 2020-03-04 2020-05-29 中国石油大学(北京) KT model-based hole seam parameter inversion method and device and storage medium
CN111208566A (en) * 2020-03-04 2020-05-29 中国石油大学(北京) Hole seam parameter inversion method and device based on SCA model and storage medium

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
李东庆: "致密砂岩储层岩石物理及AVO特征的实验研究", 《中国优秀博硕士学位论文全文数据库(博士)基础科学辑》 *
欧阳芳 等: "基于微观孔隙结构特征的速度频散和衰减模拟", 《地球物理学报》 *
欧阳芳 等: "基于等效介质理论的孔隙纵横比分布反演", 《地球物理学报》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114924318A (en) * 2022-05-06 2022-08-19 中国石油大学(华东) Seismic rock physical modeling method for stable prediction of mineral modulus

Similar Documents

Publication Publication Date Title
CN111425193B (en) Reservoir compressibility evaluation method based on clustering analysis logging rock physical facies division
CN109374497B (en) Rock micro-pore structure testing method
CN113009561A (en) DEM model-based seismic wave velocity parameter determination method, device and equipment
CN109490965B (en) Method and device for quantitatively evaluating formation heterogeneity
CN113009565A (en) Method, device and equipment for determining seismic wave velocity parameters based on SCA model
CN103293563B (en) Method for determining rock fracture development degree and fluid property of oil and gas reservoir
CN108181654A (en) AVAF analogy methods and device based on multi-scale rock physical model
US9892321B2 (en) Using maximal inscribed spheres for image-based rock property estimation
CN112946783B (en) Hydrate saturation determination method, device and equipment
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
EA027440B1 (en) Method of predicting the pressure sensitivity of seismic velocity within reservoir rocks
CN105931125B (en) Method for predicting yield of compact oil staged multi-cluster volume fracturing horizontal well
Aquino-López et al. Modeling and inversion of elastic wave velocities and electrical conductivity in clastic formations with structural and dispersed shales
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
CN110261898A (en) Well logging and seismic velocity matching process based on the analysis of earthquake petrophysics experiment
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
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

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

Application publication date: 20210622

RJ01 Rejection of invention patent application after publication