Prediction method for deformation risk grade of hydraulic fracturing casing

Abstract

The invention discloses a method for predicting the deformation risk level of a hydraulic fracturing casing pipe. The method is specifically implemented according to the following steps that 1, fracture attributes of a coherent body, a curvature body and an ant body are calculated; 2, lithologic attributes of amplitude, texture and phase are calculated; 3, calculating an attribute correlation coefficient Cij according to the fracture attributes and the lithologic attributes in the steps 1 and 2; step 4, risk grade classification: a first-level casing deformation risk: Cij belongs to [0.75, 1], a second-level casing deformation risk: Cij belongs to [0.5, 0.75], a third-level casing deformation risk: Cij belongs to [0.25, 0.5], and a fourth-level casing deformation risk: Cij belongs to [0, 0.25]; and step 5, determining the deformation risk level of the casing at each depth of the target layer according to the step 3 and the step 4, and taking corresponding measures. According to the method for predicting the deformation risk level of the hydraulic fracturing casing pipe, the integrity of a shaft is guaranteed by formulating the risk level and an early warning mechanism, and an oil and gas well is smoothly put into production after fracturing is finished.

Description

Technical Field

[0001] The invention belongs to the technical field of oil and gas field development and relates to a method for predicting the risk level of hydraulic fracturing casing deformation. Background Art

[0002] With the continuous development of seismic technology and the demand for oilfield development and production, seismic technology is increasingly being used in the development and production of oil reservoirs. In recent years, microseismic monitoring has played an important role in reservoir transformation. Unconventional oil and gas development practices at home and abroad have shown that hydraulic fracturing plays an important role in the migration and accumulation of oil and gas, and also has an important influence on the physical properties of oil, gas, water and rocks. During the hydraulic fracturing process, geological parameters, engineering parameters and seismic parameters have a great influence on the deformation of the casing. Therefore, studying the deformation law of the casing and the risk level prediction method before hydraulic fracturing can solve many important problems in the hydraulic fracturing process.

[0003] At present, the biggest problem facing the effective development of low-yield and low-efficiency oil and gas resources is what development technology to adopt to maximize the reserve utilization rate. Hydraulic fracturing is the most effective technical means to achieve the effective development of unconventional oil and gas reservoirs. Large-scale fracturing has become normalized. In the process of hydraulic fracturing, how to ensure the integrity of the wellbore, avoid casing deformation during large-scale fracturing or reduce the probability of casing deformation is a problem that troubles fracturing engineers. Therefore, predicting the position of casing deformation before fracturing is an effective technical means to solve the risk of casing deformation. The patent "A method and device for determining the deformation area of ​​fracturing casing" (publication number: CN106199712B, publication date: 2018.6.1) discloses a method for determining the deformation area of ​​fracturing casing, including: 1) establishing a demonstration physical model that is compatible with the geological factors of a preset study area; 2) determining the maximum curvature attribute in the seismic data of the preset area; 3) the area corresponding to the change rate of the three-dimensional rock mechanics parameters in the preset study area is greater than the preset change rate threshold and the maximum curvature attribute is greater than the preset curvature attribute is determined as the deformation area of ​​the fracturing casing. This method uses the curvature attribute in seismic data without adding sensitive attributes such as coherence and ant volume, and has a single targeted application. Summary of the invention

[0004] The purpose of the present invention is to provide a method for predicting the risk level of hydraulic fracturing casing deformation, which has the characteristic of ensuring the integrity of the wellbore.

[0005] The technical solution adopted by the present invention is a method for predicting the risk level of hydraulic fracturing casing deformation, which is specifically implemented according to the following steps:

[0006] Step 1, calculate the fracture properties of the coherence body, curvature body and ant body;

[0007] Step 2, calculating the lithological attributes of amplitude, texture and phase;

[0008] Step 3: Calculate the attribute correlation coefficient c based on the fracture attributes and lithology attributes in steps 1 and 2 ij ;

[0009] Step 4: Risk level classification: Level 1 risk: C ij ∈[0.75, 1], secondary set change risk: C ij ∈[0.5, 0.75], risk of three-level set change: C ij ∈[0.25, 0.5], the risk of the fourth-level set change: C ij ∈[0, 0.25];

[0010] Step 5: Determine the casing deformation risk level at each depth of the target layer based on steps 3 and 4 and take corresponding measures.

[0011] The present invention is also characterized in that:

[0012] The calculation of the coherent volume properties in step 1 is specifically as follows:

[0013] Step 1.11. The three-dimensional gradient vector at the point (x, y, t) can be expressed as but: Where u is the three-dimensional seismic data, is the gradient;

[0014] Step 1.12: Smooth the three-dimensional gradient vector in step 1.11 to obtain the structure tensor:

[0015] Among them, S is the structure tensor,

[0016] Introducing the weighting factor into the above formula, we get the structure tensor matrix:

[0017] Among them, S S is the structure tensor matrix, ω is the weighting factor, ω=(S x 2 +S y 2 ) 2 ;

[0018] Step 1.13: Calculate the coherent volume properties according to the structure tensor matrix in step 1.12:

[0019] Where p is the apparent dip angle of the seismic data of the main survey line in the x direction, q is the apparent dip angle of seismic data in the y direction, λ is the matrix eigenvalue of the structure tensor matrix.

[0020] The calculation of the curvature body properties in step 1 is as follows:

[0021] Step 1.21, the two-dimensional surface fitting formula is: T(x, y) = ax 2 +by 2 +cxy+dx+ey+f, where a, b, c, d, e, f are coefficients;

[0022] Step 1.22: Calculate the fitting coefficients a, b, c, d, e, f according to step 1.21:

[0023]

[0024]

[0025]

[0026]

[0027]

[0028]

[0029] Where, t 1 , t 2 , t 3 , t 4 , t 5 , t 6 , t 7 , t 8 , t 9 is t 5 The 3×3 sample value at the center;

[0030] Step 1.23, calculate the curvature body properties according to step 1.22:

[0031] Maximum curvature: K max =K m +(K m 2 -K g ) 1 / 2 ,

[0032] Minimum curvature: K min =K m -(K m 2 -K g ) 1 / 2 ,

[0033] Where: K g is the Gaussian curvature,

[0034] K m is the mean curvature,

[0035] but:

[0036] Maximum positive curvature: K maxpc =(a+b)+[(ab) 2 +c 2 ] 1 / 2 ,

[0037] Maximum negative curvature: K minpc =(a+b)-[(ab) 2 +c 2 ] 1 / 2 .

[0038] The calculation of ant body attributes in step 1 is as follows:

[0039] The ant algorithm is used to calculate the transfer probability of an ant seeking the next node:

[0040]

[0041] When ants pass by, the pheromones on the path will increase and will also evaporate over time, so:

[0042]

[0043] Where: i is the current location of the ant; j is the location of the ant's next path; τ ij (t) is the pheromone remaining on the line connecting nodes i and j at time t; η ij is the expected degree from position i to position j; α and β represent the relative importance of the visibility of the pheromone trace and route respectively; τ ij (t+n) is the pheromone on path ij after n migrations; ρ is the residual degree of pheromone; 1-ρ is the volatility of pheromone; Δτ ij is the pheromone remaining on path ij in one cycle.

[0044] The calculation of the amplitude attribute in step 2 is as follows:

[0045] Where: N is the number of sample points in the time window, x i is the amplitude value of the i-th sample point in the time window.

[0046] The calculation of texture attributes in step 2 is as follows:

[0047] The gray-level co-occurrence matrix is: [p(i, j, α, β)] N×N , where α and β represent the directions of the grayscale pairs in three-dimensional space, N represents the grayscale level, i represents the i-th row of the grayscale co-occurrence matrix, and j represents the j-th column of the grayscale co-occurrence matrix.

[0048] According to the gray-level co-occurrence matrix, four texture attributes, namely entropy S, energy E, contrast C and homogeneity H, which are sensitive to cracks, are calculated:

[0049] S=-∑ i Σ j p ij logp ij , E=Σ i Σ j p ij 2 ,

[0050] C=Σ i ∑ j (ij) 2 p ij ,

[0051] Where: S is entropy, E is energy, C is contrast, H is homogeneity, i represents the i-th row of the gray-level co-occurrence matrix, j represents the j-th column of the gray-level co-occurrence matrix, and p is probability.

[0052] The calculation of the phase attribute in step 2 is specifically:

[0053] Where: ω is the circular frequency, the time domain sampling interval is 1, a is a constant, |a|<1.

[0054] Attribute correlation coefficient c ij The calculation is:

[0055] According to the six attributes calculated in steps 1 and 2, each attribute has m variables, X ij Represents the jth index of the i-th attribute after data transformation, and the data matrix is:

[0056]

[0057] but:

[0058] The distance coefficient between the i-th point and the j-th point:

[0059] The i-th attribute [x i1 , x i2 , …, x im ] and the jth attribute [x j1 , x j2 , …, x jm ]:

[0060]

[0061] Where: C ij is the attribute correlation coefficient, C ij ∈[0, 1].

[0062] Step 5 is as follows:

[0063] From steps 3 and 4, we can know that:

[0064] When C ij When ∈[0.75, 1], it belongs to the first-level casing change risk, and no fracturing measures are taken;

[0065] When C ij When ∈[0.5, 0.75], it belongs to the secondary casing change risk, and the fracturing scale is reduced;

[0066] When C ij When ∈[0.25, 0.5], it belongs to the third-level casing change risk, and the fracturing scale is adjusted randomly;

[0067] When C ij When ∈[0, 0.25], the risk of level 4 set change is high and the construction is normal.

[0068] The beneficial effects of the present invention are:

[0069] (1) The method for predicting the deformation risk level of hydraulic fracturing casing of the present invention ensures the integrity of the wellbore, and the oil and gas well is smoothly put into production after the fracturing is completed.

[0070] (2) The method for predicting the risk level of casing deformation during hydraulic fracturing of the present invention uses multi-parameter seismic attributes of fracture attributes and lithology attributes to predict the casing deformation risk section, effectively avoiding the risk of casing deformation during hydraulic fracturing, ensuring the smooth production of oil and gas wells after hydraulic fracturing, and providing technical support for oil and gas field development. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 It is the risk section of casing change predicted by the coherent body attribute in the prediction method of the present invention;

[0072] Figure 2 It is the well section with casing change risk predicted by curvature body attribute in the prediction method of the present invention;

[0073] Figure 3 It is the well section with risk of casing change predicted by the ant body attribute in the prediction method of the present invention;

[0074] Figure 4 It is the well section with casing change risk predicted by amplitude attribute in the prediction method of the present invention;

[0075] Figure 5 It is the well section with risk of casing change predicted by texture attribute in the prediction method of the present invention;

[0076] Figure 6 It is the well section with casing change risk predicted by phase attribute in the prediction method of the present invention;

[0077] Figure 7 It is a prediction diagram of risk level of known well casing change in the prediction method of the present invention. DETAILED DESCRIPTION

[0078] The present invention is described in detail below in conjunction with specific implementation modes.

[0079] Example 1

[0080] The method for predicting the deformation risk level of hydraulic fracturing casing of the present invention is specifically implemented according to the following steps:

[0081] Step 1, calculate the fracture properties of the coherence body, curvature body and ant body;

[0082] Step 2, calculating the lithological attributes of amplitude, texture and phase;

[0083] Step 3: Calculate the attribute correlation coefficient c based on the fracture attributes and lithology attributes in steps 1 and 2 ij ;

[0084] Step 4: Risk level classification: Level 1 risk: C ij ∈[0.75, 1], secondary set change risk: C ij ∈[0.5, 0.75], risk of three-level set change: C ij ∈[0.25, 0.5], the risk of the fourth-level set change: C ij ∈[0, 0.25];

[0085] Step 5: Determine the casing deformation risk level at each depth of the target layer based on steps 3 and 4 and take corresponding measures.

[0086] Example 2

[0087] The method for predicting the deformation risk level of hydraulic fracturing casing of the present invention is specifically implemented according to the following steps:

[0088] Step 1, calculate the fracture properties of the coherence body, curvature body and ant body;

[0089] The calculation of the coherent volume properties in step 1 is specifically as follows:

[0090] Step 1.11. The three-dimensional gradient vector at the point (x, y, t) can be expressed as but: Where u is the three-dimensional seismic data, is the gradient;

[0091] Step 1.12: Smooth the three-dimensional gradient vector in step 1.11 to obtain the structure tensor:

[0092] Among them, S is the structure tensor,

[0093] Introducing the weighting factor into the above formula, we get the structure tensor matrix:

[0094] Among them, S S is the structure tensor matrix, ω is the weighting factor, ω=(S x 2 +S y 2 ) 2 ;

[0095] Step 1.13: Calculate the coherent volume properties according to the structure tensor matrix in step 1.12:

[0096] Where p is the apparent dip angle of the seismic data of the main survey line in the x direction, q is the apparent dip angle of seismic data in the y direction, λ is the matrix eigenvalue of the structure tensor matrix;

[0097] The calculation of the curvature body properties in step 1 is as follows:

[0098] Step 1.21, the two-dimensional surface fitting formula is: T(x, y) = ax 2 +by 2 +cxy+dx+ey+f, where a, b, c, d, e, f are coefficients;

[0099] Step 1.22: Calculate the fitting coefficients a, b, c, d, e, f according to step 1.21:

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106] Where, t 1 , t 2 , t 3 , t 4 , t 5 , t 6 , t 7 , t 8 , t 9 is t 5 The 3×3 sample value at the center;

[0107] Step 1.23, calculate the curvature body properties according to step 1.22:

[0108] Maximum curvature: K max =K m +(K m 2 -K g ) 12 ,

[0109] Minimum curvature: K min =K m -(K m 2 -K g ) 12 ,

[0110] Where: K g is the Gaussian curvature,

[0111] K m is the mean curvature,

[0112] but:

[0113] Maximum positive curvature: K maxpc =(a+b)+[(ab) 2 +c 2 ] 12 ,

[0114] Maximum negative curvature: K minpc =(a+b)-[(ab) 2 +c 2 ] 12 ;

[0115] The calculation of ant body attributes in step 1 is as follows:

[0116] The ant algorithm is used to calculate the transfer probability of an ant seeking the next node:

[0117]

[0118] When ants pass by, the pheromones on the path will increase and will also evaporate over time, so:

[0119]

[0120] Where: i is the current location of the ant; j is the location of the ant's next path; τ ij (t) is the pheromone remaining on the line connecting nodes i and j at time t; η ij is the expected degree from position i to position j; α and β represent the relative importance of the visibility of the pheromone trace and route respectively; τij (t+n) is the pheromone on path ij after n migrations; ρ is the residual degree of pheromone; 1-ρ is the volatility of pheromone; Δτ ij is the pheromone left on path ij in one cycle;

[0121] Step 2, calculating the lithological attributes of amplitude, texture and phase;

[0122] The calculation of the amplitude attribute in step 2 is as follows:

[0123] Where: N is the number of sample points in the time window, x i is the amplitude value of the i-th sample point in the time window;

[0124] The calculation of texture attributes in step 2 is as follows:

[0125] The gray-level co-occurrence matrix is: [p(i, j, α, β)] N×N , where α and β represent the directions of the grayscale pairs in three-dimensional space, N represents the grayscale level, i represents the i-th row of the grayscale co-occurrence matrix, and j represents the j-th column of the grayscale co-occurrence matrix.

[0126] According to the gray-level co-occurrence matrix, four texture attributes, namely entropy S, energy E, contrast C and homogeneity H, which are sensitive to cracks, are calculated:

[0127] S=-∑ i ∑ j p ij logp ij , E=∑ i ∑ j p ij 2 ,

[0128] C=∑ i ∑ j (ij) 2 p ij ,

[0129] Where: S is entropy, E is energy, C is contrast, H is homogeneity, i represents the i-th row of the gray-level co-occurrence matrix, j represents the j-th column of the gray-level co-occurrence matrix, and p is probability;

[0130] The calculation of the phase attribute in step 2 is specifically:

[0131] Where: ω is the circular frequency, the time domain sampling interval is 1, a is a constant, |a|<1;

[0132] Step 3: Calculate the attribute correlation coefficient C based on the fracture attributes and lithology attributes in steps 1 and 2 ij ;

[0133] Step 4: Risk level classification: Level 1 risk: C ij ∈[0.75, 1], secondary set change risk: C ij ∈[0.5, 0.75], risk of three-level set change: C ij ∈[0.25, 0.5], the risk of the fourth-level set change: C ij ∈[0, 0.25];

[0134] Step 5: Determine the casing deformation risk level at each depth of the target layer based on steps 3 and 4 and take corresponding measures.

[0135] Example 3

[0136] The method for predicting the deformation risk level of hydraulic fracturing casing of the present invention is specifically implemented according to the following steps:

[0137] Step 1, calculate the fracture properties of the coherence body, curvature body and ant body;

[0138] The calculation of the coherent volume properties in step 1 is specifically as follows:

[0139] Step 1.11. The three-dimensional gradient vector at the point (x, y, t) can be expressed as but: Where u is the three-dimensional seismic data, is the gradient;

[0140] Step 1.12: Smooth the three-dimensional gradient vector in step 1.11 to obtain the structure tensor:

[0141] Among them, S is the structure tensor,

[0142] Introducing the weighting factor into the above formula, we get the structure tensor matrix:

[0143] Among them, S S is the structure tensor matrix, ω is the weighting factor, ω=(S x 2 +S y 2 ) 2 ;

[0144] Step 1.13: Calculate the coherent volume properties according to the structure tensor matrix in step 1.12:

[0145] Where p is the apparent dip angle of the seismic data of the main survey line in the x direction, q is the apparent dip angle of seismic data in the y direction, λ is the matrix eigenvalue of the structure tensor matrix;

[0146] Calculating 3D gradient vectors using seismic data The gradient vector of each point in the three-dimensional seismic data consists of three elements, representing the x, y and t directions respectively. At the location point (x, y, t), the three-dimensional gradient vector describing the dip and azimuth of the seismic event axis can be expressed as Construct the structure tensor matrix S of seismic data S , because seismic data is affected by noise and other interferences, the gradient calculation needs to be smoothed, and the obtained structure tensor is a 3×3 semi-positive symmetric matrix with an eigenvalue greater than or equal to 0. Therefore, the eigenvalue decomposition method is easy to obtain the eigenvector, and the eigenvector corresponding to the maximum eigenvalue is perpendicular to the reflection interface and represents the normal direction of the seismic event axis. In order to improve the signal-to-noise ratio of seismic data calculation, the weighting factor is introduced to overcome the problem of instability of the structure tensor in the above formula.

[0147] The weighting factor is calculated using complex seismic traces, and the seismic data structure tensor is smoothed; the weighting factor ω is related to the actual data and the direction of the structure tensor vector, and is usually calculated using the complex trace analysis method: ω = (S x 2 +S y 2 ) 2 , using the instantaneous phase Φ as input data, calculate the data structure tensor: Where: u H is the Hilbert transform of the time axis t. The instantaneous wave number sum can be calculated by the following formula:

[0148] The calculation of the curvature body properties in step 1 is as follows:

[0149] Step 1.21, the two-dimensional surface fitting formula is: T(x, y) = ax 2 +by 2 +cxy+dx+ey+f, where a, b, c, d, e, f are coefficients;

[0150] Step 1.22: Calculate the fitting coefficients a, b, c, d, e, f according to step 1.21:

[0151]

[0152]

[0153]

[0154]

[0155]

[0156]

[0157] Where, t 1 , t 2 , t 3 , t 4 , t 5 , t 6 , t 7 , t 8 , t 9 is t 5 The 3×3 sample value at the center;

[0158] Step 1.23, calculate the curvature body properties according to step 1.22:

[0159] Maximum curvature: K max =K m +(K m 2 -K g ) 1 / 2 ,

[0160] Minimum curvature: K min =K m -(K m 2 -K g ) 1 / 2 ,

[0161] Where: K g is the Gaussian curvature,

[0162] K m is the mean curvature,

[0163] but:

[0164] Maximum positive curvature: K maxpc =(a+b)+[(ab) 2 +c 2 ] 1 / 2 ,

[0165] Maximum negative curvature: K minpc =(a+b)-[(ab) 2 +c 2 ] 1 / 2 ;

[0166] The curvature attribute is extracted from the seismic data using the difference method, and the derivatives of the curves and surfaces composed of the interpreted discrete data are calculated. Among the infinite normal curvatures, the one with the largest absolute value is the maximum curvature Kmax, and the one orthogonal to it is the minimum curvature Kmin. These two curvatures represent the two extreme values ​​of curvature. They are the measurements of the maximum and minimum curvatures at each point on the interface. In 2003, Sigismundi and Soldo believed that the relative movement of the fault blocks can be easily seen on the time slice of the maximum curvature data. Faults appear as a critical point of positive and negative curvature. The positive curvature represents the upper wall of the fault, and the negative curvature represents the lower wall of the fault. Among all possible normal curvatures, the most positive curvature is the most positive curvature, and the most negative curvature is the most negative curvature. This curvature highlights the boundary and is mainly used to display faults.

[0167] The calculation of ant body attributes in step 1 is as follows:

[0168] The ant algorithm is used to calculate the transfer probability of an ant seeking the next node:

[0169]

[0170] When ants pass by, the pheromones on the path will increase and will also evaporate over time, so:

[0171]

[0172] Where: i is the current location of the ant; j is the location of the ant's next path; τ ij (t) is the pheromone remaining on the line connecting nodes i and j at time t; η ij is the expected degree from position i to position j; α and β represent the relative importance of the visibility of the pheromone trace and route respectively; τ ij (t+n) is the pheromone on path ij after n migrations; ρ is the residual degree of pheromone; 1-ρ is the volatility of pheromone; Δτ ij is the pheromone left on path ij in one cycle;

[0173] The goal is to achieve global optimization through intelligent information transmission between ant groups. The identification method of ant tracking cracks is based on the basic principles of ant algorithm. Sufficient ant populations are released in the seismic data body. When ants that meet the preset fracture conditions are found, they will "release" pheromones and summon other ants to concentrate on the fracture to track it until the fracture is tracked and identified.

[0174] After each ant migrates n times, it obtains a path, the length of which is L. kIf the path meets the requirement of the shortest path, stop. Otherwise, use the above formula to recalculate the pheromone concentration of each path and conduct a second round of search. Through the continuous updating of pheromones, it will eventually converge to the optimal path.

[0175] Coherence body, curvature body and ant body attributes describe faults or natural fractures of different scales. The coherence body attributes represent the location of large-scale faults or natural fractures, and the change of the coherence body attributes of the target layer represents the development of the faults or natural fractures; the curvature attributes represent the location of medium-scale natural fractures, and the change of the curvature body attributes of the target layer represents the development of the faults or natural fractures; the ant body attributes represent the location of small-scale natural fractures, and the change of the ant body attributes of the target layer represents the development of the faults or natural fractures. Usually, when it is difficult to describe fractures with coherence body attributes and curvature body attributes, ant body attributes are used to finely describe the target layer.

[0176] Step 2, calculating the lithological attributes of amplitude, texture and phase;

[0177] The calculation of the amplitude attribute in step 2 is as follows:

[0178] Where: N is the number of sample points in the time window, x i is the amplitude value of the i-th sample point in the time window;

[0179] The calculation of texture attributes in step 2 is as follows:

[0180] The gray-level co-occurrence matrix is: [p(i, j, α, β)] N×N , where α and β represent the directions of the grayscale pairs in three-dimensional space, N represents the grayscale level, i represents the i-th row of the grayscale co-occurrence matrix, and j represents the j-th column of the grayscale co-occurrence matrix.

[0181] According to the gray-level co-occurrence matrix, four texture attributes, namely entropy S, energy E, contrast C and homogeneity H, which are sensitive to cracks, are calculated:

[0182] S=-∑ i ∑ j p ij logp ij , E=∑ i ∑ j p ij 2 ,

[0183] C=∑ i ∑ j (ij) 2 p ij ,

[0184] Where: S is entropy, E is energy, C is contrast, H is homogeneity, i represents the i-th row of the gray-level co-occurrence matrix, j represents the j-th column of the gray-level co-occurrence matrix, and p is probability;

[0185] A variety of texture attributes are extracted using the gray-level co-occurrence matrix, which is a matrix composed of the combination probability p(i, j, α, β) between the gray levels of the seismic 3D data volume [p(i, j, α, β)] N×N . Where α and β represent the direction of the grayscale pair in three-dimensional space, and N represents the grayscale level. i represents the i-th row of the grayscale co-occurrence matrix, and j represents the j-th column of the grayscale co-occurrence matrix. The probability p is the probability that the point G(x+Δx, y+Δy, z+Δz) with a grayscale value of j appears when the distance is 1 and the directions are α and β, taking the point G(x, y, z) with a grayscale value of i as the starting position. Put the seismic data volume into a three-dimensional coordinate system, then:

[0186] When α=90°, β=0°, the x-axis is the main survey line direction, and p(i, j, α, β) is calculated as: p(i, j, α=90°, β=0°)=∑{[(x, y, z), (x+Δx, y+Δy, z+Δz)]|G(x, y, z)=i, G(x+Δx, y+Δy, z+Δz)=j, Δx=1, Δy=0, Δz=0} (1)

[0187] When α=0°, β=0°, the y-axis is the direction of the contact line, and p(i, j, α, β) is calculated as: p(i, j, α=0 / , β=0°)=∑{[(x, y, z), (x+Δx, y+Δy, z+Δz)]|G(x, y, z)=i, G(x+Δx, y+Δy, z+Δz)=j, Δx=0, Δy=1, Δz=0} (2)

[0188] When β = 0°, the y-axis is the direction of the contact line, and p(i, j, α, β) is calculated as:

[0189] p(i, j, α, β=90°)=∑{[(x, y, z), (x+Δx, y+Δy, z+Δz)]|G(x, y, z)=i, G(x+Δx, y+Δy, z+Δz)=j, Δx=0, Δy=0, Δz=1} (3)

[0190] Formula (1) searches for adjacent grayscales (i, j) along the main survey line, formula (2) searches for adjacent grayscales (i, j) along the connecting survey line, and formula (3) searches for adjacent grayscales (i, j) along the time direction. Using the grayscale co-occurrence matrix, 14 texture attributes can be calculated, and four texture attributes, entropy S, energy E, contrast C, and homogeneity H, which are more sensitive to cracks, are preferably calculated. Entropy reflects the complexity of the texture. The larger the entropy value, the more complex the distribution of the texture in the image; energy reflects the stability of texture changes. The larger the energy, the more stable the texture changes; contrast reflects the clarity of the texture. The larger the contrast, the clearer the texture; homogeneity reflects the uniformity of texture distribution. The larger the homogeneity, the more uniform the texture distribution. According to the representative meaning of these four texture attributes, it can be seen that the larger the entropy value, the smaller the energy, the larger the contrast, the smaller the homogeneity, and the more developed the cracks.

[0191] The calculation of the phase attribute in step 2 is specifically:

[0192] Where: ω is the circular frequency, the time domain sampling interval is 1, a is a constant, |a|<1;

[0193] Z-transformation transforms the time domain signal into an expression in the complex frequency domain. In actual phase calculation, the principal value of the inverse tangent is not directly taken, but the numerator and denominator in the brackets are used as the vertical and horizontal coordinates to determine the quadrant of the phase.

[0194] Amplitude, texture and phase attributes characterize the changes in reservoir lithology. Amplitude attributes are extracted based on seismic data. Amplitude attributes represent changes in oil, gas and lithology. Large differences in amplitude attributes of the target layer indicate large changes in reservoir physical properties. Texture attributes are extracted based on seismic data. Texture attributes mainly describe changes in special lithology bodies. Large differences in texture attributes of the target layer indicate large changes in reservoir lithology. Phase attributes are extracted based on seismic data. Phase attributes mainly reflect faults and large sedimentary features. Large differences in phase attributes indicate large changes in reservoir physical properties.

[0195] Step 3: Calculate the attribute correlation coefficient C based on the fracture attributes and lithology attributes in steps 1 and 2 ij ;

[0196] Attribute correlation coefficient C ij The calculation is:

[0197] According to the six attributes calculated in steps 1 and 2, each attribute has m variables, X ij Represents the jth index of the i-th attribute after data transformation, and the data matrix is:

[0198]

[0199] The relationship between the rows of the matrix is ​​equivalent to the correlation between the attributes. If the six attributes are regarded as six attribute points in the m-dimensional space, then:

[0200] The distance coefficient between the i-th point and the j-th point: d ij The smaller the value, the smaller the attribute x i and x j The closer they are, the more they are grouped together;

[0201] The i-th attribute [x i1 , x i2 , …, x im ] and the jth attribute [x j1 , x j2 , …, x jm ]:

[0202]

[0203] Where: C ij is the attribute correlation coefficient, C ij ∈[0, 1], C ij The closer it is to 1, the higher the attribute correlation obtained;

[0204] Step 4: Risk level classification: Level 1 risk: C ij ∈[0.75, 1], secondary set change risk: C ij ∈[0.5, 0.75], risk of three-level set change: C ij ∈[0.25, 0.5], the risk of the fourth-level set change: C ij ∈[0, 0.25];

[0205] Step 5: Determine the casing deformation risk level at each depth of the target layer based on steps 3 and 4 and take corresponding measures;

[0206] Step 5 is as follows:

[0207] From steps 3 and 4, we can know that:

[0208] When C ij When ∈[0.75, 1], it belongs to the first-level casing change risk, and no fracturing measures are taken;

[0209] When C ij When ∈[0.5, 0.75], it belongs to the secondary casing change risk, and the fracturing scale is reduced;

[0210] When C ij When ∈[0.25, 0.5], it belongs to the third-level casing change risk, and the fracturing scale is adjusted randomly;

[0211] When C ijWhen ∈[0, 0.25], the risk of level 4 set change is high and the construction is normal.

[0212] Coherence body attributes, curvature body attributes, and ant body attributes represent fracture attributes, while amplitude attributes, texture attributes, and phase attributes represent lithology attributes. Cluster analysis is performed on fracture attributes and lithology attributes to select sensitive attributes to casing change during hydraulic fracturing of the target layer. The higher the correlation coefficient, the greater the probability of casing change.

[0213] The method for predicting the deformation risk level of hydraulic fracturing casing of the present invention is as follows: the well coherence body attribute, curvature body attribute and ant body attribute are known to predict the casing deformation risk well sections respectively. Figure 1-Figure 3 As shown in the figure, the amplitude attribute, texture attribute and phase attribute predict the casing change risk section respectively as follows Figure 4-Figure 6 As shown, the horizontal coordinates are all in the east-west direction, and the vertical coordinates are all in the north-south direction. The coherence body attributes predict the casing change risk sections, which are the 1st and 5th sections; the curvature body attributes predict the casing change risk sections, which are the 1st, 5th and 6th sections; the ant body attributes predict the casing change risk sections, which are the 1st, 2nd, 5th and 7th sections; the amplitude attributes predict the casing change risk sections, which are the 2nd and 7th sections; the texture attributes predict the casing change risk sections, which are the 2nd and 3rd sections; the phase attributes predict the casing change risk sections, which are the 1st, 5th and 6th sections.

[0214] The prediction diagram of casing change risk level of known wells by the prediction method of the present invention is as follows: Figure 7 As shown, gray to black corresponds to the risk level from low to high.

[0215] The correlation coefficients of the coherence volume attributes, curvature volume attributes, ant volume attributes, amplitude attributes, texture attributes, and phase attributes of each depth segment of the target layer are calculated respectively; the correlation coefficient of the first segment is C ij =0.92, the correlation coefficient of the second segment is C ij =0.45, the correlation coefficient of the third segment is C ij =0.18, the correlation coefficient of the fourth segment is C ij =0.16, the correlation coefficient of the 5th segment is C ij =0.88, the correlation coefficient of the sixth segment is C ij =0.68, the correlation coefficient of the 7th segment is C ij =0.54, then the casing change risk level prediction table of the known wells is shown in Table 1:

[0216] Table 1

[0217] Predict casing change risk sections Correlation coefficient Risk level of the change Paragraph 1 0.92 Level 1 risk section Paragraph 2 0.45 Level 3 risk section Paragraph 3 0.18 Level 4 risk section Paragraph 4 0.16 Level 4 risk section Paragraph 5 0.88 Level 1 risk section Paragraph 6 0.68 Secondary risk section Paragraph 7 0.54 Level 3 risk section

[0218] The method for predicting the risk level of casing deformation in hydraulic fracturing of the present invention aims to classify the influencing factors that cause casing deformation, formulate risk levels and early warning mechanisms, optimize hydraulic fracturing engineering parameters in well sections with higher risk levels, including: liquid volume, sand volume, displacement, temporary plugging agent addition time, etc., to ensure that the oil and gas wells are smoothly put into production after the hydraulic fracturing is completed, avoid engineering safety accidents, and provide technical support for oil and gas field development; through this application, during the hydraulic fracturing process, in the first-level casing change risk well section, measures of not fracturing are taken, in the second-level casing change risk well section, measures of reducing the fracturing scale are taken, in the third-level casing change risk well section, the fracturing scale is randomly adjusted during the fracturing construction process; in the fourth-level casing change risk well section, normal construction is carried out; this method ensures the integrity of the wellbore, and the oil and gas wells are smoothly put into production after the fracturing is completed.

Claims

1. A method for predicting the risk level of hydraulic fracturing casing deformation, characterized in that: Follow the steps below to implement it: Step 1, calculate the fracture properties of the coherence body, curvature body and ant body; Step 2, calculating the lithological attributes of amplitude, texture and phase; Step 3: Calculate the attribute correlation coefficient C based on the fracture attributes and lithology attributes in steps 1 and 2 ij ; Step 4: Risk level classification: Level 1 risk: C ij ∈[0.75, 1], secondary set change risk: C ij ∈[0.5, 0.75], risk of three-level set change: C ij ∈[0.25, 0.5], the risk of the fourth-level set change: C ij ∈[0, 0.25]; Step 5: Determine the casing deformation risk level at each depth of the target layer based on steps 3 and 4 and take corresponding measures.

2. The method for predicting casing deformation risk level according to claim 1, characterized in that: The calculation of the coherent volume attributes in step 1 is specifically as follows: Step 1.

11. The three-dimensional gradient vector at the point (x, y, t) can be expressed as but: Where u is the three-dimensional seismic data, is the gradient; Step 1.12: Smooth the three-dimensional gradient vector in step 1.11 to obtain the structure tensor: Among them, S is the structure tensor, Introducing the weighting factor into the above formula, we get the structure tensor matrix: Among them, S S is the structure tensor matrix, ω is the weighting factor, ω=(S x 2 +S y 2 ) 2 ; Step 1.13: Calculate the coherent volume properties according to the structure tensor matrix in step 1.12: Where p is the apparent dip angle of the seismic data of the main survey line in the x direction, q is the apparent dip angle of seismic data in the y direction, λ is the matrix eigenvalue of the structure tensor matrix.

3. The method for predicting casing deformation risk level according to claim 1, characterized in that: The calculation of the curvature body attribute in step 1 is specifically as follows: Step 1.21, the two-dimensional surface fitting formula is: T(x, y) = ax 2 +by 2 +cxy+dx+ey+f, where a, b, c, d, e, f are coefficients; Step 1.22: Calculate the fitting coefficients a, b, c, d, e, f according to step 1.21: Where, t1, t2, t3, t4, t5, t6, t7, t8, t9 are 3×3 sample values ​​centered on t5; Step 1.23, calculate the curvature body properties according to step 1.22: Maximum curvature: K max =K m +(K m 2 -K g ) 1 / 2 , Minimum curvature: K min =K m -(K m 2 -K g ) 1 / 2 , Where: K g is the Gaussian curvature, K m is the mean curvature, but: Maximum positive curvature: K maxpc =(a+b)+[(ab) 2 +c 2 ] 1 / 2 , Maximum negative curvature: K minpc =(a+b)-[(ab) 2 +c 2 ] 1 / 2 .

4. The method for predicting casing deformation risk level according to claim 1, characterized in that: The calculation of the ant body attributes in step 1 is specifically as follows: The ant algorithm is used to calculate the transfer probability of an ant seeking the next node: When ants pass by, the pheromones on the path will increase and will also evaporate over time, so: Where: i is the current location of the ant; j is the location of the next path of the ant; τ ij (t) is the pheromone remaining on the line connecting nodes i and j at time t; η ij is the expected degree from position i to position j; α and β represent the relative importance of the visibility of the pheromone trace and route respectively; τ ij (t+n) is the pheromone on path ij after n migrations; ρ is the residual degree of pheromone; 1-ρ is the volatility of pheromone; Δτ ij is the pheromone remaining on path ij in one cycle.

5. The method for predicting casing deformation risk level according to claim 1, characterized in that: The calculation of the amplitude attribute in step 2 is specifically as follows: Where: N is the number of sample points in the time window, x i is the amplitude value of the i-th sample point in the time window.

6. The method for predicting casing deformation risk level according to claim 1, characterized in that: The calculation of the texture attribute in step 2 is specifically as follows: The gray-level co-occurrence matrix is: [p(i, j, α, β)] N×N , where α and β represent the directions of the grayscale pairs in three-dimensional space, N represents the grayscale level, i represents the i-th row of the grayscale co-occurrence matrix, and j represents the j-th column of the grayscale co-occurrence matrix. According to the gray-level co-occurrence matrix, four texture attributes, namely entropy S, energy E, contrast C and homogeneity H, which are sensitive to cracks, are calculated: S=-∑ i ∑ j p ij logp ij ,E=∑ i ∑ j p ij 2 , C=∑ i ∑ j (i-j) 2 p ij , Where: S is entropy, E is energy, C is contrast, H is homogeneity, i represents the i-th row of the gray-level co-occurrence matrix, j represents the j-th column of the gray-level co-occurrence matrix, and p is probability.

7. The method for predicting casing deformation risk level according to claim 1, characterized in that: The calculation of the phase attribute in step 2 is specifically as follows: Where: ω is the circular frequency, the time domain sampling interval is 1, a is a constant, |a|<1.

8. The method for predicting casing deformation risk level according to claim 1, characterized in that: The attribute correlation coefficient C ij The calculation is: According to the six attributes calculated in steps 1 and 2, each attribute has m variables, X ij Represents the jth index of the i-th attribute after data transformation, and the data matrix is: but: The distance coefficient between the i-th point and the j-th point: The i-th attribute [x i1 , x i2 , …, x im ] and the jth attribute [x j1 , x j2 , …, x jm ]: Where: C ij is the attribute correlation coefficient, C ij ∈[0, 1].

9. The method for predicting casing deformation risk level according to claim 1, characterized in that: The step 5 is specifically as follows: From steps 3 and 4, we can know that: When C ij When ∈[0.75, 1], it belongs to the first-level casing change risk, and no fracturing measures are taken; When C ij When ∈[0.5, 0.75], it belongs to the secondary casing change risk, and the fracturing scale is reduced; When C ij When ∈[0.25, 0.5], it belongs to the third-level casing change risk, and the fracturing scale is adjusted randomly; When C ij ∈[0, 0.25], it belongs to the fourth level of set change risk and normal construction.

Images