Fluid identification method based on GBP equation and array sound wave data time-frequency analysis

Abstract

The invention discloses a fluid identification method based on a GBP equation and array sound wave data time-frequency analysis, and solves the problem that a traditional saturation formula cannot accurately calculate water saturation. The method comprises the steps that data are preprocessed, borehole correction is conducted, and a G-B-P equation is established to describe the relation between igneous rock fluid saturation and acoustic parameters; carrying out finite element numerical simulation on the in-well array acoustic logging information based on a G-B-P equation to obtain an acoustic full wave train signal; performing time-frequency analysis on the full-wave data to obtain IMF component information; determining accumulated energy of an IMF modal component, and performing porosity correction on the accumulated energy of the IMF modal component; and proposing a fluid identification factor, establishing a fluid discrimination graph, and identifying fluid properties according to the fluid discrimination graph. According to the method, fluid identification is carried out by using non-resistivity data, the problem of alteration of complex conductivity of igneous rock is avoided, and research on igneous rock reservoir prediction and CO2 storage monitoring is facilitated.

Description

Technical Field

[0001] The invention relates to the technical field of natural gas exploration, and in particular to a fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis. Background Art

[0002] In recent years, igneous rock reservoirs have been the focus of natural gas exploration and development. Igneous rock reservoirs can help solve energy needs. In addition, when oil and gas reservoirs are exhausted, volcanic rock reservoirs provide good storage space, transportation channels and storage capacity, which is very beneficial for the storage of CO2. Fluid identification technology is of great significance to reservoir prediction and carbon dioxide reserves monitoring.

[0003] The lithology and mineral composition, porosity and pore structure characteristics of igneous rocks will affect the distribution and saturation of fluids. In altered igneous reservoirs, the complex mechanism of rock conductivity makes it impossible for traditional saturation formulas to accurately calculate water saturation. Resistivity data is affected by the combined effects of lithology, physical properties and fluid properties, and high-resistance water layers and low-resistance gas layers will affect fluid identification in altered igneous reservoirs.

[0004] This method is a new fluid identification method based on the Gassmann-Brie-Patchchy (GBP) saturation calculation formula and array acoustic logging. The results of practical application show that this method has a good application effect in fluid identification of igneous gas reservoirs. Summary of the invention

[0005] The purpose of the present invention is to provide a fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis to solve the technical problem that the traditional saturation formula cannot accurately calculate the water saturation in altered igneous rock reservoirs.

[0006] To achieve the above object, the present invention provides the following technical solutions: The present invention provides a fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis, which specifically includes the following steps: S1: Preprocess array acoustic wave and conventional logging data, perform borehole correction, and establish GBP equation based on core experiments to describe the relationship between igneous rock fluid saturation and acoustic parameters; S2: Based on the GBP equation, finite element numerical simulation of the array acoustic logging data in the well was performed to obtain the full acoustic wave train signals under different lithology, porosity and water saturation conditions; S3: The ICEEMDAN modal decomposition method is used in combination with SPWVD to perform time-frequency analysis on the full-wave data measured by the array acoustic wave to obtain the IMF component information; S4: determining the cumulative energy of the IMF modal component, and performing porosity correction on the cumulative energy of the IMF modal component; S5: Considering the irreducible water saturation, mineral composition and physical parameters of core analysis, as well as the influence of porosity and lithology, the fluid identification factor FII is proposed and the fluid discrimination diagram is established. S6: Identify fluid properties based on the fluid discrimination diagram, including calculating the FII critical value of standard gas and water formations from the density logging value at each depth point.

[0007] Furthermore, in step S1, the establishment of the GBP equation is specifically shown in the following formula: (1), In the formula, K pat is the rock bulk modulus of the patch model; K sat is the rock bulk modulus of the uniform model; S w is water saturation; K is the rock bulk modulus, which tends to be patchy in the high saturation stage and to be uniform in the medium and low saturation stages; f is the fitting coefficient, f≥1, and the calculated value of K is adjusted to be between the uniform fluid distribution model and the patch model; In the Brie equation, e is an empirical coefficient related to the connectivity and distribution uniformity of pore fluids, e ≥ 1, (2), In the formula, K fl is the volume weighted average of the air-water modulus, K g and K w are the bulk moduli of gas and water, respectively.

[0008] Furthermore, in step S2, the device used in the finite element numerical simulation is set to establish an isotropic formation model under an axisymmetric configuration, and then an acoustic logging tool with a transmitter and two receivers is placed in the wellbore.

[0009] Furthermore, in step S3, the ICEEMDAN modal decomposition method is specifically implemented as follows: S31: adding Gaussian white noise to the full-wave signal of the array sound wave to be decomposed to obtain a noisy signal; S32: Calculate each noisy signal The difference between the first IMF and the average value is determined to obtain the residual r1 of the first decomposition; S33: subtract the first residual r1 from the original signal x(n) to obtain a first component M1 of the original signal; S34: when i is greater than or equal to 2, constructing an i-th signal N containing controllable noise; S35: Calculate each The difference between the i-th IMF and the average value is then calculated to obtain the residual r of the i-th decomposition i ; S36: Using residual r i Subtract the residual r i-1 , get the i-th component M of the original signal i ; S37: Set i=i+1, and return to step S34 to calculate the next value until all IMF components are obtained.

[0010] Furthermore, in step S4, the boundary conditions are as follows: The frequency boundaries of a time-frequency distribution are defined as follows: (14), Where M t (f) is the frequency edge, and its physical meaning is the energy spectrum density of the signal. SPWV(t,f) represents SPWVD. The time boundary is defined as: (15), Where M f (t) represents the time edge of the time-frequency distribution. The physical meaning of the time edge of the time-frequency distribution represents the instantaneous energy of the signal. The cumulative energy of the time-frequency distribution is defined as: (16), Where E tf It is the cumulative energy of the time-frequency distribution, which represents the total energy of the signal.

[0011] Furthermore, the porosity of IMF has an influence on the cumulative energy E tf The specific formula for correcting the influence of is: IMF1:E tfIMF1C =E tfIMF1 • (17), IMF2: (18), In the formula, E tfIMF1C and E tfIMF2C are the corrected accumulated energies of IMF1 and IMF2, is the porosity.

[0012] Furthermore, in step 5, the specific expression of the fluid identification factor FII is as follows: (19), Where FII is a dimensionless parameter and is positively correlated with the cumulative energy of IMF1 after removing the influence of porosity, and negatively correlated with the cumulative energy of IMF2 after removing the influence of porosity.

[0013] Based on the above technical solution, the embodiments of the present invention can at least produce the following technical effects: (1) The fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis provided by the present invention is more accurate in describing the relationship between mixed fluid and sound velocity value for igneous rock cores by using GBP equation. Finite element simulation of in-well array acoustic wave logging data can obtain acoustic full wave logging data of igneous rocks and pyroclastic rock layers with different water saturation and porosity. Then, after comparing the effects of different modal decomposition methods, it was determined that the ICEEMDAN method is most effective in eliminating modal mixing. Therefore, the ICEEMDAN-SPWVD time-frequency analysis method is used to process the full wave data. Comparing the time-frequency analysis results, it is found that the influence of porosity and lithology density on the cumulative energy amplitude after modal decomposition is significantly greater than the influence of fluid saturation. These results show that before making a fluid discrimination diagram based on the difference in fluid properties, the influence of porosity and lithology should be considered. After correcting the influence of porosity, the cumulative energy values ​​at different porosity levels are basically within one order of magnitude, which provides the necessary conditions for subsequent research on the influence of saturation on cumulative energy.

[0014] (2) The fluid identification method provided by the present invention based on the GBP equation and the time-frequency analysis of array acoustic wave data proposes the fluid identification factor FII and constructs the FII-Den fluid identification map. Through practical application, it is found that the value of the water layer is significantly higher than that of the gas layer, so better fluid identification can be performed through the image value. In addition, the identification map can be used to effectively identify the properties of the fluid.

[0015] (3) The fluid identification method based on the GBP equation and array acoustic wave data time-frequency analysis provided by the present invention can be seen from the results of digital core electrical simulation that the complex conductivity mechanism in altered igneous rocks makes it difficult to accurately calculate the water saturation using the traditional saturation calculation formula. Therefore, the present invention uses non-resistivity data for fluid identification, bypassing the problem of complex conductivity of altered igneous rocks. This is a new method for fluid identification in igneous reservoirs. The method was applied to igneous reservoirs in the southwest region, and the application results were consistent with the experimental conclusions, verifying the effectiveness of the proposed method. This method is helpful for the research on igneous reservoir prediction and CO2 storage monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the structures shown in these drawings without paying creative work.

[0017] Figure 1 The figures are the calculation results of different fluid distribution models of the present invention; (a) is the calculation result of GBP equation under different f values; (b) is the comparison of GB formula, Patchy formula, GBP equation and experimental data of A-1 sample; (c) is the comparison of GB formula, Patchy formula, GB-P equation and experimental data of A-5 sample; Figure 2 is a relationship diagram between the bulk modulus of dry rock and the porosity of different lithologies in the present invention; Figure 3 It is the FII-DEN fluid discrimination diagram of the igneous rock layer of the present invention; Figure 4 It is the result of the identification of igneous reservoir fluids in the southwest region of the present invention. DETAILED DESCRIPTION

[0018] The technical solutions in the embodiments of the present invention will be described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. In addition, the technical solutions between the various embodiments can be combined with each other, but they must be based on the ability of ordinary technicians in this field to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0019] This embodiment provides a fluid identification method based on the GBP equation and array acoustic wave data time-frequency analysis, which specifically includes the following steps: S1: Preprocess array acoustic and conventional logging data, perform borehole correction, and establish the GBP equation based on core experiments to describe the relationship between igneous rock fluid saturation and acoustic parameters.

[0020] The establishment of GBP equation is as follows: (1), In the formula, K pat is the rock bulk modulus of the patch model; K sat is the rock bulk modulus of the homogeneous model (GB equation); Sw is water saturation; K is the rock bulk modulus, which tends to be patchy in the high saturation stage and to be uniform in the medium and low saturation stages; f is the fitting coefficient, f≥1, and the calculated value of K is adjusted to be between the uniform fluid distribution model and the patch model; In the Brie equation, e is an empirical coefficient related to the connectivity and distribution uniformity of pore fluid, e≥1.

[0021] (2), In the formula, K fl is the volume weighted average of the air-water modulus, K g and K w are the bulk moduli of gas and water, respectively.

[0022] like Figure 1 As shown in the figure, the calculation results of different fluid distribution models in this embodiment include the calculation results of GBP equation (f=2, 5, 10), GB equation (e=8); and Patchy equation. Obviously, the calculation results of different models are very different. In the GBP equation, the value of f will affect the calculation results. Figure 1 (a) and Figure 1 (b) It can be seen that the experimental results are in good agreement with the calculation results of the GBP equation (f=5). Compared with the GB equation, the GBP equation is more consistent with the experimental results in the high water saturation range. These findings make the GBP equation more suitable for tight reservoirs with high water saturation, such as igneous rocks and tight sandstone reservoirs. It provides a reliable basis for simulating the full wave train data of different water saturations.

[0023] S2: Based on the GBP equation, finite element numerical simulation of the in-well array acoustic logging data was performed to obtain the full acoustic wave train signals under different lithology, porosity and water saturation conditions.

[0024] The finite element numerical simulation device is set to establish an isotropic formation model in an axisymmetric configuration, and then place an acoustic logging tool with a transmitter and two receivers in the wellbore; the distance between the transmitter and the first receiver is set to 2m, and the distance between the receivers is 1m. The borehole radius is 0.1m, and the borehole is filled with water. The formation around the borehole is 6 meters high and 4 meters wide. Due to the complex mineral composition of igneous rocks and pyroclastic rocks, the effective average modulus and density are used to represent the modulus and density of the rock skeleton, which specifically includes the following steps: S21: Measurement of sound velocity in basic igneous and pyroclastic rocks at different saturation states; S22: Based on the sound velocity value, the bulk modulus of dry rock is calculated using the following formula: (3), In the formula, K dry is the bulk modulus of dry rock, Vp and Vs represent the longitudinal wave velocity and shear wave velocity, respectively, and ρdry refers to the density of dry rock.

[0025] like Figure 2 As shown, when fitting the relationship between the two parameters of the porosity and bulk modulus information of the rock sample, it is assumed that the rock is completely composed of pores, which means that the porosity is equal to 1 and the bulk modulus of the dry rock is equal to the bulk modulus of the natural gas. It is observed that the bulk modulus of the dry sample has a good correlation with the porosity, which can be obtained: Basic igneous rock dry rock bulk modulus: K dry =0.3492 -1.056 (4), Pyroclastic dry rock bulk modulus: K dry =0.3359 -1.056 (5), In the formula, is the porosity of the sample. On the fluid-structure-acoustic coupling boundary, the continuous conversion of fluid pressure energy to structural deformation is considered, including fluid load and structural acceleration borne by the fluid. Therefore, from a mathematical point of view, its boundary condition can be expressed as: (6), Among them, u tt is the solid acceleration, n is the surface normal direction, p t is the total sound pressure, F A Represents the load acting on the structure.

[0026] S23: Set the sound source As for the sound source, a circuit is provided in this embodiment, which is connected to the piezoelectric device, drives the device and receives the signal by using the piezoelectric effect. The excitation voltage input is determined as follows: (7), Wherein, t is time, T0 is the reciprocal of frequency f0, and frequency f0 is 10kHz.

[0027] S3: The ICEEMDAN modal decomposition method is used in combination with SPWVD to perform time-frequency analysis on the full-wave data measured by the array acoustic wave to obtain the IMF1 and IMF2 component information.

[0028] The specific implementation steps of the ICEEMDAN algorithm are as follows: S31: Let x(n) be the full-wave signal of the array acoustic wave to be decomposed, define Ek(·) as the kth IMF component obtained by EMD, and then add it to the original signal as special noise, as shown below: (8), In the formula, w (n) is the nth Gaussian white noise to be added, α1 is the expected signal-to-noise ratio of the first decomposition signal; is the signal after adding noise; S32: Calculate each The difference between the first IMF and the average value is determined, and the residual r1 of the first decomposition is obtained as follows: (9), Where <·> is the operator for calculating the average value of N signals; S33: Subtract the first residual r1 from the original signal x(n) to obtain the first component of the original signal, denoted as M1: M1=x(n)-r1(10) S34: When i is greater than or equal to 2, construct the i-th signal N containing controllable noise as shown below: (11), S35: Calculate each The difference between the i-th IMF and the i-th IMF is then averaged to obtain the residual r of the i-th decomposition i as follows: (12), S36: Then use the residual r i Subtract the residual r i-1 , we get the i-th component of the original signal, denoted as M i : M i =r i-1 -r i (13), S37: Set i=i+1, and return to step S34 to calculate the next value until all IMF components are obtained.

[0029] S4. Determine the cumulative energy of the modal components of IMF1 and IMF2, and perform porosity correction on the cumulative energy of the modal components of IMF1 and IMF2.

[0030] The frequency boundaries of a time-frequency distribution are defined as follows: (14), Where M t (f) is the frequency edge, and its physical meaning is the energy spectrum density of the signal. SPWV(t,f) represents SPWVD.

[0031] The time boundary is defined as: (15), Where M f (t) represents the time edge of the time-frequency distribution. The physical meaning of the time edge of the time-frequency distribution represents the instantaneous energy of the signal.

[0032] The cumulative energy of the time-frequency distribution is defined as: (16), Where E tf It is the cumulative energy of the time-frequency distribution, which represents the total energy of the signal.

[0033] The porosity of IMF1 and IMF2 has an influence on the accumulated energy E tf The specific formula for correcting the influence of is: IMF1:E tfIMF1C =E tfIMF1 • (17), IMF2: (18), In the formula, E tfIMF1C and E tfIMF2C are the corrected accumulated energies of IMF1 and IMF2, is the porosity.

[0034] S4. Considering the irreducible water saturation, mineral composition and physical parameters of core analysis, as well as the influence of porosity and lithology, the fluid identification factor FII is proposed and the fluid discrimination diagram is established.

[0035] In this embodiment, the FII-DEN fluid discrimination diagram for igneous rock layers is established. Figure 3 shown.

[0036] The specific expression of fluid identification factor FII is as follows: (19), Where FII is a dimensionless parameter and is positively correlated with the cumulative energy of IMF1 after removing the influence of porosity, and negatively correlated with the cumulative energy of IMF2 after removing the influence of porosity.

[0037] S5. Identify fluid properties based on the fluid discrimination diagram, including calculating the FII critical value of standard gas and water formations from the density logging value at each depth point.

[0038] In this embodiment, a fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis of the present invention is actually applied to the fluid identification results of igneous reservoirs in the southwest region, such as Figure 4As shown in the figure, if the FII value calculated by the array acoustic wave is greater than or equal to the water line (blue) in the figure, it can be judged that the layer is a water layer. However, if the value is less than the natural gas line (orange) in the figure, the layer is judged to be a natural gas layer. In addition, if the calculated FII value is between the two lines, it is estimated that the layer is a gas-water co-layer.

[0039] The above shows and describes the basic principles and main features of the present invention and the advantages of the present invention. It should be understood by those skilled in the art that the present invention is not limited to the above embodiments, and the above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, and these changes and improvements fall within the scope of the present invention to be protected. The scope of the present invention to be protected is defined by the attached claims and their equivalents.

Claims

1. A fluid identification method based on GBP equation and time-frequency analysis of array acoustic wave data, characterized in that: The specific steps include: S1: Preprocess array acoustic wave and conventional logging data, perform borehole correction, and establish GBP equation based on core experiments to describe the relationship between igneous rock fluid saturation and acoustic parameters; S2: Based on the GBP equation, finite element numerical simulation of the array acoustic logging data in the well was performed to obtain the full acoustic wave train signals under different lithology, porosity and water saturation conditions; S3: The ICEEMDAN modal decomposition method is used in combination with SPWVD to perform time-frequency analysis on the full-wave data measured by the array acoustic wave to obtain the IMF component information; S4: determining the cumulative energy of the IMF modal component, and performing porosity correction on the cumulative energy of the IMF modal component; S5: Considering the irreducible water saturation, mineral composition and physical parameters of core analysis, as well as the influence of porosity and lithology, the fluid identification factor FII is proposed and the fluid discrimination diagram is established. S6: Identify fluid properties based on the fluid discrimination diagram, including calculating the FII critical value of standard gas and water formations from the density logging value at each depth point.

2. The fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis according to claim 1 is characterized in that: In step S1, the establishment of the GBP equation is specifically shown in the following formula: （1）， In the formula, K pat is the rock bulk modulus of the patch model; K sat is the rock bulk modulus of the uniform model; S w is water saturation; K is the rock bulk modulus, which tends to be patchy in the high saturation stage and to be uniform in the medium and low saturation stages; f is the fitting coefficient, f≥1, and the calculated value of K is adjusted to be between the uniform fluid distribution model and the patch model; In the Brie equation, e is an empirical coefficient related to the connectivity and distribution uniformity of pore fluids, e≥1, （2）， In the formula, K fl is the volume weighted average of the air-water modulus, K g and K w are the bulk moduli of gas and water, respectively.

3. The fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis according to claim 1 is characterized in that: In step S2, the device used in the finite element numerical simulation is set to establish an isotropic formation model under an axisymmetric configuration, and then an acoustic logging tool with a transmitter and two receivers is placed in the wellbore.

4. The fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis according to claim 1 is characterized in that: In step S3, the specific implementation steps of the ICEEMDAN modal decomposition method are as follows: S31: adding Gaussian white noise to the full-wave signal of the array sound wave to be decomposed to obtain a noisy signal; S32: Calculate each noisy signal The difference between the first IMF and the average value is determined to obtain the residual r1 of the first decomposition; S33: subtract the first residual r1 from the original signal x(n) to obtain a first component M1 of the original signal; S34: when i is greater than or equal to 2, constructing an i-th signal N containing controllable noise; S35: Calculate each The difference between the i-th IMF and the average value is then calculated to obtain the residual r of the i-th decomposition i ; S36: Using residual r i Subtract the residual r i-1 , get the i-th component M of the original signal i ; S37: Set i=i+1, and return to step S34 to calculate the next value until all IMF components are obtained.

5. The fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis according to claim 1 is characterized in that: In step S4, the boundary conditions are as follows: The frequency boundaries of the time-frequency distribution are defined as follows: （14）， Where M t (f) is the frequency edge, and its physical meaning is the energy spectrum density of the signal. SPWV(t,f) represents SPWVD. The time boundary is defined as: （15）， Where M f (t) represents the time edge of the time-frequency distribution. The physical meaning of the time edge of the time-frequency distribution represents the instantaneous energy of the signal. The cumulative energy of the time-frequency distribution is defined as: （16）， Where E tf It is the cumulative energy of the time-frequency distribution, which represents the total energy of the signal.

6. The fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis according to claim 5 is characterized in that: Effect of porosity on the accumulated energy E of IMF tf The specific formula for correcting the influence of is: IMF1:E tfIMF1C =AND tfIMF1 • (17), IMF2: (18), In the formula, E tfIMF1C and E tfIMF2C are the corrected accumulated energies of IMF1 and IMF2, is the porosity.

7. The fluid identification method based on GBP equation and array acoustic wave data time-frequency analysis according to claim 1 is characterized in that: In step 5, the specific expression of the fluid identification factor FII is as follows: （19）， Where FII is a dimensionless parameter and is positively correlated with the cumulative energy of IMF1 after removing the influence of porosity, and negatively correlated with the cumulative energy of IMF2 after removing the influence of porosity.

Images