A method for predicting pre-stack fractures of lacustrine carbonate rocks based on fractional order total variation

Abstract

The present application relates to the field of geophysical inversion and oil and gas reservoir prediction, and particularly relates to a method for predicting fractures of lacustrine carbonate rock based on fractional order total variation. The method comprises the following steps: obtaining an initial model of a to-be-inverted parameter by preprocessing seismic data; constructing a forward model of prestack elastic parameters and an objective function based on fractional order total variation; updating the inversion of the log matrix L by combining the alternating direction multiplier method, the gradient being zero and the extreme value, and eigenvalue decomposition; obtaining Young's modulus, Poisson's ratio and density according to the relationship between Young's modulus, Poisson's ratio and the elastic three-parameter log matrix L i+1 ; and obtaining the fracture prediction result by slicing the inversion result of Young's modulus and Poisson's ratio according to the horizon. The present application uses fractional order total variation to establish a sparse regularization term and construct an objective function of prestack elastic three-parameter seismic inversion, thereby improving the accuracy and stability of the inversion result and the accuracy of fracture prediction.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of geophysical inversion and oil and gas reservoir prediction, and in particular to a lake carbonate pre-stack fracture prediction method based on fractional order total variation. BACKGROUND

[0002] Seismic inversion is an important means of oil and gas reservoir prediction. It is a process of establishing a forward model according to the mathematical relationship between the detected seismic record and the physical quantity to be inverted, and using an optimization method to solve the optimal estimation of the forward model. Seismic inversion technology based on sparse regularization is an important method of seismic inversion. It introduces sparse regularization constraints into the forward model, and successfully improves the resolution and robustness of the inversion result by using sparse information.

[0003] With the increasing difficulty of obtaining conventional oil and gas resources, global oil and gas exploration has entered an era of simultaneous development of conventional and unconventional oil and gas. Rock elastic parameters such as Young's modulus and Poisson's ratio are proportional to the porosity of the lake carbonate reservoir, so pre-stack elastic parameter inversion is an effective way to predict fractures from seismic data. For example, the Bayesian inversion method, Alemie and Sacchi introduced the multivariate Cauchy distribution into Bayesian inversion. This prior distribution can combine the correlation between model parameters; Huber distribution, t distribution, etc. have also been introduced into Bayesian inversion. The above prior distributions all assume that the reflectivity is a random variable, and its distribution satisfies a certain distribution characterized by a certain parameter. When the reflectivity distribution does not satisfy the specific distribution, the estimation result of the reflectivity may have errors, and the accuracy and resolution of the pre-stack seismic inversion of elastic parameters are low.

[0004] CN114428301A discloses a pre-stack elastic parameter inversion method, which comprises obtaining pre-stack seismic data; obtaining the reflectivity of the elastic parameter; constructing a forward matrix of the pre-stack seismic data; obtaining a likelihood function of the seismic data according to the forward matrix; obtaining a target function of the elastic parameter inversion according to the likelihood function and prior information; based on the maximum likelihood estimation method, the likelihood function is converted into a maximum marginal likelihood function; the maximum marginal likelihood function is solved to obtain an estimation factor; and the elastic parameter is inverted according to the estimation factor and the target function. By introducing automatic correlation discriminant prior information, the prior information of the elastic parameter is fused into the inversion.

[0005] CN110857997A discloses a step-by-step pre-stack elastic parameter inversion method based on lateral constraint, which comprises the following steps:

[0006] 1) Based on the PP wave and PS wave angle gathers matched with the same phase axis, the PP wave and PS wave angle domain reflectivity is obtained;

[0007] 2) Using the Stewart reflection coefficient formula for multi-wave joint inversion, and using the Kalman filtering algorithm to constrain the inversion equation horizontally, the longitudinal wave velocity variation rate and the transverse wave velocity variation rate are obtained;

[0008] 3) Based on the longitudinal wave velocity variation rate and the transverse wave velocity variation rate, the Aki-Richards three-parameter reflection coefficient formula is substituted, and the Kalman filtering algorithm is used to constrain the inversion equation horizontally to obtain the density variation rate;

[0009] 4) The inversion results obtained in steps 2) and 3) are compensated for low frequency to obtain the absolute value of the final prestack elastic parameter.

[0010] However, the existing prestack elastic parameter inversion method has low accuracy of Young's modulus and Poisson's ratio, resulting in poor crack prediction effect. Therefore, an inversion method is needed to overcome the above problems. SUMMARY

[0011] The main purpose of the present application is to provide a lake carbonate rock prestack crack prediction method based on fractional order total variation. The method can improve the accuracy and stability of the inversion results of Young's modulus and Poisson's ratio, and further obtain high-precision crack prediction results through high-precision inversion results of Young's modulus and Poisson's ratio; overcome the shortcomings of the prior art.

[0012] To achieve the above purpose, the technical scheme adopted by the present application is as follows:

[0013] The present application provides a lake carbonate rock prestack crack prediction method based on fractional order total variation, which comprises the following steps:

[0014] Step 1: pre-processing the seismic data to obtain the initial model of the parameter to be inverted;

[0015] Step 2: according to the processed seismic data, the Young's modulus, Poisson's ratio and density are combined to form a logarithmic matrix L0, and a prestack elastic parameter forward model and an objective function based on fractional order total variation are constructed;

[0016] Step 3: combining the alternating direction multiplier method, the gradient being zero and the eigenvalue decomposition, the logarithmic matrix is inverted and updated to obtain the updated logarithmic matrix L i ;

[0017] Step 4: judge whether the values before and after updating satisfy ||L i+1 -L i ||2 / ||L i ||2>tol, if yes, return to step 3 for loop; if no, according to the Young's modulus, Poisson's ratio and density and the elastic three-parameter logarithmic matrix L i+1The Young's modulus, Poisson's ratio and density are obtained from the relationship of the seismic data and the fracture prediction results are obtained by slicing the Young's modulus and Poisson's ratio according to the layer.

[0018] Step 5: The Young's modulus and Poisson's ratio inversion results are sliced according to the layer to obtain the fracture prediction results.

[0019] Further, in step 1, the seismic data includes seismic records, wavelet data and logging data; the layer information is extracted from the seismic records and the logging data is filtered by interpolation to obtain the initial model of the Young's modulus, Poisson's ratio and density.

[0020] Further, in step 2, the prestack elastic parameter forward model and the objective function based on the fractional order total variation are constructed, which specifically includes the following steps:

[0021] The longitudinal and transverse fractional order difference matrix and the longitudinal first order difference matrix are constructed by using the longitudinal and transverse fractional order difference vector,

[0022]

[0023]

[0024] D=d*E (3i-1)×(3i-1)

[0025] Wherein, * represents convolution operation, represents the transverse fractional order difference matrix, represents the longitudinal fractional order difference matrix, D represents the longitudinal first order difference matrix, d = [-11] represents the longitudinal first order difference vector, i represents the total number of rows of a single parameter to be inverted, and j represents the total number of columns of a single parameter to be inverted;

[0026] The forward model of the prestack elastic three parameters based on the fractional order total variation is constructed:

[0027]

[0028] Wherein, |||1 represents the L1 norm, |||2 represents the L2 norm, * represents convolution operation; S is a seismic record matrix, W is a wavelet matrix, and G represents a coefficient matrix of the elastic three parameters to be inverted;

[0029] On the basis of the forward model, the Lagrange multiplier term R y , R x and the dual term C y , C x are introduced to obtain the prestack elastic three parameter seismic inversion objective function based on the fractional order total variation:

[0030]

[0031] Wherein, η is a dual weight coefficient.

[0032] Further, the relationship between the elastic three-parameter logarithmic matrix and the Young's modulus, Poisson's ratio and density is expressed by the following formula:

[0033] L = [lnE lnσ lnρ] T

[0034] wherein L represents the elastic three-parameter logarithmic matrix, E represents the Young's modulus, σ represents the Poisson's ratio, ρ represents the density, ln represents the logarithmic operation, and T represents the transpose of the matrix;

[0035] Further, the longitudinal and transverse fractional order difference vectors are expressed by the following formula:

[0036]

[0037] wherein ψ a (k) = (-1) k Γ(a+1) / [Γ(k+1)Γ(a-k+1)];

[0038] Further, S is composed of the seismic records S θ1 , S θ2 and S θ3 from three angles, and is expressed by the following formula:

[0039] S = [S θ1 , S θ2 , S θ3 ] T

[0040] W is expressed by the following formula:

[0041]

[0042] wherein w q represents the qth data of the seismic wavelet, and q represents the length of the wavelet;

[0043] G is expressed by the following formula:

[0044]

[0045] wherein G E represents the coefficient matrix of the Young's modulus, G σ represents the coefficient matrix of the Poisson's ratio, and G ρ represents the coefficient matrix of the density.

[0046] Further, the step 3 includes the following specific steps:

[0047] The initial elastic three-parameter logarithmic matrix L0 is taken as the initial model L 1 of the fractional order full variation pre-stack elastic three-parameter seismic inversion, and the initial value R of the Lagrange multiplier term is introducedx = R y = 0 and its dual C x = C y = 0.

[0048] The elastic three-parameter log matrix L is updated by using the alternating direction multiplier algorithm, the rule of taking extreme value when gradient is zero and eigenvalue decomposition i+1 :

[0049]

[0050] where U and V are square matrices and The orthogonal matrix obtained after eigenvalue decomposition:

[0051]

[0052] where SVD represents eigenvalue decomposition, and Λ and Σ are the eigenvalue diagonal matrices obtained after decomposition;

[0053] The update formula is as follows:

[0054]

[0055] where a represents the number of eigenvalues in the eigenvalue diagonal matrix Λ, and b represents the number of eigenvalues in the eigenvalue diagonal matrix Σ,

[0056] The Lagrange multiplier item R is updated according to the alternating direction multiplier algorithm and the soft threshold shrinkage algorithm x i+1 , R y i+1 ;

[0057] The dual item C is updated according to the alternating direction multiplier algorithm and Fermat's lemma

[0058] Further, the updated Lagrange multiplier item R x i+1 , R y i+1 is represented by the following formula:

[0059]

[0060] Further, the updated dual item C is represented by the following formula:

[0061]

[0062] Further, in step 4, the parameters to be inverted and the elastic three-parameter log matrix Li+1 The relationship is expressed by the following formula:

[0063]

[0064] Compared with the prior art, the present invention has the following advantages:

[0065] This invention employs fractional total variation to establish sparse regularization terms and constructs an objective function for pre-stack elastic three-parameter seismic inversion, thereby improving the accuracy and stability of the inversion results and thus enhancing the accuracy of crack prediction. Attached Figure Description

[0066] Figure 1 This is a flowchart of the method described in an embodiment of the present invention.

[0067] Figure 2 For the noisy seismic profile of the method described in the embodiments of the present invention: (1) represents S θ1 (2) indicates S θ2 (3) indicates S θ3 .

[0068] Figure 3 The following is a schematic diagram of the initial models of Young's modulus, Poisson's ratio and density in the method described in the embodiments of the present invention: (1) represents the initial model of Young's modulus, (2) represents the initial model of Poisson's ratio, and (3) represents the initial model of density.

[0069] Figure 4 The following is a schematic cross-sectional view of the inversion results of Young's modulus, Poisson's ratio and density in the method described in the embodiments of the present invention: (1) represents the inversion result of Young's modulus, (2) represents the inversion result of Poisson's ratio, and (3) represents the inversion result of density.

[0070] Figure 5 The method described in this embodiment of the invention provides the result of predicting cracks along the layer slice using Young's modulus.

[0071] Figure 6 This refers to the Poisson's ratio crack prediction result along the slice in the method described in the embodiments of the present invention. Detailed Implementation

[0072] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0073] It is to be noted that the terms used herein are only intended to describe specific embodiments and are not intended to limit the exemplary embodiments according to the present application. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise, and it is further understood that the terms "comprising" and / or "including" when used in this specification, specify the presence of stated features, steps, operations, and / or combinations thereof.

[0074] In order to enable a person skilled in the art to more clearly understand the technical solutions of the present application, the technical solutions of the present application will be described in detail below in combination with specific embodiments.

[0075] Embodiment 1

[0076] Step 1: Pretreatment of lacustrine carbonate seismic data and obtaining initial model of parameters to be inverted:

[0077] Step 1.1: Collecting seismic records S, wavelet data w and well logging data;

[0078] Step 1.2: Extracting horizon information from seismic records and performing interpolation filtering on well logging data to obtain initial models of Young's modulus, Poisson's ratio and density.

[0079] Step 2: Combining Young's modulus, Poisson's ratio and density into an elastic three-parameter matrix according to the processed seismic data, and calculating the logarithmic matrix L0 of the elastic three-parameter matrix, calculating the fractional order difference matrix to construct a prestack elastic parameter forward model based on fractional order total variation and an objective function, and the step is specifically:

[0080] Step 2.1: Combining an elastic three-parameter matrix based on the initial models of Young's modulus, Poisson's ratio and density, and obtaining an initial elastic three-parameter logarithmic matrix L0 according to the recursive relationship between the elastic three-parameter matrix and its logarithm; the relationship between the elastic three-parameter logarithmic matrix and Young's modulus, Poisson's ratio and density is shown in formula 1:

[0081] L = [lnE lnσ lnρ] T Formula (1)

[0082] Wherein, L represents the elastic three-parameter logarithmic matrix, E represents Young's modulus, σ represents Poisson's ratio, ρ represents density, ln represents logarithmic operation, and T represents the transpose of the matrix;

[0083] Step 2.2: Calculating the vertical and horizontal fractional order difference vectors:

[0084] The calculation formula is shown in formula 2 and 3:

[0085] ψ a (k) = (-1) kΓ(a+1) / [Γ(k+1)Γ(a-k+1)] Formula (2)

[0086]

[0087] k, a are fractional order difference coefficients.

[0088] Step 2.3: Constructing longitudinal and transverse fractional order difference matrix and longitudinal first order difference matrix by using longitudinal and transverse fractional order difference vector, as shown in formulas 4, 5 and 6

[0089]

[0090]

[0091] D=d*E (3i-1)×(3i-1) Formula (6)

[0092] wherein * represents convolution operation, represents transverse fractional order difference matrix, represents longitudinal fractional order difference matrix, D represents longitudinal first order difference matrix, d=[-1 1] represents longitudinal first order difference vector, E represents unit matrix, i represents total row number of a single parameter to be inverted, and j represents total column number of a single parameter to be inverted.

[0093] Step 2.4: Constructing forward elastic three-parameter model based on fractional order total variation, as shown in formula 7:

[0094]

[0095] wherein |||1 represents L1 norm, |||2 represents L2 norm, * represents convolution operation, λ is difference regularization factor, and μ is fidelity term weight coefficient. S is composed of near, middle and far three angle seismic records S θ1 , S θ2 and S θ3 , as shown in formula 8

[0096] S=[S θ1 , S θ2 , S θ3 ] T Formula (8)

[0097] W is wavelet matrix, as shown in formula 9

[0098]

[0099] wherein w q represents the qth data of seismic wavelet, and q represents length of wavelet.

[0100] G represents coefficient matrix of elastic three parameters to be inverted, as shown in formula 10

[0101]

[0102] where G E is the coefficient matrix of Young's modulus, G σ is the coefficient matrix of Poisson's ratio, and G ρ is the coefficient matrix of density.

[0103] Step 2.5: Introduce the Lagrange multiplier term R y , R x and the dual term C y , C x to obtain the prestack elastic three-parameter seismic inversion objective function based on fractional order total variation as shown in Equation 11:

[0104]

[0105] where η is the dual weight coefficient.

[0106] Step 3: Update L0 by combining the alternating direction multiplier method and the forward model to obtain the updated elastic three-parameter logarithm matrix L i+1 ; the step is specifically:

[0107] Step 3.1: Take the initial elastic three-parameter logarithm matrix L0 as the initial model L 1 of the prestack elastic three-parameter seismic inversion based on fractional order total variation, and introduce the initial value R x of the Lagrange multiplier term R y = 0 and its dual term C x = C y = 0;

[0108] Step 3.2: Update the elastic three-parameter logarithm matrix L i+1 using the alternating direction multiplier algorithm, the criterion of taking the extreme value when the gradient is zero, and eigenvalue decomposition, and the calculation principle is shown in Equation 12:

[0109]

[0110] where U and V are the orthogonal matrices obtained after eigenvalue decomposition of the square matrices and , as shown in Equation 13

[0111]

[0112] where SVD represents eigenvalue decomposition, and Λ and Σ are the diagonal matrices of eigenvalues obtained after decomposition.

[0113] In Equation 12, the update formula is shown in Equation 14

[0114]

[0115] Where a represents the number of eigenvalues ​​in the eigenvalue diagonal matrix Λ, and b represents the number of eigenvalues ​​in the eigenvalue diagonal matrix Σ.

[0116] Step 3.3: Update the Lagrange multiplier term R according to the alternating direction multiplier algorithm and the soft threshold shrinkage algorithm. x i+1 R y i+1 The calculation formula is shown in Formula 15:

[0117]

[0118] Step 3.4: Update the dual term according to the alternating direction multiplier algorithm and Fermat's Lemma. The calculation formula is shown in Formula 16:

[0119]

[0120] Step 4: Determine if the values ​​before and after the update satisfy ||L i+1 -L i ||2 / ||L i If ||2>tol, then return to step 3 and loop; otherwise, based on Young's modulus, Poisson's ratio, and the logarithmic matrix of the three parameters of density and elasticity, L... i+1 The relationship between Young's modulus, Poisson's ratio, and density is obtained as follows:

[0121] Based on the parameters to be inverted and the logarithmic matrix of the three elastic parameters L i+1 The parameters to be inverted are obtained from the relationship shown in Formula 17:

[0122]

[0123] Step 5: Slice the Young's modulus and Poisson's ratio inversion results by layer to obtain the crack prediction results.

[0124] like Figure 4 As shown, the inversion results of the pre-stack elastic three parameters better conform to the trend of the seismic record compared to the initial model, proving the correctness of the method. Figure 5The left side represents the crack prediction result corresponding to the Young's modulus inversion result, and the right side represents the crack prediction result corresponding to the Poisson's ratio inversion result. The application adopts a fractional order total variation to establish a sparse regularization term and construct a target function of prestack elastic three-parameter seismic inversion, further solves the target function through an alternating direction multiplier method to obtain the Young's modulus and Poisson's ratio inversion result, and obtains the crack prediction result through layer slicing along the inversion result, and proposes a prestack crack prediction method for lacustrine carbonate rocks based on fractional order total variation.

[0125] The above embodiment is a preferred embodiment of the present application, but the embodiments of the present application are not limited by the above embodiment, and any change, modification, substitution, combination, simplification made without departing from the spirit and principle of the present application should be an equivalent replacement mode, and all are included in the protection scope of the present application.

Claims

1. A method for predicting pre-stack fractures of lacustrine carbonate rocks based on fractional order total variation, characterized in that, The method comprises the following steps: Step 1: pre-processing seismic data to obtain an initial model of parameters to be inverted; Step 2: constructing a log matrix L0 of Young's modulus, Poisson's ratio and density according to the processed seismic data, and constructing a pre-stack elastic parameter forward model and an objective function based on fractional order total variation; Step 3: Inverse update the log matrix L by combining the alternating direction method of multipliers, taking extreme value when the gradient is zero, and eigenvalue decomposition i ; Step 4: judge whether the values before and after updating satisfy ||L i+1 -L i ||2 / ||L i ||2>tol, if yes, return to step 3 for a loop; if no, obtain Young's modulus, Poisson's ratio and density according to the relationship between the Young's modulus, Poisson's ratio and density and the elastic three-parameter logarithmic matrix L i+1 ​ Step 5: slicing the Young's modulus and Poisson's ratio inversion results according to horizons to obtain fracture prediction results; The step 2 of constructing a pre-stack elastic parameter forward model and an objective function based on fractional order total variation comprises the following steps: a longitudinal and transverse fractional order difference matrix and a longitudinal first order difference matrix are constructed by using a longitudinal and transverse fractional order difference vector, D = d * E (3m-1)×(3m-1) wherein * represents a convolution operation, represents a transverse fractional-order difference matrix, represents a longitudinal fractional-order difference matrix, D represents a longitudinal first-order difference matrix, d = [-1 1] represents a longitudinal first-order difference vector, d x a represents a transverse fractional-order difference vector, d y a represents a longitudinal fractional-order difference vector, m represents the total number of rows of a single parameter to be inverted; j represents the total number of columns of a single parameter to be inverted; a pre-stack elastic three-parameter forward model based on fractional order total variation is constructed: where ‖ ‖1 represents an L1 norm, ‖ ‖2 represents an L2 norm, * represents a convolution operation; S is a seismic record matrix, W is a wavelet matrix, and G represents a coefficient matrix of elastic three parameters to be inverted; On the basis of the forward model, the Lagrange multiplier term R y is introduced x and the dual term C y , C x is obtained. The prestack elastic three-parameter seismic inversion objective function based on fractional order total variation is obtained: where η is a dual item weight coefficient.

2. The method of claim 1, wherein, In step 1, the seismic data comprises seismic records, wavelet data and logging data; horizon information is extracted from the seismic records, and the logging data is subjected to interpolation filtering to obtain an initial model of Young's modulus, Poisson's ratio and density.

3. The method of claim 1, wherein, The relationship between the elastic three-parameter log matrix and Young's modulus, Poisson's ratio and density is represented by the following formula: L = [ln E ln σ ln p] T where L represents the elastic three-parameter log matrix, E represents Young's modulus, σ represents Poisson's ratio, ρ represents density, ln represents a logarithm operation, and T represents a transpose of a matrix.

4. The method of claim 1, wherein, The longitudinal and transverse fractional order difference vector is represented by the following formula: where ψ a (k) = (-1) k Γ(a + 1) / [Γ(k + 1)Γ(a - k + 1)]; 5. The method of claim 1, wherein, S is a seismic record S from three angles, mesial distal θ1 , S θ2 and S θ3 consists of, expressed by the following formula: S = [S θ1 ,S θ2 ,S θ3 ] T W is represented by the following formula: where w q represents the qth data of the seismic wavelet, q represents the length of the wavelet; G is represented by the following formula: where G E represents a coefficient matrix of the Young's modulus, G σ represents a coefficient matrix of the Poisson's ratio, G ρ represents a coefficient matrix of the density.

6. The method of claim 1, wherein, The step 3 comprises the following steps: Taking the initial elastic three-parameter log matrix L0 as the initial model L of the fractional order total variation prestack elastic three-parameter seismic inversion 1 And introducing the initial value R of the Lagrange multiplier term x = R y = 0 and its dual C x = C y = 0; The elastic three-parameter log-matrix L is updated by using an alternating direction multiplier algorithm, a gradient zero-takes-extreme-value criterion and eigenvalue decomposition i+1 : where U and V are square matrices and orthogonal matrix obtained after eigenvalue decomposition: where SVD represents eigenvalue decomposition, and Λ and Σ are eigenvalue diagonal matrices obtained after decomposition; The update formula is as follows: where a represents the number of eigenvalues in the eigenvalue diagonal matrix Λ, b represents the number of eigenvalues in the eigenvalue diagonal matrix Σ, The Lagrange multiplier term R is updated according to the alternating direction multiplier algorithm and the soft threshold shrinkage algorithm x i+1 , R y i+1 ; updating the dual according to the alternating direction method of multipliers and fermat's little theorem 7. The method of claim 6, wherein, Updated Lagrange multiplier term R x i+1 , R y i+1 is expressed by the following equation:

8. The method of claim 6, wherein, updated pair is shown by the following equation:

9. The method of claim 1, wherein, In step 4, the relationship between the parameters to be inverted and the elastic three-parameter log matrix L i+1 is expressed by the following equation:

Images