A multi-wave joint AVO inversion method, device and electronic equipment for fractured reservoirs

Abstract

The present invention provides a multi-wave joint AVO inversion method, device and electronic equipment for fractured reservoirs, which can compress VSP seismic observation records corresponding to PS waves, PS1 waves and PS2 waves of fractured reservoirs into the time domain of the first arrival of zero-bias VSP longitudinal waves for AVA gather conversion, thereby obtaining AVA seismic observation records with time matching of each wave without losing fracture information. Afterwards, based on the inversion equation corresponding to the initial model for deriving fracture parameters, the anisotropic parameters and isotropic parameters of the fractured reservoir are inverted for imaging. The inversion equation simulates the AVA seismic synthetic record according to the precise reflection coefficient based on the initial model, and calculates the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model. This gradient calculation method effectively solves the problem that there is no linear relationship between the precise reflection coefficient and the objective function, and can improve the inversion performance based on the precise reflection coefficient.

Description

Technical Field

[0001] This document relates to the field of seismic exploration technology, and in particular to a multi-wave combined AVO inversion method, device and electronic equipment for fractured reservoirs. Background Art

[0002] Due to tectonic movements such as folding and squeezing, deep strata uplift, etc., vertical fractures can develop in the strata. Such fractures can not only serve as storage space for the reservoir, but also improve the permeability of the reservoir, thereby increasing the economic recoverability of the reservoir. In geophysics, the medium in which a group of vertical fractures develop is approximately equivalent to a transversely isotropic (HTI) medium with a horizontal symmetry axis, which becomes the basic model used for fracture prediction in various reservoirs such as tight sandstone, carbonate rock, shale, and hydrate.

[0003] With the development of multi-component seismic exploration, the joint inversion method of multiple reflection waves has developed rapidly, greatly improving the accuracy of reservoir prediction. In HTI media, shear waves can be divided into fast shear waves (S1 waves) and slow shear waves (S2 waves). Since the polarization direction of the S1 wave is parallel to the isotropic plane, the reflection coefficient of the converted fast shear wave (PS1 wave) is only affected by the isotropic background parameters. However, the polarization of the converted slow shear wave (PS2 wave) is perpendicular to the fracture surface and contains information on isotropic and anisotropic parameters. Therefore, in the amplitude variation with incident angle (AVA) response, split shear waves can provide more inversion parameter information. In addition, the time delay between the S1 wave and the S2 wave caused by the difference in propagation speed also characterizes the degree of fracture development.

[0004] At present, due to the relative difficulty in matching multi-wave data, there is no joint inversion application for PP waves, PS1 waves and PS2 waves. Summary of the invention

[0005] The embodiment of the present invention aims to provide a multi-wave joint AVO inversion method, device and electronic equipment based on fractured reservoirs, which can perform AVO inversion on fractured VSP reservoir data combining PP waves, PS1 waves and PS2 waves.

[0006] In order to achieve the above purpose, the embodiments of this specification are implemented as follows:

[0007] In a first aspect, a multi-wave joint AVO inversion method for a fractured reservoir is provided, comprising:

[0008] The VSP seismic observation records corresponding to PS waves, PS1 waves and PS2 waves of fractured reservoirs are compressed into the time domain of the first arrival of zero-bias VSP longitudinal waves to perform AVA gather conversion and obtain AVA seismic observation records.

[0009] Constructing an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropy parameters in AVA gather data, and the inversion objective function is a Bayesian theory function;

[0010] Constructing an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate AVA seismic synthetic records based on the initial model according to the precise reflection coefficient, and calculating the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model;

[0011] Substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir;

[0012] Based on the anisotropic parameters of the fractured reservoir, the fracture state of the fractured reservoir is determined.

[0013] In a second aspect, a multi-wave joint AVO inversion device for fractured reservoirs is provided, comprising:

[0014] The multi-wave matching module is used to compress the VSP seismic observation records corresponding to the PS wave, PS1 wave and PS2 wave of the fractured reservoir into the time domain of the first arrival of the zero-bias VSP longitudinal wave to perform AVA gather conversion to obtain the AVA seismic observation records;

[0015] A model function construction module, used to construct an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropy parameters in AVA gather data, and the inversion objective function is a Bayesian theory function;

[0016] An inversion equation construction model is used to construct an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate AVA seismic synthetic records based on the initial model according to the precise reflection coefficient, and calculate the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model;

[0017] An inversion execution module, used for substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir;

[0018] The fracture state determination module is used to determine the fracture state of the fracture reservoir based on the anisotropic parameters of the fracture reservoir.

[0019] In a third aspect, an electronic device is provided, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the computer program is executed by the processor:

[0020] The VSP seismic observation records corresponding to PS waves, PS1 waves and PS2 waves of fractured reservoirs are compressed into the time domain of the first arrival of zero-bias VSP longitudinal waves to perform AVA gather conversion and obtain AVA seismic observation records.

[0021] Constructing an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropy parameters in AVA gather data, and the inversion objective function is a Bayesian theory function;

[0022] Constructing an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate AVA seismic synthetic records based on the initial model according to the precise reflection coefficient, and calculating the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model;

[0023] Substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir;

[0024] Based on the anisotropic parameters of the fractured reservoir, the fracture state of the fractured reservoir is determined.

[0025] In a fourth aspect, a computer readable storage medium is provided, wherein a computer program is stored on the computer readable storage medium, and when the computer program is executed by a processor, the following steps are implemented:

[0026] The VSP seismic observation records corresponding to PS waves, PS1 waves and PS2 waves of fractured reservoirs are compressed into the time domain of the first arrival of zero-bias VSP longitudinal waves to perform AVA gather conversion and obtain AVA seismic observation records.

[0027] Constructing an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropy parameters in AVA gather data, and the inversion objective function is a Bayesian theory function;

[0028] Constructing an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate AVA seismic synthetic records based on the initial model according to the precise reflection coefficient, and calculating the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model;

[0029] Substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir;

[0030] Based on the anisotropic parameters of the fractured reservoir, the fracture state of the fractured reservoir is determined.

[0031] The present invention can compress the VSP seismic observation records of PS waves, PS1 waves and PS2 waves corresponding to the fractured reservoir into the time domain of the first arrival of the zero-bias VSP longitudinal wave for AVA gather conversion, thereby obtaining the AVA seismic observation records of the PS waves, PS1 waves and PS2 waves with time matching without losing the fracture information; then, construct an initial model for deriving fracture parameters, and invert the anisotropic parameters of the fractured reservoir based on the inversion equation corresponding to the initial model. Among them, the inversion equation can simulate the AVA seismic synthetic record according to the precise reflection coefficient based on the initial model, and calculate the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model. This gradient calculation method can effectively solve the problem that there is no linear relationship between the precise reflection coefficient and the objective function, and can improve the inversion performance based on the precise reflection coefficient, and finally obtain more accurate anisotropic parameters for imaging the structure of the fractured reservoir. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in the embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0033] Figure 1 A schematic diagram of a first flow chart of a multi-wave joint inversion method for anisotropic parameters of fractured reservoirs provided in an embodiment of the present invention.

[0034] Figure 2 A second flow chart of the multi-wave joint inversion method for anisotropic parameters of fractured reservoirs provided in an embodiment of the present invention.

[0035] Figure 3Schematic diagram of the P-wave first arrival of zero-bias VSP data (a) and the P-wave and S-wave layer velocities obtained from the zero-bias VSP data (b).

[0036] Figure 4 It is the shot record of Z, R and T components after preprocessing.

[0037] Figure 5 It is the shot record of separated Z, R and T components.

[0038] Figure 6 To match to t PFB Shot records of the Z, R, and T components of the domain.

[0039] Figure 7 To match to t0 PFB Shot records of the Z, R, and T components of the domain.

[0040] Figure 8 Fracture orientation calculated for Alford rotation Schematic diagram of seismic records of fast and slow shear wave delays τ.

[0041] Fig. 9 To match to t0 PFB Schematic diagram of seismic records of PP waves, PS1 waves and PS2 waves in the domain.

[0042] Fig.10 Schematic diagram of the seismic records of the PP wave and PS wave incident angle gathers (a, b), reflection point X coordinate gathers (c, d), and reflection point Z coordinate gathers (e) matched to the t0PFB domain.

[0043] Fig.11 Schematic diagram of the initial model for the inversion parameters λ, μ, and ρ.

[0044] Fig.12 The crack direction beside the well Schematic diagram of the time difference τ between the fast full wave and the fast full wave.

[0045] Fig.13 Schematic diagram of the fracture density e near the well.

[0046] Fig.14 Schematic diagram of the depth-domain AVA gathers of PP, PS1, and PS2 waves generated for the actual data (a, c, e) and the initial model (b, d, f).

[0047] Fig.15 Schematic diagram of the depth domain superposition recording of PP, PS1 and PS2 waves.

[0048] Fig.16 Schematic diagram of the angle-independent wavelets of PP wave, PS1 wave, and PS2 wave in the depth domain.

[0049] Fig.17 Schematic diagram of the joint inversion results.

[0050] Fig.18 A schematic diagram of the structure of a multi-wave joint inversion device for anisotropic parameters of fractured reservoirs provided in an embodiment of the present invention.

[0051] Fig.19 A schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0052] In order to enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of this specification, not all of the embodiments. Based on the embodiments in this specification, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of this specification.

[0053] As mentioned above, multi-component seismic exploration can greatly improve the accuracy of reservoir prediction and is an important development direction in the current seismic exploration field. Due to the difficulty in data matching of multi-wave components, there is currently no practical inversion application of combined PP waves, PS1 waves and PS2 waves.

[0054] To this end, the present invention proposes a technical solution for AVO inversion by combining PP wave, PS1 wave and PS2 wave for fractured reservoirs, so as to provide more accurate fracture parameters for fractured reservoirs, such as fracture density of lithology and fluid.

[0055] Figure 1 1 is a flowchart of a multi-wave joint AVO inversion method for a fractured reservoir provided by an embodiment of the present invention, comprising the following steps:

[0056] S102, compressing the VSP seismic observation records corresponding to the PS wave, PS1 wave and PS2 wave of the fractured reservoir into the time domain of the first arrival of the zero-bias VSP longitudinal wave to perform AVA gather conversion to obtain the AVA seismic observation records.

[0057] Specifically, based on the relationship between the first arrival of the longitudinal wave and the depth of the detector, this step matches and compresses the VSP seismic observation records of the PP wave, PS1 wave and PS2 wave to the time domain of the first arrival of the zero-bias VSP longitudinal wave, and obtains the fast and slow shear wave delays and the crack orientation through Alford rotation; then, the VSP seismic observation records compressed to the time domain of the first arrival of the zero-bias VSP longitudinal wave are subjected to layer stripping of the R component and the T component to obtain the VSP seismic observation records of the PS1 wave and the PS2 wave respectively; the fast and slow shear wave delays τ are used to correct the time difference of the PS2 wave relative to the PS1 wave in the VSP seismic observation record; based on the P wave and S wave in the VSP seismic observation record, the VSP seismic observation record is corrected to obtain the time difference of the PS2 wave relative to the PS1 wave. The layer velocity of the wave is used to perform ray tracing on the PP wave, PS1 wave and PS2 wave (the input is the isotropic layer velocity obtained from the zero-bias VSP seismic observation record) to obtain the coordinate gathers of the incident angle and reflection point of the PP wave, PS1 wave and PS2 wave of each shot; the coordinate gathers of the incident angle and reflection point of the PP wave, PS1 wave and PS2 wave of each shot are converted into AVA gathers to obtain AVA seismic observation records.

[0058] S104, constructing an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropic parameters in AVA gather data, and the inversion objective function is a Bayesian theory function.

[0059] It should be understood that the initial model also includes isotropic parameters, that is, the fracture parameters are calculated based on the logging data of isotropic parameters and anisotropic parameters, which is not limited in this article.

[0060] Optionally, the initial model is used to calculate the fracture density based on the layer velocity and fast and slow shear wave delays in the AVA gather data.

[0061] The derivation formula of the initial model for the crack density is:

[0062] Where e represents the crack density, L represents the depth window converted from the initial arrival velocity during the Alford rotation scan, and V S1 represents the shear wave velocity of the S1 wave, τ represents the time delay between fast and slow shear waves, g=μ / (λ+2μ), μ and λ are Lamé inversion parameters.

[0063] Correspondingly, the formula of the objective function is: J(M)=[A prs -A syn (M ite ,i)] T [A prs - A syn (M ite ,i)] T +K(Mite ), where J(M) represents the objective function, i and k represent the incident angle and sampling point respectively; M represents the inversion parameter vector corresponding to the sampling point k. k ,μ k ,ρ k ,e k ] T , ρ represents the rock density, A prs Indicates the AVA record input in the inversion, M ite represents the iteratively updated inversion parameter, K(M ite ) represents the regularized Cauchy distribution term, A syn Indicates the AVA record used to update the initial model. syn (M ite ,i)=W*R(W ite ,i), W PP , W PS1 , W PS2 Represent the wavelet matrices of PP wave, PS1 and PS2 wave respectively;

[0064] The formula for the initial model to simulate the AVA seismic synthetic record according to the accurate reflection coefficient is R(M ite ,i)=[R PP (M ite ,i),R PS1 (M ite ,i),R PS2 (M ite ,i)] T , R PP(Mite,i) , R PS1 (M ite ,i) and R PS2 (M ite ,i) respectively represent the use of M corresponding to PP wave, PS1 and PS2 wave ite Calculate the reflection coefficient matrix at incident angle i.

[0065] S106, constructing an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate the AVA seismic synthetic record according to the precise reflection coefficient based on the initial model, and calculate the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model.

[0066] Optionally, the inversion equation for iterating the initial model is: ite =M0+ΔM; where M0 represents the initial model, ΔM represents the iterative increment, ΔM=-(H+λI)-1 G T κ, H represents the Hessian matrix, I represents the diagonal matrix, and κ represents the Jacobian matrix of the objective function G represents the second-order form of the Jacobian matrix of the objective function, They represent the Jacobian matrices of PP wave, PS1 wave and PS2 wave at incident angle i respectively.

[0067] Optionally, when the L2-norm of the objective function J(M) is less than a preset threshold, the inversion equation ends the iteration of the initial model.

[0068] S108, substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir.

[0069] It should be understood that the inversion parameters in this step also include isotropic parameters, and this article does not give examples and elaborate on the inversion of isotropic parameters.

[0070] S110, determining the fracture development state of the fracture reservoir based on the anisotropic parameters of the fracture reservoir.

[0071] The fracture development state may include, but is not limited to, information such as fracture density, which is not specifically limited herein.

[0072] The present invention can compress the VSP seismic observation records of the PS wave, PS1 wave and PS2 wave corresponding to the fractured reservoir into the time domain of the zero-bias VSP longitudinal wave first arrival for AVA gather conversion, thereby obtaining the AVA seismic observation records of the PS wave, PS1 wave and PS2 wave time matching without losing the fracture information; then, construct an initial model for deriving fracture parameters, and invert the anisotropic parameters of the fractured reservoir based on the inversion equation corresponding to the initial model. Among them, the inversion equation can simulate the AVA seismic synthetic record according to the precise reflection coefficient based on the initial model, and calculate the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model. This gradient calculation method can effectively solve the problem that there is no linear relationship between the precise reflection coefficient and the objective function, and can improve the inversion performance based on the precise reflection coefficient, and finally obtain more accurate anisotropic parameters for imaging the structure of the fractured reservoir.

[0073] The method of the embodiment of the present invention is described in detail below.

[0074] According to Bakulin's modeling theory, the embodiment of the present invention simplifies the inversion parameters into Lamé parameters λ, μ, rock density ρ and fracture density e. The use of accurate reflection coefficients is suitable for strong contrast and strong anisotropic formations. In order to solve the nonlinear problem caused by accurate reflection coefficients, a second-order Levenberg–Marquardt (LM) algorithm for Bayesian objective function is proposed to perform AVA inversion of HTI medium depth domain joint PP wave, PS1 wave and PS2 wave.

[0075] In order to protect the fracture information contained in the wave field, the VSP seismic observation records were preprocessed by using the processing method of protecting the data vector characteristics. Then, the data of the three components were matched and compressed to the zero-bias VSP longitudinal wave first arrival time domain based on the relationship between the first arrival of the longitudinal wave and the depth of the detector. The R component and the T component were layer-stripped by Alford rotation to obtain the VSP seismic observation records of the PS1 wave and the PS2 wave, and the fracture azimuth was calculated. After using τ to eliminate the time difference of the VSP seismic observation record of the PS2 wave, the PP wave, PS1 wave and PS2 wave matched to the zero-bias VSP P-wave first arrival time domain were obtained.

[0076] After eliminating the anisotropic effect, the incident angle and reflection point coordinates were obtained through isotropic ray tracing to convert the depth domain AVA gathers. A calculation formula for the initial fracture density model was proposed using the layer velocity obtained from zero-bias VSP seismic observation records and τ obtained from Alford rotation to achieve inversion.

[0077] Among them, the principle of the inversion equation of nonlinear least squares under Bayesian constraints is as follows:

[0078] The calculation of the reflection coefficient requires the stiffness coefficient c ij . The inversion parameters λ, μ, ρ and e are used to calculate c by Bakulin's modeling method. ij Here, the normal and tangential weaknesses of the crack are defined as Δ N and Δ T , then in dry cracks or cracks filled with fluid there are:

[0079]

[0080]

[0081] Where g = μ / (λ+2μ). In a fluid-filled fracture, there are:

[0082] Δ N =0 (3)

[0083] According to the linear slip theory, the crack orientation is The HTI medium stiffness coefficient matrix can be written as:

[0084] c 11 =(λ+2μ)(1-Δ N ),

[0085] c 13 =λ(1-Δ N ),

[0086]

[0087]

[0088] c 44 =μ,

[0089] c 55 =μ(1-Δ T ),

[0090] c 66 =c 55 . (4)

[0091] Establish the objective function based on the maximum probability of the posterior distribution of the Cauchy distribution: J(M) = [A prs -A syn (M ite ,i)] T [A prs -A syn (M ite ,i)] T +K(M ite ) (5)

[0092] Where J(M) represents the objective function, i and k represent the incident angle and sampling point respectively; M represents the inversion parameter vector corresponding to the sampling point k. k ,μ k ,ρ k ,e k ] T , ρ represents the rock density, A syn Indicates the AVA record input in the inversion, M ite represents the iteratively updated inversion parameter, K(M ite ) represents the regularized Cauchy distribution term, A syn Represents an AVA record used to update the initial model.

[0093] A syn It can be expressed as:

[0094] A syn (M ite ,i)=W*R(Mite ,i) (6)

[0095]

[0096] Among them, W PP , W PS1 , W PS2 Represent the wavelet matrices of PP wave, PS1 and PS2 wave respectively;

[0097] Formula (6) defines the formula for the initial model to simulate the AVA seismic synthetic record according to the accurate reflection coefficient forward modeling:

[0098] R(M ite ,i)=[R PP (M ite ,i),R PS1 (M ite ,i),R PS2 (M ite ,i)] T (8)

[0099] Among them, R PP (M ite ,i),R PS1 (M ite ,i) and R PS2 (M ite ,i) respectively represent the use of M corresponding to PP wave, PS1 and PS2 wave ite Calculate the reflection coefficient matrix at incident angle i.

[0100] For HTI media, the reflection coefficient is generally approximated by a linear equation. However, in order to make the inversion equation applicable to a wider range of conditions, including strong contrast and strong anisotropy, the exact reflection coefficient M is used in formula (8): ite Calculated by the iterative increment ΔM in each iteration and the initial model M0:

[0101] M ite =M0+ΔM (9)

[0102] Where ΔM is the correction amount of all inversion parameters at sampling point k.

[0103] ΔM=[Δλ k ,Δμ k ,Δρ k ,Δ e k ] T (10)

[0104] Since there is a nonlinear relationship between the inversion parameter vector M and the reflection coefficient matrix R, the LM algorithm can be used here to solve formula (10). Usually, the iteration parameter ΔM is solved by obtaining the first-order Taylor expansion of the objective function. Under weak contrast and weak anisotropy conditions, the inversion results under this approximation are reliable. However, in the case of strong anisotropy or strong contrast, the linear relationship between the anisotropy parameter and the reflection coefficient is poor.

[0105] Therefore, a second-order approximation of the objective function is proposed:

[0106]

[0107] The partial derivative with respect to ΔM is:

[0108]

[0109] For the partial derivative in formula (12), the Jacobian matrix of the objective function is defined as:

[0110]

[0111] And the second-order form:

[0112]

[0113] In the calculation of the exact reflection coefficient, it is difficult to obtain an explicit derivative expression due to matrix nesting. Therefore, the central difference method is used to obtain a second-order approximate solution of the first and second derivatives of the objective function with respect to the inversion parameters.

[0114]

[0115]

[0116] Where p represents four inversion parameters, represents the difference step size. A small value is usually given based on experience. However, if If the given value is not within a reasonable order of magnitude, it will lead to non-convergence of multi-parameter inversion. Therefore, in order to improve the stability of differential calculation, the standard deviation p / σ of the parameters is used p To limit the difference step size:

[0117]

[0118]

[0119] Substituting formula (13) and (14) into formula (12), for the lth iteration, we have:

[0120] κ(M0+ΔMl )+G(M0+ΔM l )ΔM (i) =0, (19)

[0121] The above formula can be written as:

[0122] HΔM (i) =-G T κ, (20)

[0123] Where H is the Hessian matrix, H = G T G. Adding the damping factor to formula (20) yields the model update using the second-order LM iterative algorithm:

[0124] ΔM=-(H+λI) -1 G T κ (21)

[0125] In each iteration, ΔM is added to the initial model M0 to update M ite , and then update A using formula (6) syn , until the L2-norm of J(M) is small enough to stop the inversion.

[0126] Correspondingly, under the layer matching principle of PP wave, PS1 wave and PS2 wave:

[0127] Step 1, such as Figure 3 As shown in the figure, the first arrival time t0 of the zero-bias P-wave is picked up by using the first arrival wave recorded by zero-bias VSP seismic observation. PFB The relationship between the detector number and Figure 3 a), and the layer velocity models of P and S waves ( Figure 3 b). This is crucial for the subsequent layer matching, AVA gather conversion and the establishment of the initial model. Because the incident angle of the zero-bias VSP seismic observation record is approximately 0, the seismic wave mainly exhibits the isotropic properties of the medium when propagating vertically. Therefore, the layer velocity model obtained can be approximated to the isotropic background velocity.

[0128] Step 2, such as Figure 4 As shown in the figure, taking a shot at an offset distance of 1012 m as an example, the VSP seismic observation records of the walkaround three-component (PS wave, PS1 wave and PS2 wave) are preprocessed by amplitude preservation methods such as horizontal rotation, noise suppression and spherical diffusion compensation while preserving the amplitude characteristic relationship between the three components.

[0129] Step 3, such as Figure 5 As shown, the uplink and downlink waves are separated for the preprocessed data in step 2 to obtain the separated direct wave (downlink wave) and reflected wave (uplink wave).

[0130] Step 4: From the direct wave record obtained in step 3, the first arrival time t of the longitudinal wave corresponding to different detector numbers can be picked up. PFB , through which the detector number -t is obtained PFB It should be noted that the detector number -t obtained for different offset shots PFB The relationship is different.

[0131] Step 5, such as Figure 6 As shown, using the detector number obtained in step 4 -t PFB Relationship, the reflected wave field obtained in step 3 is compressed to t according to the detector number. PFB The well-matched strong reflection axes on the three components represent the geological landmark layers T1j1 (for example: the interface between the Xiangxi Formation and the Jialingjiang Formation of the Triassic System), the bottom interface of the target layer T1j1 and T1f1 (the interface between the Permian System and the Triassic System).

[0132] Step 6 At this point, the layer matching between the Z, R, and T components has been completed, but due to the different detector numbers of the guns -t PFB The relationships are different, and layer matching has not been performed between different guns. Figure 7 As shown, replace t in step 5 PFB The Z, R, and T components of the domain are recorded according to the detector number in step 1 - t0 PFB Relationship, match the Z, R, T components of all guns to t0 PFB domain. At this time, the layer matching of the Z, R, and T components of different shots is completed. It is found that due to the development of cracks, the split S1 wave and S2 wave mix at different speeds. In the range surrounded by the curve, the number of reflection phase axes is greater in the R component than in the Z component. In addition, since the polarization direction of the split shear wave is affected by the crack azimuth and is not parallel to the ZR plane, Figure 7 The T component records the projection of the split shear wave. All these phenomena are distributed near the target formation T1j1, that is, near the actual fracture reservoir.

[0133] Step 7: The Z component data can be considered as PP waves, but the R and T components are a mixture of PS1 and PS2 waves, and wave field separation is required. Figure 8 As shown, replace t0 in step 6 PFB The R and T components of the domain are scanned by Alford rotation to obtain the fast and slow shear wave delay τ and the crack orientation In order to eliminate the influence of the crack development in the overburden on the Alford rotation results of the underlying layer, The R component and the T component are layer-stripped to obtain the PS1 and PS2 wave records, and the PS2 wave time difference is eliminated by the fast and slow shear wave delay τ to match the PS1 wave record. Fig. 9 As shown in Figure 3, within the range of crack development, the wave field mixing phenomenon of split shear waves is eliminated, and the reflection phase axes of PP waves and separated PS1 waves are better matched than those of Z and R components.

[0134] So far, we have obtained the matching t0 PFB PP, PS1, PS2 wave cannon records in the domain.

[0135] Correspondingly, refer to Figure 3 As shown in (b), the principle of AVA gather generation is as follows:

[0136] Step 8: In traditional AVA gather generation, the incident angle is calculated by ray tracing based on the layer velocity. However, in pre-stack data processing, the S2 wave layer velocity is difficult to obtain, which makes it difficult to generate AVA gathers. Through the above step 7, the wave fields of PP wave, PS1 wave, and PS2 wave have been flattened, and the PS2 wave time difference caused by anisotropy has been eliminated, matching it to t0. PFB For PS2 waves, after time difference correction, the S-wave velocity of the isotropic background can be used for simulation. Figure 3 As shown in B, the isotropic P-wave and S-wave layer velocities obtained in step 1 are used for ray tracing. The incident angle and reflection point coordinates of the PP and PS waves of each shot obtained by ray tracing are obtained as follows after natural interpolation: Fig.10 As shown in Figure 2, it should be noted that the PS wave represents the upgoing converted shear wave before splitting. Because the PP wave and the PS wave have the same geological reflection point, the Z coordinates of the PP wave and the PS wave reflection point are the same.

[0137] Step 9: Based on the incident angle and reflection point coordinate gathers obtained in step 8, the shot records of all matching PP waves, PS1 waves, and PS2 waves obtained in step 7 are converted into depth domain AVA gathers.

[0138] Correspondingly, refer to Figure 2 c shows the principle of establishing the initial model:

[0139] Step 10: Since the P-wave and S-wave in zero-offset VSP data propagate in isotropic planes and are not affected by vertical cracks when the incident angle is zero, the interval velocity ( Figure 3 b) P-wave velocity and S-wave velocity V as isotropic background P and V S The corresponding initial model of inversion parameters λ, μ, and ρ can be calculated by the following formula:

[0140] ρ Gardner =-0.0296V P 2 +0.461VP +0.963. (22)

[0141]

[0142]

[0143] in, Fig.11 Schematic diagram of the initial model for the inversion parameters λ, μ, and ρ.

[0144] Step 11: Substitute the fast and slow shear wave delays τ and crack orientations in step 7 ( Figure 8 ) Using a method similar to step 9, the incident angle and reflection point coordinate gathers are shifted to Fig.12 The spatial coordinates of the imaging range. The results beside the well are as follows Fig.12 Crack orientation It will be used to calculate the reflection coefficient later, but for the initial model e, it needs to be converted with τ.

[0145] Step 12, derive the initial model formula of crack density e.

[0146] For the anisotropy parameter γ defined by Thomsen (Ruger, 2002):

[0147]

[0148] Combining equation (24) with equation (2), we get:

[0149]

[0150] The travel time of PS1 wave is known to be tPS1:

[0151]

[0152] Where L is the time window scanned during Alford rotation, which is converted into a depth window using the initial arrival velocity. Then the time difference between the fast and slow waves is:

[0153]

[0154] Using formula (2, 25, 27), the crack density formula can be calculated:

[0155]

[0156] Substitute the interval velocity obtained in step 1 and the wellside τ obtained in step 11 into formula (28) to calculate e as Fig.13 shown.

[0157] Actual data inversion results

[0158] Test the VSP data near Luzhou, Sichuan (the above figure also uses this data). From step 9, we get Fig.14 The depth domain AVA gathers shown ( Fig.14 ), the actual data AVA gathers generated match well with the synthetic AVA gathers calculated using the initial model. The results after AVA gather stacking are as follows Fig.16 As shown in Figure 2, the inversion range is the velocity ratio. The inversion uses the angle-independent depth-domain statistical wavelet extracted from the AVA gather. Using the inversion, we get Fig.17 The inversion results are shown.

[0159] Fig.17 The inversion results of isotropic background parameters representing lithological characteristics are consistent with the geological interface. The high fracture density range is 2420–2490 m, which is consistent with the depth of the fractured reservoir verified by geological background and production data.

[0160] Fig.18 1 is a schematic diagram of the structure of a multi-wave joint AVO inversion device for a fractured reservoir provided by an embodiment of the present invention, comprising:

[0161] The multi-wave matching module 1810 is used to compress the VSP seismic observation records corresponding to the PS wave, PS1 wave and PS2 wave of the fractured reservoir into the time domain of the first arrival of the zero-bias VSP longitudinal wave to perform AVA gather conversion to obtain the AVA seismic observation records.

[0162] The model function construction module 1820 is used to construct an initial model and an inversion objective function corresponding to the initial model based on AVA seismic observation records, wherein the initial model is used to calculate fracture parameters based on logging data containing anisotropy parameters in the AVA data set, and the inversion objective function is a Bayesian theoretical function.

[0163] The inversion equation construction model 1830 is used to construct the inversion equation corresponding to the initial model, and the inversion equation is used to forward simulate the AVA seismic synthetic record according to the precise reflection coefficient based on the initial model, and calculate the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative results of the initial model.

[0164] The inversion execution module 1840 is used to substitute the initial model into the inversion equation for inversion to obtain the anisotropic parameters of the fractured reservoir.

[0165] The fracture state determination module 1850 is used to determine the fracture state of the fractured reservoir based on the anisotropic parameters of the fractured reservoir.

[0166] Optionally, the multi-wave matching module 1810 compresses the VSP seismic observation records of the PS wave, PS1 wave and PS2 wave corresponding to the fractured reservoir into the time domain of the zero-bias VSP P-wave first arrival for AVA gather conversion to obtain the AVA seismic observation records, including: based on the relationship between the P-wave first arrival and the detector depth, matching and compressing the VSP seismic observation records of the PP wave, PS1 wave and PS2 wave into the time domain of the zero-bias VSP P-wave first arrival, and obtaining the fast and slow S-wave delays and the fracture orientation through Alford rotation; performing layer stripping of the R component and the T component on the VSP seismic observation records compressed into the time domain of the zero-bias VSP P-wave first arrival to obtain the VSP seismic observation records of the PS1 wave and the PS2 wave respectively; using the fast and slow S-wave delays to correct the time difference of the PS2 wave relative to the PS1 wave in the VSP seismic observation record; performing ray tracing on the PP wave, PS1 wave and PS2 wave based on the layer velocities of the P wave and the S wave in the VSP seismic observation record to obtain the PP wave, PS1 wave and PS2 wave of each shot. The coordinate gathers of the incident angles and reflection points of the PP waves, PS1 waves and PS2 waves of each shot are converted into AVA gathers to obtain AVA seismic observation records.

[0167] Optionally, the initial model is used to calculate the fracture density based on the layer velocity and fast and slow shear wave delays in the AVA gather data.

[0168] Optionally, the derivation formula of the initial model for the crack density is:

[0169] Where e represents the crack density, L represents the depth window converted from the initial arrival velocity during the Alford rotation scan, and V S1 represents the shear wave velocity of the S1 wave, τ represents the time delay between fast and slow shear waves, g=μ / (λ+2μ), μ and λ are Lamé inversion parameters.

[0170] Optionally, the objective function is formulated as follows: J(M)=[A prs -A syn (M ite ,i)] T [A prs - A syn (M ite ,i)] T +K(M ite ), where J(M) represents the objective function, i and k represent the incident angle and sampling point respectively; M represents the inversion parameter vector corresponding to the sampling point k. k ,μ k ,ρ k ,e k ]T , ρ represents the rock density, A prs Indicates the AVA record input in the inversion, M ite represents the iteratively updated inversion parameter, K(M ite ) represents the regularized Cauchy distribution term, A syn Indicates the AVA record used to update the initial model. syn (M ite ,i)=W*R(M ite ,i), W PP , W PS1 , W PS2 Represent the wavelet matrices of PP wave, PS1 and PS2 wave respectively;

[0171] The formula for the initial model to simulate the AVA seismic synthetic record according to the accurate reflection coefficient is R(M ite ,i)=[R PP (M ite ,i),R PS1 (M ite ,i),R PS2 (M ite ,i)] T , R PP (M ite ,i), R PS1 (M ite ,i) and R PS2 (M ite ,i) respectively represent the use of M corresponding to PP wave, PS1 and PS2 wave ite Calculate the reflection coefficient matrix at incident angle i.

[0172] Optionally, the inversion equation for iterating the initial model is: ite =M0+ΔM; where M0 represents the initial model, ΔM represents the iterative increment, ΔM=-(H+λI) -1 G T κ, H represents the Hessian matrix, I represents the diagonal matrix, and κ represents the Jacobian matrix of the objective function G represents the second-order form of the Jacobian matrix of the objective function, They represent the Jacobian matrices of PP wave, PS1 wave and PS2 wave at incident angle i respectively.

[0173] Optionally, when the L2-norm of J(M) of the objective function is less than a preset threshold, the inversion equation ends the iteration of the initial model.

[0174] Fig.18 This is a schematic diagram of the structure of an electronic device according to an embodiment of this specification. Fig.18 At the hardware level, the electronic device includes a processor, and optionally also includes an internal bus, a network interface, and a memory. The memory may include a memory, such as a high-speed random access memory (RAM), and may also include a non-volatile memory (non-volatile memory), such as at least one disk storage. Of course, the electronic device may also include hardware required for other services.

[0175] The processor, network interface and memory can be interconnected through an internal bus, which can be an ISA (Industry Standard Architecture) bus, a PCI (Peripheral Component Interconnect) bus or an EISA (Extended Industry Standard Architecture) bus. The bus can be divided into an address bus, a data bus, a control bus, etc. For ease of representation, Fig.18 Only one bidirectional arrow is used in the diagram, but this does not mean that there is only one bus or only one type of bus.

[0176] The memory is used to store the program. Specifically, the program may include a program code, and the program code includes a computer operation instruction. The memory may include a memory and a non-volatile memory, and provides instructions and data to the processor.

[0177] The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs it, forming the above-mentioned multi-wave joint AVO inversion device at the logical level. Correspondingly, the processor executes the program stored in the memory and is specifically used to perform the following operations:

[0178] The VSP seismic observation records corresponding to PS waves, PS1 waves and PS2 waves of fractured reservoirs are compressed into the time domain of the first arrival of zero-bias VSP longitudinal waves to perform AVA gather conversion and obtain AVA seismic observation records.

[0179] Constructing an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropy parameters in AVA gather data, and the inversion objective function is a Bayesian theory function;

[0180] Constructing an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate AVA seismic synthetic records based on the initial model according to the precise reflection coefficient, and calculating the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model;

[0181] Substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir;

[0182] Based on the anisotropic parameters of the fractured reservoir, the fracture state of the fractured reservoir is determined.

[0183] The above is as in this manual Figure 1 The method disclosed in the illustrated embodiment can be applied to a processor and implemented by the processor. The processor may be an integrated circuit chip with signal processing capabilities. In the implementation process, each step of the above method can be completed by an integrated logic circuit of hardware in the processor or an instruction in the form of software. The above processor can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components. The methods, steps and logic block diagrams disclosed in the embodiments of the present invention can be implemented or executed. The general-purpose processor can be a microprocessor or the processor can also be any conventional processor, etc. The steps of the method disclosed in conjunction with the embodiments of the present invention can be directly embodied as a hardware decoding processor for execution, or can be executed by a combination of hardware and software modules in the decoding processor. The software module can be located in a mature storage medium in the field such as random access memory, flash memory, read-only memory, programmable read-only memory or electrically erasable programmable memory, register, etc. The storage medium is located in the memory, and the processor reads the information in the memory and completes the steps of the above method in combination with its hardware.

[0184] It should be understood that the electronic device of the embodiment of the present invention can enable the service processing device to implement the corresponding Figure 1 The steps and functions in the method shown in the figure are not repeated in this article because the principles are the same.

[0185] Of course, in addition to software implementation, the electronic device of this specification does not exclude other implementation methods, such as logic devices or a combination of software and hardware, etc., that is to say, the execution subject of the following processing flow is not limited to each logic unit, but can also be hardware or logic devices.

[0186] In addition, an embodiment of the present invention further provides a computer-readable storage medium, which stores one or more programs, and the one or more programs include instructions.

[0187] When the above instructions are executed by a portable electronic device including multiple application programs, the portable electronic device can execute Figure 1 The steps of the method shown include:

[0188] The VSP seismic observation records corresponding to PS waves, PS1 waves and PS2 waves of fractured reservoirs are compressed into the time domain of the first arrival of zero-bias VSP longitudinal waves and converted into AVA gathers to obtain AVA seismic observation records.

[0189] An initial model and an inversion objective function corresponding to the initial model are constructed, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropic parameters in AVA gather data, and the inversion objective function is a Bayesian theory function.

[0190] An inversion equation corresponding to the initial model is constructed, and the inversion equation is used to forward simulate the AVA seismic synthetic record according to the precise reflection coefficient based on the initial model, and calculate the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model.

[0191] The initial model is substituted into the inversion equation for inversion to obtain the anisotropic parameters of the fractured reservoir.

[0192] Based on the anisotropic parameters of the fractured reservoir, the fracture state of the fractured reservoir is determined.

[0193] Those skilled in the art will appreciate that the embodiments of this specification may be provided as methods, systems or computer program products. Therefore, this specification may take the form of a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware. Moreover, this specification may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0194] The above is a description of a specific embodiment of the specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in an order different from that in the embodiments and still achieve the desired results. In addition, the processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0195] The above are only embodiments of this specification and are not intended to limit this specification. For those skilled in the art, this specification may be subject to various changes and modifications. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of this specification shall be included in the scope of the claims of this specification. In addition, all other embodiments obtained by ordinary technicians in this field without creative work shall fall within the scope of protection of this document.

Claims

1. A multi-wave joint inversion method for anisotropic parameters of fractured reservoirs, characterized in that: include: The VSP seismic observation records of the fractured reservoir corresponding to the PS wave, PS1 wave and PS2 wave are compressed into the time domain of the first arrival of the zero-bias VSP longitudinal wave and converted into AVA gathers to obtain the AVA seismic observation records. Constructing an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropy parameters in AVA gather data, and the inversion objective function is a Bayesian theory function; Constructing an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate AVA seismic synthetic records based on the initial model according to the precise reflection coefficient, and calculating the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model; Substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir; Determining the fracture development state of the fracture reservoir based on the anisotropic parameters of the fracture reservoir; The initial model is used to calculate the fracture density based on the layer velocity and fast and slow shear wave delays in the AVA gather data.

2. The method according to claim 1, characterized in that The VSP seismic observation records of the fractured reservoir corresponding to the PS wave, PS1 wave and PS2 wave are compressed into the time domain of the zero-bias VSP longitudinal wave first arrival and converted into AVA gathers to obtain the AVA seismic observation records, including: Based on the relationship between the first arrival of the P-wave and the depth of the geophone, the VSP seismic observation records of the PP wave, PS1 wave and PS2 wave are matched and compressed into the time domain of the first arrival of the zero-bias VSP P-wave, and the fast and slow S-wave delays and the fracture orientation are obtained through Alford rotation. The VSP seismic observation records compressed to the time domain of the first arrival of the zero-bias VSP longitudinal wave are subjected to layer stripping of the R component and the T component to obtain the VSP seismic observation records of the PS1 wave and the PS2 wave respectively; The time difference between PS2 wave and PS1 wave in VSP seismic observation record is corrected by using fast and slow shear wave delays; Based on the layer velocities of P and S waves in the VSP seismic observation records, PP, PS1 and PS2 waves are ray traced to obtain the incident angles and coordinate gathers of the reflection points of PP, PS1 and PS2 waves for each shot; The coordinate gathers of the incident angles and reflection points of the PP wave, PS1 wave and PS2 wave of each shot are converted into AVA gathers to obtain AVA seismic observation records.

3. The method according to claim 1, characterized in that The derivation formula of the initial model for the crack density is: Where e represents the crack density, L represents the depth window converted from the initial arrival velocity during the Alford rotation scan, and V S1 represents the shear wave velocity of S1 wave, τ represents the time delay between fast and slow shear waves, g=μ / (λ+2μ), μ and λ are Lamé inversion parameters.

4. The method according to claim 3, characterized in that The formula of the objective function is: J(M) = [A prs -A syn (M ite ,i)] T [A prs -A syn (M ite ,i)] T +K(M ite ), where J(M) represents the objective function, i and k represent the incident angle and sampling point respectively; M represents the inversion parameter vector corresponding to the sampling point k. k ,μ k ,ρ k ,e k ] T , ρ represents the rock density, A prs Indicates the AVA record input in the inversion, M ite represents the iteratively updated inversion parameter, K(M ite ) represents the regularized Cauchy distribution term, A syn Represents the AVA record used to update the initial model. W PS2 Represent the wavelet matrices of PP wave, PS1 and PS2 wave respectively; The formula for the initial model to simulate the AVA seismic synthetic record according to the accurate reflection coefficient is R(M ite ,i)=[R PP (M ite ,i),R PS1 (M ite ,i),R PS2 (M ite ,i)] T , R PP (M ite ,i),R PS1 (M ite ,i) and R PS2 (M ite ,i) respectively represent the use of M corresponding to PP wave, PS1 and PS2 wave ite Calculate the reflection coefficient matrix at incident angle i.

5. The method according to claim 4, characterized in that The iterative formula for the initial model in the inversion equation is: ite =M0+ΔM; where M0 represents the initial model, ΔM represents the iterative increment, ΔM=-(H+λI) -1 G T κ, H represents the Hessian matrix, I represents the diagonal matrix, and κ represents the Jacobian matrix of the objective function G represents the second-order form of the Jacobian matrix of the objective function, They represent the Jacobian matrices of PP wave, PS1 wave and PS2 wave at incident angle i respectively.

6. The method according to claim 5, characterized in that The inversion equation ends the iteration of the initial model when the L2-norm of J(M) of the objective function is less than a preset threshold.

7. A multi-wave joint inversion device for anisotropic parameters of fractured reservoirs, characterized in that: include: The multi-wave matching module is used to compress the VSP seismic observation records corresponding to the PS wave, PS1 wave and PS2 wave of the fractured reservoir into the time domain of the first arrival of the zero-bias VSP longitudinal wave to perform AVA gather conversion to obtain the AVA seismic observation records; A model function construction module, used to construct an initial model and an inversion objective function corresponding to the initial model, wherein the initial model is used to calculate fracture parameters according to logging data containing anisotropy parameters in AVA gather data, and the inversion objective function is a Bayesian theory function; An inversion equation construction model is used to construct an inversion equation corresponding to the initial model, wherein the inversion equation is used to forward simulate AVA seismic synthetic records based on the initial model according to the precise reflection coefficient, and calculate the gradient of the initial model according to the partial derivative of the forward operator and the second-order approximate solution of the inversion objective function to iterate the initial model, so as to invert the anisotropic parameters of the fractured reservoir according to the iterative result of the initial model; An inversion execution module, used for substituting the initial model into the inversion equation for inversion to obtain anisotropic parameters of the fractured reservoir; A fracture state determination module, used for determining the fracture state of the fracture reservoir based on anisotropic parameters of the fracture reservoir; The initial model is used to calculate the fracture density based on the layer velocity and fast and slow shear wave delays in the AVA gather data.

8. An electronic device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the computer program is executed by the processor to execute the method according to claim 1.

9. A computer-readable storage medium having a computer program stored thereon, wherein the computer program implements the method of claim 1 when executed by a processor.

Images