Prestack wide-azimuth seismic fracture prediction method based on maximum likelihood attribute constraint

Abstract

A pre-stack wide-azimuth seismic fracture prediction method based on maximum likelihood attribute constraint belongs to the technical field of pre-stack fracture prediction, and comprises the following steps: firstly, using post-stack maximum likelihood attributes to construct a fracture probability distribution model containing anisotropy information of an underground medium, the model being used for simulating large-scale underground fractures and faults; then, the post-stack maximum likelihood attribute is used as a fracture development prior information constraint to be added into a target function, so that the reliability and the stability of pre-stack fracture prediction are improved; and finally, adding Gaussian distribution and a smooth background model into a target function by using wide-azimuth seismic data by adopting a pre-stack anisotropy inversion method under a Bayesian framework so as to improve the rationality and stability of inversion. Compared with crack prediction carried out by directly utilizing wide-azimuth seismic data, the crack prediction technical effect in the invention has the advantages that multiple solution lines of pre-stack crack prediction are greatly reduced, and the accuracy of a crack prediction result is improved.

Description

Technical Field

[0001] The invention belongs to the technical field of pre-stack fracture prediction, and in particular relates to a maximum likelihood attribute constrained pre-stack wide-azimuth seismic fracture prediction method. Background Art

[0002] Fractures are crucial in oil and gas exploration and development because they serve as crucial storage and migration pathways within reservoirs. Fractures can form naturally within rocks or through geological processes such as tectonic movement, folding, and faulting. By analyzing fracture distribution, strike, density, width, and other characteristics, optimal drilling locations and production methods can be determined, improving the success rate and efficiency of oil and gas exploration and development. Effective fracture prediction from pre-stack wide-azimuth seismic data has become a key technical challenge in oil and gas exploration. By leveraging the azimuth amplitude variations in pre-stack wide-azimuth seismic data, direct prediction of fracture parameters can be achieved.

[0003] Typically, fractures can be characterized using a variety of methods, including geology, core drilling, well logging, and seismic data. Core drilling and well logging, while the most intuitive means of fracture identification, only indicate fracture development around the wellbore. Fracture prediction using seismic data, however, can predict fracture patterns across the entire region. There are relatively numerous methods for seismic fracture prediction, the most common of which are post-stack attributes such as coherence, curvature, and maximum likelihood. These attributes are typically effective for predicting large and medium-scale fractures where seismic events are significantly misaligned or deflected. They can provide insights into the macroscopic development patterns and stress states of the entire region, but they offer limited insights into fracture details. Pre-stack wide-azimuth seismic inversion, based on the variation of seismic amplitude with offset and azimuth, reveals how different fracture orientations and densities vary with amplitude and azimuth. This allows the inversion of fracture physical parameters from pre-stack wide-azimuth seismic data. This method can predict detailed fracture parameters based on actual wide-azimuth seismic data and is therefore widely used in both theoretical research and practical data interpretation.

[0004] Due to the strong heterogeneity and anisotropy of fractured reservoirs, when the signal-to-noise ratio of pre-stack wide-azimuth seismic data is low or the fracture development pattern in the study area is relatively complex, the reliability of pre-stack fracture prediction using the variation of amplitude with offset and azimuth is relatively poor. In order to better achieve fine characterization of fractures in low signal-to-noise ratio data or in areas with complex fracture development, the prior information of large-scale fractures based on geometric attributes can be integrated on the basis of wide-azimuth pre-stack inversion. Under the constraint of this prior information, the accuracy of fracture parameter inversion can be comprehensively improved. Therefore, the present invention carries out wide-azimuth pre-stack fracture prediction under the constraint of post-stack maximum likelihood geometric attributes to achieve high-precision pre-stack fracture detection in areas with low signal-to-noise ratio or complex fracture development. Summary of the Invention

[0005] In order to solve the above problems, the present invention proposes a maximum likelihood attribute constrained pre-stack wide-azimuth earthquake fracture prediction method, comprising the following steps:

[0006] Step A: Using post-stack seismic data as a driver, perform maximum likelihood attribute detection on the post-stack seismic data to determine the development pattern of post-stack large-scale faults;

[0007] Step A1: Based on the seismic data, the gradient structure tensor algorithm is used to obtain the apparent dip angle of the formation along the main survey line and the connecting lateral line;

[0008] Step A2: Under the constraints of the apparent dip angles of the main survey line and the connecting lateral line, seismic similarity detection is performed based on seismic multi-channel cross-correlation;

[0009] Step A3: Calculate the maximum likelihood attribute through seismic multi-channel similarity to clarify the development pattern of large-scale faults;

[0010] Step B: Using wide-azimuth seismic data to obtain elastic impedance at different orientations, and constructing a direct characterization relationship between elastic impedance structure and fracture parameters;

[0011] Step C: Under the constraints of the Bayesian framework, the maximum likelihood attributes in step A are used as prior information to establish the inversion target functional for prestack anisotropic fracture prediction, thereby realizing prestack fracture prediction under the maximum likelihood attribute prior constraints.

[0012] Furthermore, in step A1, based on the seismic data, the gradient structure tensor algorithm is used to obtain the apparent dip angle of the strata along the main survey line and the tie lateral line;

[0013] The gradient structure tensor algorithm is used to obtain the amplitude vector gradient g from the seismic data. It is characterized as the derivative of the seismic amplitude u in three directions and is represented as follows:

[0014]

[0015] Where g represents the amplitude vector gradient, x, y and z represent the three-dimensional spatial directions of the seismic data, u represents the seismic amplitude, g x 、g y and g z are the derivatives of the earthquake amplitude along the x, y and z directions respectively;

[0016] The GST matrix is ​​constructed using directional derivatives as follows:

[0017]

[0018] The gradient structure tensor matrix T is expressed as follows using eigenvectors and eigenvalues:

[0019]

[0020] Where λ1, λ2 and λ3 represent the three non-negative eigenvalues ​​of the matrix T, v1, v2 and v3 represent the three eigenvectors corresponding to the three eigenvalues ​​λ1, λ2 and λ3, and λ1 ≥ λ2 ≥ λ3 ≥ 0. The apparent dip angle of the formation along the main survey line and the connecting survey line is expressed as:

[0021]

[0022] Where p represents the apparent dip of the formation along the main survey line, q represents the apparent dip of the formation along the connecting survey line, and v1(x), v1(y) and v1(z) represent the components of the first eigenvector v1 along the x, y and z directions respectively.

[0023] Furthermore, in step A2, under the constraints of the apparent dip angles of the main survey line and the contact lateral line, seismic similarity detection is performed based on seismic multi-channel cross-correlation;

[0024] Under the constraint of the apparent dip angle of the formation, driven by seismic attributes, multi-channel cross-correlation is performed on earthquakes to obtain seismic similarity. To obtain the similarity of a central seismic trace, it is necessary to analyze the window near the trace. A rectangular analysis window is selected. Assuming that the analysis window contains a total of N traces, the similarity C of the trace is expressed as:

[0025]

[0026] The subscript n represents the nth seismic data in the analysis window, N represents the total number of channels in the analysis window, and x represents the total number of channels in the analysis window. n and y n represents the distance of the nth data from the center point in the x and y directions, p represents the apparent dip of the formation along the main survey line, q represents the apparent dip of the formation along the connecting survey line, u represents the seismic amplitude, and the superscript H represents the Hilbert transformation of the seismic amplitude;

[0027] By increasing the longitudinal time window analysis window of multiple similarities to avoid the minimum value of the denominator in formula (5), a more robust similarity algorithm with the longitudinal analysis time window is characterized as follows:

[0028]

[0029] Where K represents the number of sample points above and below the central seismic trace, and the total number of sample points in the analysis time window is 2K+1.

[0030] Furthermore, in step A3, the maximum likelihood attribute is calculated by seismic multi-channel similarity to clarify the development pattern of large-scale faults;

[0031] The similarity of seismic traces is obtained by using apparent dip angle and earthquake calculation to further enhance the difference between earthquake faults and non-faults. The likelihood attribute L is further introduced as follows:

[0032] L=1-C m (7)

[0033] Where C represents the multi-channel similarity of the analysis center seismic trace, and the superscript m adjusts the difference between the maximum and minimum values ​​of the similarity by using an exponent;

[0034] Obtain the likelihood attribute surface distributed along the cross section, use the fault inclination and dip to determine the fracture inclination and dip to scan, assuming that the fracture inclination κ∈[κ min ,κ max ],inclination The dip and inclination are denoted by Δκ and The scanning interval is scanned, and by comparing the likelihood attributes L under different inclinations and inclinations, the maximum likelihood attribute L is selected. max That is the maximum likelihood attribute, expressed as:

[0035]

[0036] Among them L max is the maximum likelihood attribute, κ is the fracture tendency, is the fracture inclination angle.

[0037] Furthermore, in step B, the elastic impedance at different orientations is obtained using wide-azimuth seismic data, and a direct characterization relationship between the elastic impedance structure and the fracture parameters is constructed;

[0038] Under the action of ground stress and overlying stratum pressure, underground rocks exhibit directional cracks. The vertical cracks with directional arrangement are approximately equivalent to a transversely isotropic HTI medium with a horizontal symmetry axis. Based on the anisotropy assumption, the reflection coefficient of the double-layer HTI medium is as follows:

[0039]

[0040] Among them, ΔI P Indicates the difference in longitudinal wave impedance between the upper and lower layers of media. represents the mean longitudinal wave impedance of the upper and lower layers of media, Δα represents the difference in longitudinal wave velocity between the upper and lower layers of media, represents the mean longitudinal wave velocity of the upper and lower layers of the medium, Δβ represents the difference in shear wave velocity between the upper and lower layers of the medium, represents the mean shear wave velocity of the upper and lower layers, θ and φ represent the incident angle and azimuth respectively, and ΔG represents the difference in shear modulus between the upper and lower layers. represents the mean shear modulus of the upper and lower layers, ε (v) , δ(v) and γ (v) is the Thomsen anisotropy parameter of the HTI medium, Δε (v) , Δδ (v) and Δγ (v) Indicates the difference in Thomsen anisotropy parameters between the upper and lower media;

[0041] According to Schoenberg linear sliding theory and Hudson theory, the Thomsen anisotropy parameter ε (v) , δ (v) and γ (v) and the crack normal weakness Δ N and tangential weakness Δ T The connections are as follows:

[0042]

[0043] Where g represents the square of the inverse mean of the ratio of the longitudinal and transverse wave velocities of the upper and lower layers. Combining Equation (9) with Equation (10), the reflection coefficient of the HTI medium is obtained as the normal weakness Δ N and tangential weakness Δ T Representational form:

[0044]

[0045] Among them, ΔI P Indicates the difference in longitudinal wave impedance between the upper and lower layers of media. Indicates the average longitudinal wave impedance of the upper and lower layers, ΔI S Indicates the difference in shear wave impedance between the upper and lower layers of the medium. represents the mean shear wave impedance of the upper and lower layers, θ and φ represent the incident angle and azimuth respectively, g represents the square of the inverse mean of the ratio of the longitudinal and shear wave velocities of the upper and lower layers, Δ N and Δ T They represent the normal weakness and tangential weakness of the crack respectively.

[0046] Furthermore, according to Schoenberg linear sliding theory, the crack density is characterized as:

[0047]

[0048] Following Connolly's idea of ​​elastic impedance, the elastic impedance is expressed as a reflection coefficient as follows:

[0049]

[0050] Substituting equation (13) into equation (11), the reflection coefficient is expressed in the form of elastic impedance:

[0051]

[0052] Using Equation (14), we can obtain a direct relationship between elastic impedance and crack parameters. By solving the elastic impedance at different azimuth angles, we can obtain a direct representation of the crack parameters. The inversion problem can then be expressed as:

[0053] d=Gm (15)

[0054] Where d represents the observation parameter, i.e., the elastic impedance at different azimuths and incident angles, G represents the coefficient matrix, and m represents the fracture parameters and elastic parameters to be inverted.

[0055] Furthermore, in step C, under the constraints of the Bayesian framework, the maximum likelihood attribute in step A is used as prior information to establish an inversion target functional for prestack anisotropic fracture prediction, thereby achieving prestack fracture prediction under the maximum likelihood attribute prior constraint;

[0056] Bayesian theory uses the Bayesian formula, combined with the overall observation sample and prior information, to obtain the probability density of model parameters under incomplete statistics, which is specifically expressed as:

[0057]

[0058] Where m represents the fracture parameters and elastic parameters to be inverted, d represents the observed parameters, p(md) represents the posterior probability density of the parameters to be inverted, p(dm) represents the likelihood function of the observed data, p(d) represents the probability density function of the known observed data, and p(m) represents the prior information of the parameters to be inverted.

[0059] Furthermore, the earthquake maximum likelihood attribute L obtained in step A is max The prior information of the parameters to be inverted is added as a constraint to the inversion objective function, and the inversion objective function is expressed as:

[0060]

[0061] Where m represents the fracture parameters and elastic parameters to be inverted, d represents the observed parameters, p(md) represents the posterior probability density of the parameters to be inverted, η1 and η2 represent the inversion weight coefficients, and M l The initial model of the parameters to be inverted is represented by the inversion objective functional of fracture prediction under the maximum likelihood attribute constraint:

[0062]

[0063] The inversion objective functional in Equation (18) is solved to finally realize the pre-stack wide azimuthal fracture prediction under the maximum likelihood attribute constraint.

[0064] The beneficial effects of the present invention are as follows: the present invention solves the problem of instability in fracture prediction inversion. First, the post-stack maximum likelihood attribute is used to construct a fracture probability distribution model containing anisotropic information of the underground medium. The model can be used to simulate large-scale underground fractures and faults. Then, the post-stack maximum likelihood attribute is added to the objective function as a priori information constraint on fracture development to improve the reliability and stability of pre-stack fracture prediction. Finally, the pre-stack anisotropic inversion method under the Bayesian framework is used to utilize wide-azimuth seismic data, and the Gaussian distribution and smooth background model are added to the objective function to improve the rationality and stability of the inversion. Compared with directly using wide-azimuth seismic data for fracture prediction, the fracture prediction technology in this invention greatly reduces the multi-solution of pre-stack fracture prediction and improves the accuracy of fracture prediction results. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 This is the flow chart of the pre-stack wide-azimuth earthquake fracture prediction technology under maximum likelihood attribute constraints;

[0066] Figure 2 This is the post-stack seismic profile through Well A;

[0067] Figure 3 is the maximum likelihood attribute profile through well A;

[0068] Figure 4 The superimposed cross-section diagrams of the angle section at different azimuths and incident angles through Well A;

[0069] Figure 5 This is a fracture density profile predicted by conventional pre-stack anisotropic inversion of the present invention;

[0070] Figure 6 Predict fracture density profiles for prestack anisotropy inversion under maximum likelihood attribute constraints. DETAILED DESCRIPTION

[0071] To make the technical means and objectives of the present invention easier to understand, the present invention is further described below in conjunction with specific embodiments. A maximum likelihood attribute-constrained pre-stack wide-azimuth seismic fracture prediction method mainly includes the following steps:

[0072] Step A: Using post-stack seismic data as a driver, perform maximum likelihood attribute detection on the post-stack seismic data to determine the development pattern of post-stack large-scale faults;

[0073] Cracks in underground rocks serve as the main migration channels and storage spaces for oil and gas. Under the action of underground stress, the development of cracks usually has a high correlation with the main faults. Seismic coherence is a traditional method for detecting earthquake faults. It mainly uses the similarity of adjacent seismic traces in seismic data to characterize the discontinuities in the strata. However, when the seismic response characteristics are weak or the seismic signal-to-noise ratio is low, the discontinuity of the earthquake is also relatively strong. At this time, the coherence attribute has a relatively weak ability to characterize the fault. Even if the fault really exists, when the fault distance is less than the wavelength of the seismic wave, the seismic phase axis is relatively continuous. Therefore, simply using the continuity of the seismic phase axis cannot efficiently identify the fault. Therefore, in order to solve the problem that continuity cannot efficiently identify faults, the present invention proposes a new maximum likelihood attribute. The main technical ideas are as follows:

[0074] Step A1: Based on the seismic data, the gradient structure tensor (GST) algorithm is used to obtain the apparent dip angle of the seismic formation along the main survey line and the connecting lateral line;

[0075] The Gradient Structure Tensor (GST) algorithm is used to obtain the amplitude vector gradient g from seismic data. It can be represented as the derivative of the seismic amplitude u in three directions. The main characteristics are as follows:

[0076]

[0077] Where g represents the amplitude vector gradient, x, y and z represent the three-dimensional spatial directions of the seismic data, u represents the seismic amplitude, g x 、g y and g z are the derivatives of the earthquake amplitude along the x, y and z directions respectively.

[0078] The GST matrix is ​​constructed using directional derivatives as follows:

[0079]

[0080] The gradient structure tensor matrix T is expressed as follows using eigenvectors and eigenvalues:

[0081]

[0082] Where λ1, λ2 and λ3 represent the three non-negative eigenvalues ​​of the matrix T, v1, v2 and v3 represent the three eigenvectors corresponding to the three eigenvalues ​​λ1, λ2 and λ3, and λ1 ≥ λ2 ≥ λ3 ≥ 0. The apparent dip of the formation along the main survey line and the connecting survey line can be expressed as:

[0083]

[0084] Where p represents the apparent dip of the formation along the main survey line, q represents the apparent dip of the formation along the connecting survey line, and v1(x), v1(y) and v1(z) represent the components of the first eigenvector v1 along the x, y and z directions respectively.

[0085] Step A2: Under the constraints of the apparent dip angles of the main survey line and the connecting lateral line, seismic similarity detection is performed based on seismic multi-channel cross-correlation;

[0086] Under the constraint of the apparent dip angle of the formation, driven by seismic attributes, multi-channel cross-correlation is performed on earthquakes to obtain seismic similarity. To obtain the similarity of a central seismic trace, it is necessary to analyze the window near the trace. Usually, a rectangular analysis window is selected. Assuming that the analysis window contains a total of N traces, the similarity C of the trace can be expressed as:

[0087]

[0088] The subscript n represents the nth seismic data in the analysis window, N represents the total number of channels in the analysis window, and x represents the total number of channels in the analysis window. n and y n They represent the distance of the nth data from the center point in the x and y directions, p represents the apparent dip of the stratum along the main survey line, q represents the apparent dip of the stratum along the connecting survey line, u represents the seismic amplitude, and the superscript H represents the Hilbert transform of the seismic amplitude. This can effectively avoid the problem of gradient instability caused by slow changes in seismic amplitude when the earthquake is at a strong peak or trough.

[0089] In order to improve the stability of the similarity algorithm, the longitudinal time window analysis window of multiple similarities can be added to avoid the minimum value of the denominator in formula (5). Therefore, a more robust similarity algorithm with the introduction of the longitudinal analysis window can be characterized as follows:

[0090]

[0091] K represents the number of sample points above and below the central seismic trace, so the total number of sample points in the analysis time window is 2K+1.

[0092] Step A3: Calculate the maximum likelihood attribute through seismic multi-channel similarity to clarify the development pattern of large-scale faults;

[0093] The similarity of seismic traces is obtained by using apparent dip angle and earthquake calculation. In order to further enhance the difference between earthquake faults and non-faults, the likelihood attribute L is further introduced as follows:

[0094] L=1-C m (7)

[0095] Where C represents the multi-channel similarity of the analysis center seismic trace, and the superscript m adjusts the difference between the maximum and minimum similarity values ​​by using an exponent. The larger the exponent, the greater the difference between faults and non-faults.

[0096] In order to obtain the likelihood attribute surface distributed along the cross section, the fault dip and inclination are usually used to determine it. However, the dip and inclination of the fault surface are unknown, so the dip and inclination of the fault need to be scanned. Assume that the dip of the fault κ∈[κ min ,κ max ],inclination The dip and inclination are denoted by Δκ and The scanning interval is scanned, and by comparing the likelihood attributes L under different inclinations and inclinations, the maximum likelihood attribute L is selected. max That is the maximum likelihood attribute, which can be expressed as:

[0097]

[0098] Among them L max is the maximum likelihood attribute, κ is the fracture tendency, is the fracture inclination angle.

[0099] Step B: Using wide-azimuth seismic data to obtain elastic impedance at different orientations, and constructing a direct characterization relationship between elastic impedance structure and fracture parameters;

[0100] Under the influence of geostress and overlying strata, underground rocks typically exhibit directional fractures. We equate these directional vertical fractures to a transversely isotropic (HTI) medium with a horizontal axis of symmetry. Based on the anisotropy assumption, Ruger (1996) derives the reflection coefficient of a two-layer HTI medium as follows:

[0101]

[0102] Among them, ΔI P Indicates the difference in longitudinal wave impedance between the upper and lower layers of media. represents the mean longitudinal wave impedance of the upper and lower layers of media, Δα represents the difference in longitudinal wave velocity between the upper and lower layers of media, represents the mean longitudinal wave velocity of the upper and lower layers of the medium, Δβ represents the difference in shear wave velocity between the upper and lower layers of the medium, represents the mean shear wave velocity of the upper and lower layers, θ and φ represent the incident angle and azimuth respectively, and ΔG represents the difference in shear modulus between the upper and lower layers. represents the mean shear modulus of the upper and lower layers, ε (v) , δ (v) and γ (v) is the Thomsen anisotropy parameter of the HTI medium, Δε (v) , Δδ(v) and Δγ (v) Indicates the difference in Thomsen anisotropy parameters between the upper and lower media.

[0103] According to Schoenberg linear sliding theory and Hudson theory, the Thomsen anisotropy parameter ε can be (v) , δ (v) and γ (v) and the crack normal weakness Δ N and tangential weakness Δ T The connections are as follows:

[0104]

[0105] Where g represents the square of the inverse mean of the ratio of the longitudinal and transverse wave velocities of the upper and lower layers. Combining Equation (9) with Equation (10), we can obtain the reflection coefficient of the HTI medium in terms of the normal weakness Δ N and tangential weakness Δ T Representational form:

[0106]

[0107] Among them, ΔI P Indicates the difference in longitudinal wave impedance between the upper and lower layers of media. Indicates the average longitudinal wave impedance of the upper and lower layers, ΔI S Indicates the difference in shear wave impedance between the upper and lower layers of the medium. represents the mean shear wave impedance of the upper and lower layers, θ and φ represent the incident angle and azimuth respectively, g represents the square of the inverse mean of the ratio of the longitudinal and shear wave velocities of the upper and lower layers, Δ N and Δ T They represent the normal weakness and tangential weakness of the crack respectively.

[0108] According to Schoenberg linear sliding theory, the crack density can be characterized as:

[0109]

[0110] Therefore, it is only necessary to change the normal weakness of the crack Δ N and tangential weakness Δ T By finding this, we can calculate the fracture density in the formation. Following Connolly's elastic impedance concept, the elastic impedance can be expressed as a reflection coefficient as follows:

[0111]

[0112] Substituting equation (13) into equation (11), the reflection coefficient is expressed in the form of elastic impedance:

[0113]

[0114] Equation (14) can be used to obtain a direct relationship between elastic impedance and crack parameters. Therefore, by solving the elastic impedance at different azimuth angles, a direct representation of the crack parameters can be obtained, and the inversion problem can be expressed as:

[0115] d=Gm (15)

[0116] Where d represents the observation parameter, i.e., the elastic impedance at different azimuths and incident angles, G represents the coefficient matrix, and m represents the fracture parameters and elastic parameters to be inverted.

[0117] Step C: Under the constraints of the Bayesian framework, the maximum likelihood attribute in step A is used as prior information to establish the inversion target functional for prestack anisotropic fracture prediction, thus realizing prestack fracture prediction under the maximum likelihood attribute prior constraint.

[0118] Bayesian theory uses the Bayesian formula, combined with the overall observation sample and prior information, to obtain the probability density of model parameters under incomplete statistics, which is specifically expressed as:

[0119]

[0120] Where m represents the fracture parameters and elastic parameters to be inverted, d represents the observed parameters, p(md) represents the posterior probability density of the parameters to be inverted, p(dm) represents the likelihood function of the observed data, p(d) represents the probability density function of the known observed data, and p(m) represents the prior information of the parameters to be inverted.

[0121] Therefore, the earthquake maximum likelihood attribute L obtained in step A is max The prior information of the parameters to be inverted is added as a constraint to the inversion objective function, and the inversion objective function can be expressed as:

[0122]

[0123] Where m represents the fracture parameters and elastic parameters to be inverted, d represents the observed parameters, p(md) represents the posterior probability density of the parameters to be inverted, η1 and η2 represent the inversion weight coefficients, and M l represents the initial model of the parameters to be inverted. Therefore, the inversion objective functional for fracture prediction under the maximum likelihood attribute constraint can be expressed as:

[0124]

[0125] The inversion objective functional in equation (18) is solved to ultimately achieve pre-stack wide-azimuth fracture prediction under the maximum likelihood attribute constraint. Compared with directly using wide-azimuth seismic data for fracture prediction, the fracture prediction technology in this invention greatly reduces the multi-solution problem of pre-stack fracture prediction and improves the accuracy of fracture prediction results.

[0126] The above invention technical process (attached Figure 1 ), applied to a certain eastern actual work area, firstly the post-stack seismic data (attached Figure 2 ) to obtain the maximum likelihood attribute (attached Figure 3 ), and then the direct relationship between the crack parameters and elastic impedance is obtained by deduction, and then the maximum likelihood attribute (attached) is combined with the Bayesian framework. Figure 3 ) as a priori information, added to the target functional of prestack anisotropic inversion, and the variation of azimuth and incident angle with offset distance (Appendix Figure 4 ), pre-stack anisotropy inversion was performed to obtain parameters such as fracture density. By comparing the conventional pre-stack anisotropy inversion results (Appendix Figure 5 ) and the prestack anisotropy inversion results using maximum likelihood attribute constraints (Appendix Figure 6 ) shows that the fracture density inversion results under the maximum likelihood attribute constraint are more consistent with the well measured results. The technology in the patent of this invention can effectively reduce the multi-solution of pre-stack anisotropic inversion and improve the accuracy of pre-stack fracture prediction inversion.

[0127] To address the problem of fracture prediction in areas with low signal-to-noise ratios or complex fracture development, this paper uses prestack wide-azimuth seismic data as a basis and utilizes the variation of amplitude with azimuth and offset to perform prestack anisotropic inversion. Prestack anisotropic inversion not only considers elastic parameters such as P- and S-wave velocities and density, but also anisotropic parameters caused by fractures. Different elastic and anisotropic parameters contribute significantly to the reflection coefficient, resulting in a high degree of ambiguity in prestack fracture prediction results.

[0128] To avoid this ambiguity, the present invention conducts fracture detection based on post-stack maximum likelihood geometric attributes based on post-stack seismic data, and then conducts wide-azimuth pre-stack fracture prediction under the constraints of post-stack attributes, thereby improving the accuracy of fracture detection and reducing the ambiguity of fracture prediction results.

[0129] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed in the present invention, who makes equivalent replacements or changes based on the technical solutions and concepts of the present invention, should be covered by the scope of protection of the present invention.

Claims

1. A maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method, characterized in that: The following steps are involved: Step A: Using post-stack seismic data as a driver, perform maximum likelihood attribute detection on the post-stack seismic data to determine the development pattern of post-stack large-scale faults; Step A1: Based on the seismic data, the gradient structure tensor algorithm is used to obtain the apparent dip angle of the formation along the main survey line and the connecting lateral line; Step A2: Under the constraints of the apparent dip angles of the main survey line and the connecting lateral line, seismic similarity detection is performed based on seismic multi-channel cross-correlation; Step A3: Calculate the maximum likelihood attribute through seismic multi-channel similarity to clarify the development pattern of large-scale faults; Step B: Using wide-azimuth seismic data to obtain elastic impedance at different orientations, and constructing a direct characterization relationship between elastic impedance structure and fracture parameters; Step C: Under the constraints of the Bayesian framework, the maximum likelihood attributes in step A are used as prior information to establish the inversion target functional for prestack anisotropic fracture prediction, thereby realizing prestack fracture prediction under the maximum likelihood attribute prior constraints.

2. The maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method according to claim 1, characterized in that: In the step A1, based on the seismic data, the gradient structure tensor algorithm is used to obtain the apparent dip angle of the strata along the main survey line and the connecting lateral line; The gradient structure tensor algorithm is used to obtain the amplitude vector gradient g from the seismic data. It is characterized as the derivative of the seismic amplitude u in three directions and is represented as follows: Where g represents the amplitude vector gradient, x, y and z represent the three-dimensional spatial directions of the seismic data, u represents the seismic amplitude, g x 、g y and g z are the derivatives of the earthquake amplitude along the x, y and z directions respectively; The GST matrix is ​​constructed using directional derivatives as follows: The gradient structure tensor matrix T is expressed as follows using eigenvectors and eigenvalues: Where λ1, λ2 and λ3 represent the three non-negative eigenvalues ​​of the matrix T, v1, v2 and v3 represent the three eigenvectors corresponding to the three eigenvalues ​​of λ1, λ2 and λ3, and λ1 ≥ λ2 ≥ λ3 ≥ 0, then the apparent dip angle of the formation along the main survey line and the connecting survey line is expressed as: Where p represents the apparent dip of the formation along the main survey line, q represents the apparent dip of the formation along the connecting survey line, and v1(x), v1(y) and v1(z) represent the components of the first eigenvector v1 along the x, y and z directions respectively.

3. The maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method according to claim 1, characterized in that: In step A2, under the constraints of the apparent dip angles of the main survey line and the contact lateral line, seismic similarity detection is performed based on seismic multi-channel cross-correlation; Under the constraint of the apparent dip angle of the formation, driven by seismic attributes, multi-channel cross-correlation is performed on earthquakes to obtain seismic similarity. To obtain the similarity of a central seismic trace, it is necessary to analyze the window near the trace. A rectangular analysis window is selected. Assuming that the analysis window contains a total of N traces, the similarity C of the trace is expressed as: The subscript n represents the nth seismic data in the analysis window, N represents the total number of channels in the analysis window, and x represents the total number of channels in the analysis window. n and y n represents the distance of the nth data from the center point in the x and y directions, p represents the apparent dip of the formation along the main survey line, q represents the apparent dip of the formation along the connecting survey line, u represents the seismic amplitude, and the superscript H represents the Hilbert transformation of the seismic amplitude; By increasing the longitudinal time window analysis window of multiple similarities to avoid the minimum value of the denominator in formula (5), a more robust similarity algorithm with the longitudinal analysis time window is characterized as follows: Where K represents the number of sample points above and below the central seismic trace, and the total number of sample points in the analysis time window is 2K+1.

4. The maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method according to claim 1, characterized in that: In step A3, the maximum likelihood attribute is calculated by seismic multi-channel similarity to clarify the development pattern of large-scale faults; The similarity of seismic traces is obtained by using apparent dip angle and earthquake calculation to further enhance the difference between earthquake faults and non-faults. The likelihood attribute L is further introduced as follows: L=1-C m (7) Where C represents the multi-channel similarity of the analysis center seismic trace, and the superscript m adjusts the difference between the maximum and minimum values ​​of the similarity by using an exponent; Obtain the likelihood attribute surface distributed along the cross section, use the fault inclination and dip to determine the fracture inclination and dip to scan, assuming that the fracture inclination κ∈[κ min ,κ max ], the inclination angle θ∈[θ min ,θ max ], the inclination and dip are scanned at the scanning intervals of Δκ and Δθ respectively, and the maximum likelihood attribute L is selected by comparing the likelihood attributes L under different inclinations and dips. max That is the maximum likelihood attribute, expressed as: L max =max[1-C(κ,θ) m ] (8) Among them L max is the maximum likelihood attribute, κ is the fault tendency, and θ is the fault dip angle.

5. The maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method according to claim 1, characterized in that: In step B, the elastic impedance at different orientations is obtained using wide-azimuth seismic data, and a direct characterization relationship between the elastic impedance structure and the fracture parameters is constructed; Under the action of ground stress and overlying stratum pressure, underground rocks exhibit directional cracks. The vertical cracks with directional arrangement are approximately equivalent to a transversely isotropic HTI medium with a horizontal symmetry axis. Based on the anisotropy assumption, the reflection coefficient of the double-layer HTI medium is as follows: Among them, ΔI P Indicates the difference in longitudinal wave impedance between the upper and lower layers of media. represents the mean longitudinal wave impedance of the upper and lower layers of media, Δα represents the difference in longitudinal wave velocity between the upper and lower layers of media, represents the mean longitudinal wave velocity of the upper and lower layers of the medium, Δβ represents the difference in shear wave velocity between the upper and lower layers of the medium, represents the mean shear wave velocity of the upper and lower layers, θ and φ represent the incident angle and azimuth respectively, and ΔG represents the difference in shear modulus between the upper and lower layers. represents the mean shear modulus of the upper and lower layers, ε (v) , δ (v) and γ (v) is the Thomsen anisotropy parameter of the HTI medium, Δε (v) , Δδ (v) and Δγ (v) Indicates the difference in Thomsen anisotropy parameters between the upper and lower media; According to Schoenberg linear sliding theory and Hudson theory, the Thomsen anisotropy parameter ε (v) , δ (v) and γ (v) and the crack normal weakness Δ N and tangential weakness Δ T The connections are as follows: Where g represents the square of the inverse mean of the ratio of the longitudinal and transverse wave velocities of the upper and lower layers. Combining Equation (9) with Equation (10), the reflection coefficient of the HTI medium is obtained as the normal weakness Δ N and tangential weakness Δ T Representational form: Among them, ΔI P Indicates the difference in longitudinal wave impedance between the upper and lower layers of media. Indicates the average longitudinal wave impedance of the upper and lower layers, ΔI S Indicates the difference in shear wave impedance between the upper and lower layers of the medium. represents the mean shear wave impedance of the upper and lower layers, θ and φ represent the incident angle and azimuth respectively, g represents the square of the inverse mean of the ratio of the longitudinal and shear wave velocities of the upper and lower layers, Δ N and Δ T They represent the normal weakness and tangential weakness of the crack respectively.

6. The maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method according to claim 5, characterized in that: According to Schoenberg linear sliding theory, the crack density is characterized as: Following Connolly's idea of ​​elastic impedance, the elastic impedance is expressed as a reflection coefficient as follows: Substituting equation (13) into equation (11), the reflection coefficient is expressed in the form of elastic impedance: Using Equation (14), we can obtain a direct relationship between elastic impedance and crack parameters. By solving the elastic impedance at different azimuth angles, we can obtain a direct representation of the crack parameters. The inversion problem can then be expressed as: d=Gm (15) Where d represents the observation parameter, i.e., the elastic impedance at different azimuths and incident angles, G represents the coefficient matrix, and m represents the fracture parameters and elastic parameters to be inverted.

7. The maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method according to claim 6, characterized in that: In step C, under the constraints of the Bayesian framework, the maximum likelihood attribute in step A is used as prior information to establish an inversion target functional for prestack anisotropic fracture prediction, thereby realizing prestack fracture prediction under the maximum likelihood attribute prior constraint; Bayesian theory uses the Bayesian formula, combined with the overall observation sample and prior information, to obtain the probability density of model parameters under incomplete statistics, which is specifically expressed as: Where m represents the fracture parameters and elastic parameters to be inverted, d represents the observed parameters, p(m|d) represents the posterior probability density of the parameters to be inverted, p(d|m) represents the likelihood function of the observed data, p(d) represents the probability density function of the known observed data, and p(m) represents the prior information of the parameters to be inverted.

8. The maximum likelihood attribute-constrained pre-stack wide-azimuth earthquake fracture prediction method according to claim 7, characterized in that: The earthquake maximum likelihood attribute L obtained in step A is max The prior information of the parameters to be inverted is added as a constraint to the inversion objective function, and the inversion objective function is expressed as: Where m represents the fracture parameters and elastic parameters to be inverted, d represents the observed parameters, p(m|d) represents the posterior probability density of the parameters to be inverted, η1 and η2 represent the inversion weight coefficients, and M l The initial model of the parameters to be inverted is represented by the inversion objective functional of fracture prediction under the maximum likelihood attribute constraint: The inversion objective functional in Equation (18) is solved to finally realize the pre-stack wide azimuthal fracture prediction under the maximum likelihood attribute constraint.

Images