Five-dimensional seismic data step-by-step inversion method, device and electronic equipment for shale reservoirs

Abstract

This invention provides a method, apparatus, and electronic device for step-by-step inversion of five-dimensional seismic data in shale reservoirs. The method includes: acquiring five-dimensional seismic data and acquiring fracture dip angles; inverting a target azimuth elastic impedance data set from the five-dimensional seismic data, the target azimuth elastic impedance data set including target azimuth elastic impedance data under various angle combinations; determining the target azimuth elastic impedance difference function based on the fracture dip angle; and inverting and predicting fracture weakness parameters based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function; determining a target incident angle-related elastic impedance function; and inverting and predicting anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted and predicted fracture weakness parameters. Embodiments of this invention can improve the accuracy of the inversion results.

Description

Technical Field

[0001] This invention relates to the field of oil and gas exploration technology, and in particular to a method, apparatus and electronic equipment for step-by-step inversion of five-dimensional seismic data of shale reservoirs. Background Technology

[0002] The exploration and development of shale oil and gas reservoirs is currently a research hotspot and a challenge. Accurately characterizing their anisotropy is crucial for improving the accuracy of seismic interpretation and the reliability of reservoir prediction. Fractures are a key component of shale reservoirs, controlling not only the accumulation and migration of oil and gas but also determining the recovery rate. Related technologies typically focus on shale reservoirs with vertical fractures (equivalent to orthotropic media), but inclined fractures are also commonly found in shale reservoirs. This significantly affects the propagation and reflection amplitude of seismic waves, leading to lower accuracy in inversion results. Therefore, a satisfactory solution for improving the accuracy of inversion results has not yet been found. Summary of the Invention

[0003] In view of this, embodiments of the present invention provide a step-by-step inversion method, apparatus, and electronic equipment for five-dimensional seismic data of shale reservoirs to solve the problem of low accuracy of inversion results caused by related technologies. That is, embodiments of the present invention can determine the target azimuth elastic impedance difference function through fracture dip angle, etc., to invert fracture weakness parameters, and invert anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance through the elastic impedance function related to the target incident angle. By fully considering the influence of tilted fractures, the accuracy of inversion results can be effectively improved. In other words, embodiments of the present invention can effectively improve the accuracy of inversion results through step-by-step inversion, thereby effectively improving the accuracy of shale reservoir prediction.

[0004] According to one aspect of the present invention, a step-by-step inversion method for five-dimensional seismic data of shale reservoirs is provided, the method comprising: Five-dimensional seismic data and fracture dip angles are acquired. The five-dimensional seismic data includes multiple seismic traces, with each seismic trace corresponding to an incident angle and a relative azimuth angle. The target azimuth elastic impedance data set is derived from the five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data under each angle combination in multiple angle combinations. An angle combination includes any incident angle among J incident angles and any relative azimuth angle among K relative azimuth angles, where J and K are both positive integers. Based on the crack dip angle, the target azimuth elastic impedance difference function is determined, and the parameters to be inverted in the target azimuth elastic impedance difference function include crack weakness parameters; and based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, the crack weakness parameters are inverted and predicted. The target incident angle-related elastic impedance function is determined, and the parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance; and based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance are inverted and predicted.

[0005] According to another aspect of the present invention, a step-by-step inversion device for five-dimensional seismic data of shale reservoirs is provided, the device comprising: The acquisition unit is used to acquire five-dimensional seismic data and to acquire fracture dip angle. The five-dimensional seismic data includes multiple seismic traces, and each seismic trace corresponds to an incident angle and a relative azimuth angle. The processing unit is used to invert the target azimuth elastic impedance data set from the five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data under each angle combination in multiple angle combinations. An angle combination includes any incident angle among J incident angles and any relative azimuth angle among K relative azimuth angles, where J and K are both positive integers. The processing unit is further configured to determine the target azimuth elastic impedance difference function based on the crack dip angle, wherein the parameters to be inverted in the target azimuth elastic impedance difference function include crack weakness parameters; and to invert and predict the crack weakness parameters based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function. The processing unit is further configured to determine the target incident angle-related elastic impedance function, wherein the parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance; and based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, to invert and predict the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance.

[0006] According to another aspect of the present invention, an electronic device is provided, the electronic device including a processor and a memory storing a program, wherein the program includes instructions that, when executed by the processor, cause the processor to perform the methods mentioned above.

[0007] According to another aspect of the present invention, a non-transitory computer-readable storage medium storing computer instructions for causing a computer to perform the methods mentioned above is provided.

[0008] According to another aspect of the present invention, a computer program product is provided, comprising a computer program, wherein the computer program, when executed by a processor, is used to cause the computer to perform the methods mentioned above.

[0009] This invention can acquire five-dimensional seismic data and fracture dip angles. The five-dimensional seismic data includes multiple seismic traces, each corresponding to an incident angle and a relative azimuth. A target azimuth elastic impedance data set is then derived from the five-dimensional seismic data. This set includes target azimuth elastic impedance data for each of multiple angle combinations, where each angle combination includes any one of J incident angles and any one of K relative azimuth angles. Based on the fracture dip angle, a target azimuth elastic impedance difference function is determined. The parameters to be inverted in this function include fracture weakness parameters. Finally, based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, fracture weakness parameters are predicted and inverted. Furthermore, the target incident angle-related elastic impedance function can be determined. The parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance. Based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted fracture weakness parameters, the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance are inverted and predicted. It can be seen that the embodiments of the present invention can determine the target azimuth elastic impedance difference function through fracture dip angle, etc., to invert fracture weakness parameters, and invert the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance through the target incident angle-related elastic impedance function. This fully considers the influence of tilted fractures and can effectively improve the accuracy of the inversion results. That is, the embodiments of the present invention can effectively improve the accuracy of the inversion results through step-by-step inversion, thereby effectively improving the accuracy of shale reservoir prediction, etc. Attached Figure Description

[0010] Further details, features, and advantages of the invention are disclosed in the following description of exemplary embodiments in conjunction with the accompanying drawings, in which: Figure 1 A flowchart illustrating a step-by-step inversion method for five-dimensional seismic data of shale reservoirs according to an exemplary embodiment of the present invention is shown. Figure 2 A schematic diagram of a monoclinic medium model according to an exemplary embodiment of the present invention is shown; Figure 3 A schematic diagram illustrating the effect of variations in isotropic background parameters and anisotropic parameters on the reflectance coefficient according to an exemplary embodiment of the present invention is shown. Figure 4 A schematic diagram illustrating the effect of a change in crack weakness parameters on the reflection coefficient according to an exemplary embodiment of the present invention is shown. Figure 5 A schematic diagram illustrating a comparison of the accuracy of reflection coefficient curves according to an exemplary embodiment of the present invention is shown; Figure 6 A flowchart illustrating another step-by-step inversion method for five-dimensional seismic data of shale reservoirs according to an exemplary embodiment of the present invention is shown; Figure 7 A schematic diagram of a synthetic azimuth seismic gather according to an exemplary embodiment of the present invention is shown; Figure 8 A schematic diagram of an orientation elastic impedance inversion result according to an exemplary embodiment of the present invention is shown; Figure 9 A schematic diagram of a model parameter inversion result according to an exemplary embodiment of the present invention is shown; Figure 10 A schematic diagram of a profile of parameter inversion results for a fractured shale oil reservoir model according to an exemplary embodiment of the present invention is shown. Figure 11 A schematic diagram of a well logging cross-plot analysis according to an exemplary embodiment of the present invention is shown; Figure 12 A schematic block diagram of a step-by-step inversion apparatus for five-dimensional seismic data of shale reservoirs according to an exemplary embodiment of the present invention is shown. Figure 13 A structural block diagram of an exemplary electronic device that can be used to implement embodiments of the present invention is shown. Detailed Implementation

[0011] Embodiments of the present invention will now be described in more detail with reference to the accompanying drawings. While some embodiments of the invention are shown in the drawings, it should be understood that the invention can be implemented in various forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided to provide a more thorough and complete understanding of the invention. It should be understood that the accompanying drawings and embodiments are for illustrative purposes only and are not intended to limit the scope of protection of the invention.

[0012] It should be understood that the various steps described in the method embodiments of the present invention may be performed in different orders and / or in parallel. Furthermore, the method embodiments may include additional steps and / or omit the steps shown. The scope of the present invention is not limited in this respect.

[0013] The term "comprising" and its variations as used herein are open-ended, meaning "including but not limited to". The term "based on" means "at least partially based on". The term "one embodiment" means "at least one embodiment"; the term "another embodiment" means "at least one additional embodiment"; the term "some embodiments" means "at least some embodiments". Definitions of other terms will be given in the following description. It should be noted that the concepts of "first", "second", etc., mentioned in this invention are used only to distinguish different devices, modules, or units, and are not intended to limit the order of functions performed by these devices, modules, or units or their interdependencies.

[0014] It should be noted that the terms "a" and "a plurality of" used in this invention are illustrative rather than restrictive. Those skilled in the art should understand that, unless otherwise expressly indicated in the context, they should be understood as "one or more".

[0015] The names of the messages or information exchanged between the multiple devices in the embodiments of the present invention are for illustrative purposes only and are not intended to limit the scope of these messages or information.

[0016] It should be noted that the execution subject of the step-by-step inversion method for shale reservoir five-dimensional seismic data provided in this embodiment of the invention (also referred to as the step-by-step inversion method based on shale reservoir five-dimensional seismic data or the step-by-step inversion method based on five-dimensional seismic data, etc.) can be one or more electronic devices, and this invention does not limit this. The electronic device can be a terminal (i.e., a client) or a server. When the execution subject includes multiple electronic devices, and these multiple electronic devices include at least one terminal and at least one server, the step-by-step inversion method for shale reservoir five-dimensional seismic data provided in this embodiment of the invention can be jointly executed by the terminal and the server. Accordingly, the terminal mentioned here may include, but is not limited to: smartphones, laptops, desktop computers, intelligent voice interaction devices, etc. The server mentioned here can be an independent physical server, a server cluster or distributed system composed of multiple physical servers, or a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communication, middleware services, domain name services, security services, CDN (Content Delivery Network), and big data and artificial intelligence platforms, etc.

[0017] Based on the above description, this invention proposes a step-by-step inversion method for five-dimensional seismic data of shale reservoirs (also referred to as a reservoir prediction method). This step-by-step inversion method for five-dimensional seismic data of shale reservoirs can be executed by the aforementioned electronic device (terminal or server); or, this step-by-step inversion method for five-dimensional seismic data of shale reservoirs can be executed jointly by a terminal and a server. For ease of explanation, the following description will use the execution of this step-by-step inversion method for five-dimensional seismic data of shale reservoirs by an electronic device as an example; such as Figure 1 As shown, the step-by-step inversion method for five-dimensional seismic data of shale reservoirs may include the following steps S101-S104: S101, acquire five-dimensional seismic data and acquire fracture dip angle. The five-dimensional seismic data includes multiple seismic traces, and each seismic trace corresponds to an incident angle and a relative azimuth angle.

[0018] Optionally, the aforementioned five-dimensional seismic data (also referred to as observed five-dimensional seismic data or target five-dimensional seismic data, etc.) can be five-dimensional seismic data of the target area. Optionally, the target area can be any area, and this embodiment of the invention does not limit it. The five-dimensional seismic data may include three spatial dimensions, offset distance (i.e., incident angle, which can be simply referred to as incident angle), and azimuth (i.e., azimuth angle); based on this, the five-dimensional seismic data can be a set of seismic data (i.e., seismic trace data) covering three spatial dimensions + incident angle + azimuth angle information.

[0019] In this embodiment of the invention, the methods for acquiring five-dimensional seismic data may include, but are not limited to, the following: The first acquisition method: The electronic device stores five-dimensional seismic data in its own storage space. In this case, the electronic device can acquire five-dimensional seismic data from its own storage space.

[0020] The second method of acquisition: electronic devices can obtain earthquake data download links and use those links to download five-dimensional earthquake data, etc.

[0021] Optionally, the fracture dip angle can be the dip angle of a set of inclined fractures in the target area; in this embodiment of the invention, the fracture dip angle can be prior information for step-by-step inversion; optionally, the above-mentioned fracture dip angle can also be called prior information on fracture dip angle. Optionally, under the assumptions of long-wavelength seismic approximation and weak anisotropy, a shale reservoir containing a set of directionally arranged inclined fractures can be equivalent to a monoclinic medium, such as... Figure 2 As shown; where, Figure 2 A set of diagonal lines (i.e., solid diagonal lines) arranged along the z-axis section can represent a set of inclined fractures. In this case, a shale reservoir containing a set of inclined fractures can be equivalent to a monoclinic medium, and Figure 2 In this context, v can represent the crack dip angle (i.e., the angle between the axis of symmetry of the inclined crack and the z-axis of the coordinate system); and... Figure 2In the horizontal cross-section, the dashed lines other than the coordinate axes can represent ray projections (i.e., the projection of the plane formed by the incident P-wave (also called the incident P-wave) and the reflected P-wave (also called the reflected P-wave), θ can represent the incident angle, and φ can represent the seismic observation azimuth (also called the seismic observation azimuth angle or absolute azimuth angle). It should be understood that a relative azimuth angle can be the angle between a seismic observation azimuth and the azimuth of the crack's axis of symmetry, such as |φ-φ syn |,φ syn The azimuth axis of the crack can be represented as the orientation of the crack's axis of symmetry. Based on this, without considering the orientation of the crack's axis of symmetry, a relative azimuth angle can be a seismic observation orientation. Alternatively, the orientation of the crack's axis of symmetry can also be represented as the orientation of the crack's normal direction.

[0022] Optionally, the electronic device may store the fracture dip angle of any region in its own storage space, in which case the electronic device can obtain the fracture dip angle from its own storage space; or, the electronic device may obtain a fracture dip angle download link to download the fracture dip angle according to the fracture dip angle download link, thereby obtaining the fracture dip angle, etc.; the embodiments of the present invention do not limit this. The fracture dip angle may be prior information (i.e., prior information on the fracture dip angle of an inclined fracture); for example, the fracture dip angle may be obtained through imaging logging data, or it may be obtained through core data observation, etc.

[0023] S102, retrieves the target azimuth elastic impedance data set from the five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data under each angle combination. An angle combination includes any incident angle from J incident angles and any relative azimuth angle from K relative azimuth angles, where J and K are both positive integers.

[0024] Optionally, the five-dimensional seismic data can be five-dimensional seismic data with partial angular superposition, such that all incident angles are superimposed into J incident angles, and all relative azimuths are superimposed into K relative azimuths. The five-dimensional seismic data can include seismic trace data for various angle combinations. The seismic trace data for one angle combination corresponds to the incident angle and relative azimuth in the corresponding angle combination. The seismic trace data for one angle combination can also be referred to as the seismic trace data for one incident angle and one relative azimuth. Based on this, the data for one incident angle and one relative azimuth can also be represented as the data for the angle combination formed by the corresponding incident angle and the corresponding relative azimuth, or it can be represented as the data corresponding to the corresponding incident angle and the corresponding relative azimuth.

[0025] Optionally, for the j-th incident angle among J incident angles and the k-th relative azimuth angle among K relative azimuth angles (i.e., for the angle combination formed by the j-th incident angle and the k-th relative azimuth angle, j∈[1,J], k∈[1,K]), the relationship between the seismic trace data (also referred to as seismic data) and the target azimuth elastic impedance data (also referred to as elastic impedance data or azimuth elastic impedance data) under the j-th incident angle and the k-th relative azimuth angle can be expressed in matrix form as shown in Equation 1.1: Formula 1.1 Among them, S PP G can represent seismic trace data at the j-th incident angle and the k-th relative azimuth angle. p X can represent the forward and inverse operators of azimuth elastic impedance at the j-th incident angle and the k-th relative azimuth angle, and X can represent the target azimuth elastic impedance data at the j-th incident angle and the k-th relative azimuth angle, as shown in Formula 1.2: Equation 1.2 Among them, s i It can represent the j-th incident angle (i.e. θ j and the kth relative azimuth angle (i.e. φ k In the seismic trace data of the i-th sampling point (where a seismic trace data can include i sampling points, i is a positive integer), W can represent the wavelet matrix (here, it can be the wavelet matrix at the j-th incident angle and the k-th relative azimuth angle), D0 can represent the difference operation matrix (i.e., the difference operator, here, it can be the difference operator corresponding to the seismic trace data at the j-th incident angle and the k-th relative azimuth angle), InAEI can represent the logarithm of the azimuth elastic impedance AEI, AEI i This can represent the value of the i-th sampling point in the target azimuth elastic impedance data at the j-th incident angle and the k-th relative azimuth angle (i.e., the azimuth elastic impedance of the i-th sampling point). Optionally, the difference operator can be set according to experience or actual needs, and this embodiment of the invention does not limit this; optionally, the difference operator can be used to calculate the difference between the parameters of the upper and lower medium model. Optionally, a target azimuth elastic impedance data can be an inverted azimuth elastic impedance data, which can be used to calculate the inversion target value in the subsequent inversion prediction process, so that the model simulation value approximates the corresponding inversion target value.

[0026] Optionally, in embodiments of the present invention, each seismic trace data in multiple seismic trace data can also be represented in S. PP In the middle, that is, at this time S PP It can include data from each seismic trace (i.e., the seismic data to be inverted), G pIt may include forward and inverse operators for azimuth elastic impedance under various angle combinations, where X can represent the target azimuth elastic impedance data under various angle combinations (i.e., the target azimuth elastic impedance data set), and then the target azimuth elastic impedance data under various angle combinations can be inverted in one step, and so on; the embodiments of the present invention do not limit this.

[0027] Furthermore, given the initial model X of the azimuth elastic impedance pri (For example, it may include the initial model of azimuth elastic impedance under the j-th incident angle and the k-th relative azimuth angle, or the initial model of azimuth elastic impedance under various angle combinations.) After that, the damped least squares inversion algorithm based on model constraints can be used to solve Equation 1.1 (also known as the seismic elastic relationship function), thereby obtaining the target azimuth elastic impedance data (such as the target azimuth elastic impedance data under the j-th incident angle and the k-th relative azimuth angle, or the target azimuth elastic impedance data set); for example, as shown in Equation 1.3, the target azimuth elastic impedance data can be expressed as: Equation 1.3 Where σ can represent the damping factor related to the signal-to-noise ratio, C m This can represent the azimuth elastic impedance covariance matrix (i.e., the model parameter covariance matrix). Optionally, the damping factor can be set according to experience or actual needs, and this embodiment of the invention does not limit this. It should be noted that this embodiment of the invention does not limit the specific method for determining the model parameter covariance matrix; for example, a well logging statistical estimation algorithm can be used to determine the model parameter covariance matrix, such as C. m =P T In P / V, V represents the total number of logging samples in the target logging data (i.e., the target logging data may include V logging samples), and P represents the mean-reduced logging parameter sample matrix. For example, assuming each row of the target logging data represents one logging sample, the mean of each column of the target logging data can be calculated, and the mean of its column can be subtracted from each element of the target logging data to eliminate the overall parameter offset, resulting in a mean-reduced logging parameter sample matrix (also known as a zero-mean matrix), and so on.

[0028] Based on this, when retrieving the target azimuth elastic impedance data set from five-dimensional seismic data, the electronic device can determine multiple wavelet matrices from the five-dimensional seismic data. These multiple wavelet matrices include wavelet matrices for various angle combinations, meaning that the wavelet matrix for any angle combination can be determined from the seismic trace data for any angle combination. For example, for any angle combination among multiple angle combinations (including all angle combinations composed of J incident angles and K relative azimuth angles), the azimuth seismic wavelet (also called a seismic wavelet) for that angle combination can be extracted from the seismic trace data for that angle combination. The wavelet matrix for that angle combination can then be determined using the azimuth seismic wavelet, thus achieving the determination of the wavelet matrix for any angle combination from the seismic trace data for any angle combination. Optionally, a wavelet matrix can be a convolution matrix constructed from an azimuth seismic wavelet, or when an azimuth seismic wavelet is directly represented in matrix form, a wavelet matrix can be an azimuth seismic wavelet, etc. This embodiment of the invention does not limit this.

[0029] Then, based on multiple wavelet matrices, azimuth elastic impedance forward and inverse operators can be constructed (including azimuth elastic impedance forward and inverse operators under various angle combinations). For example, the azimuth elastic impedance forward and inverse operators under the j-th incident angle and the k-th relative azimuth angle can be constructed using the wavelet matrices under the j-th incident angle and the k-th relative azimuth angle (i.e., azimuth elastic impedance forward and inverse operators under various angle combinations can be constructed separately). Alternatively, multiple wavelet matrices can be used to construct an azimuth elastic impedance forward and inverse operator that includes azimuth elastic impedance forward and inverse operators under various angle combinations (i.e., the azimuth elastic impedance forward and inverse operators under various angle combinations can be treated as a whole). For example, multiple wavelet matrices and difference operators (which can include difference operators corresponding to each seismic trace data, or difference operators corresponding to each angle combination) can be used to construct the azimuth elastic impedance forward and inverse operators. Optionally, the difference operators corresponding to different angle combinations can be the same or different, and this embodiment of the invention does not limit this.

[0030] Furthermore, an azimuth elastic impedance covariance matrix can be constructed based on the target logging data; that is, the mean-reduced logging parameter sample matrix can be determined using the target logging data, and the azimuth elastic impedance covariance matrix can be determined using the mean-reduced logging parameter sample matrix. Optionally, the electronic device can obtain the target logging data from its own storage space, or download the target logging data according to a target logging data download link, etc.; this embodiment of the invention does not limit this.

[0031] Accordingly, the electronic device can determine the initial azimuth elastic impedance model (which may include the initial azimuth elastic impedance model under various angle combinations), and based on the initial azimuth elastic impedance model, the forward and inverse azimuth elastic impedance operators, and the azimuth elastic impedance covariance matrix, it can inversely derive the target azimuth elastic impedance data set from the five-dimensional seismic data. For example, based on the initial azimuth elastic impedance model, the forward and inverse azimuth elastic impedance operators, and the azimuth elastic impedance covariance matrix at the j-th incident angle and the k-th relative azimuth angle, the target azimuth elastic impedance data at the j-th incident angle and the k-th relative azimuth angle can be calculated, thereby realizing the inversion of the target azimuth elastic impedance data at the j-th incident angle and the k-th relative azimuth angle from the seismic trace data at the j-th incident angle and the k-th relative azimuth angle. Alternatively, it can inversely derive the target azimuth elastic impedance data set in one step, and so on. Optionally, the initial azimuth elastic impedance model under various angle combinations can be set according to experience or actual needs, or it can be determined based on well logging data, etc.; this embodiment of the invention does not limit this. Optionally, the initial model of azimuth elastic impedance under different angle combinations can be the same or different, and this embodiment of the invention does not limit this. Optionally, the target logging data can be logging data without angles, that is, the azimuth elastic impedance covariance matrix under different angle combinations can be the same.

[0032] As can be seen, the embodiments of the present invention can use Formula 1.3 to inversely derive the target azimuth elastic impedance data set from five-dimensional seismic data based on the initial model of azimuth elastic impedance, the forward and inverse azimuth elastic impedance operators, and the azimuth elastic impedance covariance matrix; in other words, the embodiments of the present invention can use Formula 1.3 to inversely derive the target azimuth elastic impedance data under different incident angles and different relative azimuth angles from five-dimensional seismic data.

[0033] S103. Based on the crack dip angle, determine the target azimuth elastic impedance difference function. The parameters to be inverted in the target azimuth elastic impedance difference function include crack weakness parameters. Based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, invert and predict the crack weakness parameters.

[0034] Optionally, the target azimuth elastic impedance difference function can be determined by the target monoclinic medium azimuth elastic impedance function, which supports its use for inversion prediction.

[0035] Optionally, the electronic device can also acquire the initial monoclinic medium azimuth seismic reflection coefficient equation. The parameters in the initial monoclinic medium azimuth seismic reflection coefficient equation (i.e., the parameters to be inverted) include the P-wave modulus of the isotropic background, the shear modulus of the isotropic background, the P-wave impedance, the anisotropy parameter, and the fracture weakness parameter. A monoclinic medium azimuth seismic reflection coefficient equation can be a azimuth seismic reflection coefficient equation describing the monoclinic medium (i.e., the monoclinic anisotropic medium). For example, the initial monoclinic medium azimuth seismic reflection coefficient equation can be as shown in Equation 1.4: Equation 1.4 Among them, R pp It can represent the azimuth earthquake reflection coefficient. θ It can represent the angle of incidence. φ The relative azimuth angle can be represented by ν, the crack dip angle by v, M and μ by M and μ by μ, respectively, representing the P-wave modulus and shear modulus of an isotropic background. Z represents the P-wave impedance (i.e., the P-wave impedance of an isotropic background, Z = sqrt(ν). ρ M)= ρ α, ρ α can represent density, and α can represent P-wave velocity; ε and δ can represent anisotropy parameters (i.e., Thomsen anisotropy parameters describing the vertical and horizontal isotropy of shale), ε can be used to measure P-wave anisotropy (also called P-wave anisotropy parameter), and δ has no clear physical meaning; δ N and δ T The normal and tangential weakness parameters of the crack can be represented separately. That is, the crack weakness parameters can include normal weakness parameters (also called normal crack weakness parameters) and tangential weakness parameters (also called tangential crack weakness parameters). Furthermore, the symbol "-" at the top of the denominator indicates that the model parameters of the upper and lower media are averaged (i.e., the average value between the model parameters of the upper and lower media). This can represent the difference between the model parameters of the upper and lower media (i.e., the difference between the model parameters of the upper and lower media). It can be seen that the initial monoclinic medium azimuth seismic reflection coefficient equation includes M, μ, Z, ε, δ, δ N and δ T These are the 7 model parameters.

[0036] Furthermore, the expressions for each coefficient in Formula 1.4 can be shown in Formula 1.5: Formula 1.5 Where g = μ / M.

[0037] For example, such as Figure 3As shown, the effects of variations in isotropic background parameters (i.e., M, μ, Z) and anisotropic parameters (i.e., ε and δ) on the reflection coefficient (i.e., azimuth seismic reflection coefficient) are related to the incident angle, and this variation is most significant at large incident angles. Figure 3 The curves in the graph can represent the changes of the parameter under different values, such as Figure 3 The curves in subgraph (a) can be represented from bottom to top (i.e., the direction indicated by the arrows) as follows: The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp With the change of incident angle; Figure 3 Subgraph (b) contains curves that, from top to bottom (i.e., the arrows indicate the direction), can be represented as follows: The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp With the change of incident angle; Figure 3 The curves included in subgraph (c) can be represented from bottom to top (i.e., the direction indicated by the arrows) as follows: The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp With the change of incident angle; Figure 3 The curves included in subgraph (d) can be represented from bottom to top (i.e., the direction indicated by the arrows) as follows: The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp With the change of incident angle; Figure 3 The curves included in subgraph (e) can be represented from bottom to top (i.e., the direction indicated by the arrows) as follows: The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp The reflection coefficient varies with the incident angle; the direction of the arrow can be used to indicate correlation, such as an upward arrow indicating a positive correlation between the reflection coefficient and the parameter, and a downward arrow indicating a negative correlation. Furthermore, as... Figure 4 As shown, the influence of crack weakness parameters on the reflection coefficient varies with the crack dip angle; among them, Figure 4 The neutron plot (a) represents a crack dip angle of v = 60°, and the five surfaces (from top to bottom) are respectively... δ N The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp With changes in the angle of incidence and relative azimuth; Figure 4 The neutron diagram (b) can represent a crack dip angle of v = 90°, with the five surfaces (from top to bottom) representing... δ N The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. ppWith changes in the angle of incidence and relative azimuth; Figure 4 The neutron diagram (c) represents a crack dip angle of v = 60°, and the five surfaces (from bottom to top) are respectively... δ T The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp With changes in the angle of incidence and relative azimuth; Figure 4 The neutron diagram (d) represents a crack dip angle of v = 90°, and the five surfaces (from bottom to top) are respectively... δ T The reflection coefficient R is calculated when the values ​​are -0.2, -0.1, 0, 0.1, and 0.2. pp It varies with the angle of incidence and the relative azimuth.

[0038] It is evident that the influence of fracture dip angle on the reflection coefficient is significant when predicting fractured reservoirs. Furthermore, the influence of anisotropic parameters and fracture weakness parameters on the reflection coefficient is far less than that of elastic parameters (such as M and Z) in an isotropic background (e.g., the curvature of curves related to parameters like M is greater, or the range of variation for different reflection coefficient curves corresponding to parameters like M is larger). Moreover, the influence of anisotropic parameters on the reflection coefficient is similar to that of the P-wave modulus (e.g., the range of variation is similar and both are positively correlated), indicating that it is difficult to simultaneously invert the seven unknown parameters in the equation. It should be noted that... Figure 3 The horizontal axis can be in degrees (°, i.e., the angle of incidence can be in degrees), and the vertical axis is dimensionless (i.e., the reflection coefficient is dimensionless, i.e., it has no unit); correspondingly, Figure 4 Intermediate angles (such as θ, φ The units for all of these (etc.) can be degrees, and the color bars (i.e., color scales) can represent the azimuth seismic reflection coefficient (dimensionless).

[0039] Based on this, the electronic device can determine the expressions for the anisotropic P-wave modulus and the anisotropic shear modulus, and adjust the initial monoclinic medium azimuth seismic reflection coefficient equation to a monoclinic medium azimuth seismic reflection coefficient equation based on these expressions. The parameters in the monoclinic medium azimuth seismic reflection coefficient equation include fracture weakness parameters, anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance. Specifically, the expression for the anisotropic P-wave modulus A can be Mexp(ε) (i.e., A = Mexp(ε)), and the expression for the anisotropic shear modulus B can be μexp[(ε-δ) / (4g)] (i.e., B = μexp[(ε-δ) / (4g)]). In this embodiment of the invention, sin 2 θ tan 2 θ =tan 2 θ-sin 2 θ ,as well as Therefore, the initial monoclinic medium azimuth seismic reflection coefficient equation can be substituted into it, and the monoclinic medium azimuth seismic reflection coefficient equation can be re-expressed based on the expressions for the anisotropic P-wave modulus and the anisotropic shear modulus. This allows the initial monoclinic medium azimuth seismic reflection coefficient equation to be adjusted into the monoclinic medium azimuth seismic reflection coefficient equation.

[0040] For example, the equation for the azimuth seismic reflection coefficient of a monoclinic medium can be shown in Equation 1.6: Equation 1.6 Among them, a A ( θ )=a M ( θ ), a B ( θ )=a μ ( θ As can be seen, the azimuth reflection coefficient equation for monoclinic media can include five model parameters: anisotropic P-wave modulus, anisotropic shear modulus, P-wave impedance, normal weakness parameter, and tangential weakness parameter, and each model parameter has a clear physical meaning. For example, such as... Figure 5 As shown, in this embodiment of the invention, Equation 1.6 (also known as the approximate equation of the invention), Equation 1.4 (also known as the original equation), and Pšen are used. The accuracy of the equations was verified by comparing them with the reflection coefficient equations for arbitrarily weakly anisotropic media proposed by Ik and Martins (which can be simply referred to as the P&M equations). It can be seen that under different fracture dip angles (e.g., 50° and 80°), the azimuth seismic reflection coefficient equation for monoclinic media proposed in this embodiment (also called the approximate equation for the azimuth seismic reflection coefficient of monoclinic media PP wave) matches well with the curves calculated by the original equations (i.e., the reflection coefficient curves). This indicates that using the azimuth seismic reflection coefficient equation for monoclinic media proposed in this embodiment to predict shale reservoirs with inclined fractures is reasonable. In other words, the azimuth seismic reflection coefficient equation for PP wave (i.e., the azimuth seismic reflection coefficient equation for monoclinic media) and the azimuth elastic impedance equation (i.e., the azimuth elastic impedance function of the target monoclinic medium) proposed in this embodiment are more suitable for describing shale reservoirs with inclined fractures. Furthermore, the equations derived in this embodiment contain only 5 model parameters, each with a clear physical meaning. Figure 5 The neutron plot (a) can represent the reflection coefficient curves (i.e., the reflection coefficient as a function of the incident angle) under various calculation methods when the crack dip angle is 50°. Figure 5Neutron plot (b) can represent the reflection coefficient curves under various calculation methods when the crack inclination angle is 80°. It should be noted that the unit of each angle (such as incident angle, relative azimuth angle, crack inclination angle, etc.) in this embodiment of the invention is degrees. Optionally, a function can also be called an equation or expression, etc., and this embodiment of the invention does not limit this.

[0041] Furthermore, the electronic equipment can determine the relationship between the seismic reflection coefficient (i.e., azimuth seismic reflection coefficient) and the elastic impedance (i.e., azimuth elastic impedance), and based on the equation and relationship of the azimuth seismic reflection coefficient of the monocline medium, construct the azimuth elastic impedance function of the target monocline medium (also called the azimuth elastic impedance equation of the target monocline medium). It should be understood that the seismic reflection coefficient reflects the information of the stratigraphic interface (i.e., the stratigraphic interface), while the elastic impedance reflects the information of the stratigraphic layers (i.e., the stratigraphic layers). The relationship between the seismic reflection coefficient and the elastic impedance can be expressed as shown in Equation 1.7: Equation 1.7 Accordingly, when constructing the azimuth elastic impedance function of the target monoclinic medium based on the equation and relational expression of the azimuth seismic reflection coefficient of the monoclinic medium, mathematical operations can be performed using Equations 1.6 and 1.7 to obtain the azimuth elastic impedance function of the target monoclinic medium (i.e., the normalized expression for the azimuth elastic impedance of the monoclinic medium). The azimuth elastic impedance function of the target monoclinic medium can be shown in Equation 1.8: Formula 1.8 Wherein, AEI represents azimuth elastic impedance (i.e., azimuth elastic impedance data), AEI0 = sqrt(ρ0M0). The subscript "0" in Formula 1.8 can represent reference values ​​for normalized model parameters, such as AEI0 representing the azimuth elastic impedance reference value, ρ0 representing the density reference value, M0 representing the P-wave modulus reference value for isotropic backgrounds, A0 representing the anisotropic P-wave modulus reference value, B0 representing the anisotropic shear modulus reference value, Z0 representing the P-wave impedance reference value, and so on. Optionally, each reference value can be set according to experience or actual needs, and this embodiment of the invention does not limit this; for example, it can be obtained by averaging well logging data, or by consulting literature, etc.

[0042] S104, determine the target incident angle related elastic impedance function. The parameters to be inverted in the target incident angle related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance. Based on the target azimuth elastic impedance data set, the target incident angle related elastic impedance function, and the inverted predicted crack weakness parameters, invert and predict the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance.

[0043] Based on this, embodiments of the present invention can use the seismic step-by-step inversion method (i.e., the step-by-step inversion method) to estimate five model parameters in the azimuth elastic impedance function of the target monoclinic medium.

[0044] This invention can acquire five-dimensional seismic data and fracture dip angles. The five-dimensional seismic data includes multiple seismic traces, each corresponding to an incident angle and a relative azimuth. A target azimuth elastic impedance data set is then derived from the five-dimensional seismic data. This set includes target azimuth elastic impedance data for each of multiple angle combinations, where each angle combination includes any one of J incident angles and any one of K relative azimuth angles. Based on the fracture dip angle, a target azimuth elastic impedance difference function is determined. The parameters to be inverted in this function include fracture weakness parameters. Finally, based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, fracture weakness parameters are predicted and inverted. Furthermore, the target incident angle-related elastic impedance function can be determined. The parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance. Based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted fracture weakness parameters, the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance are inverted and predicted. It can be seen that the embodiments of the present invention can determine the target azimuth elastic impedance difference function through fracture dip angle, etc., to invert fracture weakness parameters, and invert the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance through the target incident angle-related elastic impedance function. This fully considers the influence of tilted fractures and can effectively improve the accuracy of the inversion results. That is, the embodiments of the present invention can effectively improve the accuracy of the inversion results through step-by-step inversion, thereby effectively improving the accuracy of shale reservoir prediction, etc.

[0045] Based on the above description, this embodiment of the invention also proposes a more specific method for step-by-step inversion of five-dimensional seismic data of shale reservoirs. Accordingly, this method for step-by-step inversion of five-dimensional seismic data of shale reservoirs can be executed by the aforementioned electronic device (terminal or server); or, this method can be executed jointly by a terminal and a server. For ease of explanation, the following description will use the execution of this method for step-by-step inversion of five-dimensional seismic data of shale reservoirs by an electronic device as an example; please refer to [link to relevant documentation]. Figure 6 The step-by-step inversion method for five-dimensional seismic data of shale reservoirs may include the following steps S601-S605: S601, acquire five-dimensional seismic data and acquire fracture dip angle. The five-dimensional seismic data includes multiple seismic traces, and each seismic trace corresponds to an incident angle and a relative azimuth angle.

[0046] S602, the target azimuth elastic impedance data set is derived from five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data under each angle combination. An angle combination includes any incident angle from J incident angles and any relative azimuth angle from K relative azimuth angles, where J and K are both positive integers.

[0047] Since subsequent inversion predictions all use Newton's method for nonlinear inversion solutions, the principle of Newton's method is introduced here first, as follows: In geophysical exploration, the nonlinear functional relationship between observation data and model parameters can be represented by a matrix as shown in Equation 2.1: Equation 2.1 Where d represents the observed data vector, G represents the nonlinear forward modeling operator, and m represents the model parameter vector. Equation 2.1 can be solved using Newton's method as shown in Equation 2.2: Equation 2.2 in, β The step size for each iteration can be represented by m0, which can represent the initial value of the model parameters (such as the initial value of fracture weakness parameters, anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance, etc.) or the model parameters output from the previous iteration (i.e., the m obtained from the previous iteration can be used as m0 in the current iteration). Optionally, the initial value of the model parameters can be set according to experience or actual needs, and this embodiment of the invention does not limit this; for example, it can be determined by well logging data, or obtained by seismic interpretation of layer constraints, etc. Optionally, the step size for each iteration can be set according to experience or actual needs, and this embodiment of the invention does not limit this. Correspondingly, the perturbation update amount of the model parameters in each iteration m can be represented as shown in Formula 2.3: Equation 2.3 Where H can be the Hessian matrix. Correspondingly, H and Y can be represented as shown in Equation 2.4: Equation 2.4 in, d mod / m can represent the first derivative with respect to the model parameters (which can be substituted into m0 for calculation). 2 d mod / m 2 It can represent the second derivative with respect to the model parameters, d modThis can represent the model forward modeling composite data (i.e., model simulation data, which takes the current iteration stratigraphic model m0 as input, i.e., inputting m0 under the current iteration to obtain model simulation data, such as obtaining the azimuth elastic impedance difference value simulated by the model (which can be called the azimuth elastic impedance difference simulation value), etc.), d input It can represent the input measured observation data (such as the azimuth elastic impedance difference value obtained through the target azimuth elastic impedance data set (which can be called the azimuth elastic impedance difference target value) etc.).

[0048] S603, determine the first relative azimuth and the second relative azimuth from K relative azimuths, and determine the target incident angle from J incident angles.

[0049] Optionally, both the first relative azimuth and the second relative azimuth can be set according to experience or actual needs; or, the electronic device can randomly select any two different relative azimuths from the K relative azimuths to serve as the first relative azimuth and the second relative azimuth, respectively; or, the first relative azimuth can be randomly selected from the K relative azimuths, and a second relative azimuth can be determined from the K relative azimuths whose difference from the first relative azimuth is greater than a preset relative azimuth difference threshold, etc.; this embodiment of the invention does not limit this. Optionally, the preset relative azimuth difference threshold can be set according to experience or actual needs, and this embodiment of the invention does not limit this.

[0050] Optionally, when determining the target angle of incidence from the J angles of incidence, the maximum angle of incidence can be determined from the J angles of incidence and used as the target angle of incidence; alternatively, an angle of incidence can be randomly selected from all angles of incidence greater than a preset angle of incidence threshold from the J angles of incidence, and the selected angle of incidence can be used as the target angle of incidence, etc.; this embodiment of the invention does not limit this. Optionally, the preset angle of incidence threshold can be set according to experience or actual needs, and this embodiment of the invention does not limit this.

[0051] Optionally, the number of target incident angles can be one or more, and this embodiment of the invention does not limit this. For example, assuming the number of target incident angles is O (O is a positive integer), the top O largest incident angles among the J incident angles can all be used as target incident angles, or O incident angles can be randomly selected from all incident angles greater than a preset incident angle threshold among the J incident angles, and the selected incident angles can be used as target incident angles, etc. This embodiment of the invention does not limit this. For ease of explanation, the following description will use a single target incident angle as an example.

[0052] S604. Based on the first relative azimuth angle, the second relative azimuth angle, the target incident angle, and the crack dip angle, the target azimuth elastic impedance difference function is determined. The parameters to be inverted in the target azimuth elastic impedance difference function include the crack weakness parameters. Based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, the crack weakness parameters are inverted and predicted.

[0053] Optionally, when determining the target azimuth elastic impedance difference function based on the first relative azimuth angle, the second relative azimuth angle, the target incident angle, and the crack dip angle, the electronic device can determine the target monoclinic medium azimuth elastic impedance function. The parameters to be inverted in the target monoclinic medium azimuth elastic impedance function include the crack weakness parameter, the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance. The target incident angle, the first relative azimuth angle, and the crack dip angle can be substituted into the target monoclinic medium azimuth elastic impedance function to obtain the azimuth elastic impedance substitution formula under the first relative azimuth angle. The target incident angle, the second relative azimuth angle, and the crack dip angle can also be substituted into the target monoclinic medium azimuth elastic impedance function to obtain the azimuth elastic impedance substitution formula under the second relative azimuth angle. Then, the ratio between the azimuth elastic impedance substitution formula under the first relative azimuth angle and the azimuth elastic impedance substitution formula under the second relative azimuth angle can be calculated to obtain the target azimuth elastic impedance difference function. Based on this, embodiments of the present invention can subtract the azimuth elastic impedance data under the first relative azimuth angle and the second relative azimuth angle from the azimuth elastic impedance data under the second relative azimuth angle to obtain the target azimuth elastic impedance difference function. For example, the target azimuth elastic impedance difference function can be shown in Formula 2.5: Equation 2.5 in, φ p It can represent the first relative azimuth angle. φ q It can represent the second relative azimuth angle, DEI ( θ , φ p , φ q The azimuth elastic impedance difference of a monoclinic medium can be represented by (v), that is, the target azimuth elastic impedance difference function can be DEI( θ , φ p , φ q The expression for v); and, at this time, the incident angle can be the target incident angle.

[0054] Correspondingly, the target azimuth elastic impedance difference function can be expressed in matrix form as shown in Equation 2.6: Equation 2.6 Here, we take the example of the target incident angle being the j-th incident angle, and the seismic trace data containing i sampling points, as an example, d1=[DEI1( θ j , φ p , φ q ,v) , DEI i ( θ j , φ p , φ q ,v)] T m1=[δ N ,δ T ] T , , G1 can represent a nonlinear forward modeling operator related to the incident angle, azimuth angle, and crack dip angle. In this case, the model parameters may include normal crack weakness parameters and tangential crack weakness parameters, and the initial values ​​of the model parameters may include the initial values ​​of the normal crack weakness parameters and the initial values ​​of the tangential crack weakness parameters; DEI i ( θ j , φ p , φ q ,v) can represent the value of the i-th sampling point in the directional elastic impedance difference of a monoclinic medium.

[0055] Based on this, when inverting and predicting crack weakness parameters based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, the electronic device can use Newton's method to solve the target azimuth elastic impedance difference function based on the target azimuth elastic impedance data set, thereby obtaining the inverted and predicted crack weakness parameters (i.e., the inversion result of the crack weakness parameters). In other words, Newton's method can be used to invert and predict normal crack weakness parameters and tangential crack weakness parameters. For example, the first and second derivatives of the target azimuth elastic impedance difference function with respect to crack weakness parameters (including normal crack weakness parameters and tangential crack weakness parameters) can be calculated using formulas 2.5 and 2.6, respectively, as shown in formula 2.7. Equation 2.7 in, d1 / m1 can represent the first derivative with respect to the crack weakness parameter (which can be substituted into m0 for calculation), and can also represent the first and second derivatives with respect to the normal crack weakness parameter and the tangential crack weakness parameter, respectively; correspondingly, the incident angle in Formula 2.7 can be the target incident angle.

[0056] Based on this, by substituting Formula 2.7 into Formulas 2.2-2.4, the normal crack weakness parameters and tangential crack weakness parameters can be obtained by inversion using the Newton method, thus obtaining the inverted predicted crack weakness parameters.

[0057] Correspondingly, when using Newton's method to solve the target azimuth elastic impedance difference function based on the target azimuth elastic impedance data set to obtain the inverted predicted crack weakness parameters (i.e., when inverting and predicting crack weakness parameters based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function), the target azimuth elastic impedance data at the target incident angle and the first relative azimuth angle, as well as the target azimuth elastic impedance data at the target incident angle and the second relative azimuth angle, can be determined from the target azimuth elastic impedance data set. Then, the target azimuth elastic impedance difference target value (also called the azimuth elastic impedance difference observation value, used as d) can be calculated using the target azimuth elastic impedance data at the target incident angle and the first relative azimuth angle and the target azimuth elastic impedance data at the target incident angle and the second relative azimuth angle. input ); and can use the crack weakness parameters under the current iteration (i.e., m0, which may include the normal crack weakness parameters and tangential crack weakness parameters under the current iteration) to calculate the simulated value of the azimuth elastic impedance difference under the current iteration (e.g., the normal crack weakness parameters and tangential crack weakness parameters under the current iteration can be substituted into Formula 2.5 to calculate DEI, which can be used as the simulated value of the azimuth elastic impedance difference under the current iteration, and thus as d modFurthermore, the first and second derivatives of the target azimuth elastic impedance difference function with respect to the crack weakness parameter can be calculated (also known as the first and second derivatives under the current iteration, i.e., obtained by substituting the crack weakness parameter under the current iteration). Using these first and second derivatives, the simulated value of the azimuth elastic impedance difference, and the simulated value of the azimuth elastic impedance difference, the perturbation update amount of the crack weakness parameter is calculated. This perturbation update amount is then used to update the crack weakness parameter (i.e., the perturbation update amount of the crack weakness parameter is substituted into Formula 2.2 to calculate the crack weakness parameter under the next iteration). Based on this, the electronic device can use the crack weakness parameter under the next iteration as the crack weakness parameter under the current iteration, and iteratively execute the above process of using the crack weakness parameter under the current iteration to calculate the simulated value of the azimuth elastic impedance difference under the current iteration until the iteration stopping condition is met, thereby obtaining the inverted predicted crack weakness parameter (i.e., the crack weakness parameter at the time of reaching the iteration stopping condition is used as the inverted predicted crack weakness parameter). Optionally, the iteration stopping condition can be set according to experience or actual needs, and this embodiment of the invention does not limit this. For example, the iteration stopping condition can be that the difference between the simulated value of the azimuth elastic impedance difference and the target value of the azimuth elastic impedance difference in the current iteration is less than a preset difference threshold, or that the change between the simulated values ​​of the azimuth elastic impedance difference in consecutive iterations is less than a preset change threshold, and so on. Optionally, both the preset difference threshold and the preset change threshold can be set according to experience or actual needs, and this embodiment of the invention does not limit this.

[0058] Optionally, when calculating the target value of the difference in azimuth elastic impedance using the target incident angle and the first relative azimuth angle and the target incident angle and the second relative azimuth angle, the ratio between the target azimuth elastic impedance data under the target incident angle and the first relative azimuth angle and the target azimuth elastic impedance data under the target incident angle and the second relative azimuth angle can be used as the target value of the difference in azimuth elastic impedance.

[0059] S605, determine the target incident angle related elastic impedance function, the parameters to be inverted in the target incident angle related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus and P-wave impedance; and based on the target azimuth elastic impedance data set, the target incident angle related elastic impedance function and the inverted predicted crack weakness parameters, invert and predict the anisotropic P-wave modulus, anisotropic shear modulus and P-wave impedance.

[0060] Optionally, when determining the target incident angle-related elastic impedance function, the electronic device can determine the target monoclinic medium azimuth elastic impedance function, which includes the crack weakness parameter response impedance expression. This expression indicates the anisotropic impedance caused by the crack weakness parameter. The ratio between the target monoclinic medium azimuth elastic impedance function and the crack weakness parameter response impedance expression is calculated to obtain the target incident angle-related elastic impedance function. Based on this, embodiments of the present invention can use azimuth elastic impedance data minus the anisotropic impedance caused by the crack weakness parameter to obtain elastic impedance data that is only related to the incident angle information. The crack weakness parameter response impedance can be the anisotropic impedance caused by the crack weakness parameter.

[0061] For example, according to Equation 1.8, the elastic impedance data that is only related to the incident angle information can be expressed as shown in Equation 2.8, that is, the elastic impedance function related to the target incident angle can be expressed as shown in Equation 2.8: Equation 2.8 Where EI can represent the elastic impedance value related to the incident angle, and the denominator exp[.] in Formula 2.8 can represent the expression for the response impedance of the crack weakness parameter.

[0062] Optionally, Equation 2.8 can be expressed in matrix form as shown in Equation 2.9: Equation 2.9 Here, taking the j-th incident angle as an example, d2=[EI1(θ)] j ), …, EI i (θ j )] T m2=[A,B,Z] T A=[A1, …,A i ] T B = [B1, …, B i ] T Z = [Z1, …, Z i ] T G2 can represent a nonlinear forward modeling operator that depends only on the incident angle information, EI i (θ j ) can represent the value of the i-th sampling point in the incident angle-dependent elastic impedance.

[0063] According to Equations 2.8 and 2.9, the first and second derivatives of the elastic impedance function related to the target incident angle with respect to the model parameters (which can be anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance) can be calculated respectively, as shown in Equation 2.10: Equation 2.10 in, d2 / m2 represents the first derivative with respect to the anisotropic longitudinal wave modulus, the anisotropic shear modulus, and the longitudinal wave impedance. 2 EI( θ ) / A 2 It can represent the second derivative with respect to the anisotropic longitudinal wave modulus. 2 EI( θ ) / B 2 It can represent the second derivative with respect to the anisotropic shear modulus. 2 EI( θ ) / Z 2 It can represent the second derivative with respect to the longitudinal wave impedance, and so on.

[0064] Correspondingly, analogous to the method for solving crack weakness parameters, Equation 2.10 can be substituted into Equations 2.2-2.4 to use Newton's method to invert and solve for the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance. In other words, when inverting and predicting the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, Newton's method can be used to invert and solve for the target incident angle-related elastic impedance function based on the target azimuth elastic impedance data set and the inverted predicted crack weakness parameters, obtaining the inverted predicted anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance (i.e., the inversion results of the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance).

[0065] Based on this, when inverting and predicting the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance using the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, the target value of the incident angle-related elastic impedance at each of the J incident angles (also called the incident angle-related elastic impedance observation value, which serves as the d value at each incident angle) can be determined based on the target azimuth elastic impedance data set and the inverted predicted crack weakness parameters. input Based on the incident angle-related elastic impedance target value and the target incident angle-related elastic impedance function at each incident angle, the anisotropic longitudinal wave modulus, anisotropic shear modulus, and longitudinal wave impedance are inverted and predicted.

[0066] Optionally, when inverting and predicting the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance based on the target value of the incident angle-related elastic impedance and the target incident angle-related elastic impedance function at each incident angle, the model parameters under the current iteration (i.e., m0, which can include the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance under the current iteration) can be used to calculate the simulated value of the incident angle-related elastic impedance at each incident angle under the current iteration (e.g., any incident angle and the model parameters under the current iteration can be substituted into Formula 2.8 to calculate EI, which can be used as the simulated value of the incident angle-related elastic impedance at any incident angle under the current iteration, and thus as the d value at any incident angle under the current iteration). mod Correspondingly, the first and second derivatives of the target incident angle-related elastic impedance function with respect to the model parameters can be calculated, thereby obtaining the first and second derivatives of each incident angle under the current iteration (e.g., substituting any incident angle and the model parameters under the current iteration into the derivative expression to obtain the first and second derivatives of any incident angle under the current iteration). Then, using the first and second derivatives of each incident angle under the current iteration, the target value of the incident angle-related elastic impedance under each incident angle, and the simulated value of the incident angle-related elastic impedance under the current iteration, the model parameter perturbation update amount under each incident angle can be calculated, so as to determine the model parameters under the next iteration (i.e., the anisotropic longitudinal wave modulus, anisotropic shear modulus, and longitudinal wave impedance under the next iteration) based on the model parameter perturbation update amount under each incident angle. Based on this, the electronic device can use the model parameters of the next iteration as the model parameters of the current iteration, and iteratively execute the above process using the model parameters of the current iteration to calculate the simulated values ​​of the angle-related elastic impedance for each incident angle in the current iteration, until the iteration stopping condition is met, thereby obtaining the inverted predicted anisotropic longitudinal wave modulus, anisotropic shear modulus, and longitudinal wave impedance. Optionally, when calculating the model parameter perturbation update amount for each incident angle using the first and second derivatives of each incident angle in the current iteration, the target value of the angle-related elastic impedance for each incident angle, and the simulated value of the angle-related elastic impedance for each incident angle in the current iteration, for any incident angle among the J incident angles, the model parameter perturbation update amount for any incident angle can be calculated using formulas 2.3 and 2.4, using the first and second derivatives of any incident angle in the current iteration, the target value of the angle-related elastic impedance for any incident angle, and the simulated value of the angle-related elastic impedance for any incident angle in the current iteration.

[0067] Optionally, when determining the model parameters for the next iteration based on the model parameter perturbation update amounts at each incident angle, the average of the model parameter perturbation update amounts at each incident angle can be used as the integrated model parameter perturbation update amount. The integrated model parameter perturbation update amount and the model parameters for the current iteration are then used to determine the model parameters for the next iteration. (For example, the integrated model parameter perturbation update amount (i.e., the average of the model parameter perturbation update amounts at each incident angle) can be used as the integrated model parameter perturbation update amount. Substitute m) and the model parameters under the current iteration (i.e., m0) into Formula 2.2 to calculate the model parameters under the next iteration; or, the median of the model parameter perturbation update amounts under each incident angle can be used as the integrated model parameter perturbation update amount, and the integrated model parameter perturbation update amount and the model parameters under the current iteration can be used to determine the model parameters under the next iteration, etc.; the embodiments of the present invention do not limit this.

[0068] Optionally, when determining the target value of angle-related elastic impedance for each of the J incident angles based on the target azimuth elastic impedance data set and the inverted predicted crack weakness parameters, for the j-th incident angle, the target azimuth elastic impedance data for each angle combination (i.e., including the j-th incident angle) of the j-th incident angle can be determined from the target azimuth elastic impedance data set. The inverted predicted crack weakness parameters are then substituted into the crack weakness parameter response impedance expression to calculate the crack weakness parameter response impedance for each angle combination of the j-th incident angle. Then, the ratio between the target azimuth elastic impedance data and the crack weakness parameter response impedance for each angle combination of the j-th incident angle can be calculated to obtain the target value of angle-related elastic impedance for each angle combination of the j-th incident angle. Thus, the target value of angle-related elastic impedance for the j-th incident angle can be determined using the target value of angle-related elastic impedance for each angle combination of the j-th incident angle. Optionally, when determining the target value of the angle-related elastic impedance under each angle combination of the j-th incident angle, the mean of the target values ​​of the angle-related elastic impedance under each angle combination of the j-th incident angle can be used as the target value of the angle-related elastic impedance under the j-th incident angle; or, the median of the target values ​​of the angle-related elastic impedance under each angle combination of the j-th incident angle can be used as the target value of the angle-related elastic impedance under the j-th incident angle, and so on; the embodiments of the present invention do not limit this.

[0069] As can be seen, the embodiments of the present invention use azimuth elastic impedance data to subtract the anisotropic impedance caused by the crack weakness parameter to obtain elastic impedance data that is only related to the incident angle information, and use Newton's method to invert and predict the anisotropic longitudinal wave modulus, anisotropic shear modulus and P-wave impedance parameters.

[0070] In this embodiment of the invention, in order to further verify the feasibility and effectiveness of the step-by-step inversion method for five-dimensional seismic data of shale reservoirs proposed in this embodiment of the invention, synthetic seismic gather tests and actual field tests were conducted respectively.

[0071] On one hand, time-domain logging curves of a fractured shale reservoir were selected for synthetic testing. First, the fracture dip angle was fixed at 70°. A azimuth seismic gather was synthesized by convolving a Ricker wavelet with a dominant frequency of 24 Hz and the azimuth seismic reflection coefficient equation of the monoclinic medium (as shown in Equation 1.6). Then, Gaussian random noise with a signal-to-noise ratio of 10:1 was added to the synthesized azimuth seismic gather (also called the synthetic azimuth seismic gather). This allowed the noisy synthetic azimuth seismic gather to be used as seismic data for inversion. Figure 7 As shown, the incident angles of the synthesized azimuth seismic gathers are 8°, 18°, and 28°, respectively, and the relative azimuth angles (also referred to as azimuth angles) are 0°, 45°, and 90°, respectively. Accordingly, the step-by-step inversion method proposed in this embodiment of the invention can be used to estimate the model parameters in the equations; firstly, the target azimuth elastic impedance dataset can be estimated using a model-constrained damped least squares inversion algorithm. For example, as... Figure 8 As shown, with incident angle θ1=8°, incident angle θ2=18°, and relative azimuth angle... φ 1=0° and relative azimuth angle φ Taking each angle combination formed by 2=45° as an example, the target azimuth elastic impedance data (i.e., the azimuth elastic impedance data indicated by the inversion result) obtained by the embodiment of the present invention matches the true value well. Figure 8 The unit of AEI (i.e., target azimuth elastic impedance data, also known as azimuth elastic impedance inversion result) can be kg / m³ × m / s (i.e., kg / m³·m / s); where AEI(θ1, φ 1) Can represent the incident angle θ1 and the relative azimuth angle φ Target azimuth elastic impedance data under 1, AEI(θ2, φ 1) It can represent the incident angle θ2 and the relative azimuth angle. φ Target azimuth elastic impedance data under 1, AEI(θ1, φ 2) It can represent the incident angle θ1 and the relative azimuth angle. φ Target azimuth elastic impedance data under θ2, AEI(θ2, φ 2) It can represent the incident angle θ2 and the relative azimuth angle. φ2. Target azimuth elastic impedance data. Then, the difference data of azimuth elastic impedance corresponding to a large incident angle can be selected to predict the crack weakness parameters. For example, an incident angle of 28° can be selected as the target incident angle to determine the target azimuth elastic impedance difference function for inversion prediction of crack weakness parameters. Finally, elastic impedance data containing only incident angle information (i.e., the elastic impedance function related to the target incident angle) is used to estimate the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance. The inversion results (i.e., model parameter inversion results) are as follows: Figure 9 As shown, when the signal-to-noise ratio is 10:1, the inversion results of the method proposed in this embodiment of the invention match the true values ​​well, indicating that the method proposed in this embodiment of the invention is reasonable and feasible. Figure 9 The unit of anisotropic longitudinal wave modulus can be gigapascal, the unit of anisotropic shear modulus can be gigapascal, the unit of longitudinal wave impedance can be kilogram per cubic meter × meter per second (also known as meter per second kilogram per cubic meter), and both normal crack weakness (i.e., normal crack weakness parameter) and tangential crack weakness (i.e., tangential crack weakness parameter) are dimensionless.

[0072] On the other hand, field data collected from a fractured shale oil reservoir in eastern my country were used to further test the method proposed in this embodiment. The working area has a burial depth exceeding 3400 meters, and the target layer is extensively developed with organic-rich shale. The mineral composition is mainly quartz and carbonate minerals, followed by clay content. Based on mineral composition, it is mainly divided into organic-rich carbonate-rich shale facies and organic-rich mixed shale facies. Well logging and core testing data revealed that the target layer mainly develops inclined fractures, with the degree of fracture development gradually decreasing with increasing depth. Furthermore, imaging logging data showed that the fracture dip is 170° (i.e., the azimuth of the fracture axis), and the fracture dip angle is mainly concentrated around 50°; therefore, the target reservoir is equivalent to a monoclinic medium. The average seismic incident angles used for inversion were 8°, 18°, and 28°, with corresponding relative azimuth angles of 45°, 105°, and 165°. The model parameters were estimated using the step-by-step inversion method proposed in this embodiment of the invention. The inversion results (i.e., the inversion result profile of the fractured shale oil reservoir model parameters) are as follows: Figure 10 As shown; where, Figure 10 The neutron diagram (a) can represent the inversion results of the anisotropic P-wave modulus. Figure 10 Neutron plot (b) can represent the inversion results of anisotropic shear modulus. Figure 10 The neutron diagram (c) can represent the inversion result of the longitudinal wave impedance. Figure 10 The neutron diagram (d) can represent the inversion results of the normal crack weakness parameters. Figure 10The neutron plot (e) represents the inversion results of the tangential fracture weakness parameters. It can be seen that the inverted predicted anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance increase with depth in the longitudinal direction, and have good overall continuity in the transverse direction (i.e., the transverse direction is basically consistent), but retain the gradual characteristics of formation variation (i.e., there are local variations in the transverse direction). The inverted predicted fracture weakness parameters are mainly concentrated in the shallow to mid-depth region where the P-wave impedance is relatively low (i.e., the fracture weakness parameters are larger in the location of lower impedance), which is consistent with the prior knowledge of imaging logging. Figure 10 The curves in the figure represent the high-shear filtering display of the model parameters corresponding to the logging data at the well location. It can be seen that the inversion results of the above model parameters match well with the well curve (i.e., the curve peak to the right indicates a higher value, and the peak to the left indicates a lower value), further confirming the feasibility and effectiveness of the step-by-step inversion method proposed in this embodiment of the invention. Figure 10 The unit of anisotropic longitudinal wave modulus can be gigapascal, the unit of anisotropic shear modulus can be gigapascal, the unit of longitudinal wave impedance can be kilograms per cubic meter × meters per second, and both normal crack weakness and tangential crack weakness are dimensionless.

[0073] On the other hand, through cross-plot analysis of limited logging data from fractured shale oil reservoirs, this embodiment of the invention reveals that the cross-plot analysis results of the ratio (A / B) between the anisotropic P-wave modulus A and the anisotropic shear modulus B, and the anisotropic shear modulus (B), help distinguish between carbonate-rich shale facies and mixed shale facies in shale oil reservoirs. Figure 11 As shown. This means that, in practical applications, the lithofacies of shale reservoirs can be identified through well logging cross-analysis results of the inversion results obtained in this embodiment of the invention. Among them, Figure 11 GPa in this context can be represented as Gigapascal (or simply GPa).

[0074] This invention, after acquiring five-dimensional seismic data and the fracture dip angle, can invert a target azimuth elastic impedance data set from the five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data for each of multiple angle combinations. An angle combination includes any one of J incident angles and any one of K relative azimuth angles. A first relative azimuth angle and a second relative azimuth angle are determined from the K relative azimuth angles, and a target incident angle is determined from the J incident angles. Then, based on the first relative azimuth angle, the second relative azimuth angle, the target incident angle, and the fracture dip angle, a target azimuth elastic impedance difference function can be determined. The parameters to be inverted in the target azimuth elastic impedance difference function include fracture weakness parameters. Finally, based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, the fracture weakness parameters are predicted through inversion. Furthermore, the target incident angle-related elastic impedance function can be determined. The parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance. Based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted fracture weakness parameters, the anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance are inverted and predicted. It is evident that this embodiment of the invention considers the influence of fracture dip angle on the inversion results. It can directly calculate the second derivative with respect to the model parameters using the elastic impedance equations derived according to this embodiment (such as the target azimuth elastic impedance difference function, the target incident angle-related elastic impedance function, etc.), resulting in a concise calculation process. Moreover, this embodiment of the invention can reduce the number of model parameters in each inversion by introducing a seismic step-by-step inversion strategy, thereby reducing the ill-conditionedness and uncertainty of multi-parameter inversion of shale reservoirs with dipped fractures, and thus improving the reliability of the model parameter predictions in the equations, effectively improving the accuracy of the inversion results.

[0075] Based on the description of the relevant embodiments of the above-mentioned step-by-step inversion method for five-dimensional seismic data of shale reservoirs, this invention also proposes a step-by-step inversion device for five-dimensional seismic data of shale reservoirs. This device can be a computer program (including program code) running on an electronic device; such as... Figure 12 As shown, the shale reservoir five-dimensional seismic data step-by-step inversion device may include an acquisition unit 1201 and a processing unit 1202. This shale reservoir five-dimensional seismic data step-by-step inversion device can perform... Figure 1 or Figure 6 The step-by-step inversion method for 5D seismic data of shale reservoirs shown herein, i.e., the step-by-step inversion device for 5D seismic data of shale reservoirs can operate the above-mentioned unit: The acquisition unit 1201 is used to acquire five-dimensional seismic data and acquire fracture dip angle. The five-dimensional seismic data includes multiple seismic traces, and one seismic trace corresponds to one incident angle and one relative azimuth angle. Processing unit 1202 is used to retrieve a target azimuth elastic impedance data set from the five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data under each angle combination in multiple angle combinations. An angle combination includes any one of J incident angles and any one of K relative azimuth angles, where J and K are both positive integers. The processing unit 1202 is further configured to determine the target azimuth elastic impedance difference function based on the crack dip angle, wherein the parameters to be inverted in the target azimuth elastic impedance difference function include crack weakness parameters; and to invert and predict the crack weakness parameters based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function. The processing unit 1202 is further configured to determine the target incident angle-related elastic impedance function, wherein the parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance; and based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, to invert and predict the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance.

[0076] In one embodiment, when processing unit 1202 retrieves the target azimuth elastic impedance data set from the five-dimensional seismic data, it may specifically be used to: Multiple wavelet matrices are determined from the five-dimensional seismic data, including wavelet matrices under various angle combinations; Based on the multiple wavelet matrices, azimuth elastic impedance forward and inverse operators are constructed; and based on the target logging data, azimuth elastic impedance covariance matrix is ​​constructed. An initial model of azimuth elastic impedance is determined, and based on the initial model of azimuth elastic impedance, the forward and inverse azimuth elastic impedance operators, and the azimuth elastic impedance covariance matrix, the target azimuth elastic impedance data set is inverted from the five-dimensional seismic data.

[0077] In another embodiment, when determining the target azimuth elastic impedance difference function based on the crack dip angle, the processing unit 1202 may specifically be used for: The first relative azimuth angle and the second relative azimuth angle are determined from the K relative azimuth angles, and the target incident angle is determined from the J incident angles; Based on the first relative azimuth angle, the second relative azimuth angle, the target incident angle, and the crack inclination angle, the target azimuth elastic impedance difference function is determined.

[0078] In another embodiment, when determining the target azimuth elastic impedance difference function based on the first relative azimuth angle, the second relative azimuth angle, the target incident angle, and the crack dip angle, the processing unit 1202 may specifically be used for: The azimuth elastic impedance function of the target monoclinic medium is determined, wherein the parameters to be inverted in the azimuth elastic impedance function of the target monoclinic medium include the crack weakness parameter, the anisotropic longitudinal wave modulus, the anisotropic shear modulus, and the longitudinal wave impedance. Substituting the target incident angle, the first relative azimuth angle, and the crack inclination angle into the azimuth elastic impedance function of the target monoclinic medium, we obtain the azimuth elastic impedance substitution formula under the first relative azimuth angle. Substituting the target incident angle, the second relative azimuth angle, and the crack inclination angle into the azimuth elastic impedance function of the target monoclinic medium, we obtain the azimuth elastic impedance under the second relative azimuth angle; The ratio between the azimuth elastic impedance under the first relative azimuth angle and the azimuth elastic impedance under the second relative azimuth angle is calculated to obtain the target azimuth elastic impedance difference function.

[0079] In another embodiment, when determining the elastic impedance function related to the target incident angle, the processing unit 1202 may specifically be used for: The azimuth elastic impedance function of the target monoclinic medium is determined, wherein the azimuth elastic impedance function of the target monoclinic medium includes the crack weakness parameter response impedance expression, which is used to indicate the anisotropic impedance caused by the crack weakness parameter. The ratio between the azimuth elastic impedance function of the target monoclinic medium and the response impedance expression of the crack weakness parameter is calculated to obtain the target incident angle-related elastic impedance function.

[0080] In another embodiment, when processing unit 1202 inverts and predicts the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, it can specifically be used for: Based on the target orientation elastic impedance data set and the inverted predicted crack weakness parameters, the incident angle-related elastic impedance target values ​​for each of the J incident angles are determined. Based on the incident angle-related elastic impedance target value and the target incident angle-related elastic impedance function at each incident angle, the anisotropic longitudinal wave modulus, the anisotropic shear modulus, and the longitudinal wave impedance are inverted and predicted.

[0081] In another embodiment, the acquisition unit 1201 can also be used for: The initial monoclinic medium azimuth seismic reflection coefficient equation is obtained. The parameters in the initial monoclinic medium azimuth seismic reflection coefficient equation include the P-wave modulus of the isotropic background, the shear modulus of the isotropic background, the P-wave impedance, the anisotropic parameter, and the fracture weakness parameter. Processing unit 1202 can also be used for: The expressions for the anisotropic P-wave modulus and the anisotropic shear modulus are determined, and based on the expressions for the anisotropic P-wave modulus and the anisotropic shear modulus, the initial monoclinic medium azimuth seismic reflection coefficient equation is adjusted to a monoclinic medium azimuth seismic reflection coefficient equation. The parameters in the monoclinic medium azimuth seismic reflection coefficient equation include the fracture weakness parameter, the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance. The relationship between the seismic reflection coefficient and the elastic impedance is determined, and based on the azimuth seismic reflection coefficient equation of the monoclinic medium and the relationship expression, the azimuth elastic impedance function of the target monoclinic medium is constructed.

[0082] According to one embodiment of the present invention, Figure 12 Each unit in the illustrated shale reservoir five-dimensional seismic data step-by-step inversion device can be individually or entirely merged into one or more other units, or one or more of the units can be further divided into multiple functionally smaller units. This achieves the same operation without affecting the technical effects of the embodiments of the present invention. The above units are based on logical function division. In practical applications, the function of one unit can also be implemented by multiple units, or the function of multiple units can be implemented by one unit. In other embodiments of the present invention, any shale reservoir five-dimensional seismic data step-by-step inversion device may also include other units. In practical applications, these functions can also be implemented with the assistance of other units, and can be implemented collaboratively by multiple units.

[0083] According to another embodiment of the present invention, it is possible to perform operations such as those described above by running on a general-purpose electronic device, such as a computer, which includes processing elements and storage elements such as a central processing unit (CPU), random access memory (RAM), and read-only memory (ROM). Figure 1 or Figure 6 The computer program (including program code) involved in each step of the corresponding method shown, to construct such... Figure 12 The diagram illustrates a step-by-step inversion device for five-dimensional seismic data of shale reservoirs, and a method for implementing this embodiment of the invention. The computer program can be stored on, for example, a computer storage medium, loaded onto the aforementioned electronic device via the computer storage medium, and run therein.

[0084] Based on the description of the method and apparatus embodiments above, an exemplary embodiment of the present invention also provides an electronic device, including: at least one processor; and a memory communicatively connected to the at least one processor. The memory stores a computer program executable by the at least one processor, which, when executed by the at least one processor, causes the electronic device to perform the method according to an embodiment of the present invention.

[0085] An exemplary embodiment of the present invention also provides a non-transitory computer-readable storage medium storing a computer program, wherein the computer program, when executed by a computer's processor, is used to cause the computer to perform a method according to an embodiment of the present invention.

[0086] An exemplary embodiment of the present invention also provides a computer program product, including a computer program, wherein, when executed by a computer's processor, the computer program is used to cause the computer to perform a method according to an embodiment of the present invention.

[0087] refer to Figure 13 The present invention will now be described in the form of a structural block diagram of an electronic device 1300 that can serve as a server or client of the present invention, which is an example of a hardware device that can be applied to various aspects of the present invention. The electronic device is intended to represent various forms of digital electronic computer devices, such as laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. The electronic device can also represent various forms of mobile devices, such as personal digital processors, cellular phones, smartphones, wearable devices, and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely illustrative and are not intended to limit the implementation of the invention described and / or claimed herein.

[0088] like Figure 13 As shown, the electronic device 1300 includes a computing unit 1301, which can perform various appropriate actions and processes according to a computer program stored in a read-only memory (ROM) 1302 or a computer program loaded from a storage unit 1308 into a random access memory (RAM) 1303. The RAM 1303 may also store various programs and data required for the operation of the electronic device 1300. The computing unit 1301, ROM 1302, and RAM 1303 are interconnected via a bus 1304. An input / output (I / O) interface 1305 is also connected to the bus 1304.

[0089] Multiple components in electronic device 1300 are connected to I / O interface 1305, including: input unit 1306, output unit 1307, storage unit 1308, and communication unit 1309. Input unit 1306 can be any type of device capable of inputting information to electronic device 1300. Input unit 1306 can receive input digital or character information and generate key signal inputs related to user settings and / or function control of electronic device. Output unit 1307 can be any type of device capable of presenting information and may include, but is not limited to, a display, speaker, video / audio output terminal, vibrator, and / or printer. Storage unit 1308 may include, but is not limited to, disk and optical disk. Communication unit 1309 allows electronic device 1300 to exchange information / data with other devices through computer networks such as the Internet and / or various telecommunications networks, and may include, but is not limited to, modems, network cards, infrared communication devices, wireless communication transceivers, and / or chipsets, such as Bluetooth™ devices, WiFi devices, WiMax devices, cellular communication devices, and / or the like.

[0090] The computing unit 1301 can be a variety of general-purpose and / or special-purpose processing components with processing and computing capabilities. Some examples of the computing unit 1301 include, but are not limited to, a central processing unit (CPU), a graphics processing unit (GPU), various special-purpose artificial intelligence (AI) computing chips, various computing units running machine learning model algorithms, a digital signal processor (DSP), and any suitable processor, controller, microcontroller, etc. The computing unit 1301 performs the various methods and processes described above. For example, in some embodiments, the step-by-step inversion method for five-dimensional seismic data of shale reservoirs can be implemented as a computer software program tangibly contained in a machine-readable medium, such as storage unit 1308. In some embodiments, part or all of the computer program can be loaded and / or installed on electronic device 1300 via ROM 1302 and / or communication unit 1309. In some embodiments, the computing unit 1301 can be configured to perform the step-by-step inversion method for five-dimensional seismic data of shale reservoirs by any other suitable means (e.g., by means of firmware).

[0091] The program code used to implement the methods of the present invention can be written in any combination of one or more programming languages. This program code can be provided to a processor or controller of a general-purpose computer, special-purpose computer, or other programmable data processing device, such that when executed by the processor or controller, the program code causes the functions / operations specified in the flowcharts and / or block diagrams to be implemented. The program code can be executed entirely on the machine, partially on the machine, as a standalone software package partially on the machine and partially on a remote machine, or entirely on a remote machine or server.

[0092] In the context of this invention, a machine-readable medium can be a tangible medium that may contain or store a program for use by or in conjunction with an instruction execution system, apparatus, or device. A machine-readable medium can be a machine-readable signal medium or a machine-readable storage medium. Machine-readable media can include, but are not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatus, or devices, or any suitable combination of the foregoing. More specific examples of machine-readable storage media include electrical connections based on one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fibers, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination of the foregoing.

[0093] As used herein, the terms "machine-readable medium" and "computer-readable medium" refer to any computer program product, device, and / or apparatus (e.g., disk, optical disk, memory, programmable logic device (PLD)) for providing machine instructions and / or data to a programmable processor, including machine-readable media that receive machine instructions as machine-readable signals. The term "machine-readable signal" refers to any signal for providing machine instructions and / or data to a programmable processor.

[0094] To provide interaction with a user, the systems and techniques described herein can be implemented on a computer having: a display device for displaying information to the user (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor); and a keyboard and pointing device (e.g., a mouse or trackball) through which the user provides input to the computer. Other types of devices can also be used to provide interaction with the user; for example, feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and input from the user can be received in any form (including sound input, voice input, or tactile input).

[0095] The systems and technologies described herein can be implemented in computing systems that include backend components (e.g., as a data server), or computing systems that include middleware components (e.g., an application server), or computing systems that include frontend components (e.g., a user computer with a graphical user interface or web browser through which a user can interact with implementations of the systems and technologies described herein), or any combination of such backend, middleware, or frontend components. The components of the system can be interconnected via digital data communication of any form or medium (e.g., a communication network). Examples of communication networks include local area networks (LANs), wide area networks (WANs), and the Internet.

[0096] Furthermore, it should be understood that the above-disclosed embodiments are merely preferred embodiments of the present invention and should not be construed as limiting the scope of the present invention. Therefore, any equivalent variations made in accordance with the claims of the present invention are still within the scope of the present invention.

Claims

1. A method of stepwise inversion of five-dimensional seismic data for shale reservoirs, characterized in that, include: Five-dimensional seismic data and fracture dip angles are acquired. The five-dimensional seismic data includes multiple seismic traces, with each seismic trace corresponding to an incident angle and a relative azimuth angle. The target azimuth elastic impedance data set is derived from the five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data under each angle combination in multiple angle combinations. An angle combination includes any incident angle among J incident angles and any relative azimuth angle among K relative azimuth angles, where J and K are both positive integers. Based on the crack dip angle, the target azimuth elastic impedance difference function is determined, and the parameters to be inverted in the target azimuth elastic impedance difference function include crack weakness parameters; Based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function, the crack weakness parameters are inverted and predicted; The target incident angle-related elastic impedance function is determined, and the parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance; and based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance are inverted and predicted.

2. The method of claim 1, wherein, The target azimuth elastic impedance data set retrieved from the five-dimensional seismic data includes: Multiple wavelet matrices are determined from the five-dimensional seismic data, including wavelet matrices under various angle combinations; Based on the multiple wavelet matrices, azimuth elastic impedance forward and inverse operators are constructed; and based on the target logging data, azimuth elastic impedance covariance matrix is ​​constructed. An initial model of azimuth elastic impedance is determined, and based on the initial model of azimuth elastic impedance, the forward and inverse azimuth elastic impedance operators, and the azimuth elastic impedance covariance matrix, the target azimuth elastic impedance data set is inverted from the five-dimensional seismic data.

3. The method according to claim 1 or 2, characterized in that, The determination of the target azimuth elastic impedance difference function based on the crack dip angle includes: The first relative azimuth angle and the second relative azimuth angle are determined from the K relative azimuth angles, and the target incident angle is determined from the J incident angles; Based on the first relative azimuth angle, the second relative azimuth angle, the target incident angle, and the crack inclination angle, the target azimuth elastic impedance difference function is determined.

4. The method of claim 3, wherein, The step of determining the target azimuth elastic impedance difference function based on the first relative azimuth angle, the second relative azimuth angle, the target incident angle, and the crack dip angle includes: The azimuth elastic impedance function of the target monoclinic medium is determined, wherein the parameters to be inverted in the azimuth elastic impedance function of the target monoclinic medium include the crack weakness parameter, the anisotropic longitudinal wave modulus, the anisotropic shear modulus, and the longitudinal wave impedance. Substituting the target incident angle, the first relative azimuth angle, and the crack inclination angle into the azimuth elastic impedance function of the target monoclinic medium, we obtain the azimuth elastic impedance substitution formula under the first relative azimuth angle. Substituting the target incident angle, the second relative azimuth angle, and the crack inclination angle into the azimuth elastic impedance function of the target monoclinic medium, we obtain the azimuth elastic impedance under the second relative azimuth angle; The ratio between the azimuth elastic impedance under the first relative azimuth angle and the azimuth elastic impedance under the second relative azimuth angle is calculated to obtain the target azimuth elastic impedance difference function.

5. The method according to claim 1 or 2, characterized in that, The determination of the target incident angle-related elastic impedance function includes: The azimuth elastic impedance function of the target monoclinic medium is determined, wherein the azimuth elastic impedance function of the target monoclinic medium includes the crack weakness parameter response impedance expression, which is used to indicate the anisotropic impedance caused by the crack weakness parameter. The ratio between the azimuth elastic impedance function of the target monoclinic medium and the response impedance expression of the crack weakness parameter is calculated to obtain the target incident angle-related elastic impedance function.

6. The method of claim 1 or 2, wherein, The inversion prediction of the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters includes: Based on the target orientation elastic impedance data set and the inverted predicted crack weakness parameters, the incident angle-related elastic impedance target values ​​for each of the J incident angles are determined. Based on the incident angle-related elastic impedance target value and the target incident angle-related elastic impedance function at each incident angle, the anisotropic longitudinal wave modulus, the anisotropic shear modulus, and the longitudinal wave impedance are inverted and predicted.

7. The method according to claim 1 or 2, characterized in that, The method further includes: The initial monoclinic medium azimuth seismic reflection coefficient equation is obtained. The parameters in the initial monoclinic medium azimuth seismic reflection coefficient equation include the P-wave modulus of the isotropic background, the shear modulus of the isotropic background, the P-wave impedance, the anisotropic parameter, and the fracture weakness parameter. The expressions for the anisotropic P-wave modulus and the anisotropic shear modulus are determined, and based on the expressions for the anisotropic P-wave modulus and the anisotropic shear modulus, the initial monoclinic medium azimuth seismic reflection coefficient equation is adjusted to a monoclinic medium azimuth seismic reflection coefficient equation. The parameters in the monoclinic medium azimuth seismic reflection coefficient equation include the fracture weakness parameter, the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance. The relationship between the seismic reflection coefficient and the elastic impedance is determined, and based on the azimuth seismic reflection coefficient equation of the monoclinic medium and the relationship expression, the azimuth elastic impedance function of the target monoclinic medium is constructed.

8. An apparatus for stepwise inversion of five-dimensional seismic data of shale reservoirs, characterized in that, The device includes: The acquisition unit is used to acquire five-dimensional seismic data and to acquire fracture dip angle. The five-dimensional seismic data includes multiple seismic traces, and each seismic trace corresponds to an incident angle and a relative azimuth angle. The processing unit is used to invert the target azimuth elastic impedance data set from the five-dimensional seismic data. The target azimuth elastic impedance data set includes target azimuth elastic impedance data under each angle combination in multiple angle combinations. An angle combination includes any incident angle among J incident angles and any relative azimuth angle among K relative azimuth angles, where J and K are both positive integers. The processing unit is further configured to determine the target azimuth elastic impedance difference function based on the crack dip angle, wherein the parameters to be inverted in the target azimuth elastic impedance difference function include crack weakness parameters; and to invert and predict the crack weakness parameters based on the target azimuth elastic impedance data set and the target azimuth elastic impedance difference function. The processing unit is further configured to determine the target incident angle-related elastic impedance function, wherein the parameters to be inverted in the target incident angle-related elastic impedance function include anisotropic P-wave modulus, anisotropic shear modulus, and P-wave impedance; and based on the target azimuth elastic impedance data set, the target incident angle-related elastic impedance function, and the inverted predicted crack weakness parameters, to invert and predict the anisotropic P-wave modulus, the anisotropic shear modulus, and the P-wave impedance.

9. An electronic device, comprising: include: processor; as well as Stored program memory, The program includes instructions that, when executed by the processor, cause the processor to perform the method according to any one of claims 1-7.

10. A non-transitory computer-readable storage medium storing computer instructions, wherein, The computer instructions are used to cause the computer to perform the method according to any one of claims 1-7.

Images