Well shock frequency dispersion correction method suitable for carbonate rock fracture-vuggy reservoir

Abstract

A well seismic velocity frequency dispersion correction method suitable for a fracture-vuggy carbonate reservoir comprises the following steps: (1) on the basis of fracture-vuggy depicted by imaging logging FMI, carrying out conventional logging curve fracture prediction and fracture characteristic parameter calculation based on a support vector machine; (2) selecting a target fracture-vug type carbonate reservoir core to carry out lithology, physical property and low-frequency rock physical parameter measurement; (3) constructing a full-band fractured-vuggy carbonate rock physical model-pore elastic equivalent medium model; (4) calibrating a pore elasticity equivalent medium model by experimental data; and (5) well seismic velocity dispersion correction. According to the method, on the basis that FMI imaging data depicts fracture characteristics, the fracture recognition method based on a conventional logging curve is formed, the fracture recognition precision reaches 85%, and a data foundation is laid for describing an elastic wave propagation rule of an igneous rock fracture reservoir; a pore elasticity equivalent medium model more suitable for the fracture-vug type carbonate reservoir is constructed, and the calculated seismic frequency band speed is higher in precision.

Description

Technical Field

[0001] The invention relates to the technical field of reservoir rock physical property research, and in particular to a well seismic dispersion correction method. Background Art

[0002] When seismic waves propagate in the formation, the velocity of seismic waves of different frequencies will change, and energy will attenuate in the frequency band where the velocity changes dramatically. This characteristic of velocity changing with frequency is called velocity dispersion. In actual production, due to the inconsistency of well seismic velocity data, there is always a problem of mismatching well seismic data when making synthetic seismic records from logging data. Stretching or compression is required (only when the cable is stretched or compressed, it makes sense to do so, but generally the cable cannot be detected to be stretched or compressed) to improve the accuracy of synthetic record production. In addition to scale problems, wellbore quality problems, anisotropy and other logging data quality problems, the difference in measurement frequency cannot be ignored. In actual engineering, the measurement frequency of acoustic logging is generally in the range of 2KHz to 20KHz, which is much higher than the frequency of seismic waves in conventional surface exploration. The significant frequency difference makes the propagation speed of logging acoustic waves and surface exploration seismic waves different in heterogeneous formations containing fluids. In terms of intrinsic physical mechanism, there is a mismatch between logging synthetic records and seismic traces.

[0003] In recent years, some scholars have carried out research on well-seismic dispersion correction methods. Zhao Jianguo et al. measured the full-band P- and S-wave velocities of ancient Ordovician-Cambrian carbonate rocks using low-frequency measurement equipment, and used the seismic band velocities at different pressures to fit the relationship between formation pressure and seismic band velocity. Then, the logging band velocities at different pressures were used to fit the relationship between formation pressure and logging band velocity, and then the logging and seismic velocities were fitted. This method only considered the relationship between formation pressure and well-seismic velocity, and did not consider the influence of formation mineral composition, microscopic pores, fluid saturation and viscosity. Therefore, the prediction effect is less accurate. Deng Jixin et al. used the resonant Q model to correct well-seismic dispersion, in which the value of Q is a constant. However, it is known from experiments that if there is dispersion between well-seismic measurements, Q must not be a constant, so the accuracy of well-seismic dispersion correction is also low. Zhang Yuanzhong et al. extrapolated ultrasonic experimental data to achieve dispersion correction between ultrasonic and logging data, but did not involve seismic bands.

[0004] P-wave and S-wave velocity data are irreplaceable in oil and gas exploration, and play a very important role in time-depth conversion, synthetic record production, seismic data inversion, reservoir modeling, etc. However, in some formations, such as fracture-cavity carbonate reservoirs, due to the presence of cracks, which then produce non-negligible velocity dispersion, it is risky to directly use the P-wave and S-wave velocities of the logging frequency band instead of the P-wave and S-wave velocities of the seismic frequency band to predict reservoir fluids. Summary of the invention

[0005] To solve the problem of the influence of well-seismic dispersion, the present invention provides a well-seismic velocity dispersion correction method suitable for fractured-vuggy carbonate reservoirs.

[0006] The technical solution of the present invention is as follows:

[0007] A well-seismic velocity dispersion correction method suitable for fractured-vuggy carbonate reservoirs, characterized by comprising the following steps:

[0008] (1) Based on the characterization of fractures and vugs by imaging logging FMI, carry out fracture prediction of conventional logging curves based on support vector machines and calculate their fracture characteristic parameters;

[0009] (2) Select cores of the target fractured-vuggy carbonate reservoir to carry out measurements of lithology, physical properties, and low-frequency rock physics parameters;

[0010] (3) Construct a full-frequency carbonate fractured-vuggy rock physics model - a pore-elastic equivalent medium model;

[0011] (4) Calibrate the pore-elastic equivalent medium model with experimental data;

[0012] (5) Well-seismic velocity dispersion correction.

[0013] Preferably, step (1) specifically includes the following steps:

[0014] ① Select wells in the work area containing FMI imaging logging data, count the characteristics of fractures or holes and the response relationships between non-fractures and conventional curves, and establish training samples;

[0015] ② Normalize the sample data;

[0016] ③ Select the kernel function and penalty factor;

[0017] ④ Discriminate high-conductivity fractures;

[0018] ⑤ Calculate fracture characteristic parameters.

[0019] More preferably, in step ①, count no less than 60 conventional curve response values of high-conductivity fractures and no less than 40 conventional curve response values of non-high-conductivity fractures.

[0020] More preferably, in step ②, the normalization interval is [-1, 1]. The normalization calculation formula is as follows:

[0021]

[0022] where x i 、x imax 、x iminrespectively represent the original logging value of the i-th conventional curve, the maximum original logging value of the i-th conventional curve, the minimum original logging value of the i-th conventional curve, and y i represents the normalized value of the i-th conventional curve.

[0023] Further preferably, in step ③, a Gaussian kernel is selected as the kernel function of the support vector machine method, as shown in the following formula:

[0024]

[0025] where m is the imaging logging eigenvalue, d is the order of the polynomial, σ is the width of the Gaussian distribution, and the value of the penalty factor C is 800.

[0026] Further preferably, in step ④, through the selection of samples, feature normalization, kernel function, and penalty factor, the support vector machine is trained using the training sample set to obtain a discriminant function for high-conductivity fracture identification based on the support vector machine method, as shown in the following formula:

[0027] n = sgn{∑α i n i k(y i , m) + b}

[0028] where α i is the gradient coefficient, n and n i are the outputs of the predicted high-conductivity fracture samples and the support vector samples respectively, k(y i , m) is the Gaussian kernel function, and b is the intercept.

[0029] Further preferably, in step ⑤,

[0030] Calculation of fracture angle, fracture width, and fracture density:

[0031]

[0032] Y < 0, θ is less than 30 degrees; 0 < Y < 0.1, θ is between 30 - 60 degrees; Y > 0.1, θ is greater than 60 degrees;

[0033] where Y is the discriminant factor, θ is the fracture dip angle, and qualitative discrimination is carried out in different regions. R d is the deep lateral resistivity, and R s is the shallow lateral resistivity;

[0034] Fracture porosity Φ f is obtained using the dual laterolog resistivity calculation formula: When Rd > Rs:

[0035]

[0036] When Rd < Rs:

[0037]

[0038] The crack width (ε) is in micrometers and is obtained by the dual laterolog resistivity calculation formula: When Rd > Rs:

[0039]

[0040] When Rd > Rs:

[0041]

[0042] Where: Rd and Rs are the deep and shallow laterolog resistivities, unit: Ω·m; ε is the crack width, unit: micrometer, Rmf and Rw are the mud filtrate and formation water resistivities, unit: Ω·m; Φ f is the crack porosity, unit: %; m f is the crack cementation index, which is 1.5; For the Ordovician, Rw is taken as 0.012 Ω·m, R mf = 2.8; R b is the resistivity of the tight limestone layer, with a value of 40000 Ω·m.

[0043] Preferably, in step (2), select the core of the fracture-vuggy carbonate rock in the target formation to carry out parameter measurements including porosity, permeability, density, and mineral composition, and then carry out low-frequency rock physical parameter measurements under formation conditions. The measurement frequency is 12 - 200 Hz, and the saturation states are dry, water-saturated, and oil-saturated, to obtain the porosity, permeability, density, mineral composition data of the core, and the low-frequency P-wave and S-wave velocity data of the core under formation conditions.

[0044] Preferably, in step (3),

[0045] The bulk modulus K eff and shear modulus u eff of the equivalent medium of the poroelastic equivalent medium model are expressed as:

[0046]

[0047]

[0048] Among them, ω is the angular frequency, K is the bulk modulus of the skeleton, μ is the shear modulus, ε f is the microcrack density, a is the crack radius, r is the microcrack radius, ε c is the crack density, θ is the dip angle, is the porosity, λ is the Lame coefficient, Kc, Kp, A, and B are all intermediate parameters without specific physical meanings. A and B are respectively

[0049]

[0050]

[0051] Among them,

[0052]

[0053]

[0054] k f is the fluid bulk modulus, τ is the time scale, which controls the frequency band range where the dispersion characteristics appear, is proportional to the viscosity η of the fluid, inversely proportional to the permeability k, and is related to the fracture radius a. Its reciprocal is called the characteristic frequency;

[0055] For a smaller aspect ratio:

[0056]

[0057] is the particle scale, ν is the Poisson's ratio of the solid particles; according to the formula proposed by Connell and Budiansky, the P-wave velocity V P and the S-wave velocity V S are respectively:

[0058]

[0059]

[0060] Among them, Re is to take the real part, and ρ is the density.

[0061] Preferably, in step (5),

[0062] First, obtain the lithology component data through well logging interpretation result data, then calculate the frame modulus based on the Voigt-Reuss-Hill model. On this basis, combine the well logging interpreted porosity and predicted fracture parameters to calculate the dry rock elastic modulus based on the Berrymann model. Then, further combine the saturation data interpreted by well logging to calculate the P-wave velocity at well logging frequency based on the poroelastic equivalent medium model. If the absolute difference between the calculated P-wave velocity and the measured P-wave velocity is less than the error, then directly calculate the P- and S-wave velocities in the seismic frequency band. If the difference between the calculated P-wave velocity and the measured P-wave velocity is greater than the error, then adjust the microcrack density and calculate again according to the above process, and repeat the cycle until the absolute difference between the calculated P-wave velocity and the measured P-wave velocity is less than the error;

[0063] ① Frame modulus calculation - Voigt-Reuss-Hill average theory:

[0064] Given the relative content and elastic modulus of each phase of the known rock medium, the upper and lower limits of the effective elastic modulus of the rock medium are calculated using the Voigt and Reuss formulas. Among them, the upper limit Voigt formula represents the same strain state, where for each phase constituting the rock medium under the same strain, the ratio of stress to strain of the rock medium; the Reuss lower limit represents the same stress state, where for each phase constituting the rock medium under the same stress, the ratio of stress to strain of the rock medium. The Hill formula is the arithmetic mean result of the Voigt upper limit and the Reuss lower limit.

[0065]

[0066] Where M V 、M R 、M VRH represent the values of the bulk modulus K, shear modulus μ, and Young's modulus E obtained by the Voigt, Reuss, and Hill methods; f i is the volume content of the i-th component constituting the rock medium; M i is the elastic modulus of the i-th component;

[0067] ② Calculation of the elastic modulus of dry rock:

[0068] The Berrymann model is the general form of the self-consistent approximation for N-phase composites:

[0069]

[0070]

[0071] In the i-th material, x i is its volume fraction, and P and Q represent geometric factors. The superscript i on P and Q indicates that these factors are in the background medium with self-consistent effective moduli and that contain material i; k i , u i are the bulk modulus and shear modulus of the i-th mineral respectively. The modulus of the inclusions is set to zero to simulate dry pores;

[0072] ③ Calculation of the full-frequency P-wave and S-wave velocities of fluid-bearing rocks:

[0073] After calculating the elastic modulus of dry rock, combined with the well logging interpretation result data such as porosity, permeability, fluid saturation and other parameters, as well as the predicted fracture characteristic parameters, the calculation of the full-frequency P-wave and S-wave velocities can be realized based on the poroelastic equivalent medium model.

[0074] The beneficial technical effects of the present invention are as follows:

[0075] Based on the fracture and cave data characterized by FMI imaging logging, this paper first forms a method for quantitatively characterizing fracture characteristic parameters of conventional logging curves. Secondly, starting from the petrophysical experiments on carbonate fracture and cave cores, the rationality of the constructed petrophysical model is verified by finely calibrating the petrophysical experimental data and the petrophysical model. Then, based on the petrophysical model, by adjusting the microfracture parameters to make them coincide with or be less than the error of the calculated logging data and the actual logging data, the calculation of the P-wave and S-wave velocities in the seismic frequency band of fracture and cave carbonate reservoirs is realized.

[0076] A well-seismic dispersion correction method suitable for large-scale fracture reservoirs in carbonate rocks of the present invention has advantages that other technologies do not have. Its specific advantages and characteristics are shown in the following aspects:

[0077] First, based on characterizing fracture characteristics with FMI imaging data, a fracture identification method based on conventional logging curves is formed, and the fracture identification accuracy reaches 85%, laying a data foundation for describing the elastic wave propagation law of igneous rock fracture reservoirs.

[0078] Second, a pore-elastic equivalent medium model more suitable for fracture and cave carbonate reservoirs that simultaneously considers the microscale and mesoscale is constructed, and the seismic frequency band velocity calculated by using logging velocity as a constraint for cyclic iteration has higher accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 Schematic diagram of the actual pore structure of the rock in the embodiment of the present invention;

[0080] Figure 2 Schematic diagram of the equivalent pore structure of the rock;

[0081] Figure 3 Comparison between the measured P-wave velocity data of the core of Well Tuofu 27 and the model;

[0082] Figure 4 Comparison between the measured S-wave velocity data of the core of Well Tuofu 27 and the model;

[0083] Figure 5 Flow chart for identifying fractures with conventional logging curves;

[0084] Figure 6 Flow chart for well-seismic dispersion correction;

[0085] Figure 7 FMI imaging interpretation fracture indication diagram of Well Aiding 6;

[0086] Figure 8 Effect diagram of dispersion correction of Well No. 1;

[0087] Figure 9 Effect diagram of dispersion correction of Well No. 2;

[0088] Figure 10 It is the frequency dispersion correction effect diagram of Well No. 3. Specific implementation manners

[0089] The present invention will be further described in detail below through specific embodiments, but it does not limit the technical solution of the present invention. All changes or equivalent replacements based on the present invention shall fall within the protection scope of the present invention.

[0090] Embodiment

[0091] The present invention provides a well-seismic velocity dispersion correction method suitable for fractured-vuggy carbonate reservoirs, realizing high-precision well-seismic dispersion correction for fractured-vuggy reservoirs, including the following steps:

[0092] (1) Based on the characterization of fractures and vugs by imaging logging (FMI), carry out fracture prediction of conventional logging curves based on support vector machines and calculation of their fracture characteristic parameters.

[0093] Fullbore Formation MicroImager (FMI) is a high-resolution formation imaging logging tool. It conducts micro-resistivity scanning measurements through button electrodes attached to the wellbore wall and converts the differences in formation resistivity near the wellbore wall into clear image brightness changes. When geological bodies such as fractures and solution pores appear in the formation, the resistivity measurement values change. Therefore, FMI can clearly characterize the formation fracture characteristic parameters. However, the cost of FMI logging is expensive, and the single-well measurement price is more than a hundred times that of conventional curve logging. Therefore, in most cases, there are at most 1 to 2 wells in a work area that contain imaging logging data. For conventional logging curves, such as dual laterolog resistivity, microspherical focused resistivity, acoustic velocity, density, well diameter, natural gamma, compensated neutron and other logging data, when fractures appear in the formation, different responses will also occur on these logging curves. However, due to the existence of multiple solutions, the fracture identification accuracy is relatively low. Therefore, a method for predicting fractures (vugs) in conventional logging curves based on support vector machines and calculating their fracture characteristic parameters is developed based on the characterization of fractures and holes by FMI data. This method takes six curves of dual laterolog resistivity, acoustic velocity, density, neutron, and well diameter as inputs, uses the (Support Vector Machine) SVM method to identify high-conductivity fractures and solution holes in carbonate rocks in the study area, and then calculates the fracture characteristic parameters according to empirical formulas. The calculation process is as Figure 5 .

[0094] ① Select wells in the work area that contain FMI imaging logging data, count the fracture (vug) characteristics and the response relationships between non-fractures and conventional curves, and establish training samples:

[0095] The six conventional curves with high responsiveness to fractures in conventional logging curves are Vp, CNL, RD, Cal, RS, and DEN. From the measurable results of the aforementioned low-frequency experiments, it can be seen that only highly conductive fractures, i.e., fractures that are not completely filled, have a greater impact on well-seismic dispersion. Therefore, it is only necessary to statistically analyze the response relationship between the characteristics of highly conductive fractures and conventional curves. As many sample data as possible are the prerequisite for improving the prediction accuracy. Generally, the sample data should not be less than 100. Here, it is necessary to statistically analyze no less than 60 response values of conventional curves for highly conductive fractures and no less than 40 response values of conventional curves for non-highly conductive fractures.

[0096] ② Normalization processing of sample data:

[0097] To prevent a certain sample feature from being too large or too small, resulting in an unbalanced role in training, and in kernel calculations, inner product operations or exponential operations will be used, and unbalanced data may cause calculation difficulties. Generally, it is necessary to normalize the sample data, and the normalization interval is [-1, 1]. The normalization calculation formula is as follows:

[0098]

[0099] Among them, x i , x imax , x imin respectively represent the original logging value of the i-th conventional curve, the maximum original logging value of the i-th conventional curve, and the minimum original logging value of the i-th conventional curve, and y i represents the normalized value of the i-th conventional curve.

[0100] ③ Selection of kernel function and penalty factor:

[0101] The ultimate goal of the support vector machine is to find a suitable classification function to predict unknown samples. The determination of the classification function mainly involves the selection of the kernel function and the determination of the penalty factor C. These two parameters have a great impact on the accuracy and applicability of the model. The most commonly used kernel functions are linear kernel, polynomial kernel, and Gaussian kernel. In the present invention, the Gaussian kernel is selected as the kernel function of the support vector machine method, as follows:

[0102]

[0103] Among them, m is the imaging logging eigenvalue, d is the order of the polynomial, σ is the width of the Gaussian distribution, and the value of the penalty factor C is 800.

[0104] ④ Discrimination of highly conductive fractures:

[0105] Through the selection of samples, feature normalization, kernel function, and penalty factor, the support vector machine is trained using the training sample set to obtain a discrimination function for identifying highly conductive fractures based on the support vector machine method, as follows:

[0106] n = sgn{∑α i n i k(y i , m) + b} (3)

[0107] where α i is the gradient coefficient, n and n i are the output of the predicted high-conductivity fracture sample and the support vector sample respectively, k(y i , m) is the Gaussian kernel function, and b is the intercept.

[0108] ⑤ Calculate the fracture characteristic parameters:

[0109] After identifying the high-conductivity fractures, the calculation of fracture angle, fracture width, and fracture density can be achieved according to formulas (4, 5, 6, 7, 8):

[0110]

[0111] Y < 0, θ is less than 30 degrees; 0 < Y < 0.1, θ is between 30 - 60 degrees; Y > 0.1, θ is greater than 60 degrees (4)

[0112] where Y is the discrimination factor, θ is the fracture dip angle, and currently only qualitative discrimination can be performed in different regions. R d is the deep lateral resistivity, and R s is the shallow lateral resistivity.

[0113] The fracture porosity (Φ f ) is obtained using the dual lateral resistivity calculation formula: When Rd > Rs:

[0114]

[0115] When Rd < Rs:

[0116]

[0117] The fracture width (ε) in microns is also obtained using the dual lateral resistivity calculation formula: When Rd > Rs:

[0118]

[0119] When Rd > Rs:

[0120]

[0121] In the formula: Rd and Rs are the deep and shallow lateral resistivities, unit: Ω·m; ε is the fracture width, unit: micron, Rmf and Rw are the mud filtrate and formation water resistivities, unit: Ω·m; Φ f is the fracture porosity, unit: %; m fis the fracture cementation index, which is 1.5; for the Ordovician system, Rw is taken as 0.012 Ω·m, R mf = 2.8; R b is the resistivity of the tight limestone formation, with a value of 40000 Ω·m.

[0122] (2) Select core samples of the target fractured-vuggy carbonate reservoir to conduct measurements of lithology, physical properties, and low-frequency rock physics parameters:

[0123] Select core samples of the fractured-vuggy carbonate in the target formation to conduct measurements of parameters such as porosity, permeability, density, and mineral composition. Then, conduct measurements of low-frequency rock physics parameters under formation conditions. The measurement frequency is 12 - 200 Hz, and the saturation states are dry, water-saturated, and oil-saturated. Finally, obtain the data of porosity, permeability, density, mineral composition of the core samples, as well as the low-frequency P-wave and S-wave velocity data of the core samples under formation conditions, providing data support for subsequent verification of the rock physics model.

[0124] (3) Construct a full-frequency carbonate fractured-vuggy rock physics model - the poroelastic equivalent medium model:

[0125] The rock of the fractured-vuggy carbonate reservoir is actually a fluid-containing fractured porous medium. The presence of pores and fluids and the flow between the fluids in the pores and the fluids in the fractures have an important impact on the propagation of seismic waves. At this time, the classic rock physics model, the Gassmann equation, which guides seismic exploration, is no longer valid.

[0126] The poroelastic equivalent medium model is proposed based on the squirt flow mechanism of fractured porous rocks. It assumes that the pore space of the rock consists of randomly isotropically distributed oblate microcracks, spherical equidimensional pores (intercrystalline pores, solution pores, and solution caves can all be equivalent to spherical pores), and oblate fractures arranged in a certain direction. Among them, the radii of the microcracks and pores are the same as the rock grain size, while the radius of the fractures can be much larger than the grain size but smaller than the seismic wave wavelength. The microcracks, pores, and the microcracks and pores, holes can communicate with each other. Each fracture can communicate with multiple microcracks or pores, but each microcrack and each pore can communicate with at most one fracture, and the fractures do not communicate with each other. Figure 1 shows a schematic diagram of the relationship between fractures, microcracks, and pores (caves) in the model. Among them, fractures, microcracks, and pores can be regarded as different inclusions. The random isotropic distribution of microcracks and spherical equidimensional pores and the directional arrangement of fractures make the fractured porous medium described by this model have hexagonal symmetry.

[0127] When there is a wave-induced fluid pressure gradient in a fractured porous medium, this rock physics model considers two scales of wave-induced fluid flow to reach a new fluid pressure equilibrium state. They are the mesoscopic scale (larger than the pore size but smaller than the wavelength, with a typical size of dozens of centimeters) fluid flow occurring between fractures and microfractures or spherical isometric pores, and the particle scale fluid flow occurring between microfractures and spherical pores or between microfractures in different directions. Without considering mesoscopic fractures, this model degenerates into a microscopic jet flow model. The mathematical calculation process of deriving the model formula is omitted, and only the results are recorded.

[0128] The bulk modulus K of the equivalent medium of the poroelastic equivalent medium model eff and the shear modulus μ eff can be expressed as:

[0129]

[0130]

[0131] where ω is the angular frequency, K is the bulk modulus of the skeleton, μ is the shear modulus, ε f is the microfracture density, a is the fracture radius, r is the microfracture radius, ε c is the fracture density, θ is the dip angle, is the porosity, λ is the Lame coefficient, Kc, Kp, A, and B are all intermediate parameters without specific physical meanings. A and B are respectively

[0132]

[0133]

[0134] where,

[0135]

[0136]

[0137] k f is the fluid bulk modulus, τ is the time scale, which controls the frequency band range where the dispersion characteristics appear. It is proportional to the fluid viscosity η, inversely proportional to the permeability k, and related to the fracture radius a. Its reciprocal is called the characteristic frequency. For a smaller aspect ratio:

[0138]

[0139] is the particle scale, ν is the Poisson's ratio of solid particles. According to the formula proposed by Connell and Budiansky, the P-wave velocity V P and the S-wave velocity V S are respectively:

[0140]

[0141]

[0142] Among them, Re is to obtain the real part, and ρ is the density.

[0143] (4) Calibrating the poroelastic equivalent medium model with experimental data:

[0144] Select the core of Well Tuopu 27 for verification, and the test state is water-saturated. From the measured data of lithology and physical properties, it can be known that the calcite content of the core of Well Tuopu 27 is 91%, the clay content is 9%, there are a small number of unfilled fractures visible on the core scale, the angle is about 60 degrees, the width is about 0.5 mm, the density is 2.675 g / cm3, the porosity is 1.26%, the permeability is about 6 mD, and there are microfractures, but they cannot be quantitatively evaluated and are used as fine-tuning parameters.

[0145] Figure 3 、 Figure 4 shows the comparison results of the longitudinal and transverse wave velocities of the core of Well Tuopu 27 with the poroelastic equivalent medium model. From Figure 3 、 4 it can be seen that the measured longitudinal and transverse wave velocities are in good agreement with the longitudinal and transverse wave velocities calculated by the poroelastic equivalent model. The coincidence degree of the longitudinal wave velocity exceeds 90%, and the coincidence degree of the transverse wave velocity exceeds 85%; it is verified that the poroelastic equivalent medium model can better describe the propagation of elastic waves in the full-frequency fracture-vuggy carbonate reservoir, thus laying a foundation for subsequent well-seismic dispersion correction.

[0146] (5) Well-seismic velocity dispersion correction:

[0147] First, obtain the lithology component data through the well logging interpretation result data, and then calculate the frame modulus based on the Voigt-Reuss-Hill model. On this basis, combine the well logging interpretation porosity and the predicted fracture parameters to calculate the dry rock elastic modulus based on the Berrymann model. Then, further combine the saturation data obtained through well logging interpretation, and calculate the longitudinal wave velocity at the well logging frequency based on the poroelastic equivalent medium model. If the absolute difference between the calculated longitudinal wave velocity and the measured longitudinal wave velocity is less than the error, then the longitudinal and transverse wave velocities in the seismic frequency band can be directly calculated at this time. If the difference between the calculated longitudinal wave velocity and the measured longitudinal wave velocity is greater than the error, then after adjusting the microfracture density, calculate again according to the above process, and repeat the cycle until the absolute difference between the calculated longitudinal wave velocity and the measured longitudinal wave velocity is less than the error. The flow chart is as Figure 6 。

[0148] ① Calculation of frame modulus - Voigt-Reuss-Hill average theory:

[0149] Given the relative content and elastic modulus of each phase of the known rock medium, the upper and lower limits of the effective elastic modulus of the rock medium can be calculated using the Voigt and Reuss formulas (Equation 13). Among them, the upper limit Voigt formula represents the same strain state, where for each phase composing the rock medium under the same strain, the ratio of stress to strain of the rock medium; the Reuss lower limit represents the same stress state, where for each phase composing the rock medium under the same stress, the ratio of stress to strain of the rock medium. The Hill formula is the arithmetic mean result of the Voigt upper limit and the Reuss lower limit.

[0150]

[0151] In the formula, M V , M R , M VRH represent the values of the bulk modulus K, shear modulus μ, and Young's modulus E obtained by the Voigt, Reuss, and Hill methods; f i is the volume content of the i-th component of the rock medium; M i is the elastic modulus of the i-th component.

[0152] ② Calculation of the elastic modulus of dry rock:

[0153] Berryman (1980b, 1995) gave the general form of the self-consistent approximation for N-phase composites:

[0154]

[0155]

[0156] In the i-th material, x i is its volume fraction, and P and Q represent geometric factors. The superscript i on P and Q indicates that these factors are in the background medium with self-consistent effective moduli and that contains material i. k i , u i are the bulk modulus and shear modulus of the i-th mineral respectively. The summation represents considering all phases, including minerals and pores. Dry pores can be simulated by setting the inclusion modulus to zero. Fluid-saturated pores can be simulated by setting the inclusion shear modulus to zero. At the same time, these equations are coupled and must be solved by simultaneous iteration.

[0157] ③ Calculation of the full-frequency compressional and shear wave velocities of fluid-containing rocks:

[0158] After calculating the elastic modulus of dry rocks, combining with logging interpretation result data such as porosity, permeability, fluid saturation and other parameters as well as predicting fracture characteristic parameters, the calculation of full-frequency P-wave and S-wave velocities can be realized based on the poroelastic equivalent medium model.

[0159] Taking the western contiguous work area of Tahe Oilfield as an example, the well-seismic dispersion correction effect is analyzed in detail. There is an FMI imaging data well, Aiding 6 Well, in the western contiguous area of Tahe Oilfield. Figure 7 Figure 5 is the FMI imaging interpretation fracture indication diagram of Aiding 6 Well. The red tadpoles indicate the development of highly conductive fractures, the yellow tadpoles represent high-resistance fractures, and the green tadpoles represent the formation interfaces. The highly conductive fractures are the fractures that need to be identified. Therefore, the response relationships between highly conductive fractures, non-highly conductive fractures and conventional curves of this well are selected for statistics to establish training samples. Then, based on the support vector machine method, part of the above samples are normalized and pattern recognition training is carried out, and the remaining part of the samples are used as test samples for verification, and the correct rate is 85%, as shown in Table 1 below.

[0160] Table 1 Test of Remaining Samples

[0161] DEN CNL RD RS AC CAL Fracture Mark Prediction Mark 2.689 0.578 2579.08 1789.65 48.791 6.196 1 1 2.686 0.617 2433.83 1708.07 48.986 6.195 1 1 2.675 0.831 1765.15 1309.97 49.831 6.19 1 1 2.674 1.o63 1346.1 1015.56 50.376 6.188 1 1 2.673 1.216 1166.11 874.172 50.824 6.189 1 1 2.669 1.349 1112.96 826.947 50.982 6.19 1 1 2.653 1.881 952.944 681.173 51.57 6.194 1 1 2.644 2.114 1008.11 721.554 51.735 6.193 1 1 2.64 2.171 1162.53 846.59 51.764 6.194 1 1 2.639 2.179 1285.23 928.449 51.71 6.194 1 0 2.704 0.138 2571.76 2205.02 51.181 6.065 1 1 2.683 0.156 2035.04 1796.7 50.905 6.081 1 1 2.67 0.168 1790.62 1595.12 50.685 6.096 1 1 2.663 0.165 1675.45 1475.49 50.528 6.105 1 1 2.666 o.177 1585.36 1403.27 50.47 6.109 1 1 2.681 0.231 1238.63 1122.33 50.226 6.132 1 0 2.691 0.266 1107.29 1005.04 50.059 6.153 1 1 2.69 0.317 1000.28 899.722 49.877 6.179 1 1 2.688 0.338 999.646 894.626 49.735 6.19 1 0 2.684 0.429 1028.75 910.01 49.15 6.198 1 1 2.693 0.477 1120.89 1004.7 48.853 6.116 1 1 2.698 0.493 1223.97 1114.31 48.644 6.04 1 1 2.698 0.496 1454.88 1353.08 48.589 6.04 1 1 2.696 0.511 2437.22 2357.98 48.412 6.04 1 1 2.695 0.52 3062.19 2975.33 48.435 6.042 1 0 2.697 0.525 3820.59 3718.05 48.527 6.044 1 1 2.695 0.517 4033.55 3936.17 48.594 6.046 1 1 2.686 0.497 4515.94 4474.83 48.872 6.055 1 1 2.684 0.527 3871.68 3765.75 49.016 6.052 1 0

[0162] Then, based on the fracture prediction results, three wells in the work area are selected to carry out well-seismic dispersion correction. Figure 8 Figure 16 is the dispersion correction effect diagram of Well No. 1. Figure 9 Figure 18 is the dispersion correction effect diagram of Well No. 2. Figure 10 Figure 20 is the dispersion correction effect diagram of Well No. 3. Observe Figure 8 、 9, 10. Each figure includes 6 sets of well data. The first set is the lithological volume content, where the black-filled area is the shale content and the red-filled area is the limestone content; the second set is the total porosity; the third set is the predicted fracture area, and the red-filled part is the fracture porosity content; the fourth set is the water saturation curve, where the black-filled area is the water saturation and the red-filled area is the oil saturation; the fifth set, the red curve represents the calculated shear wave velocity curve at the seismic frequency band (30 Hz), and the black curve represents the calculated shear wave velocity curve at the logging frequency band (10,000 Hz); the sixth set is the calculated compressional wave velocity curve at the logging frequency band and the compressional wave velocity curve at the seismic frequency band. It can be seen from the figure that when there are fractures in the pores, there is dispersion in the compressional wave velocity and shear wave velocity at the logging frequency band, and the magnitude of dispersion is closely related to the proportion of fractures and the size of the total porosity. For example, at the depth of 6204 - 6209 meters in Well No. 1, although the proportion of fractures is relatively large, the total porosity is very small, less than 1%, so the dispersion amplitude is small. The compressional wave dispersion amplitude is about 1%, and the shear wave velocity dispersion amplitude is about 0.6%; at the depth of 6660 - 6680 meters in Well No. 2, the proportion of fractures is relatively small, but due to the large total porosity, with an average greater than 3%, the compressional wave velocity dispersion amplitude is also close to 1%, and the shear wave velocity dispersion amplitude is 0.6%; at the depth of 6650 - 6670 meters in Well No. 3, the proportion of fractures is relatively high, and at the same time the total porosity is also relatively large, about 10%, so the compressional wave velocity dispersion amplitude in this section is also relatively large, reaching about 3%, and the shear wave dispersion amplitude is about 1.3%. Through the analysis of the above three wells, it can be obtained that the dispersion amplitude is proportional to the proportion of fractures and the total porosity. Generally speaking, the Ordovician carbonate reservoir is buried relatively deep, and the formation pressure is relatively large. Large-scale fractures are either filled with minerals or the fracture diameter is very small. Therefore, overall, the dispersion amplitude is small. The maximum compressional wave velocity dispersion amplitude is less than 3%, and the maximum shear wave velocity dispersion amplitude is less than 1.3%.

[0163] The above description is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any change or replacement that can be thought of without creative work should be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope defined by the claims.

Claims

1. A well-seismic velocity dispersion correction method suitable for fractured-vuggy carbonate reservoirs, characterized in that It includes the following steps: (1) Based on the characterization of fractures and vugs by imaging logging FMI, carry out fracture prediction of conventional logging curves based on support vector machines and calculate their fracture characteristic parameters; (2) Select cores of the target fracture-vuggy carbonate reservoir to carry out measurements of lithology, physical properties, and low-frequency rock physical parameters; (3) Construct a full-frequency carbonate fracture-vuggy rock physical model - a poroelastic equivalent medium model; (4) Calibrate the poroelastic equivalent medium model with experimental data; (5) Well-seismic velocity dispersion correction.

2. The method according to claim 1, wherein Step (1) specifically includes the following steps: ① Select wells in the work area containing FMI imaging logging data, statistically analyze the fracture or vug characteristics and the response relationships between non-fractures and conventional curves, and establish training samples; ② Normalize the sample data; ③ Select the kernel function and penalty factor; ④ Identify high-conductivity fractures; ⑤ Calculate fracture characteristic parameters.

3. The method according to claim 2, wherein In step ①, statistically analyze no less than 60 conventional curve response values of high-conductivity fractures and no less than 40 conventional curve response values of non-high-conductivity fractures.

4. The method according to claim 2, wherein In step ②, the normalization interval is [-1, 1]; the normalization calculation formula is as follows: where x i , x i max , x i min respectively represent the original logging value, the original maximum logging value, and the original minimum logging value of the i-th conventional curve, and y i represents the normalized value of the i-th conventional curve.

5. The method according to claim 2, wherein In step ③, select the Gaussian kernel as the kernel function of the support vector machine method, as follows: where m is the imaging logging eigenvalue, d is the order of the polynomial, σ is the width of the Gaussian distribution, and the value of the penalty factor C is 800.

6. The method according to claim 2, wherein In step ④, through the selection of samples, feature normalization, kernel function, and penalty factor, use the training sample set to train the support vector machine to obtain a discriminant function for high-conductivity fracture identification based on the support vector machine method, as follows: n = sgn{∑α i n i k(y i , m) + b} where α i is the gradient coefficient, n, n i are respectively the output of the predicted highly conductive fracture sample and the support vector sample, k(y i , m) is the Gaussian kernel function, and b is the intercept.

7. The method according to claim 2, characterized in that In step ⑤, Calculation of fracture angle, fracture width, and fracture density: When Y < 0, θ is less than 30 degrees; when 0 < Y < 0.1, θ is between 30 - 60 degrees; when Y > 0.1, θ is greater than 60 degrees; where Y is the discrimination factor, θ is the fracture dip angle, and qualitative discrimination is carried out in different regions. R d is the deep lateral resistivity, and R s is the shallow lateral resistivity; Fracture porosity Φ f Obtained from the dual laterolog resistivity calculation formula: When Rd > Rs: When Rd < Rs: The fracture width (ε) in microns is obtained using the dual laterolog resistivity calculation formula: when Rd > Rs: When Rd > Rs: Where: Rd and Rs are the deep and shallow lateral resistivity, unit: Ω·m; ε is the fracture width, unit: micron, Rmf and Rw are the resistivity of mud filtrate and formation water, unit: Ω·m; Φ f is the fracture porosity, unit: %; m f is the fracture cementation index, which is 1.5; for the Ordovician system, Rw is taken as 0.012 Ω·m, R mf = 2.8; R b is the resistivity of the tight limestone layer, with a value of 40000 Ω·m.

8. The method according to claim 1, wherein In step (2), select cores of the fracture-vuggy carbonate rock in the target formation to carry out parameter measurements including porosity, permeability, density, and mineral composition, and then carry out low-frequency rock physical parameter measurements under formation conditions. The measurement frequency is 12 - 200 Hz, and the saturation states are dry, water-saturated, and oil-saturated, to obtain the porosity, permeability, density, mineral composition data of the cores, and the low-frequency P-wave and S-wave velocity data of the cores under formation conditions.

9. The method according to claim 1, wherein In step (3), The bulk modulus K of the equivalent medium of the poroelastic equivalent medium model eff and the shear modulus μ eff are expressed as: where ω is the angular frequency, K is the bulk modulus of the skeleton, μ is the shear modulus, ε f is the microcrack density, a is the crack radius, r is the microcrack radius, ε c is the crack density, θ is the dip angle, is the porosity, λ is the Lame coefficient, Kc, Kp, A, and B are all intermediate parameters without specific physical meanings. A and B are respectively where, k f where \(k\) is the fluid bulk modulus, \(\tau\) is the time scale that controls the frequency band where the dispersion characteristics appear, is proportional to the viscosity \(\eta\) of the fluid, inversely proportional to the permeability \(k\), and related to the fracture radius \(a\). Its reciprocal is called the characteristic frequency; For a relatively small aspect ratio: where \(l\) is the particle scale and \(\nu\) is the Poisson's ratio of the solid particles; according to the formula proposed by Connell and Budiansky, the P-wave velocity \(V_{P}\) P and the S-wave velocity \(V_{S}\) S are respectively: where Re is to find the real part and ρ is the density.

10. The method according to claim 1, wherein In step (5), First, obtain lithology component data from well logging interpretation result data. Then, calculate the frame modulus based on the Voigt-Reuss-Hill model. On this basis, combine the well logging interpreted porosity and predicted fracture parameters to calculate the dry rock elastic modulus based on the Berrymann model. Then, further combine the saturation data interpreted by well logging, and calculate the longitudinal wave velocity at well logging frequency based on the poroelastic equivalent medium model. If the absolute difference between the calculated longitudinal wave velocity and the measured longitudinal wave velocity is less than the error, then directly calculate the longitudinal and transverse wave velocities in the seismic frequency band. If the difference between the calculated longitudinal wave velocity and the measured longitudinal wave velocity is greater than the error, then adjust the microfracture density and calculate again according to the above process, repeating the cycle until the absolute difference between the calculated longitudinal wave velocity and the measured longitudinal wave velocity is less than the error; ① Frame modulus calculation - Voigt-Reuss-Hill averaging theory: Given the relative content and elastic modulus of each phase of the rock medium, use the Voigt and Reuss formulas to calculate the upper and lower limits of the effective elastic modulus of the rock medium. Among them, the upper limit Voigt formula represents the same strain state, where each phase of the rock medium has the same strain, and the ratio of stress to strain of the rock medium; The Reuss lower limit represents the same stress state, where each phase of the rock medium has the same stress, and the ratio of stress to strain of the rock medium. The Hill formula is the arithmetic mean result of the Voigt upper limit and the Reuss lower limit; where M V , M R , M VRH represent the values of the bulk modulus K, shear modulus μ, and Young's modulus E obtained by the Voigt, Reuss, and Hill methods; f i is the volume content of the i-th component that makes up the rock medium; M i is the elastic modulus of the i-th component; ② Elastic modulus calculation of dry rock: The Berrymann model is the general form of the self-consistent approximation of an N-phase composite material: In the i-th material, x i is its volume fraction, and P and Q represent geometric factors. The superscript i on P and Q indicates that these factors are in the background medium with self-consistent effective moduli and that contain material i; k i , u i are the bulk modulus and shear modulus of the i-th mineral, respectively. The inclusion modulus is set to zero to simulate dry pores; ③ Calculation of longitudinal and transverse wave velocities in the full frequency band of fluid-bearing rock: After calculating the elastic modulus of dry rock, combine parameters such as porosity, permeability, and fluid saturation from well logging interpretation result data and predicted fracture characteristic parameters, and the longitudinal and transverse wave velocities in the full frequency band can be calculated based on the poroelastic equivalent medium model.

Images