A method for azimuthally anisotropic seismic data inversion for dipping fracture rocks

Abstract

This invention is a method for inverting azimuth difference seismic data in inclined fractured rocks. Previous AVAZ inversion methods for fracture parameters were mostly based on the HTI medium assumption, with insufficient consideration for fracture dip angle. The fixed form of the inversion constraint terms made it difficult to adapt to changes in parameter statistical characteristics, and the regularization parameters were all manually adjusted, affecting inversion accuracy and efficiency. This invention establishes an approximate formula for the reflection coefficient of the TTI medium, characterized by the P-wave and S-wave velocities, density, and fracture elastic parameters of the background medium. It utilizes pre-stack azimuth difference data to perform AVAZ inversion of fracture parameters. Dynamic sparsity constraints are introduced, dynamically adjusting the sparsity of the constraints based on data characteristics, and a covariance matrix is ​​introduced to consider the statistical regularity between parameters. Combining iterative reweighted least squares and adaptive regularization parameter acquisition methods, the dynamic constraint inverse problem is transformed into a weighted... L A 2-norm objective function is used to achieve a stable solution. Data testing shows that this invention can effectively improve the inversion accuracy of crack elastic parameters.

Description

Technical Field

[0001] This invention belongs to the field of unconventional reservoir seismic exploration technology, specifically relating to a method for inverting seismic data on azimuth differences in inclined fractured rocks. Background Technology

[0002] Fractured reservoirs are a key type of unconventional oil and gas resource development reservoir in my country. Dipping fractures, a common fracture morphology in these reservoirs, significantly influence reservoir permeability and hydrocarbon accumulation patterns. Accurate acquisition of dipped fracture parameter information is crucial for guiding unconventional oil and gas exploration deployment and development planning. Seismic inversion technology is the core means of quantitatively characterizing subsurface fracture parameters. Wide-azimuth pre-stack seismic data, containing rich azimuth anisotropy information, serves as an important data foundation for fracture seismic prediction. Azimuth difference seismic data inversion methods effectively extract anisotropic signals caused by fractures by utilizing amplitude variations with azimuth and incident angles, and are widely used in fracture seismic prediction and inversion. Currently, most related studies construct inversion frameworks based on the assumption of a transversely isotropic (HTI) medium with a horizontal axis of symmetry, but this does not consider the influence of fracture dip angle on inversion accuracy. Therefore, in recent years, seismic inversion based on a transversely isotropic (TTI) medium with an inclined axis of symmetry has been proposed and applied to seismic prediction of rocks containing dipped fractures.

[0003] To address the ill-posed problem of inversion, regularization constraints are incorporated into the inversion process, including smoothing constraints based on the L2 norm, sparsity constraints based on the L1 norm, and Cauchy norm constraints, among others, to suppress noise interference and characterize the distribution features of reservoir parameters. In selecting regularization parameters, existing methods often rely on experience from existing work areas, setting fixed ranges for regularization parameters and then manually adjusting the parameters through experiments to determine the final values. However, existing methods for inverting tilted fracture orientation differences still face several unresolved issues in practical applications: First, the regularization constraints lack adaptability; existing methods often use fixed L1 or L2 norms, failing to dynamically adjust the sparsity of constraint terms based on the characteristics of the parameters to be inverted. Second, the selection of regularization parameters depends on experience; existing methods often set candidate parameter ranges manually, making it difficult to adapt to the differences in geological conditions across different work areas. Since the parameters of inclined fractures are not highly sensitive to seismic activity, and the azimuth difference data is easily contaminated by noise, the accuracy of AVAZ inversion for TTI media is insufficient. Therefore, it is necessary to improve the inversion optimization strategy to enhance the estimation accuracy of fracture elastic parameters and the adaptability to actual data. Thus, in order to address the difficulties in seismic prediction of fracture parameters in actual exploration, it is urgent to develop an adaptive and high-precision azimuth difference seismic data inversion method suitable for rocks with inclined fractures. This is of great significance for promoting high-precision exploration of fractured unconventional oil and gas reservoirs. Summary of the Invention

[0004] This invention addresses the problems existing in the prior art by providing a method for inverting seismic data on azimuth differences in inclined fractured rocks, which can improve the accuracy of seismic inversion of inclined fractured rocks.

[0005] To solve the above technical problems, the present invention provides the following technical solution: a method for inverting seismic data on azimuth differences in inclined fractured rocks, comprising the following steps:

[0006] Step 1: Calculate azimuth difference seismic data based on wide-azimuth pre-stack seismic data. Seismic wavelets with different azimuth and incident angles were set, and a seismic wavelet matrix was established. Establish an initial model for initial crack weakness parameters. Set the candidate regularization parameters;

[0007] Step 2: Based on the TTI medium reflection coefficient equation and azimuth difference data, construct data error terms and regularization constraint terms to establish the inversion objective function;

[0008] Step 3: Introduce a data-driven dynamic sparse constraint method to calculate the kurtosis coefficient of the current model parameters, dynamically select the exponential parameter in the constraint term, construct a regularization constraint term for the statistical characteristics of the adaptive model, and introduce the standard deviation estimate based on the sample covariance to consider the statistical regularity between parameters.

[0009] Step 4: Iterative optimization is performed using the iterative reweighted least squares framework to transform the non-L2 norm regularization problem of the inversion objective function into a weighted least squares subproblem. The weight matrix is ​​calculated in the outer iteration, and the optimal regularization parameter is obtained by using an adaptive estimation method for the regularization parameter, thereby obtaining a new weighted L2 norm inversion objective function.

[0010] Step 5: Use the pre-defined conjugate gradient method to solve the inversion objective function through inner-layer iteration;

[0011] Step 6: Return to steps 3 through 5 until the preset conditions are met, and output the final crack normal weakness. and crack tangential weakness Calculation results.

[0012] Furthermore, step 1 described above includes the following sub-steps:

[0013] Step 1.1: Calculate azimuth difference seismic data based on wide-azimuth pre-stack seismic data. :

[0014] (1),

[0015] in, Angle of incidence Let i be the azimuth angle. Let j be the j-th azimuth angle. For the incident angle is The i-th azimuth angle is The corresponding pre-stack seismic data, For the incident angle is The j-th azimuth angle is Corresponding pre-stack seismic data;

[0016] Step 1.2: Extract seismic wavelets with different azimuth and incident angles from the pre-stack azimuth seismic data, and establish the seismic wavelet matrix. ;

[0017] Step 1.3: Establish the initial crack parameter model :

[0018] (2),

[0019] Where n is the number of crack weakness parameters, For crack normal weakness, The tangential weakness of the crack; where, Crack density, For intermediate variables, where For the longitudinal wave velocity of the isotropic background medium, The transverse wave velocity is given by the isotropic background medium. and The bulk modulus and shear modulus of the fluid filling the fracture. The aspect ratio of the crack. Given the shear modulus of the background medium, the fracture normal weakness and fracture tangential weakness at the wellhead are calculated using the following formula:

[0020] (3),

[0021] (4),

[0022] Step 1.4: Based on the calculated fracture normal and tangential weaknesses at the wellhead, a low-frequency initial model at the wellhead is obtained through Backus averaging. For initial two-dimensional or three-dimensional models, hierarchical interpretation constraints are obtained through kriging interpolation.

[0023] Step 1.5: Set a set of candidate regularization parameter values. , Let be the regularization parameter for the j-th candidate. The number of candidate regularization parameters.

[0024] Furthermore, step 2 described above includes the following sub-steps:

[0025] Step 2.1: TTI medium PP wave reflection coefficient expressed by crack weakness, background medium elastic parameters, and crack dip angle. It includes the isotropic term of the background medium rock. and the anisotropy term caused by cracks The expression is as follows:

[0026] (5),

[0027] in, Angle of incidence It is the azimuth angle. The angle of inclination of the crack;

[0028] Step 2.2: Calculate the anisotropy term caused by the crack. The following expression:

[0029] (6),

[0030] In the formula, The coefficient for the change in normal crack weakness. This is the coefficient for the change in tangential crack weakness. This represents the change in the normal crack weakness between the upper and lower media. This represents the change in the tangential crack weakness between the upper and lower media. Step 2.3: When the fractured rock background is an isotropic medium, It does not change with azimuth angle, and uses azimuth difference data. and TTI medium reflectivity Establish the inversion objective function:

[0031] (7),

[0032] in,

[0033] , (8),

[0034] In the formula, G represents the forward modeling coefficient matrix of the TTI medium. This represents the forward modeling coefficient matrix of the normal crack weakness. This represents the forward modeling coefficient matrix of tangential crack weakness. Indicates the parameter to be determined. This represents the difference vector of weakness along the crack normal. This represents the vector of tangential weakness differences in the crack. Represents the seismic wavelet matrix; This indicates that the L2 norm is used to calculate the data error term to meet the characteristics of the noise distribution in seismic data; This represents a regularization constraint. This represents the regularization parameter.

[0035] Furthermore, the aforementioned expression for isotropic terms as follows:

[0036] (9),

[0037] , , (10),

[0038] In the formula, The coefficient of longitudinal wave velocity reflectivity. For longitudinal wave velocity reflectivity, This is the coefficient of the transverse wave velocity reflectivity. For transverse wave velocity reflectivity, The coefficient of density reflectivity. Density reflectance.

[0039] Furthermore, the coefficient of the aforementioned normal crack weakness variation... Coefficient of tangential crack weakness variation The expression is as follows:

[0040] (11), (12).

[0041] Furthermore, step 3 mentioned above specifically involves: calculating the statistical characteristics of the distribution pattern of the current model parameters in real time, i.e., the kurtosis coefficient, and dynamically selecting the exponential parameter β(m) of the sparse constraint norm based on the numerical value of this characteristic, so that the mathematical form of the constraint term can be flexibly adjusted according to the statistical characteristics of the model parameters themselves; introducing a scaling term composed of the estimated standard deviation values ​​of each parameter into the constraint term, which is calculated from the sample covariance matrix of the model parameters;

[0042] A data-driven elastic sparse regularization framework is adopted. The core idea of ​​this framework is to dynamically adjust the sparse constraint strength of the regularization term based on the real-time statistical characteristics of the model parameter distribution, thereby improving the solution accuracy and robustness of the inversion problem. The inversion objective function then becomes as follows:

[0043] (13),

[0044] in, It is the k-th component of the model parameter vector m. The regularization weight factor for the k-th parameter is not fixed, but is a function of the higher-order statistical properties of the model parameter vector m, which realizes the smooth switching of the regularization mode between the smooth L2 norm and the sparse L1 norm.

[0045] The dynamic parameter β(m) is adaptively determined by the following formula:

[0046] (14),

[0047] In the formula, the statistic ϑ(m) used for decision-making is the kurtosis coefficient, and its calculation formula is:

[0048] (15),

[0049] in, The second central moment of the parameter vector m. The fourth central moment of the parameter vector m, Let m be the mean of the parameter vector. The decision threshold is set to 3 if we want it to correspond to the kurtosis of the standard normal distribution.

[0050] When ϑ(m) > 3, it indicates that the parameter distribution has sharp peaks and thick tails, automatically activating the strong sparsity promotion mode, i.e., using L... 0.3 The quasi-norm; when ϑ(m)≤3, it indicates that the parameter distribution is relatively flat or close to a Gaussian distribution, and the algorithm switches to a robust sparse mode, that is, adopts the L1 norm;

[0051] A scaling parameter based on sample covariance estimation is introduced into the constraint terms. The covariance matrix of the model parameters can be calculated using the following formula:

[0052] (16),

[0053] in, Let s be the model vector of the s-th sample. The sample mean. Let C be the total number of samples, C be the covariance matrix, which characterizes the volatility and interrelationships among the parameters, and the scale parameter be... The parameter m is defined as the square root of the k-th element on the diagonal of the covariance matrix. k The sample standard deviation estimate can be obtained by the following formula:

[0054] (17),

[0055] In the formula, Let be the element in the k-th row and k-th column of the covariance matrix C.

[0056] Furthermore, step 4 described above includes the following sub-steps:

[0057] Step 4.1: In the adaptive iterative reweighting algorithm, a two-layer iterative architecture is adopted. The outer layer iterates and updates the weight matrix and regularization parameters. In the k-th iteration of the outer layer iterates, the regularization constraint term is... By introducing a weight matrix The local approximation is as follows:

[0058] (18),

[0059] Through the above approximation, the original objective function is transformed into a weighted least squares problem with respect to m in each outer iteration;

[0060] In the formula, The model parameter vector for the kth iteration; In the k-th iteration, the i-th model parameter component; β (k) Let be the dynamic regularization exponent parameter for the k-th iteration; In the k-th iteration, the weight coefficient corresponding to the i-th model parameter component. It is a diagonal weight matrix. is the dimension of the weight matrix;

[0061] Step 4.2: Determine the optimal regularization parameter using an adaptive estimation method for the regularization parameter. Based on the candidate regularization parameter values ​​set in step 1 Calculate the current The corresponding regularization parameter adaptive selection criterion function CF value:

[0062] (19),

[0063] in, The total number of single-channel data points participating in the inversion. For identity matrix, influence matrix In the k-th iteration, it takes the form:

[0064] (20),

[0065] Step 4.3: Calculate the CF value corresponding to all candidate regularization parameters, and select the one with the smallest CF value. As the optimal regularization parameter At this point, the inversion objective function is locally approximated as a weighted L2 norm regularization problem, forming a new inversion objective function:

[0066] (twenty one),

[0067] in, Let be the inversion objective function for the model parameter vector at the k-th iteration.

[0068] Furthermore, the aforementioned diagonal element The expression is as follows:

[0069] (twenty two),

[0070] In the formula, The scaling parameter is the i-th model parameter component estimated based on the parameter covariance.

[0071] Furthermore, step 5 mentioned above includes the following sub-steps:

[0072] Step 5.1: Initialize the parameters of the inner iteration. Let the initial model at the inner iteration number t=0 be... The initial residual is calculated using the following formula. :

[0073] (twenty three),

[0074] in, As an intermediate variable, This is the model parameter vector when the inner layer iteration count is 0 and the outer layer iteration count is k.

[0075] Step 5.2: Perform preconditioning by extracting intermediate variables. The inverse of the diagonal matrix formed by the diagonal elements yields the precondition matrix. The residuals are preprocessed using the following formula:

[0076] (twenty four),

[0077] in, Let be the residual of the t-th inner iteration. It is a residual Representation in the preprocessing space;

[0078] Step 5.3: Determine the search direction The update is performed, and the search direction for the t-th inner iteration is as follows:

[0079] (25),

[0080] in, Let be the conjugate gradient update coefficients of the t-th inner iteration;

[0081] Step 5.4: Update the model vector and residual vector using the following formula:

[0082] (26),

[0083] (27),

[0084] in, Let be the step size, and its calculation formula is: .

[0085] Furthermore, step 6 above specifically involves: conducting a double-layer iterative loop of outer and inner layers to update the model parameter vector. When the model parameter m of the k-th outer layer iteration satisfies the following formula, the output model parameters are the final inversion results of the crack normal weakness and tangential weakness:

[0086] (28),

[0087] in, This is the preset convergence threshold.

[0088] Compared with the prior art, the beneficial technical effects of the present invention using the above technical solution are as follows:

[0089] Previous AVAZ inversion methods for fracture parameters were mostly based on the HTI medium assumption. While achieving good application results, these methods were only applicable to vertical or high-angle fractures, failing to consider the influence of fracture dip angle and thus unsuitable for reservoirs with numerous diagonal fractures. Furthermore, the inversion algorithms used fixed constraint terms, making it difficult to adapt to changes in parameter statistical characteristics. Additionally, the selection of regularization parameters directly affects the inversion results, but in previous AVAZ inversions, regularization parameters were determined manually through repeated adjustments based on experience, affecting both accuracy and computational efficiency. This invention establishes an approximate formula for the reflection coefficient of the TTI medium, characterized by the P-wave and S-wave velocities, density, and fracture elastic parameters of the background medium. Using pre-stack positional difference data, AVAZ inversion is performed on diagonally fractured rocks. A data-driven dynamic sparsity constraint is introduced, dynamically adjusting the sparsity of the constraint based on parameter data characteristics, and a covariance matrix is ​​introduced to account for statistical regularities between parameters. Combining iterative reweighted least squares and adaptive regularization parameter acquisition methods, the dynamic constraint inverse problem is transformed into a weighted L2 norm objective function, achieving a stable solution for fracture weakness parameters. Synthetic data testing demonstrates that this invention effectively improves the inversion accuracy of fracture elastic parameters. This invention can provide crucial support for fracture seismic prediction; the predicted fracture weakness parameters can be directly used to estimate fracture density, fracture fluid information, etc., and can also support subsequent prediction of fractured reservoir in-situ stress. It has significant practical application value. Attached Figure Description

[0090] Figure 1 This is a flowchart illustrating the implementation of the algorithm of this invention.

[0091] Figure 2(a) shows the AVAZ inversion results of TTI media with manually selected regularization parameters combined with conventional L2 norm constraints, and Figure 2(b) shows the AVAZ inversion results of TTI media with the adaptive regularization parameter method combined with conventional L2 norm constraints.

[0092] Figure 3(a) shows the inversion results of the adaptive regularization parameter acquisition method combined with conventional L2 norm constraints, and Figure 3(b) shows the inversion results of the adaptive regularization parameter acquisition method combined with data-driven dynamic sparse constraints.

[0093] Figure 4 shows the noise resistance test of the TTI inversion method (adaptive regularization parameter acquisition method combined with data-driven dynamic sparse constraints) of the present invention: (a) TTI medium inversion results with no noise input data, Figure 4(b) TTI medium inversion results with a signal-to-noise ratio of 12 input data, Figure 4(c) TTI medium inversion results with a signal-to-noise ratio of 6 input data. Detailed Implementation

[0094] To better understand the technical content of the present invention, specific embodiments are described below in conjunction with the accompanying drawings.

[0095] In this invention, various aspects of the invention are described with reference to the accompanying drawings, in which numerous illustrative embodiments are shown. Embodiments of the invention are not limited to those depicted in the drawings. It should be understood that the invention is implemented through any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed herein are not limited to any particular implementation. Furthermore, some aspects of the invention disclosed may be used alone or in any suitable combination with other aspects of the invention disclosed.

[0096] like Figure 1 As shown, this invention provides a method for inverting seismic data on azimuth differences in inclined fractured rocks, with the following steps:

[0097] Step 1: Calculate azimuth difference seismic data based on wide-azimuth pre-stack seismic data. Seismic wavelets with different azimuth and incident angles were set, and a seismic wavelet matrix was established. Establish an initial model for initial crack weakness parameters. Set the candidate regularization parameters;

[0098] Step 2: Based on the TTI medium reflection coefficient equation and azimuth difference data, construct data error terms and regularization constraint terms to establish the inversion objective function;

[0099] Step 3: Introduce a data-driven dynamic sparse constraint method to calculate the kurtosis coefficient of the current model parameters, dynamically select the exponential parameter in the constraint term, construct a regularization constraint term that is adaptive to the statistical characteristics of the model, and introduce the standard deviation estimate based on the sample covariance to consider the statistical regularity between parameters.

[0100] Step 4: Iterative optimization is performed using the iterative reweighted least squares framework to transform the non-L2 norm regularization problem of the inversion objective function into a weighted least squares subproblem. The weight matrix is ​​calculated in the outer iteration, and the optimal regularization parameter is obtained by using an adaptive estimation method for the regularization parameter, thereby obtaining a new weighted L2 norm inversion objective function.

[0101] Step 5: Use the pre-defined conjugate gradient method to solve the inversion objective function through inner-layer iteration;

[0102] Step 6: Return to steps 3 through 5 until the preset conditions are met, and output the final crack normal weakness. and crack tangential weakness Calculation results.

[0103] In a preferred embodiment of the present invention, step 1 includes the following sub-steps:

[0104] Step 1.1: Calculate azimuth difference seismic data based on wide-azimuth pre-stack seismic data. :

[0105] (1),

[0106] in, Angle of incidence Let i be the azimuth angle. Let j be the j-th azimuth angle. The angle of incidence is The i-th azimuth angle The corresponding pre-stack seismic data, The angle of incidence is The j-th azimuth angle Corresponding pre-stack seismic data;

[0107] Step 1.2: Extract seismic wavelets with different azimuth and incident angles from the pre-stack azimuth seismic data, and establish the seismic wavelet matrix. ;

[0108] Step 1.3: Establish the initial crack parameter model :

[0109] (2),

[0110] Where n is the number of crack weakness parameters, For crack normal weakness, The fracture tangential weakness is given by the formula; the fracture normal weakness and fracture tangential weakness at the wellhead are calculated using the following formula:

[0111] (3),

[0112] (4),

[0113] Step 1.4: Based on the calculated fracture normal and tangential weaknesses at the wellhead, a low-frequency initial model at the wellhead is obtained through Backus averaging. For two-dimensional or three-dimensional initial models, hierarchical interpretation constraints are obtained through kriging interpolation.

[0114] Step 1.5: Set a set of candidate regularization parameter values. , Let be the regularization parameter for the j-th candidate. The number of candidate regularization parameters.

[0115] In a preferred embodiment of the present invention, step 2 includes the following sub-steps:

[0116] Step 2.1: TTI medium PP wave reflection coefficient expressed by crack weakness, background medium elastic parameters, and crack dip angle. It includes the isotropic term of the background medium rock. and the anisotropy term caused by cracks The expression is as follows:

[0117] (5),

[0118] in, Angle of incidence It is the azimuth angle. The angle of inclination of the crack;

[0119] Step 2.2: Calculate the anisotropy term caused by the crack. The following expression:

[0120] (6),

[0121] In the formula, The coefficient for the change in normal crack weakness. This is the coefficient for the change in tangential crack weakness. This represents the change in the normal crack weakness between the upper and lower media. This represents the change in the tangential crack weakness between the upper and lower media.

[0122] Step 2.3: When the fractured rock background is an isotropic medium, It does not change with azimuth angle, and uses azimuth difference data. and TTI medium reflectivity Establish the inversion objective function:

[0123] (7),

[0124] in,

[0125] , (8),

[0126] In the formula, G represents the forward modeling coefficient matrix of the TTI medium. This represents the forward modeling coefficient matrix of the normal crack weakness. This represents the forward modeling coefficient matrix of tangential crack weakness. Indicates the parameter to be determined. This represents the difference vector of weakness along the crack normal. This represents the vector of tangential weakness differences in the crack. Represents the seismic wavelet matrix; This indicates that the L2 norm is used to calculate the data error term to meet the characteristics of the noise distribution in seismic data; This represents a regularization constraint. This represents the regularization parameter.

[0127] As a preferred embodiment of the present invention, the expression for isotropic terms... as follows:

[0128] (9),

[0129] , , (10),

[0130] In the formula, The coefficient of longitudinal wave velocity reflectivity. For longitudinal wave velocity reflectivity, This is the coefficient of the transverse wave velocity reflectivity. For transverse wave velocity reflectivity, The coefficient of density reflectivity. Density reflectance.

[0131] As a preferred embodiment of the present invention, the coefficient of the normal crack weakness variation Coefficient of tangential crack weakness variation The expression is as follows: (11),

[0132] (12).

[0133] As a preferred embodiment of the present invention, step 3 specifically involves: constructing and introducing a regularization constraint term with statistical adaptive capability, specifically including: calculating the statistical characteristics of the distribution pattern of the current model parameters in real time, i.e., the kurtosis coefficient, and dynamically selecting the exponential parameter β(m) of the sparse constraint norm based on the numerical value of the characteristic, so that the mathematical form of the constraint term can be flexibly adjusted according to the statistical characteristics of the model parameters themselves.

[0134] Furthermore, to effectively describe and utilize the statistical correlation between model parameters, a scaling term consisting of the estimated standard deviations of each parameter is introduced into the constraint terms. This term is calculated from the sample covariance matrix of the model parameters.

[0135] A data-driven elastic sparse regularization framework is adopted. The core idea of ​​this framework is to dynamically adjust the sparse constraint strength of the regularization term based on the real-time statistical characteristics of the model parameter distribution, thereby improving the solution accuracy and robustness of the inversion problem. The inversion objective function then becomes as follows:

[0136] (13),

[0137] in, It is the k-th component of the model parameter vector m. The regularization weight factor for the k-th parameter is not fixed, but is a function of the higher-order statistical properties of the model parameter vector m, thus realizing the adaptive switching of the regularization mode between favoring smooth constraints and favoring sparse constraints.

[0138] The dynamic parameter β(m) is adaptively determined by the following formula:

[0139] (14),

[0140] In the formula, the statistic ϑ(m) used for decision-making is the kurtosis coefficient, and its calculation formula is:

[0141] (15),

[0142] in, The second central moment of the parameter vector m. The fourth central moment of the parameter vector m, Let m be the mean of the parameter vector. The decision threshold is set to 3 if we want it to correspond to the kurtosis of the standard normal distribution.

[0143] When ϑ(m) > 3, it indicates that the parameter distribution has sharp peaks and thick tails, automatically activating the strong sparsity promotion mode, i.e., using L... 0.3Quasi-norm; when ϑ(m)≤3, it indicates that the parameter distribution is relatively flat or close to Gaussian distribution, and the algorithm switches to robust sparse mode, that is, adopts L1 norm.

[0144] To effectively integrate the prior spatial structure information and statistical regularities of the model parameters, a scaling parameter based on sample covariance estimation was introduced into the constraint terms. The covariance matrix of the model parameters can be calculated using the following formula:

[0145] (16),

[0146] in, Let s be the model vector of the s-th sample. The sample mean. Let C be the total number of samples, C be the covariance matrix, which characterizes the volatility and interrelationships among the parameters, and the scale parameter be... The parameter m is defined as the square root of the k-th element on the diagonal of the covariance matrix. k The sample standard deviation estimate can be obtained by the following formula:

[0147] (17),

[0148] In the formula, Let be the element in the k-th row and k-th column of the covariance matrix C.

[0149] In a preferred embodiment of the present invention, step 4 includes the following sub-steps:

[0150] Step 4.1: In the adaptive iterative reweighting algorithm, a two-layer iterative architecture is adopted. The outer layer iterates and updates the weight matrix and regularization parameters. In the k-th iteration of the outer layer iterates, the regularization constraint term is... By introducing a weight matrix The local approximation is a quadratic form, as shown in the following expression:

[0151] (18),

[0152] Through the above approximation, the original objective function is transformed into a weighted least squares problem with respect to m in each outer iteration;

[0153] In the formula, The model parameter vector for the kth iteration; In the k-th iteration, the i-th model parameter component; β (k) Let be the dynamic regularization exponent parameter for the k-th iteration; In the k-th iteration, the weight coefficient corresponding to the i-th model parameter component. It is a diagonal weight matrix. is the dimension of the weight matrix;

[0154] Step 4.2: Determine the optimal regularization parameter using an adaptive estimation method for the regularization parameter. Based on the candidate regularization parameter values ​​set in step 1 Calculate the current The corresponding regularization parameter adaptive selection criterion function CF value:

[0155] (19),

[0156] in, The total number of single-channel data points participating in the inversion. For identity matrix, influence matrix In the k-th iteration, it takes the form:

[0157] (20),

[0158] Step 4.3: Calculate the CF value corresponding to all candidate regularization parameters, and select the one with the smallest CF value. As the optimal regularization parameter At this point, the inversion objective function is locally approximated as a weighted L2 norm regularization problem, forming a new inversion objective function:

[0159] (twenty one),

[0160] in, Let be the inversion objective function for the model parameter vector at the k-th iteration.

[0161] As a preferred embodiment of the present invention, the diagonal element The expression is as follows:

[0162] (twenty two),

[0163] In the formula, The scaling parameter is the i-th model parameter component estimated based on the parameter covariance.

[0164] In a preferred embodiment of the present invention, step 5 includes the following sub-steps:

[0165] Step 5.1: Initialize the parameters of the inner iteration. Let the initial model at the inner iteration number t=0 be... The initial residual is calculated using the following formula. :

[0166] (twenty three),

[0167] in, As an intermediate variable, This is the model parameter vector when the inner layer iteration count is 0 and the outer layer iteration count is k.

[0168] Step 5.2: Perform preconditioning by extracting intermediate variables. The inverse of the diagonal matrix formed by the diagonal elements yields the precondition matrix. The residuals are preprocessed using the following formula:

[0169] (twenty four),

[0170] in, Let be the residual of the t-th inner iteration. It is a residual Representation in the preprocessing space;

[0171] Step 5.3: Determine the search direction The update is performed, and the search direction for the t-th inner iteration is as follows:

[0172] (25),

[0173] in, Let be the conjugate gradient update coefficients of the t-th inner iteration;

[0174] Step 5.4: Update the model vector and residual vector using the following formula:

[0175] (26),

[0176] (27),

[0177] in, Let be the step size, and its calculation formula is: .

[0178] In a preferred embodiment of the present invention, step 6 specifically involves: conducting a double-layer iterative loop of outer and inner layers, updating the model parameter vector, and when the model parameter m of the k-th outer layer iteration satisfies the following formula, the output model parameters are the final inversion results of the crack normal weakness and tangential weakness:

[0179] (28).

[0180] in, This is the preset convergence threshold.

[0181] Figure 2(a) shows the AVAZ inversion results for TTI media using manually selected regularization parameters combined with conventional L2 norm constraints, and Figure 2(b) shows the AVAZ inversion results for TTI media using an adaptive regularization parameter method combined with conventional L2 norm constraints. In the figures, the black solid line represents the true value, the blue dashed line represents the initial model, and the pink solid line represents the predicted result. A comparison shows that the inversion results using manually selected regularization parameters have significant errors, particularly in the crack normal weakness δ. N and tangential weakness δ T The inversion accuracy is lower than that of adaptively acquiring regularization parameters.

[0182] Figure 3(a) shows the inversion result of the adaptive regularization parameter acquisition method combined with the conventional L2 norm, and Figure 3(b) shows the result of the adaptive regularization parameter acquisition method combined with data-driven dynamic sparse constraints. In the figures, the black solid line represents the true value, the blue dashed line represents the initial model, and the pink solid line represents the prediction result. It can be seen from the figures that the data-driven dynamic sparse constraints are superior to the conventional L2 norm constraints, and can obtain a higher accuracy in the crack normal weakness δ. N and tangential weakness δ T Inversion results.

[0183] Figure 4 shows the noise resistance test results of the TTI inversion (adaptive regularization parameter acquisition method combined with data-driven dynamic sparse constraints) proposed in this invention: (a) The input data is the TTI medium inversion result under noise-free conditions; Figure 4(b) The input data is the TTI medium inversion result with a signal-to-noise ratio of 12; Figure 4(c) The input data is the TTI medium inversion result with a signal-to-noise ratio of 6. In the figures, the black solid line represents the true value, the blue dashed line represents the initial model, and the pink solid line represents the prediction result. It can be seen that under noise-free conditions, δ N and δ T The inversion accuracy is relatively high; as the noise increases, δ T The impact is greater than δ N The overall inversion accuracy is reduced, but even under high noise conditions, the inversion results are still within an acceptable range, indicating that the present invention has good noise resistance.

[0184] While the present invention has been described above with reference to preferred embodiments, it is not intended to limit the invention. Those skilled in the art can make various modifications and refinements without departing from the spirit and scope of the invention. Therefore, the scope of protection of the present invention shall be determined by the claims.

Claims

1. A method for inverting seismic data on azimuth differences in inclined fractured rocks, characterized in that, The steps are as follows: Step 1: Calculate azimuth difference seismic data based on wide-azimuth pre-stack seismic data. Seismic wavelets with different azimuth and incident angles were set, and a seismic wavelet matrix was established. Establish an initial model for initial crack weakness parameters. Set the candidate regularization parameters; Step 2: Based on the TTI medium reflection coefficient equation and azimuth difference data, construct data error terms and regularization constraint terms to establish the inversion objective function; Step 3: Introduce a data-driven dynamic sparse constraint method to calculate the kurtosis coefficient of the current model parameters, dynamically select the exponential parameter in the constraint term, construct a regularization constraint term for the statistical characteristics of the adaptive model, and introduce the standard deviation estimate based on the sample covariance to consider the statistical regularity between parameters. Step 4: Iterative optimization is performed using the iterative reweighted least squares framework to transform the non-L2 norm regularization problem of the inversion objective function into a weighted least squares subproblem. The weight matrix is ​​calculated in the outer iteration, and the optimal regularization parameter is obtained by using an adaptive estimation method for the regularization parameter, thereby obtaining a new weighted L2 norm inversion objective function. Step 5: Use the pre-defined conjugate gradient method to solve the inversion objective function through inner-layer iteration; Step 6: Return to steps 3 through 5 until the preset conditions are met, and output the final crack normal weakness. and crack tangential weakness Calculation results.

2. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 1, characterized in that, Step 1 includes the following sub-steps: Step 1.1: Calculate azimuth difference seismic data based on wide-azimuth pre-stack seismic data. : (1)， in, Angle of incidence Let i be the azimuth angle. Let j be the j-th azimuth angle. For the incident angle is The i-th azimuth angle is The corresponding pre-stack seismic data, For the incident angle is The j-th azimuth angle is Corresponding pre-stack seismic data; Step 1.2: Extract seismic wavelets with different azimuth and incident angles from the pre-stack azimuth seismic data, and establish the seismic wavelet matrix. ; Step 1.3: Establish the initial crack parameter model : (2)， Where n is the number of crack weakness parameters, For crack normal weakness, The tangential weakness of the crack; where, Crack density, For intermediate variables, where For the longitudinal wave velocity of the isotropic background medium, The transverse wave velocity is given by the isotropic background medium. and The bulk modulus and shear modulus of the fluid filling the fracture. The aspect ratio of the crack. Given the shear modulus of the background medium, the fracture normal weakness and fracture tangential weakness at the wellhead are calculated using the following formula: (3)， (4)， Step 1.4: Based on the calculated fracture normal and tangential weaknesses at the wellhead, a low-frequency initial model at the wellhead is obtained through Backus averaging. For two-dimensional or three-dimensional initial models, hierarchical interpretation constraints are obtained through kriging interpolation. Step 1.5: Set a set of candidate regularization parameter values. , Let be the regularization parameter for the j-th candidate. The number of candidate regularization parameters.

3. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 1, characterized in that, Step 2 includes the following sub-steps: Step 2.1: TTI medium PP wave reflection coefficient expressed by crack weakness, background medium elastic parameters, and crack dip angle. It includes the isotropic term of the background medium rock. and the anisotropy term caused by cracks The expression is as follows: (5)， in, Angle of incidence It is the azimuth angle. The angle of inclination of the crack; Step 2.2: Calculate the anisotropy term caused by the crack. The following expression: (6)， In the formula, The coefficient for the change in normal crack weakness. This is the coefficient for the change in tangential crack weakness. This represents the change in the normal crack weakness between the upper and lower media. This represents the change in the tangential crack weakness between the upper and lower media. Step 2.3: When the fractured rock background is an isotropic medium, It does not change with azimuth angle, and uses azimuth difference data. and TTI medium reflectivity Establish the inversion objective function: (7)， in, ， (8)， In the formula, G represents the forward modeling coefficient matrix of the TTI medium. This represents the forward modeling coefficient matrix of the normal crack weakness. This represents the forward modeling coefficient matrix of tangential crack weakness. Indicates the parameter to be determined. This represents the difference vector of weakness along the crack normal. This represents the vector of tangential weakness differences in the crack. Represents the seismic wavelet matrix; This indicates that the L2 norm is used to calculate the data error term to meet the characteristics of the noise distribution in seismic data; This represents a regularization constraint. This represents the regularization parameter.

4. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 3, characterized in that, Expression of isotropic terms as follows: (9)， ， ， (10)， In the formula, The coefficient of longitudinal wave velocity reflectivity. For longitudinal wave velocity reflectivity, This is the coefficient of the transverse wave velocity reflectivity. For transverse wave velocity reflectivity, The coefficient of density reflectivity. Density reflectance.

5. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 3, characterized in that, Coefficient of normal crack weakness variation Coefficient of tangential crack weakness variation The expression is as follows: (11)， (12)。 6. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 3, characterized in that, Step 3 specifically involves: calculating the statistical characteristics of the current model parameters' distribution pattern, i.e., the kurtosis coefficient, in real time, and dynamically selecting the exponential parameter β(m) of the sparse constraint norm based on the numerical value of this characteristic, so that the mathematical form of the constraint term can be flexibly adjusted according to the statistical characteristics of the model parameters themselves; introducing a scaling term composed of the estimated standard deviation values ​​of each parameter into the constraint term, which is calculated from the sample covariance matrix of the model parameters. A data-driven elastic sparse regularization framework is adopted. The core idea of ​​this framework is to dynamically adjust the sparse constraint strength of the regularization term based on the real-time statistical characteristics of the model parameter distribution, thereby improving the solution accuracy and robustness of the inversion problem. The inversion objective function then becomes as follows: (13)， in, It is the k-th component of the model parameter vector m. The regularization weight factor for the k-th parameter is not fixed, but is a function of the higher-order statistical properties of the model parameter vector m, which realizes the smooth switching of the regularization mode between the smooth L2 norm and the sparse L1 norm. The dynamic parameter β(m) is adaptively determined by the following formula: (14)， In the formula, the statistic ϑ(m) used for decision-making is the kurtosis coefficient, and its calculation formula is: (15)， in, The second central moment of the parameter vector m. The fourth central moment of the parameter vector m, The mean of the parameter vector m; The decision threshold is set to 3 if we want it to correspond to the kurtosis of the standard normal distribution. When ϑ(m) > 3, it indicates that the parameter distribution has sharp peaks and thick tails, automatically activating the strong sparsity promotion mode, i.e., using L... 0.3 The quasi-norm; when ϑ(m)≤3, it indicates that the parameter distribution is relatively flat or close to a Gaussian distribution, and the algorithm switches to a robust sparse mode, that is, adopts the L1 norm; A scaling parameter based on sample covariance estimation is introduced into the constraint terms. The covariance matrix of the model parameters can be calculated using the following formula: (16)， in, Let s be the model vector of the s-th sample. The sample mean. Let C be the total number of samples, C be the covariance matrix, which characterizes the volatility and interrelationships among the parameters, and the scale parameter be... The parameter m is defined as the square root of the k-th element on the diagonal of the covariance matrix. k The sample standard deviation estimate can be obtained by the following formula: (17)， In the formula, Let be the element in the k-th row and k-th column of the covariance matrix C.

7. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 2, characterized in that, Step 4 includes the following sub-steps: Step 4.1: In the adaptive iterative reweighting algorithm, a two-layer iterative architecture is adopted. The outer layer iterates and updates the weight matrix and regularization parameters. In the k-th iteration of the outer layer iterates, the regularization constraint term is... By introducing a weight matrix The local approximation is as follows: (18)， Through the above approximation, the original objective function is transformed into a weighted least squares problem with respect to m in each outer iteration; In the formula, It is the model parameter vector for the k-th iteration; In the k-th iteration, the i-th model parameter component; β (k) Let be the dynamic regularization exponent parameter for the k-th iteration; In the k-th iteration, the weight coefficient corresponding to the i-th model parameter component. It is a diagonal weight matrix. is the dimension of the weight matrix; Step 4.2: Determine the optimal regularization parameter using an adaptive estimation method for the regularization parameter. Based on the candidate regularization parameter values ​​set in step 1 Calculate the current The corresponding regularization parameter adaptive selection criterion function CF value: (19)， in, The total number of single-channel data points participating in the inversion. For identity matrix, influence matrix In the k-th iteration, it takes the form: (20)， Step 4.3: Calculate the CF value corresponding to all candidate regularization parameters, and select the one with the smallest CF value. As the optimal regularization parameter At this point, the inversion objective function is locally approximated as a weighted L2 norm regularization problem, forming a new inversion objective function: (21)， in, Let be the inversion objective function for the model parameter vector at the k-th iteration.

8. The method for inverting seismic data with azimuth differences in inclined fractured rocks according to claim 7, characterized in that, diagonally... Yuan The expression is as follows: (22)， In the formula, is the scaling parameter of the i-th model parameter component based on parameter covariance estimation.

9. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 1, characterized in that, Step 5 includes the following sub-steps: Step 5.1: Initialize the parameters of the inner iteration. Let the initial model at the inner iteration number t=0 be... The initial residual is calculated using the following formula. : (23)， in, As an intermediate variable, This is the model parameter vector when the inner layer iteration count is 0 and the outer layer iteration count is k. Step 5.2: Perform preconditioning by extracting intermediate variables. The inverse of the diagonal matrix formed by the diagonal elements yields the precondition matrix. The residuals are preprocessed using the following formula: (24)， in, Let be the residual of the t-th inner iteration. It is a residual Representation in the preprocessing space; Step 5.3: Determine the search direction The update is performed, and the search direction for the t-th inner iteration is as follows: (25)， in, Let be the conjugate gradient update coefficients of the t-th inner iteration; Step 5.4: Update the model vector and residual vector using the following formula: (26)， (27)， in, Let be the step size, and its calculation formula is: .

10. The method for inverting azimuth difference seismic data in inclined fractured rocks according to claim 1, characterized in that, Step 6 specifically involves: conducting a double-layer iterative loop of outer and inner layers to update the model parameter vector. When the model parameter m of the k-th outer layer iteration satisfies the following formula, the output model parameters are the final inversion results of the crack normal weakness and tangential weakness: (28)， in, This is the preset convergence threshold.

Images