A method for constructing a dual-lithology comprehensive anisotropic equivalent medium model

Abstract

The present invention discloses a method for constructing a dual-lithology comprehensive anisotropic equivalent medium model, comprising the following steps: S1. constructing a sandstone skeleton; S2. constructing a shale skeleton; S3. combining sandstone and shale skeletons and supplementing fluids. Considering that the previous anisotropic models were very purposeful in analyzing the elasticity and anisotropic characteristics of rocks based on the microstructure of fractures or mudstone, the actual underground medium is often mixed with muddy bands, muddy fillings, and other conditions with both shale and sandstone, the accuracy of the rock physics model is mainly related to equivalent equivalent conditions, and a single model is limited to study the elasticity and anisotropic characteristics of the actual underground. Therefore, the present invention constructs a comprehensive anisotropic model that can characterize the anisotropy of weak rocks and the anisotropy of fractures in hard rocks in response to the above defects, and has a higher model accuracy.

Description

Technical Field

[0001] The invention relates to medium model construction, and in particular to a method for constructing a dual-lithology comprehensive anisotropic equivalent medium model. Background Art

[0002] The differences in seismic acoustic wave velocity and anisotropy are mainly caused by the physical elastic characteristics of rocks. Therefore, the study of seismic acoustic wave velocity and anisotropy characteristics can be analyzed by constructing a rock physics elasticity theoretical model. The current mainstream rock physics elasticity theoretical model uses an isotropic model. However, it is not suitable for unconventional oil and gas reservoirs. The more advanced anisotropic models currently include fracture models or theoretical models (self-consistent models, equivalent differential models, geometric distribution functions, etc.);

[0003] Fracture models (Hudson model, Shoenberg model, etc.) can construct theoretical models for studying fractures in hard rocks such as sandstone, but such models are only suitable for constructing slender fractures formed by compression of hard rocks, and the characterization of anisotropy of soft rocks is difficult and has low accuracy. Theoretical models (self-consistent model, equivalent differential model, geometric distribution function, etc.) can solve the anisotropy of softer rocks such as shale that have a specific direction shape after compression, but such models are only suitable for constructing anisotropy caused by specific shapes formed by compression of soft rocks, and the characterization of anisotropy of fractures in hard rocks is difficult and has low accuracy. Summary of the invention

[0004] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a method for constructing a dual-lithology comprehensive anisotropic equivalent medium model, which can characterize the anisotropy of soft rocks and the anisotropy of cracks in hard rocks with high accuracy.

[0005] The object of the present invention is achieved through the following technical scheme: a method for constructing a dual-lithology comprehensive anisotropic equivalent medium model, comprising the following steps:

[0006] S1. Construct the sandstone skeleton;

[0007] S2. constructing a shale skeleton;

[0008] S3. Combined sandstone and shale framework and fluid replenishment.

[0009] Furthermore, the step S1 comprises:

[0010] S101. First, based on the working area, feldspar / quartz is taken as the main component of the sandstone background matrix, and water and gas are taken as the main components of the fluid; the Wylie time average formula is used to combine the acoustic velocity and density of each component according to its volume content in the rock and fluid;

[0011]

[0012]

[0013]

[0014]

[0015]

[0016] ρf=Swρw+Sgρg

[0017] In the formula, and is the P-wave velocity, S-wave velocity and density of the sandstone matrix; Sg is the gas saturation; Sw is the water saturation; V p represents the longitudinal wave velocity, V s represents the shear wave velocity, V represents the volume fraction, ρ represents the density, and the m in the subscript sand ,qu,fs,f,w,g represent sandstone matrix, quartz, feldspar, fluid, water and gas respectively; for example, Vqu represents the volume fraction of quartz, V pqu Quartz longitudinal wave velocity, V squ Quartz shear wave velocity, ρ qu Quartz density.

[0018] S102. The velocities of the fluid and the background matrix are converted into the bulk modulus and shear modulus of the matrix using the velocity elasticity relationship:

[0019]

[0020]

[0021]

[0022]

[0023] In the formula, and are the bulk modulus, shear modulus and Lame constant of the sandstone matrix, respectively;

[0024] S103. Substitute the matrix result of step S102 into the Hudson fracture model to obtain the dry sandstone skeleton:

[0025]

[0026]

[0027]

[0028]

[0029]

[0030]

[0031]

[0032]

[0033]

[0034]

[0035]

[0036]

[0037]

[0038] in:

[0039]

[0040]

[0041]

[0042]

[0043] Finally we get:

[0044]

[0045]

[0046]

[0047]

[0048]

[0049]

[0050] In the formula, is the elastic parameter of the dry skeleton of sandstone, d represents the skeleton; When the superscript a is 0, 1, or 2, it represents the stiffness matrix of matrix, fracture correction term, and inter-fracture correction term, respectively, corresponding to the stiffness coefficient of the subscript i-th row and j-th column; Q is the characteristic parameter of the fracture, α is the aspect ratio, φ is the porosity, and ε is the fracture density; λ is the lame constant, a commonly used lithological physical and mechanical parameter; U3 and U1 are fracture correction parameters, which characterize the correction of the influence of fractures and fluids in fractures on the matrix elastic parameters.

[0051] S104. Pass Forming a matrix

[0052] Furthermore, the step S2 comprises:

[0053] S201. Shale and shale pores are combined into a shale plate with pores using SCA (anisotropic self-consistent model) and DEM (anisotropic equivalent differential model); first, a porous shale plate with a porosity of 50% is constructed using the anisotropic self-consistent equivalent medium self-consistent model:

[0054]

[0055] Where n and p are equal to 1 when they are small plates in shale, and equal to 2 when they are pores. When they are equivalent to skeletons, the elastic parameter of pores is 0, v1 = v2 = 0.5; is the geometric parameter, I is a scalar tensor, C SCA is the elastic parameter result of the equivalent calculation of the model; C 1 That is, the known shale plate stiffness matrix (can be directly obtained from reference books); C 2 That is, the stiffness matrix of the pores and their fluids, and the stiffness matrix of the fluid of the skeleton is 0;

[0056] The final C' is calculated by step S201 SCA The process is as follows:

[0057] A1. Initialize C SCA =0.5(C 1 +C 2 );

[0058] A2, C SCA Substitute the right side of the formula in step S201 to calculate C' SCA ; Then compare C' SCA With C SCA Is the error between them less than the set threshold? If so, set C' SCA Output as the final calculation result. If not, go to A3;

[0059] A3. Update C SCA , the updated C SCA Equal to C' calculated in A2SCA , then return to step A2;

[0060] S202. When the porosity content is less than 50%, the insufficient shale platelet component is supplemented using the anisotropic equivalent differential equivalent medium model to obtain:

[0061]

[0062] In the formula, v represents the content of shale mineral matrix, C 1 Represents the stiffness matrix of the shale mineral matrix; the parent rock is a uniform two-phase medium. When v = 50%, C DEM (v) = C' SCA ;

[0063] In this step, according to C' obtained in step S201, SCA Iterative calculation C DEM (1-proe), where C DEM (1-proe) The final iteration output is C a , which means the stiffness matrix of the shale platelet without crystal direction anisotropy; pore porosity; the calculation process is as follows:

[0064] The first step is to initialize C DEM (v) = C' SCA , and set the parameter dv = {(1-proe)-v} / n, where n is the set cycle parameter;

[0065] Step 2: C DEM (v) Substitute the right side of the formula in step S201 into the equation to find d(C DEM (v));

[0066] Step 3: Find C DEM (v+dv)=C DEM (v)+d(C DEM (v)); C DEM (v+dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+dv));

[0067] Step 4: Calculate C DEM (v+2dv)=C DEM (v+dv)+d(C DEM (v+dv)), C DEM (v+2dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+2dv));

[0068] Step 5: During each cycle, use the C obtained in the previous cycle DEM (v+(i-1)dv) and d(C DEM (v+(i-1)dv)) calculates C DEM (v+idv), and find d(C DEM (v+idv)), where i=1,2,…,n;

[0069] When i = n, we get C DEM (v+ndv), since v+ndv=1-proe, we get C DEM (1-proe), which is used as the final output result C a ;

[0070] S203. Use geometric functions to simulate the anisotropic characteristics of shale platelets during compaction:

[0071]

[0072]

[0073] Where φ, φ0, m 400 , n 200 are porosity, critical porosity, and model parameters respectively; for the compaction model, model parameters are: m 400 =2,n 200 =1;W 200 , W 400 W is the Legendre coefficient that describes the direction of an object in three-dimensional space in material mechanics; 200 (φ), W 400 (φ) is the porosity W 200 , W 400 ; is the W corresponding to the maximum porosity 200 , W 400 .

[0074] S204. The calculated W 200 (φ), W 400 (φ) and C a Substitute the following formula to calculate the elastic parameters of the shale skeleton:

[0075]

[0076]

[0077]

[0078]

[0079]

[0080]

[0081] In the formula

[0082]

[0083]

[0084]

[0085]

[0086] is the calculated elastic parameter of the shale platelet, a represents the matrix; a1, a2, a3 are isotropic parameters calculated by C a The calculated subscripts i and j represent the positions of the corresponding rows and columns of the stiffness matrix; represents the shale skeleton elastic parameters with directional anisotropy, i, j represent the positions of the corresponding stiffness matrix rows and columns, L, M are Lame constants;

[0087] S205.According to Constructing the Matrix

[0088] Furthermore, the step S3 comprises:

[0089] S301. Shale skeleton and sandstone skeleton are combined into a rock skeleton model using the SCA-DEM anisotropic equivalent medium model. First, a rock skeleton model with sandstone and shale accounting for 50% each is constructed using the anisotropic self-consistent equivalent medium model:

[0090]

[0091] In the formula, when n and p are equal to 1, it means sandstone, and when n and p are equal to 2, it means shale, that is, C 1 , v1 is the stiffness matrix and volume content of the sandstone skeleton; C 2 , v2 is the stiffness matrix and volume content of shale, v1 = v2 = 0.5, indicating that the volume content of each component is 50%;

[0092] Calculate the final C" through step S301 SCA The process is as follows:

[0093] B1. Calculation Initialize C SCA =0.5(C 1 +C 2 );

[0094] B2, C SCA Substitute the right side of the formula in step S201 to calculate C" SCA ; then compare C" SCA With C SCA Is the error between them less than the set threshold? If so, set C” SCA Output as the final calculation result. If not, go to B3;

[0095] B3. Update C SCA , the updated C SCA Equal to C calculated in A2" SCA , then return to step B2;

[0096] S302. The components with a volume content greater than half are supplemented using the anisotropic equivalent differential model to obtain:

[0097]

[0098] The components include sandstone and shale, wherein v is the volume fraction of the component with a volume content greater than 50% (if there is more sandstone, v is the volume fraction of sandstone, if there is more shale, v1 is the volume fraction of shale); C 1 is the stiffness matrix of the component with a volume content greater than 50%; the parent rock is the obtained equal two-phase medium. When the volume fraction of sandstone and shale obtained by logging is 50%, C DEM =C” SCA ;

[0099] In this step, according to C" obtained in step S301 SCA Iterative calculation C DEM (1-proe), where C DEM (1-proe) The final iteration output is C d , expressed as the stiffness matrix of the rock skeleton; pore porosity; the calculation process is as follows:

[0100] The first step is to initialize C DEM (v) = C” SCA , and set the parameter dv = {(1-proe)-v} / n, where n is the set cycle parameter;

[0101] Step 2: C DEM (v) Substitute the right side of the formula in step S201 into the equation to find d(C DEM (v));

[0102] Step 3: Find C DEM (v+dv)=C DEM (v)+d(C DEM (v)); C DEM (v+dv) as the new CDEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+dv));

[0103] Step 4: Calculate C DEM (v+2dv)=C DEM (v+dv)+d(C DEM (v+dv)), C DEM (v+2dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+2dv));

[0104] Step 5: During each cycle, use the C obtained in the previous cycle DEM (v+(i-1)dv) and d(C DEM (v+(i-1)dv)) calculates C DEM (v+idv), and find d(C DEM (v+idv)), where i=1,2,…,n;

[0105] When i = n, we get C DEM (v+ndv), since v+ndv=1-proe, we get C DEM (1-proe), which is used as the final output result C d ;

[0106] S303. Use Brown & Korringa anisotropic fluid substitution formula to import the elastic parameters of the fluid into the skeleton:

[0107]

[0108]

[0109] in:

[0110] S ijkl =I / C ijkl

[0111]

[0112]

[0113] Where S represents the flexibility matrix, which is equal to the unit matrix / stiffness matrix, the physical property is the stress-strain elastic matrix coefficient matrix, ijkl is the coordinate corresponding to the matrix, the superscripts d, sat, and 0 represent the dry rock skeleton, fluid-saturated rock, and matrix minerals, respectively, and eff represents the equivalent result; β flrepresents the compressibility of pore fluid; β0 represents the compressibility of minerals equal to φ represents porosity;

[0114] The rock stiffness constant obtained can be used to obtain rock acoustic and anisotropic parameters through the velocity elastic relationship and Thomson anisotropy model.

[0115]

[0116]

[0117]

[0118]

[0119]

[0120] In the formula, ρ eff , V s eff Represent the equivalent density and P- and S-wave velocities of saturated rock, respectively; ε, γ, δ are Thomson anisotropy characteristic parameters; is the elastic parameter of row i and column j corresponding to the rock saturated stiffness matrix obtained after the final fluid replacement.

[0121] The beneficial effect of the present invention is as follows: considering that the previous anisotropic models are extremely purposeful in analyzing the elasticity and anisotropic characteristics of rocks based on the microstructure of cracks or mudstone, the actual underground medium is often mixed with muddy bands, muddy fillings, etc., with both shale and sandstone. The accuracy of the rock physics model is mainly related to equivalent equivalent conditions. A single model is limited in studying the elasticity and anisotropic characteristics of the actual underground. Therefore, the present invention constructs a comprehensive anisotropic model that can characterize the anisotropy of soft rocks and the anisotropy of cracks in hard rocks in response to the above defects, and has a higher model accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0122] Figure 1 The present invention is a flow chart of the method. DETAILED DESCRIPTION

[0123] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings, but the protection scope of the present invention is not limited to the following.

[0124] like Figure 1 As shown, a method for constructing a dual-lithology comprehensive anisotropic equivalent medium model includes the following steps:

[0125] S1. Construct the sandstone skeleton;

[0126] S2. constructing a shale skeleton;

[0127] S3. Combined sandstone and shale framework and fluid replenishment.

[0128] Furthermore, the step S1 comprises:

[0129] S101. First, based on the working area, feldspar / quartz is taken as the main component of the sandstone background matrix, and water and gas are taken as the main components of the fluid; the Wylie time average formula is used to combine the acoustic velocity and density of each component according to its volume content in the rock and fluid;

[0130]

[0131]

[0132]

[0133]

[0134]

[0135] ρ f =S w ρ w +S g ρ g

[0136] In the formula, and is the longitudinal wave velocity of the sandstone matrix, S g is the gas saturation; S w is water saturation; shear wave velocity and density; V p represents the longitudinal wave velocity, V s represents the shear wave velocity, V represents the volume fraction, ρ represents the density, and the m in the subscript sand ,qu,fs,f,w,g represent sandstone matrix, quartz, feldspar, fluid, water and gas respectively; for example, V qu Represents the volume fraction of quartz, V pqu Quartz longitudinal wave velocity, V squ Quartz shear wave velocity, ρ qu Quartz density.

[0137] S102. The velocities of the fluid and the background matrix are converted into the bulk modulus and shear modulus of the matrix using the velocity elasticity relationship:

[0138]

[0139]

[0140]

[0141]

[0142] In the formula, and are the bulk modulus, shear modulus and Lame constant of the sandstone matrix, respectively;

[0143] S103. Substitute the matrix result of step S102 into the Hudson fracture model to obtain the dry sandstone skeleton:

[0144]

[0145]

[0146]

[0147]

[0148]

[0149]

[0150]

[0151]

[0152]

[0153]

[0154]

[0155]

[0156]

[0157] in:

[0158]

[0159]

[0160]

[0161]

[0162] Finally we get:

[0163]

[0164]

[0165]

[0166]

[0167]

[0168]

[0169] In the formula, are the elastic parameters of the dry skeleton of sandstone; When the superscript a is 0, 1, or 2, it represents the stiffness matrix of matrix, fracture correction term, and inter-fracture correction term, respectively, corresponding to the stiffness coefficient of the subscript i-th row and j-th column; Q is the characteristic parameter of the fracture, α is the aspect ratio, φ is the porosity, and ε is the fracture density; λ is the lame constant, a commonly used lithological physical and mechanical parameter; U3 and U1 are fracture correction parameters, which characterize the correction of the influence of fractures and fluids in fractures on the matrix elastic parameters.

[0170] In the embodiments of the present application, the Hudson model can be replaced by the Schoenberg model during the sandstone construction process.

[0171] S104. Pass Forming a matrix

[0172] Furthermore, the step S2 comprises:

[0173] S201. Shale and shale pores are combined into a shale plate with pores using SCA (anisotropic self-consistent model) and DEM (anisotropic equivalent differential model); first, a porous shale plate with a porosity of 50% is constructed using the anisotropic self-consistent equivalent medium self-consistent model:

[0174]

[0175] Where n and p are equal to 1 when they are small plates in shale, and equal to 2 when they are pores. When they are equivalent to skeletons, the elastic parameter of pores is 0, v1 = v2 = 0.5;

[0176] is the geometric parameter, I is a scalar tensor, C SCA is the elastic parameter result of the equivalent calculation of the model; C 1 That is, the known shale plate stiffness matrix (can be directly obtained from reference books); C 2 That is, the stiffness matrix of the pores and their fluids, and the stiffness matrix of the skeleton fluid is 0;

[0177] C 1 Refers to C n n is equal to 1 or Cp The corresponding parameter when p is equal to 1 is the stiffness matrix of the shale plate;

[0178] C 2 Refers to C n n is equal to 2 or C p The corresponding parameter when p is equal to 2, that is, the stiffness matrix of the pore and its fluid;

[0179] The final C' is calculated by step S201 SCA The process is as follows:

[0180] A1. Initialize C SCA =0.5(C 1 +C 2 );

[0181] A2, C SCA Substitute the right side of the formula in step S201 to calculate C' SCA ; Then compare C' SCA With C SCA Is the error between them less than the set threshold? If so, set C' SCA Output as the final calculation result. If not, go to A3;

[0182] A3. Update C SCA , the updated C SCA Equal to C' calculated in A2 SCA , then return to step A2;

[0183] S202. When the porosity content is less than 50%, the insufficient shale platelet component is supplemented using the anisotropic equivalent differential equivalent medium model to obtain:

[0184]

[0185] In the formula, v represents the content of shale mineral matrix, C 1 Represents the stiffness matrix of the shale mineral matrix; the parent rock is a uniform two-phase medium. When v = 50%, C DEM (v) = C' SCA ;

[0186] In this step, according to C' obtained in step S201, SCA Iterative calculation C DEM (1-proe), where C DEM (1-proe) The final iteration output is C a , which means the stiffness matrix of the shale platelet without crystal direction anisotropy; pore porosity; the calculation process is as follows:

[0187] The first step is to initialize C DEM(v) = C' SCA , and set the parameter dv = {(1-proe)-v} / n, where n is the set cycle parameter;

[0188] Step 2: C DEM (v) Substitute the right side of the formula in step S201 into the equation to find d(C DEM (v));

[0189] Step 3: Find C DEM (v+dv)=C DEM (v)+d(C DEM (v)); C DEM (v+dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+dv));

[0190] Step 4: Calculate C DeM (v+2dv)=C DEM (v+dv)+d(C DEM (v+dv)), C DEM (v+2dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+2dv));

[0191] Step 5: During each cycle, use the C obtained in the previous cycle DEM (v+(i-1)dv) and d(C DEM (v+(i-1)dv)) calculates C DEM (v+idv), and find d(C DEM (v+idv)), where i=1,2,…,n;

[0192] When i = n, we get C DEM (v+ndv), since v+ndv=1-proe, we get C DEM (1-proe), which is used as the final output result C a ;

[0193] S203. Use geometric functions to simulate the anisotropic characteristics of shale platelets during compaction:

[0194]

[0195]

[0196] Where φ, φ0, m 400 , n 200are porosity, critical porosity, and model parameters respectively; for the compaction model, model parameters are: m 400 =2,n 200 =1;W 200 , W 400 W is the Legendre coefficient that describes the direction of an object in three-dimensional space in material mechanics; 200 (φ), W 400 (φ) is the porosity W 200 , W 400 ; is the W corresponding to the maximum porosity 200 , W 400 .

[0197] S204. The calculated W 200 (φ), W 400 (φ) and C a Substitute the following formula to calculate the elastic parameters of the shale skeleton:

[0198]

[0199]

[0200]

[0201]

[0202]

[0203]

[0204] In the formula

[0205]

[0206]

[0207]

[0208]

[0209] are the calculated elastic parameters of the shale plate; a1, a2, a3 are isotropic parameters calculated by C a Obtained by calculation, subscripts i and j indicate the position of the corresponding stiffness matrix rows and columns; represents the shale skeleton elastic parameters with directional anisotropy, i, j represent the positions of the corresponding stiffness matrix rows and columns, L, M are Lame constants;

[0210] S205.According to Constructing the Matrix

[0211] In the embodiment of the present application, for a general anisotropic linear elastic solid, the generalized Hooke's law defines a linear relationship between its stress and strain, which is expressed as:

[0212] σ ij =c ijkl ε kl

[0213] Where σ ij ——stress tensor;

[0214] ε kl ——strain tensor;

[0215] c ijkl ——represents the elastic stiffness tensor, which is a fourth-order tensor. Since the values ​​of i, j, k, and l range from 1 to 3; each stress and strain tensor has 9 components, there are a total of 81 components, namely elastic parameters. These 81 parameters are not independent of each other, because the symmetry of stress and strain c ijkl =c jikl =c ijklk =c jilk , the number of independent elastic constants is thus reduced to 36, while the uniqueness of the strain energy potential is: c ijkl =c klij , which further reduces the number of independent elastic constants to 21, which is also the maximum number of independent elastic constants that any material can achieve.

[0216] Simplified stress, strain, stiffness and toughness tensors

[0217] σ I =c IJ e J

[0218] The elastic stiffness is expressed as a 6*6 matrix to describe anisotropy. Here, the stress vector σ and the strain vector ε are expressed as:

[0219]

[0220] There are several types of classification of medium types by stiffness matrix:

[0221] (1) Triclinic system: the most complex medium containing 21 independent elastic constants. It has no symmetry axis or symmetry plane. Its representative mineral is plagioclase.

[0222] (2) Monoclinic system: contains 13 independent elastic constants, has only one symmetry plane, and its representative mineral is hornblende;

[0223] (3) Trigonal system: It is divided into two cases containing 7 and 6 independent elastic constants, and the representative minerals are ilmenite and tourmaline respectively;

[0224] (4) Tetragonal system: It is also divided into two cases with 7 and 6 independent elastic constants, and the representative minerals are scheelite and cassiterite respectively;

[0225] The above media types are relatively rarely used. The following is a brief introduction to the commonly used stiffness matrices of several simple media:

[0226] (5) Cubic system: contains 3 independent elastic constants, has 4 cubic symmetry axes and 9 symmetry planes. Its representative mineral is garnet, and its matrix can be characterized as:

[0227]

[0228] (6) Isotropy: In an isotropic linear elastic medium, due to its maximum symmetry, there are two independent elastic constants, characterized by:

[0229]

[0230] where c 12 =c 11 -2c 44 , c 11 =λ+2μ,c 44 =2μ, λ and μ are the Lame constants of isotropic linear elastic body;

[0231] (7) Hexagonal system (hexagonal): also known as transversely isotropic medium, contains 5 independent elastic constants, the representative mineral is β-quartz, characterized by:

[0232]

[0233] where c 66 =0.5(c 11 -c 12 );

[0234] (8) Orthorhombic system: also known as orthorhombic system, contains 9 independent elastic constants, has three mutually orthogonal symmetry planes, and its representative mineral is olivine, characterized by:

[0235]

[0236] At present, geophysics mainly treats the layered structure of underground media compacted by natural deposition as equivalent to a hexagonal crystal system (transversely isotropic medium), which is also the elastic characteristic description parameter (C stiffness matrix) used in this application. It is composed of c 11 , c 12 , c13 ,c 33 ,c 44 、c 66 Characterized by.

[0237] In the embodiment of the present application, the step S104 As c 11 , c 12 , c 13 ,c 33 ,c 44 、c 66 , construct the hexagonal matrix, and you can get the matrix In step S205 As c 11 , c 12 , c 13 ,c 33 ,c 44 、c 66 , construct the hexagonal matrix, and you can get the matrix

[0238] Furthermore, the step S3 comprises:

[0239] S301. Shale skeleton and sandstone skeleton are combined into a rock skeleton model using the SCA-DEM anisotropic equivalent medium model. First, a rock skeleton model with sandstone and shale accounting for 50% each is constructed using the anisotropic self-consistent equivalent medium model:

[0240]

[0241] In the formula, when n and p are equal to 1, it means sandstone, and when n and p are equal to 2, it means shale, that is, C 1 , v1 is the stiffness matrix and volume content of the sandstone skeleton; C 2 , v2 is the stiffness matrix and volume content of shale, v1 = v2 = 0.5, indicating that the volume content of each component is 50%;

[0242] Calculate the final C" through step S301 SCA The process is as follows:

[0243] B1. Calculation Initialize C SCA =0.5(C 1 +C 2 );

[0244] B2, C SCA Substitute the right side of the formula in step S201 to calculate C" SCA ; then compare C" SCA With C SCAIs the error between them less than the set threshold? If so, set C” SCA Output as the final calculation result. If not, go to B3;

[0245] B3. Update C SCA , the updated C SCA Equal to C calculated in A2" SCA , then return to step B2;

[0246] S302. The components with a volume content greater than half are supplemented using the anisotropic equivalent differential model to obtain:

[0247]

[0248] The components include sandstone and shale, wherein v is the volume fraction of the component with a volume content greater than 50% (if there is more sandstone, v is the volume fraction of sandstone, if there is more shale, v1 is the volume fraction of shale); C 1 is the stiffness matrix of the component with a volume content greater than 50%; the parent rock is the obtained equal two-phase medium. When the volume fraction of sandstone and shale obtained by logging is 50%, C DEM =C” SCA ;

[0249] In this step, according to C" obtained in step S301 SCA Iterative calculation C DEM (1-proe), where C DEM (1-proe) The final iteration output is C d , expressed as the stiffness matrix of the rock skeleton; pore porosity; the calculation process is as follows:

[0250] The first step is to initialize C DEM (v) = C” SCA , and set the parameter dv = {(1-proe)-v} / n, where n is the set cycle parameter;

[0251] Step 2: C DEM (v) Substitute the right side of the formula in step S201 into the equation to find d(C DEM (v));

[0252] Step 3: Find C DEM (v+dv)=C DEM (v)+d(C DEM (v)); C DEM (v+dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+dv));

[0253] Step 4: Calculate CDEM (v+2dv)=C DEM (v+dv)+d(C DEM (v+dv)), C DEM (v+2dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+2dv));

[0254] Step 5: During each cycle, use the C obtained in the previous cycle DEM (v+(i-1)dv) and d(C DEM (v+(i-1)dv)) calculates C DEM (v+idv), and find d(C DEM (v+idv)), where i=1,2,…,n;

[0255] When i = n, we get C DEM (v+ndv), since v+ndv=1-proe, we get C DME (1-proe), which is used as the final output result C d ;

[0256] S303. Use Brown & Korringa anisotropic fluid substitution formula to import the elastic parameters of the fluid into the skeleton:

[0257]

[0258]

[0259] in:

[0260] S ijkl =I / C ijkl

[0261]

[0262]

[0263] Where S represents the flexibility matrix, which is equal to the unit matrix / stiffness matrix, the physical property is the stress-strain elastic matrix coefficient matrix, ijkl is the coordinate corresponding to the matrix, the superscripts d, sat, and 0 represent the dry rock skeleton, fluid-saturated rock, and matrix minerals, respectively, and eff represents the equivalent result; β fl represents the compressibility of pore fluid; β0 represents the compressibility of minerals equal to φ represents porosity;

[0264] The rock stiffness constant obtained can be used to obtain rock acoustic and anisotropic parameters through the velocity elastic relationship and Thomson anisotropy model.

[0265]

[0266]

[0267]

[0268]

[0269]

[0270] In the formula, ρ eff , V s eff Represent the equivalent density and P- and S-wave velocities of saturated rock, respectively; ε, γ, δ are Thomson anisotropy characteristic parameters; is the elastic parameter of row i and column j corresponding to the rock saturated stiffness matrix obtained after the final fluid replacement.

[0271] In the embodiment of the present application, the fracture model and the model based on shape research of weak rock are used to analyze pure sandstone and pure shale respectively. The results show that there are large differences in the equivalent results of sandstone and shale using the two models, among which the fracture model has higher accuracy for sandstone and the model based on shape research has higher accuracy for shale. This shows that the two types of rocks have typical different anisotropic characteristics, and non-corresponding model processing will produce certain errors. When the comprehensive model is used to analyze rocks with different lithologies, the results will change with the mud content in the rock, which also shows that the mud content in the rock has a strong influence on anisotropy. Finally, the actual well logging is analyzed, and the results are as analyzed above. When the mud content is high, the comprehensive model and the model based on shape research have higher accuracy, and when the mud content is low, the fracture model and the comprehensive model have higher accuracy. In addition, the comprehensive model has higher accuracy and reliability as a whole. Therefore, the comprehensive model is a model with more theoretical support, more reliable and more widely applicable.

[0272] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for constructing a dual-lithology comprehensive anisotropic equivalent medium model, characterized by: The following steps are involved: S1. Construct the sandstone skeleton; S2. constructing a shale skeleton; S3. Combined sandstone and shale skeleton and fluid replenishment: The step S3 comprises: S301. Shale skeleton and sandstone skeleton are combined into a rock skeleton model using the SCA-DEM anisotropic equivalent medium model. First, a rock skeleton model with sandstone and shale accounting for 50% each is constructed using the anisotropic self-consistent equivalent medium model: In the formula, when n and p are equal to 1, it means sandstone, and when n and p are equal to 2, it means shale, that is, C 1 , v1 is the stiffness matrix and volume content of the sandstone skeleton; C 2 , v2 is the stiffness matrix and volume content of shale, v1=v2=0.5, indicating that the volume content of each component is 50%, C SCA is the elastic parameter result of the equivalent calculation of the model, is the geometric morphological parameter, I is a scalar tensor; Calculate the final C" through step S301 SCA The process is as follows: B1. Calculation Initialize C SCA =0.5(C 1 +C 2 ); in, is the elastic parameter of the dry skeleton of sandstone, It represents the elastic parameters of shale skeleton with directional anisotropy; B2, C SCA Substitute the right side of the formula in step S201 to calculate C" SCA ; then compare C" SCA With C SCA Is the error between them less than the set threshold? If so, set C” SCA Output as the final calculation result. If not, go to B3; B3. Update C SCA , the updated C SCA Equal to C calculated in B2" SCA , then return to step B2; S302. The components with a volume content greater than half are supplemented using the anisotropic equivalent differential model to obtain: The components include sandstone and shale, wherein v is the volume fraction of the component with a volume content greater than 50%; C 1 is the stiffness matrix of the component with a volume content greater than 50%; the parent rock is the obtained equal two-phase medium. When the volume fraction of sandstone and shale obtained by logging is 50%, C DEM =C” SCA ; In this step, according to C" obtained in step S301 SCA Iterative calculation C DEM (1-proe), where C DEM (1-proe) The final iteration output is C d , expressed as the stiffness matrix of the rock skeleton; pore porosity; the calculation process is as follows: The first step is to initialize C DEM (v) = C” SCA , and set the parameter dv = {(1-proe)-v} / n, where n is the set cycle parameter; Step 2: C DEM (v) Substitute the right side of the formula in step S201 into the equation to find d(C DEM (v)); Step 3: Find C DEM (v+dv)=C DEM (v)+d(C DEM (v)); C DEM (v+dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+dv)); Step 4: Calculate C DEM (v+2dv)=C DEM (v+dv)+d(C DEM (v+dv)), C DEM (v+2dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+2dv)); Step 5: During each cycle, use the C obtained in the previous cycle DEM (v+(i-1)dv) and d(C DEM (v+(i-1)dv)) calculates C DEM (v+idv), and find d(C DEM (v+idv)), where i=1, 2, ..., n; When i = n, we get C DEM (v+ndv), since v+ndv=1-proe, we get C DEM (1-proe), and use it as the final output result C d ; S303. Import the elastic parameters of the fluid into the skeleton using the anisotropic fluid replacement formula: in: S ijkl =I / C ijkl Where S represents the flexibility matrix, which is equal to the unit matrix / stiffness matrix, the physical property is the stress-strain elastic matrix coefficient matrix, ijkl is the coordinate corresponding to the matrix, the superscripts d, sat, and 0 represent dry rock skeleton, fluid-saturated rock, and matrix minerals, respectively, and eff represents the equivalent result; β fl represents the compressibility of pore fluid; β0 represents the compressibility of minerals equal to φ represents porosity; The rock stiffness constant obtained finally is used to obtain rock acoustic and anisotropic parameters through velocity elasticity relationship and Thomson anisotropy model. In the formula, ρ eff , Represent the equivalent density and P- and S-wave velocities of saturated rock, respectively; ε, γ, δ are Thomson anisotropy characteristic parameters; is the elastic parameter of row i and column j corresponding to the rock saturated stiffness matrix obtained after the final fluid replacement.

2. The method for constructing a dual-lithology comprehensive anisotropic equivalent medium model according to claim 1, characterized in that: The step S1 comprises: S101. First, based on the working area, feldspar / quartz is taken as the main component of the sandstone background matrix, and water and gas are taken as the main components of the fluid; the Wylie time average formula is used to combine the acoustic velocity and density of each component according to its volume content in the rock and fluid; r f =S w r w +S g r g In the formula, and are the longitudinal wave velocity, shear wave velocity and density of the sandstone matrix; V p represents the longitudinal wave velocity, V s represents the shear wave velocity, V represents the volume fraction, ρ represents the density, and the m in the subscript sand , qu, fs, f, w, g represent sandstone matrix, quartz, feldspar, fluid, water and gas respectively; S g is the gas saturation; S w is the water saturation; S102. The velocities of the fluid and the background matrix are converted into the bulk modulus and shear modulus of the matrix using the velocity elasticity relationship: In the formula, and are the bulk modulus, shear modulus and Lame constant of the sandstone matrix, respectively; S103. Substitute the matrix result of step S102 into the Hudson fracture model to obtain the dry sandstone skeleton: in: Finally we get: In the formula, are the elastic parameters of the dry skeleton of sandstone; When the superscript a is 0, 1, or 2, it represents the matrix, fracture correction term, and inter-fracture correction term stiffness matrix, corresponding to the stiffness coefficient of the subscript i-th row and j-th column; q is the fracture characteristic parameter, α is the aspect ratio, φ is the porosity, and ε is the fracture density; λ is the lame constant, a commonly used rock physical and mechanical parameter; U3 and U1 are fracture correction parameters, which characterize the correction of the influence of fractures and fluids in fractures on matrix elastic parameters; S104. Pass Forming a matrix 3. The method for constructing a dual-lithology comprehensive anisotropic equivalent medium model according to claim 1, characterized in that: The step S2 comprises: S201. Shale and shale pores are combined into a shale plate with pores using an anisotropic self-consistent model and an anisotropic equivalent differential model; first, a porous shale plate with a porosity of 50% is constructed using an anisotropic self-consistent equivalent medium self-consistent model: Where n and p are equal to 1 when they are small plates in shale, and equal to 2 when they are pores. When they are equivalent to skeletons, the elastic parameters of pores are 0, v1 = v2 = 0.5; is the geometric parameter, I is a scalar tensor, C SCA is the elastic parameter result of the equivalent calculation of the model; C 1 That is the well-known shale plate stiffness matrix; C 2 That is, the stiffness matrix of the pores and their fluids, and the stiffness matrix of the fluid of the skeleton is 0; The final C' is calculated by step S201 SCA The process is as follows: A1. Initialize C SCA =0.5(C 1 +C 2 ); A2, C SCA Substitute the right side of the formula in step S201 to calculate C' SCA ; Then compare C' SCA With C SCA Is the error between them less than the set threshold? If so, set C' SCA Output as the final calculation result. If not, go to A3; A3. Update C SCA , the updated C SCA Equal to C' calculated in A2 SCA , then return to step A2; S202. When the porosity content is less than 50%, the insufficient shale platelet component is supplemented using the anisotropic equivalent differential equivalent medium model to obtain: In the formula, v represents the content of shale mineral matrix, C 1 Represents the stiffness matrix of the shale mineral matrix; the parent rock is a uniform two-phase medium. When v = 50%, C DEM (v) = C' SCA ; In this step, according to C' obtained in step S201, SCA Iterative calculation C DEM (1-proe), where C DEM (1-proe) The final iteration output is C a , which means the stiffness matrix of the shale platelet without crystal direction anisotropy; pore porosity; the calculation process is as follows: The first step is to initialize C DEM (v) = C' SCA , and set the parameter dv = {(1-proe)-v} / n, where n is the set cycle parameter; Step 2: C DEM (v) Substitute the right side of the formula in step S201 into the equation to find d(C DEM (v)); Step 3: Find C DEM (v+dv)=C DEM (v)+d(C DEM (v)); C DEM (v+dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+dv)); Step 4: Calculate C DEM (v+2dv)=C DEM (v+dv)+d(C DEM (v+dv)), C DEM (v+2dv) as the new C DEM (v) Substitute the right side of the formula in step S201 to obtain d(C DEM (v+2dv)); Step 5: During each cycle, use the C obtained in the previous cycle DEM (v+(i-1)dv) and d(C DEM (v+(i-1)dv)) calculates C DEM (v+idv), and find d(C DEM (v+idv)), where i=1, 2, ..., n; When i = n, we get C DEM (v+ndv), since v+ndv=1-proe, we get C DEM (1-proe), and use it as the final output result C a ; S203. Use geometric functions to simulate the anisotropic characteristics of shale platelets during compaction: Where φ, φ0, m 400 , n 200 are porosity, critical porosity, and model parameters respectively; for the compaction model, model parameters are: m 400 =2,n 200 =1;W 200 , W 400 W is the Legendre coefficient that describes the direction of an object in three-dimensional space in material mechanics; 200 (φ), W 400 (φ) is the porosity W 200 , W 400 ; is the W corresponding to the maximum porosity 200 , W 400 ; S204. The calculated W 200 (φ), W 400 (φ) and C a Substitute the following formula to calculate the elastic parameters of the shale skeleton: In the formula is the calculated elastic parameter of the shale plate; a1, a2, a3 are isotropic parameters calculated by C a The calculated subscripts i and j represent the positions of the corresponding rows and columns of the stiffness matrix; represents the shale skeleton elastic parameters with directional anisotropy, i, j represent the positions of the corresponding stiffness matrix rows and columns; L, M are Lame constants; S205.According to Constructing the Matrix

Images