Joint inversion method and device, storage medium and electronic equipment

Abstract

The invention provides a joint inversion method and device, a storage medium and electronic equipment, and the method comprises the steps: obtaining initial parameter data of a target region, and obtaining PP wave seismic data and PS wave seismic data of the target region, the initial parameter data comprising upper-layer medium parameter data and lower-layer medium parameter data; constructing a positive operator based on the initial parameter data; constructing a target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data; determining a target inversion result of the to-be-inverted parameter based on the target inversion function; and based on the target inversion result of the to-be-inverted parameter, calculating a target inversion result of each target parameter in a target parameter set, the target parameter set comprising an anisotropic parameter. According to the embodiment of the invention, the accuracy of seismic inversion can be improved.

Description

Technical Field

[0001] The present invention relates to the field of oil and gas seismic exploration, and in particular to a joint inversion method, device, storage medium and electronic equipment. Background Art

[0002] Prestack seismic inversion focuses on the change of amplitude with offset or incident angle, which is a key technology for describing reservoir elasticity, physical properties and anisotropy. VTI medium (anisotropic medium), as the simplest anisotropic model, usually refers to thin interlayered media or horizontal fractures in a uniform background, which plays an important role in describing seismic wave propagation. Shale, as a typical VTI medium model, has obvious anisotropic characteristics, especially in the reflection and transmission response of seismic waves, that is, shale reservoirs usually show strong anisotropy. However, related technologies usually focus on the weak anisotropy hypothesis, and it is difficult to effectively predict anisotropic parameters from actual seismic data, resulting in low accuracy of seismic inversion. Based on this, there is currently no good solution to improve the accuracy of seismic inversion. Summary of the invention

[0003] In view of this, an embodiment of the present invention provides a joint inversion method, device, storage medium and electronic device to solve the problem of low accuracy of seismic inversion in related technologies; that is, an embodiment of the present invention can implement a PP wave and PS wave joint inversion method through PP wave seismic data and PS wave seismic data to predict the anisotropic parameters of the VTI medium, thereby obtaining a target inversion result of anisotropic parameters with higher accuracy, which can effectively improve the accuracy of seismic inversion.

[0004] According to one aspect of an embodiment of the present invention, a joint inversion method is provided, the method comprising: Acquiring initial parameter data of a target area, and acquiring PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data includes upper layer medium parameter data and lower layer medium parameter data; Based on the initial parameter data, a forward operator is constructed; and based on the forward operator, the PP wave seismic data and the PS wave seismic data, a target inversion function is constructed; Based on the target inversion function, determining a target inversion result of the parameter to be inverted; Based on the target inversion results of the parameters to be inverted, target inversion results of each target parameter in a target parameter set are calculated, and the target parameter set includes anisotropic parameters.

[0005] According to another aspect of an embodiment of the present invention, a joint inversion device is provided, the device comprising: An acquisition unit, used to acquire initial parameter data of a target area, and acquire PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data includes upper layer medium parameter data and lower layer medium parameter data; A processing unit, configured to construct a forward operator based on the initial parameter data; and to construct a target inversion function based on the forward operator, the PP wave seismic data, and the PS wave seismic data; The processing unit is further used to determine a target inversion result of the parameter to be inverted based on the target inversion function; The processing unit is further used to calculate the target inversion result of each target parameter in a target parameter set based on the target inversion result of the parameter to be inverted, and the target parameter set includes anisotropic parameters.

[0006] According to another aspect of an embodiment of the present invention, an electronic device is provided, comprising a processor and a memory storing a program, wherein the program comprises instructions, which, when executed by the processor, enable the processor to perform the above-mentioned method.

[0007] According to another aspect of an embodiment of the present invention, a non-transitory computer-readable storage medium storing computer instructions is provided, wherein the computer instructions are used to enable a computer to execute the above-mentioned method.

[0008] According to another aspect of an embodiment of the present invention, a computer program product is provided, comprising a computer program, wherein when the computer program is executed by a processor, it is used to enable a computer to perform the above-mentioned method.

[0009] After obtaining the initial parameter data of the target area, as well as the PP wave seismic data and the PS wave seismic data of the target area, the embodiment of the present invention can construct a forward operator based on the initial parameter data; and construct a target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data. Further, based on the target inversion function, the target inversion result of the parameter to be inverted can be determined; and based on the target inversion result of the parameter to be inverted, the target inversion result of each target parameter in the target parameter set is calculated, and the target parameter set includes anisotropic parameters. It can be seen that the embodiment of the present invention can implement a PP wave and PS wave joint inversion method through PP wave seismic data and PS wave seismic data to predict the anisotropic parameters of the VTI medium, thereby obtaining a target inversion result of anisotropic parameters with higher accuracy, which can effectively improve the accuracy of seismic inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[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 schematic flow chart of a joint inversion method according to an exemplary embodiment of the present invention is shown; Figure 2 A schematic flow chart of another joint inversion method according to an exemplary embodiment of the present invention is shown; Figure 3 A schematic diagram showing a PP wave reflection coefficient according to an exemplary embodiment of the present invention is shown; Figure 4 A schematic diagram showing a PS wave reflection coefficient according to an exemplary embodiment of the present invention is shown; Figure 5 A schematic diagram showing an inversion result according to an exemplary embodiment of the present invention; Figure 6 A schematic diagram showing another inversion result according to an exemplary embodiment of the present invention; Figure 7 A schematic diagram showing an inversion result of an anisotropic parameter according to an exemplary embodiment of the present invention; Figure 8 A schematic diagram showing an inversion result of another anisotropic parameter according to an exemplary embodiment of the present invention; Fig. 9 A schematic diagram showing an inversion result of another anisotropy parameter according to an exemplary embodiment of the present invention; Fig.10 A schematic diagram showing an inversion result of another anisotropy parameter according to an exemplary embodiment of the present invention; Fig.11 A schematic block diagram of a joint inversion device according to an exemplary embodiment of the present invention is shown; Fig.12 A block diagram of an exemplary electronic device that can be used to implement an embodiment of the present invention is shown. DETAILED DESCRIPTION

[0011] Embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. Although certain embodiments of the present invention are shown in the accompanying drawings, it should be understood that the present invention can be implemented in various forms and should not be construed as being limited to the embodiments described herein, which are instead provided for a more thorough and complete understanding of the present invention. It should be understood that the drawings and embodiments of the present invention are only for exemplary purposes and are not intended to limit the scope of protection of the present 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. In addition, 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 "including" and its variations used in this document are open inclusions, that is, "including but not limited to". The term "based on" means "based at least in part on". The term "one embodiment" means "at least one embodiment"; the term "another embodiment" means "at least one other embodiment"; the term "some embodiments" means "at least some embodiments". Relevant definitions of other terms will be given in the following description. It should be noted that the concepts of "first", "second", etc. mentioned in the present invention are only used to distinguish different devices, modules or units, and are not used to limit the order or interdependence of the functions performed by these devices, modules or units.

[0014] It should be noted that the modifications of "one" and "plurality" mentioned in the present invention are illustrative rather than restrictive, and those skilled in the art should understand that, unless otherwise clearly indicated in the context, it should be understood as "one or more".

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

[0016] It should be noted that the execution subject of the joint inversion method provided in the embodiment of the present invention may be one or more electronic devices, which is not limited by the present invention; wherein, the electronic device may be a terminal (i.e., a client) or a server. Then, when the execution subject includes multiple electronic devices, and the multiple electronic devices include at least one terminal and at least one server, the joint inversion method provided in the embodiment of the present invention may be jointly executed by the terminal and the server. Accordingly, the terminal mentioned here may include but is not limited to: smart phones, laptops, desktop computers, intelligent voice interaction devices, and the like. The server mentioned here may be an independent physical server, or a server cluster or distributed system composed of multiple physical servers, or a cloud server that provides cloud services, cloud databases, cloud computing (cloud computing), cloud functions, cloud storage, network services, cloud communications, middleware services, domain name services, security services, CDN (Content Delivery Network), and basic cloud computing services such as big data and artificial intelligence platforms, and the like.

[0017] Based on the above description, an embodiment of the present invention proposes a joint inversion method, which can be executed by the electronic device (terminal or server) mentioned above; or, the joint inversion method can be executed by the terminal and the server together. For the sake of convenience, the following description will be given by taking the electronic device executing the joint inversion method as an example; Figure 1 As shown, the joint inversion method may include the following steps S101-S104: S101, obtaining initial parameter data of the target area, and obtaining PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data includes upper layer medium parameter data and lower layer medium parameter data.

[0018] Optionally, the target area may be any area, which is not limited in the embodiment of the present invention.

[0019] In an embodiment of the present invention, the initial parameter data may include the initial value of each target parameter in the target parameter set; optionally, the target parameter set may include but is not limited to at least one of the following: shear wave velocity (i.e., qS wave (also referred to as S wave) velocity, qS wave is a shear wave), longitudinal wave velocity (i.e., qP wave (also referred to as P wave) velocity, qP wave is a longitudinal wave), density, and anisotropy parameters, etc., which are not limited in the embodiment of the present invention. Optionally, the anisotropy parameter may include but is not limited to at least one of the following: a first anisotropy parameter (also referred to as a longitudinal anisotropy parameter or a z-direction anisotropy parameter, etc.) and a second anisotropy parameter (also referred to as a transverse anisotropy parameter or an x-direction anisotropy parameter), etc., which are not limited in the embodiment of the present invention. Optionally, the initial value of a target parameter may include the upper initial value of the corresponding target parameter in the upper medium parameter data and / or the lower initial value in the lower medium parameter data. Optionally, the initial parameter data may also include the angle between the ray and the interface normal, etc., in which case the vertical slowness described below may be determined based on the initial parameter data; or, the initial parameter data may also include longitudinal wave vertical slowness, transverse wave vertical slowness and horizontal slowness, etc.; this is not limited to the embodiments of the present invention.

[0020] Among them, PP wave can refer to that the incident wave is P wave and the reflected wave is also P wave (PP wave can also be called reflected longitudinal wave), and PS wave can refer to that the incident wave is P wave and the reflected wave is S wave (PS wave can also be called converted reflected wave).

[0021] Optionally, the initial parameter data, PP wave seismic data and PS wave seismic data may be data in a well logging model, that is, may be generated in a model test, or may be data in actual seismic data, etc.; this is not limited in the embodiment of the present invention.

[0022] Optionally, the initial parameter data, PP wave seismic data and PS wave seismic data of the target area may be obtained in the following ways, but are not limited to: The first acquisition method: the electronic device's own storage space can store initial parameter data, PP wave seismic data and PS wave seismic data for each of multiple areas, and the multiple areas include a target area. In this case, the initial parameter data, PP wave seismic data and PS wave seismic data of the target area can be acquired from its own storage space.

[0023] The second acquisition method: the electronic device can acquire the well logging data of the target area, so as to acquire the initial parameter data, PP wave seismic data and PS wave seismic data of the target area through the well logging data.

[0024] The third acquisition method: the electronic device can display the parameter initial interface. In this case, the user can perform a parameter input operation on the parameter initial interface to input the initial parameter data, PP wave seismic data and PS wave seismic data of the target area through the parameter input operation. In this case, the electronic device can respond to the parameter input operation and acquire the initial parameter data, PP wave seismic data and PS wave seismic data of the target area, etc.

[0025] S102, constructing a forward operator based on the initial parameter data; and constructing a target inversion function based on the forward operator, the PP wave seismic data, and the PS wave seismic data.

[0026] Optionally, the forward operator may include a first forward operator and a second forward operator; optionally, the first forward operator may also be called a PP wave forward operator, and the second forward operator may also be called a PS wave forward operator.

[0027] In an embodiment of the present invention, when constructing a forward operator based on initial parameter data, the electronic device may determine a stiffness matrix based on the initial parameter data, the stiffness matrix may include multiple stiffness elements (that is, a stiffness matrix may include element values ​​of each stiffness element in multiple stiffness elements; optionally, the element value of a stiffness element may be an initial element value of the corresponding stiffness element), the initial element value of a stiffness element is determined based on the initial value of each target parameter in at least one target parameter in a target parameter set, the initial parameter data include the initial values ​​of each target parameter in the target parameter set, and the initial value of a target parameter includes an upper initial value of the corresponding target parameter in the upper medium parameter data and / or a lower initial value in the lower medium parameter data. Then, based on the stiffness matrix, multiple PP wave parameter perturbation coefficients can be calculated; and based on the stiffness matrix, multiple PS wave parameter perturbation coefficients can be calculated; wherein, multiple PP wave parameter perturbation coefficients include PP wave parameter perturbation coefficients corresponding to each stiffness element in multiple stiffness elements, and multiple PS wave parameter perturbation coefficients include PS wave parameter perturbation coefficients corresponding to each stiffness element; optionally, multiple PP wave parameter perturbation coefficients can also include PP wave density parameter perturbation coefficients (i.e., the fifth PP wave parameter perturbation coefficient described below), and multiple PS wave parameter perturbation coefficients can also include PS wave density parameter perturbation coefficients (i.e., the fifth PS wave parameter perturbation coefficient described below). Based on this, the electronic device can construct a first forward operator based on multiple PP wave parameter perturbation coefficients; and construct a second forward operator based on multiple PS wave parameter perturbation coefficients, so as to realize the construction of the forward operator.

[0028] It should be noted that the elastic impedance tensor can accurately characterize the reflection and transmission coefficient at any incident angle and any elastic contrast interface. The elastic impedance tensor of the VTI medium can be expressed as shown in Formula 1.1: Formula 1.1 Where E represents the elastic impedance tensor (such as E P It can be expressed as the longitudinal elastic impedance tensor, E S The subscripts P and S can represent the qP wave and qS wave, respectively, and p and q can represent the horizontal slowness and vertical slowness, respectively (such as q P It can represent the vertical slowness of longitudinal waves, q S can represent the vertical slowness of shear waves), C can represent the stiffness matrix, C 55 etc. can represent the stiffness elements of the stiffness matrix.

[0029] In an embodiment of the present invention, the expression of each stiffness element in the stiffness matrix may be as shown in Formula 1.2, that is, the electronic device may use Formula 1.2 to calculate multiple stiffness elements in the stiffness matrix, that is, calculate the element value of each stiffness element in the multiple stiffness elements: Formula 1.2 Among them, V P It can represent the longitudinal wave velocity, V S It can represent the shear wave velocity, ρ can represent the density, δ can represent the longitudinal anisotropy parameter (i.e. the first anisotropy parameter), ε can represent the transverse anisotropy parameter (i.e. the second anisotropy parameter), and γ can represent the anisotropy parameter in the y direction (i.e. the third anisotropy parameter); f=1-C 55 / C 33 .

[0030] Based on this, the electronic device can determine the stiffness matrix based on the initial parameter data through formula 1.2, that is, determine the element value of each stiffness element; wherein the upper medium parameter data can be the parameter data of the upper medium, and the lower medium parameter data can be the parameter data of the lower medium. Optionally, when a medium is an isotropic medium, the parameter data of the medium may include the values ​​of the longitudinal wave velocity, the transverse wave velocity and the density; when a medium is a VTI medium, the parameter data of the medium may include the values ​​of each target parameter; optionally, when the parameter data of a medium does not include the value of the anisotropic parameter (that is, a medium is an isotropic medium), the value of the anisotropic parameter of the medium may be 0.

[0031] Based on this, the electronic device can determine the upper initial value (i.e., the upper medium initial value) and the lower initial value (i.e., the lower medium initial value) of each target parameter based on the initial parameter data, so as to determine the initial value of each target parameter, and perform mean operation on the upper initial value and the lower initial value of each target parameter, respectively, to obtain the parameter value of each target parameter (i.e., the parameter value of a target parameter can be initialized by the mean between the upper initial value and the lower initial value of the corresponding target parameter, and the parameter value of a target parameter can also be called the initial parameter value at this time), and use the parameter value of each target parameter to determine the initial element value of each stiffness element in multiple stiffness elements (i.e., the initial parameter value of each target parameter can be used to determine the initial element value of each stiffness element at this time), so as to determine the stiffness matrix. Exemplarily, the element value (such as the initial element value) of each stiffness element can be determined by using the parameter value of each target parameter through formula 1.2, and the parameter value of each target parameter can be substituted into formula 1.2 to determine the element value of each stiffness element. Optionally, the plurality of stiffness elements may include but are not limited to: a first stiffness element (ie, C 11 ), the second stiffness element (i.e. C 13 ), the third stiffness element (i.e. C 33 ) and the fourth stiffness element (i.e. C 55 ), etc., which is not limited in the embodiments of the present invention.

[0032] Further, when calculating multiple PP wave parameter disturbance coefficients based on the stiffness matrix, the electronic device may determine multiple intermediate coefficients based on the stiffness matrix; and calculate the PP wave diagonal parameters based on the multiple intermediate coefficients, the longitudinal wave polarization direction coefficient, and the longitudinal wave vertical slowness; based on this, multiple PP wave parameter disturbance coefficients may be calculated based on the PP wave diagonal parameters, the stiffness matrix, and the multiple intermediate coefficients. Correspondingly, when calculating multiple PS wave parameter disturbance coefficients based on the stiffness matrix, the PS wave diagonal parameters may be calculated based on the multiple intermediate coefficients, the transverse wave polarization direction coefficient, and the transverse wave vertical slowness; and multiple PS wave parameter disturbance coefficients may be calculated based on the PS wave diagonal parameters, the stiffness matrix, and the multiple intermediate coefficients.

[0033] It should be noted that in anisotropic media, the phase velocity and polarization direction change with the propagation direction. The polarization direction of VTI media can be expressed as shown in formula 1.3: Formula 1.3 Where U can represent the polarization direction (such as U P It can represent the polarization direction of longitudinal wave, U S The superscripts x and z can represent the horizontal and vertical downward directions, and G can represent the coefficient of the polarization direction (such as G P It can represent the longitudinal wave polarization direction coefficient, G SIt can represent the transverse wave polarization direction coefficient); it should be understood that the polarization direction x component of the qP wave and the qS wave are both p, that is, the same. Among them, the longitudinal wave vertical slowness can be calculated by the longitudinal wave velocity (such as using the parameter value of the longitudinal wave velocity to calculate the longitudinal wave vertical slowness (such as the ratio between the cosine value of the parameter value of the longitudinal wave incident angle and the parameter value of the longitudinal wave velocity can be used as the longitudinal wave vertical slowness to initialize the longitudinal wave vertical slowness), and the transverse wave vertical slowness can be calculated by the transverse wave velocity (such as using the parameter value of the transverse wave velocity to calculate the transverse wave vertical slowness (such as the ratio between the cosine value of the parameter value of the transverse wave incident angle and the parameter value of the transverse wave velocity can be used as the transverse wave vertical slowness to initialize the transverse wave vertical slowness); accordingly, the horizontal slowness can be determined by the parameter value of the longitudinal wave velocity or the parameter value of the transverse wave velocity to initialize the horizontal slowness, etc. Based on this, the embodiment of the present invention can initialize the stiffness matrix, the longitudinal wave vertical slowness, the transverse wave vertical slowness and the horizontal slowness, etc. through the parameter values ​​of each target parameter (that is, through the initial parameter data). Correspondingly, the electronic device can calculate the longitudinal wave polarization direction coefficient based on the parameter value of density, horizontal slowness, longitudinal wave vertical slowness and stiffness matrix (i.e., the element value of at least one stiffness element among multiple stiffness elements), and can calculate the transverse wave polarization direction coefficient based on the parameter value of density, horizontal slowness, transverse wave vertical slowness and stiffness matrix; that is, the electronic device can use the calculation formula of the polarization direction coefficient in formula 1.3 to calculate the longitudinal wave polarization direction coefficient and the transverse wave polarization direction coefficient.

[0034] Based on this, according to the elastic impedance tensor and polarization direction of the VTI medium, the propagation matrix can be expressed as shown in Formula 1.4: Formula 1.4 Under the boundary conditions of displacement continuity and stress continuity, the exact reflection coefficient can be expressed by the polarization direction and propagation matrix as shown in Formula 1.5: Formula 1.5 Where R = [R PP R SP ;R PS R SS ] is the reflection coefficient composed of qPqP (i.e. incident P wave reflects P wave), qSqP (i.e. incident S wave reflects P wave), qPqS (i.e. incident P wave reflects S wave) and qSqS (i.e. incident S wave reflects S wave). The superscripts (1) and (2) represent the upper medium and the lower medium respectively, and the subscripts "+" and "-" represent the downgoing wave and the upgoing wave respectively. Considering that the vertical slowness of the upgoing and downgoing waves is opposite, according to formulas 1.3 and 1.4, the relationship between the polarization direction of the upgoing and downgoing waves and the propagation matrix of the upgoing and downgoing waves can be shown as formula 1.6: Formula 1.6 Based on this, substituting formula 1.6 into formula 1.5, the exact reflection coefficient can be simplified to formula 1.7: Formula 1.7 Among them, matrices A and B can be expressed as shown in formula 1.8: Formula 1.8 Wherein, the superscript (m) (m=1,2) represents the upper and lower media, and I=[1 0;0 1] is the unit matrix. Based on this, the embodiment of the present invention can achieve accurate reflection coefficients for characterizing uplink and downlink waves through formula 1.7. Optionally, according to the polarization direction (i.e., formula 1.3) and the propagation matrix (i.e., formula 1.4), the specific expressions of matrices A and B can be expressed as shown in formula 1.9: Formula 1.9 Based on this, the embodiment of the present invention constructs an accurate expression of the VTI medium reflection coefficient represented by two matrices A and B (with a scale of 2×2). Accordingly, according to the specific expressions of matrices A and B (i.e., formula 1.9), formula 1.10 can be obtained: Formula 1.10 Where Λ=diag([Λ P ,Λ S ]) is a diagonal matrix, Λ P (i.e., PP wave diagonal parameters) can be calculated by substituting the first column in matrix A and the first column in matrix B into formula 1.10, Λ S (i.e., PS wave diagonal parameters) can be the result of substituting the second column in matrix A and the second column in matrix B into formula 1.10 for calculation. According to formula 1.10, the inversion operation of matrices A and B can be converted into a transposition operation, and under the weak contrast assumption (i.e., the weak contrast assumption between the upper medium and the lower medium, such as A (1) and A (2) Weak contrast, B (1) and B (2) Weak contrast, etc.) can have A=A (1) / 2+A (2) / 2, A=(A (2) -A (1) ) / 2, B=B (1) / 2+B (2) / 2, B=(B (2) -B (1) ) / 2; Based on this, the reflection coefficient can be simplified as shown in formula 1.11: Formula 1.11 It should be understood that rock physics experiments have shown that shale reservoirs usually exhibit strong anisotropy, and the applicability of conventional weak anisotropy approximation methods is limited. Based on this, the embodiment of the present invention can adopt a weak contrast hypothesis, that is, the elastic parameters and anisotropy parameters of the upper and lower layers are slightly different, which can improve the accuracy of the reflection coefficient at large angles in a medium with strong anisotropy. That is, the conventional weak anisotropy hypothesis eliminates the high-order terms of the anisotropy parameters, or eliminates the high-order terms of the anisotropy parameters that are not much different, such as retaining only the first-order terms, which causes the accuracy to be impaired, especially at large angles. The strong anisotropy parameter hypothesis mentioned in the embodiment of the present invention can effectively improve the accuracy, especially improve the accuracy of the reflection coefficient at large angles in a medium with strong anisotropy. In addition, formula 1.11 gives a general approximation of the reflection coefficient that can be applied to a strongly anisotropic VTI medium, which can be concisely represented by two 2×2 matrices and a transpose operation, thereby effectively simplifying the calculation and improving the operating efficiency.

[0035] In addition, the traditional reflection coefficient approximation equation is usually characterized by elastic parameters and anisotropic parameter disturbances, and there is a problem of large differences in contribution between the parameters, which seriously affects the stability of pre-stack seismic inversion, especially the anisotropic parameters; the general approximate equation of the reflection coefficient derived in the embodiment of the present invention (i.e., Formula 1.11) can be characterized by a stiffness matrix, which overcomes the problem of large differences in contribution of the parameters; that is, the conventional method uses five parameters for calculation, and no matter how strong the anisotropy is (i.e., no matter how large the anisotropy parameter is), it will not be very large, while the wave velocity (such as longitudinal wave velocity and shear wave velocity) is very large, and the contribution to the reflection coefficient is much greater. Parameters with large contributions are easier to invert, while parameters with small contributions are difficult to invert, resulting in difficult inversion of the anisotropic parameters and low accuracy. The embodiment of the present invention binds the anisotropic parameters and wave velocities together through the stiffness matrix, which can effectively reduce the difference, thereby facilitating seismic inversion and effectively improving the inversion accuracy of the anisotropic parameters.

[0036] Furthermore, the embodiment of the present invention can derive the approximation of the weak contrast reflection coefficients of the PP wave and PS wave in the VTI medium based on the general approximation of the reflection coefficient. Based on this, formula 1.9 can be substituted into formula 1.11 to obtain the reflection coefficient equations of the PP wave and PS wave in the VTI medium; wherein the PP wave reflection coefficient can be expressed as shown in formula 1.12: Formula 1.12 Among them, the coefficient (i.e., the intermediate coefficient) Γ 11 =C 11 / C 33 =1+2ε,Γ 55 =C 55 / C 33 =L 2 , Γ ρ =ρ / C 33 p 2 , Γ33 =C 33 / C 33 , Γ 13 = C 13 / C 33 , ; optional, Can be the difference in stiffness elements between the upper and lower layers of media (Y can be C 11 , C 13 , C 33 and C 55 ), such as using the upper initial value of each target parameter to determine the upper element value of each stiffness element, and using the lower initial value of each target parameter to determine the lower element value of each stiffness element, so that the difference between the lower element value and the upper element value of each stiffness element is used as the stiffness element difference of the corresponding stiffness element (that is, the stiffness element difference between the upper and lower media of the corresponding stiffness element). Correspondingly, It can be the difference between the density of the lower layer medium (i.e., the lower layer initial value of the density) and the density of the upper layer medium, etc.; Optionally, L can be the ratio of the transverse and longitudinal wave velocities. Based on this, the electronic device can determine multiple intermediate coefficients based on the stiffness matrix, and the multiple intermediate coefficients may include Γ 11 (also called the first intermediate coefficient), Γ 13 (also called the second intermediate coefficient), Γ 33 (also called the third intermediate coefficient), Γ 55 (also called the fourth intermediate coefficient) and Γ ρ (Also called the intermediate density coefficient). P It can be shown as formula 1.13: Formula 1.13 Based on this, the electronic device can calculate the PP wave diagonal parameter based on multiple intermediate coefficients, the longitudinal wave polarization direction coefficient and the longitudinal wave vertical slowness through formula 1.13. Correspondingly, ξ P (also known as the first reflection coefficient calculation parameter) can be shown as formula 1.14: Formula 1.14 Correspondingly, Q ψ (also known as the second reflection coefficient calculation parameter) can be shown as formula 1.15: Formula 1.15 In addition, the PS wave reflection coefficient of VTI medium can be expressed as shown in formula 1.16: Formula 1.16 Among them, Λ S (i.e., PS wave diagonal parameter) can be expressed as formula 1.17: Formula 1.17 Based on this, the electronic device can calculate the PS wave diagonal parameters based on multiple intermediate coefficients, the shear wave polarization direction coefficient and the shear wave vertical slowness through formula 1.17. (also known as the third reflection coefficient calculation parameter) can be shown as formula 1.18: Formula 1.18 Correspondingly, S (also known as the fourth reflection coefficient calculation parameter) can be shown as formula 1.19: Formula 1.19 It is worth noting that there are ratio terms of horizontal and vertical slowness in the coefficient term, and the horizontal and vertical slowness are the ratios of the sine and cosine of the incident angle to the phase velocity, respectively. Therefore, the coefficient term of the PP wave and PS wave reflection coefficient equation has nothing to do with the phase velocity, but is only related to the anisotropy parameter and the ratio of the transverse and longitudinal wave velocities, that is, the present invention implements the derivation of a PP wave and PS wave linear reflection coefficient.

[0037] Furthermore, considering that the calculation of anisotropic parameters using the stiffness matrix is ​​ill-posed, the embodiment of the present invention can construct a pre-stack seismic inversion target functional (i.e., a target inversion function) constrained by the relationship between the stiffness matrix and the anisotropic parameters. It should be understood that the stability of directly predicting anisotropic parameters from seismic data using conventional pre-stack seismic inversion methods is insufficient (this is because the contribution of anisotropic parameters is small); while the reflection coefficient equation derived by the embodiment of the present invention can be represented by the stiffness matrix, which overcomes the problem of large differences in the contribution of each parameter in the conventional reflection coefficient approximation, making it possible to predict the stiffness matrix from seismic data and indirectly and stably predict anisotropic parameters.

[0038] In the embodiment of the present invention, the PP wave and PS wave reflection coefficient approximation formula (i.e., formula 1.12 and formula 1.16) can be obtained by the stiffness coefficient (i.e., stiffness element) C 11 , C 13 , C 33 , C 55 And the disturbance term of ρ, these disturbance terms are written in the form of logarithmic differences, as shown in formula 1.20: Formula 1.20 Based on this, the forward model of PP wave and PS wave can be expressed as shown in formula 1.21: Formula 1.21 Among them, d P =[d P (θ P1 ),d P (θ P2 ),…,dP (θ Ph )] T It can be expressed by the incident angle θ P1 ,θ P2 ,…,θ Ph PP wave seismic data composed of h angle gathers, d S =[d S (θ S1 ),d P (θ S2 ),…,d P (θ Sh )] T It can be expressed by the incident angle θ S1 ,θ S2 ,…,θ Sh PS wave seismic data composed of h angle gathers, X=[X 11 ,X 13 ,X 33 ,X 55 ,ρ] T It can be the parameter to be inverted. The first forward operator G1 and the second forward operator G2 in the forward operator can be shown as formula 1.22: Formula 1.22 Among them, W P and W S may represent wavelets (such as amplitude phase wavelets) extracted from PP wave seismic data and PS wave seismic data, respectively, D is a first-order partial differential operator, v and v' may be PP wave parameter perturbation coefficients and PS wave parameter perturbation coefficients, respectively; optionally, the electronic device may use formula 1.23 to calculate the parameter perturbation coefficient (i.e., multiple PP wave parameter perturbation coefficients (including the first PP wave parameter perturbation coefficient v 11 , the second PP wave parameter perturbation coefficient v 13 、The third PP wave parameter disturbance coefficient v 33 、The fourth PP wave parameter perturbation coefficient v 55 And the fifth PP wave parameter perturbation coefficient v ρ ) and a plurality of PS wave parameter perturbation coefficients (which may include a first PS wave parameter perturbation coefficient , the second PS wave parameter perturbation coefficient , the third PS wave parameter perturbation coefficient , the fourth PS wave parameter perturbation coefficient And the fifth PS wave parameter perturbation coefficient )): Formula 1.23 Based on this, the electronic device can calculate multiple PP wave parameter disturbance coefficients based on the PP wave diagonal parameters, the stiffness matrix and the multiple intermediate coefficients through formula 1.23; specifically, the electronic device can calculate multiple second reflection coefficient calculation parameters (i.e., Q ψ ), so that the PP wave diagonal parameters, multiple intermediate coefficients and multiple second reflection coefficient calculation parameters can be used to calculate multiple PP wave parameter disturbance coefficients. Correspondingly, the electronic device can calculate multiple PS wave parameter disturbance coefficients based on the PS wave diagonal parameters, the stiffness matrix and multiple intermediate coefficients through formula 1.23; in this case, after determining the multiple intermediate coefficients based on the stiffness matrix, the electronic device can also calculate the third reflection coefficient calculation parameter and the fourth reflection coefficient calculation parameter based on the multiple intermediate coefficients and the transverse wave polarization direction coefficient through formulas 1.18 and 1.19, respectively, so that the PS wave diagonal parameters, multiple intermediate coefficients, multiple second reflection coefficient calculation parameters, the third reflection coefficient calculation parameter and the fourth reflection coefficient calculation parameter can be used to calculate multiple PS wave parameter disturbance coefficients.

[0039] In the embodiment of the present invention, when constructing the first forward operator based on multiple PP wave parameter disturbance coefficients, the first wavelet (i.e., W P ), and based on the first sub-wave and multiple PP wave parameter disturbance coefficients, construct a first forward operator, such as constructing the first forward operator by the calculation method of G1 in formula 1.22. Optionally, when extracting the first sub-wave, the seismic data of each angle channel in the PP wave seismic data can be superimposed to obtain the PP wave stacked seismic data, so as to extract the first sub-wave from the PP wave stacked seismic data; or, the sub-waves corresponding to the corresponding angle channel seismic data can be extracted from the seismic data of each angle channel in the PP wave seismic data, and the sub-waves corresponding to the seismic data of each angle channel can be averaged (i.e., the values ​​of each sampling point are averaged respectively) to obtain the first sub-wave, so as to extract the first sub-wave from the PP wave seismic data, and so on; the embodiments of the present invention are not limited to this.

[0040] Correspondingly, when constructing the second forward operator based on multiple PS wave parameter perturbation coefficients, the second wavelet (i.e., W S ), and based on the second wavelet and a plurality of PS wave parameter perturbation coefficients, construct a second forward operator, such as constructing the second forward operator by calculating G2 in formula 1.22. It should be noted that the implementation method of extracting the second wavelet may be the same as the implementation method of extracting the first wavelet, and the embodiments of the present invention will not be repeated here.

[0041] It should be understood that the coefficients in formula 1.23 are functions of the anisotropic parameters and the ratio of the transverse and longitudinal wave velocities, so the forward model is linear. In order to overcome the ill-posedness of using the predicted stiffness matrix to calculate the anisotropic parameters in prestack seismic inversion, the nonlinear relationship between the stiffness matrix and the anisotropic parameters can be considered as shown in formula 1.24: Formula 1.24 Among them, the function g in formula 1.24 can represent the function of calculating the anisotropic parameters using the stiffness matrix, that is, the terms on the right side of formula 1.24 can be used to calculate V in turn. P 、V S Based on this, according to the PP wave and PS wave forward model, the target functional (i.e., the target inversion function) constrained by the relationship between the stiffness matrix and the anisotropic parameters is constructed as shown in formula 1.25. That is to say, the electronic device can use formula 1.25 to construct the target inversion function based on the forward operator, PP wave seismic data, and PS wave seismic data: Formula 1.25 Among them, f(X) can represent the target inversion function, the first regularization term can be the PP wave seismic data error term, the second regularization term can be the PS wave seismic data error term, the third regularization term can be the stiffness matrix and anisotropic parameter relationship constraint term, the regularization parameter α is the constraint weight, the parameter α has little effect on the stiffness matrix inversion result, but has a greater impact on the final anisotropic parameter; the larger α is, the stronger the stability of the anisotropic parameter prediction. Optionally, the regularization parameter can be set according to experience, can be set according to actual needs, can be obtained according to experimental debugging, and so on; the embodiment of the present invention is not limited to this.

[0042] Based on this, when constructing a target inversion function based on the forward operator, PP wave seismic data and PS wave seismic data, the electronic device can determine the stiffness matrix and anisotropic parameter relationship constraints; and based on the first forward operator and PP wave seismic data, construct the PP wave seismic data error term; and based on the second forward operator and PS wave seismic data, construct the PS wave seismic data error term; then accordingly, the PP wave seismic data error term, the PS wave seismic data error term and the stiffness matrix and anisotropic parameter relationship constraints can be used to construct the target inversion function.

[0043] S103, determining a target inversion result of the parameter to be inverted based on the target inversion function.

[0044] In the embodiment of the present invention, the electronic device can invert the parameters to be inverted by using the target inversion function, and thus obtain the target inversion results of the parameters to be inverted.

[0045] S104, based on the target inversion result of the parameter to be inverted, calculating the target inversion result of each target parameter in the target parameter set, where the target parameter set includes anisotropic parameters.

[0046] In an embodiment of the present invention, the electronic device may determine the calculation method corresponding to each target parameter in the target parameter set based on the nonlinear relationship between the stiffness matrix and the target parameter set (as shown in Formula 1.24); and calculate the target inversion result of each target parameter according to the calculation method corresponding to each target parameter and the target inversion result of the parameter to be inverted. Optionally, at this time, the electronic device may determine the calculation method corresponding to the longitudinal wave velocity, the transverse wave velocity, the density, the first anisotropy parameter and the second anisotropy parameter, respectively, so as to calculate the target inversion results of the longitudinal wave velocity, the transverse wave velocity, the density, the first anisotropy parameter and the second anisotropy parameter respectively; that is, the target inversion result of the parameter to be inverted X may be substituted into Formula 1.24, so as to calculate the target inversion result of each target parameter according to the calculation method corresponding to each target parameter and the target inversion result of the parameter to be inverted. Exemplarily, the target inversion result of the longitudinal wave velocity may be equal to sqrt(exp(X 33 ) / exp(X ρ )), the target inversion result of shear wave velocity can be equal to sqrt(exp(X 55 ) / exp(X ρ )), the target inversion result of density can be equal to exp(X ρ ), and so on; that is, exp(X 11 ) as C 11 , exp(X 13 ) as C 13 , exp(X 33 ) as C 33 , exp(X 55 ) as C 55 , and exp(X ρ ) as ρ and substitute it into formula 1.24 to calculate the target inversion result of each target parameter respectively, that is, the target inversion result of longitudinal wave velocity, the target inversion result of shear wave velocity, the target inversion result of density, the target inversion result of the first anisotropy parameter and the target inversion result of the second anisotropy parameter can be calculated in turn.

[0047] After obtaining the initial parameter data of the target area, as well as the PP wave seismic data and the PS wave seismic data of the target area, the embodiment of the present invention can construct a forward operator based on the initial parameter data; and construct a target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data. Further, based on the target inversion function, the target inversion result of the parameter to be inverted can be determined; and based on the target inversion result of the parameter to be inverted, the target inversion result of each target parameter in the target parameter set is calculated, and the target parameter set includes anisotropic parameters. It can be seen that the embodiment of the present invention can implement a PP wave and PS wave joint inversion method through PP wave seismic data and PS wave seismic data to predict the anisotropic parameters of the VTI medium, thereby obtaining a target inversion result of anisotropic parameters with higher accuracy, which can effectively improve the accuracy of seismic inversion.

[0048] Based on the above description, the embodiment of the present invention further proposes a more specific joint inversion method. Accordingly, the joint inversion method can be executed by the electronic device (terminal or server) mentioned above; or, the joint inversion method can be executed by the terminal and the server together. For the convenience of explanation, the following description is taken as an example of an electronic device executing the joint inversion method; please refer to Figure 2 , the joint inversion method may include the following steps S201-S208: S201, obtaining initial parameter data of the target area, and obtaining PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data includes upper layer medium parameter data and lower layer medium parameter data.

[0049] S202, constructing a forward operator based on the initial parameter data; and constructing a target inversion function based on the forward operator, the PP wave seismic data, and the PS wave seismic data.

[0050] S203, determining current parameter data to be inverted of the parameter to be inverted, and calculating partial derivatives of the target inversion function based on the current parameter data to be inverted.

[0051] Optionally, in the 0th iteration (if k is used to represent the number of iterations, k=0 in this case, and k is a positive integer), the electronic device may use the initial parameter data to determine the current parameter data to be inverted (i.e., the value of X in the kth iteration) of the parameter to be inverted; that is, the parameter values ​​of each target parameter determined by the initial parameter data may be used to determine the current parameter data to be inverted so as to initialize the parameter to be inverted; illustratively, the current parameter data to be inverted may be determined based on the calculation method corresponding to each target parameter and the parameter value of each target parameter, and so on.

[0052] Optionally, when calculating the partial derivative of the target inversion function based on the current parameter data to be inverted, the electronic device may calculate the first-order partial derivative and the second-order partial derivative of the target inversion function based on the current parameter data to be inverted, so as to calculate the partial derivative of the target inversion function. Optionally, the electronic device may use formula 2.1 to calculate the first-order partial derivative and the second-order partial derivative of the target inversion function: Formula 2.1 in, and can represent the first-order partial derivative and the second-order partial derivative respectively, J and H can be the Jacobian matrix and Hessian matrix of the function g respectively, which can be specifically shown in formula 2.2: Formula 2.2 Based on this, the electronic device can calculate the Jacobian matrix and the Hessian matrix based on the current parameter data to be inverted, so as to calculate the first-order partial derivatives and the second-order partial derivatives of the target inversion function through the Jacobian matrix and the Hessian matrix.

[0053] S204, based on the partial derivative of the target inversion function, updating the parameters to be inverted to obtain updated inversion parameter data of the parameters to be inverted.

[0054] In the embodiment of the present invention, the iterative Newton method can be used for solving; according to the iterative Newton method, the updating principle of the parameter to be inverted X in each iteration can be as shown in formula 2.3, that is, the electronic device can use formula 2.3 to update the parameter to be inverted and obtain the updated inversion parameter data of the parameter to be inverted: Formula 2.3 Among them, k is the number of iterations, and the parameter μ k is the iteration step length at the kth iteration, X k It can represent the value of the parameter to be inverted at the kth iteration (i.e. the current parameter data to be inverted), X k+1 It can represent the updated inversion parameter data of the parameter to be inverted; optionally, the iteration step lengths under any two iteration times can be the same or different, which is not limited in the embodiment of the present invention; illustratively, the iteration step length under each iteration time can be set to a specified step length according to experience or actual needs, or the iteration step length under any iteration time can be randomly sampled in a preset iteration step length interval, etc. Optionally, the specified step length and the preset iteration step length interval can be set according to experience or according to actual needs, which is not limited in the embodiment of the present invention.

[0055] Optionally, for the Newton method, the iteration step size is an important parameter that affects the convergence efficiency and accuracy of the algorithm. The electronic device can also determine the iteration step size through the strong Wolfe convergence criterion. The determination condition can be written as shown in Formula 2.4: Formula 2.4 The parameter t1 usually takes a smaller value (such as t1∈(0,1)), and the parameter t2 may take a value greater than t1 (such as t2∈(t1,1)); optionally, t1 and t2 may be set according to experience or according to actual needs, which is not limited in the embodiment of the present invention; illustratively, t1 may be equal to 0.001, and t2 may be equal to 0.1. Optionally, ∆X k It can be shown as formula 2.5: Formula 2.5 S205: Based on the updated inversion parameter data and the current parameter data to be inverted, determine whether a convergence condition is met.

[0056] In one embodiment, the electronic device may use the updated inversion parameter data and the current parameter data to be inverted to calculate the inversion difference value, and when the inversion difference value is less than or equal to the preset convergence error, it is determined that the convergence condition is met; when the inversion difference value is greater than the preset convergence error, it is determined that the convergence condition is not met. Optionally, the electronic device may use formula 2.6 to calculate the inversion difference value: Formula 2.6 Wherein, Φ may represent an inversion difference value; optionally, the preset convergence error (which may be expressed as η) may be set according to experience or according to actual needs, which is not limited in the embodiment of the present invention.

[0057] In another implementation manner, the electronic device may also determine that the convergence condition is met when the number of iterations reaches a preset iteration number threshold; and / or determine that the convergence condition is not met when the number of iterations does not reach the preset iteration number threshold, etc.; the embodiment of the present invention is not limited to this. The optional preset iteration number threshold may be set according to experience or according to actual needs, and the embodiment of the present invention is not limited to this.

[0058] S206, if the convergence condition is not met, the updated inversion parameter data is used as the current parameter data to be inverted, and the partial derivative of the target inversion function is calculated based on the current parameter data to be inverted, and iterative execution is continued until the convergence condition is met.

[0059] Based on this, the electronic device can execute k=k+1 to iteratively calculate the partial derivative of the target inversion function based on the current parameter data to be inverted, and so on. Optionally, in each iteration process, the electronic device can also first use the current parameter data to be inverted to update the target inversion function, that is, use the current parameter data to be inverted to update the stiffness matrix, and then update the target inversion function based on the current stiffness matrix, so as to calculate the partial derivative of the current target inversion function based on the current parameter data to be inverted, and then continuously obtain the updated inversion parameter data under the next iteration until the convergence condition is reached. It can be seen that the embodiment of the present invention can determine the current parameter value of each target parameter through the current parameter data to be inverted to update the stiffness matrix, and so on.

[0060] S207: If the convergence condition is met, the updated inversion parameter data is used as the target inversion result of the parameter to be inverted.

[0061] In the embodiment of the present invention, the updated inversion parameter data when the convergence condition is reached can be used as the target inversion result of the parameter to be inverted.

[0062] S208, based on the target inversion result of the parameter to be inverted, calculating the target inversion result of each target parameter in the target parameter set, where the target parameter set includes anisotropic parameters.

[0063] Optionally, the electronic device may also determine the target stiffness matrix based on the target inversion results of each target parameter, that is, the target inversion results of each target parameter may be used as the parameter values ​​of the corresponding target parameters to calculate the target stiffness matrix; and the target longitudinal wave polarization direction coefficient and the target transverse wave polarization direction coefficient may be determined based on the target stiffness matrix; based on this, the electronic device may determine the target PP wave diagonal parameter based on the target longitudinal wave polarization direction coefficient, and determine the target PS wave diagonal parameter based on the target transverse wave polarization direction coefficient; thereby, the target PP wave diagonal parameter, the target stiffness matrix and the target longitudinal wave polarization direction coefficient are used to calculate the target PP wave reflection coefficient, and the target PS wave diagonal parameter, the target stiffness matrix and the target transverse wave polarization direction coefficient are used to calculate the target PS wave reflection coefficient. Optionally, the electronic device may calculate the target PP wave reflection coefficient by using the target PP wave diagonal parameter, the target stiffness matrix and the target longitudinal wave polarization direction coefficient through formula 1.12; correspondingly, the target PS wave reflection coefficient may be calculated by using the target PS wave diagonal parameter, the target stiffness matrix and the target transverse wave polarization direction coefficient through formula 1.16.

[0064] In the embodiment of the present invention, the target inversion result of each target parameter may be the target inversion result of each target parameter in the target area.

[0065] In the embodiment of the present invention, in order to further verify the feasibility of the PP wave and PS wave reflection coefficient approximation and the proposed inversion method mentioned in the embodiment of the present invention, on the one hand, the accuracy of the PP wave and PS wave reflection coefficient approximation equation is tested on a double-layer medium; illustratively, the upper layer medium can be an isotropic medium, and the lower layer medium can be a VTI medium with strong anisotropic characteristics. The embodiment of the present invention can prove that the derived approximate equation has high accuracy by comparing with the accurate reflection coefficient (such as calculated by formula 1.7), thereby effectively simplifying the calculation of the reflection coefficient while ensuring high accuracy. For example, assuming that the upper layer medium is an isotropic medium, its longitudinal and transverse wave velocities and densities are 5.5km / s (kilometers per second), 3.5km / s and 2.7g / cm 3 (g / cm3) (the anisotropy parameters can be all 0 at this time), the lower layer is a VTI medium with strong anisotropy characteristics, and its longitudinal and transverse wave velocities and densities are 4.6 km / s, 2.5 km / s and 2.65 g / cm 3 , the anisotropic parameters are δ=0.15 and ε=0.25; in this case, the embodiment of the present invention can calculate the PP wave reflection coefficient and the PS wave reflection coefficient by formula 1.12 and formula 1.16 respectively; based on this, Figure 3 As shown, the solid line can represent the exact value of the reflection coefficient, the dotted line can represent the classic weak anisotropy approximation (Rüger approximation equation), and the dotted line can represent the PP wave weak contrast reflection coefficient proposed in the embodiment of the present invention (i.e., the PP wave reflection coefficient calculated by the approximate equation proposed in the embodiment of the present invention), wherein the dotted line is substantially consistent with the solid line, which proves that the PP reflection coefficient proposed in the embodiment of the present invention is effective; accordingly, as Figure 4 As shown, the dashed line used to represent the approximate calculation of the PS wave reflection coefficient (i.e., the PS wave weak contrast reflection coefficient) is more consistent with the solid line, especially at large angles, proving that the PS wave reflection coefficient proposed in the embodiment of the present invention has a higher accuracy. Figure 3 and Figure 4 The horizontal axis is the incident angle (unit: degree (°)), and the vertical axis is the reflection coefficient.

[0066] On the other hand, the embodiments of the present invention tested the feasibility of the proposed PP wave and PS wave joint inversion method on the logging model and actual seismic data. In the model test, the noise resistance of the proposed method can be tested by adding different degrees of Gaussian white noise to the synthetic seismic data; in the actual seismic data, the stability and reliability of the proposed method can be analyzed by comparing with the conventional anisotropic seismic prediction method.

[0067] For example, the embodiment of the present invention can test the feasibility of the inversion method on a model of an actual well logging curve from a shale reservoir; wherein the P-wave velocity and density in the well logging curve can be actual well logging data, and the anisotropy parameters can be calculated by the rock physics model; optionally, the seismic data can be synthesized using the precise reflection coefficient equation (such as formula 1.7), and then different degrees of noise can be added to obtain seismic data with a signal-to-noise ratio of 8 and 4, and the anisotropy parameters can be predicted using the joint inversion method proposed in the embodiment of the present invention. Based on this, Figure 5 is the inversion result of seismic data with a signal-to-noise ratio of 8. The black solid line is the true model (i.e., the true value of each target parameter in the corresponding test area, such as calculated by the rock physics model, etc.), the gray solid line is the initial model (such as obtained based on experience or tests), and the dotted line is the target inversion result; Figure 5 It can be seen that the target inversion results are in good agreement with the true model, that is, the target inversion results of each target parameter (such as the target inversion results of the longitudinal wave velocity, etc.) are more accurate; in addition, Figure 6 The inversion results of seismic data with a signal-to-noise ratio of 4 are shown in Figure 1. The target inversion results of each target parameter are close to the true value and have high accuracy. It should be understood that compared with the inversion results with a signal-to-noise ratio of 8, the inversion accuracy of the inversion results with a signal-to-noise ratio of 4 has decreased, but it is still consistent with the real model, especially in locations with high anisotropy, which proves that the proposed method is feasible and effective. Figure 5 and Figure 6 The vertical axis is time and the horizontal axis is amplitude.

[0068] On the other hand, the embodiment of the present invention also applies the proposed joint inversion method to a certain actual shale work area, and by comparing it with the inversion method based on weak anisotropy approximation, the reliability of the proposed joint inversion method can be effectively verified. Figure 7 is the anisotropy parameter δ inverted by the conventional inversion method (such as the Rüger inversion method), and the black line is the anisotropy parameter δ calculated by the rock physics model (the black line is used to indicate that the curve on the right side (i.e., the right side of the dotted line) indicates a higher parameter value, and the curve on the left side indicates a lower parameter value, i.e., the more the black line deviates to the right, the larger the parameter value indicated). It can be seen that the result obtained by the conventional method is somewhat different from the logging curve (i.e., the black line), where the highest value on the logging curve indicates the location of the shale reservoir, and the color of the highest value on the logging curve of the conventional inversion method is darker (i.e., the inversion result is smaller), that is, the conventional inversion method fails to accurately predict the abnormal value of the anisotropy parameter in the shale reservoir; Figure 8The anisotropy parameter δ obtained by the joint inversion method proposed in the embodiment of the present invention has a good fit with the well logging curve, clearly indicating the location of the shale reservoir, that is, the target inversion result of the anisotropy parameter δ reflects that the brightness corresponding to the highest value on the black line is the brightest (that is, the parameter value is the highest), that is, the brightness of the target inversion result at the highest value on the black line is higher than the brightness of the inversion result of the conventional inversion method at the highest value on the black line, that is, the target inversion result corresponding to the highest value on the black line is greater than the inversion result of the corresponding conventional method. Correspondingly, Fig. 9 and Fig.10 The anisotropy parameter ε inverted by the conventional method and the joint inversion method proposed in the embodiment of the present invention are respectively; in comparison, the target inversion result of the anisotropy parameter ε obtained by the joint inversion method proposed in the embodiment of the present invention is more consistent with the well logging curve, proving that the inversion proposed in the embodiment of the present invention can effectively improve the stability and accuracy of anisotropy parameter prediction. Among them, Figure 7-10 The vertical axis can be time, and the horizontal axis can be the number of seismic traces.

[0069] The embodiment of the present invention can construct a forward operator based on the initial parameter data after obtaining the initial parameter data of the target area, and the PP wave seismic data and the PS wave seismic data of the target area; and construct a target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data. Then, the current parameter data to be inverted of the parameter to be inverted can be determined, and the partial derivative of the target inversion function can be calculated based on the current parameter data to be inverted; and the parameter to be inverted can be updated based on the partial derivative of the target inversion function to obtain the updated inversion parameter data of the parameter to be inverted. Based on this, it can be judged whether the convergence condition is met based on the updated inversion parameter data and the current parameter data to be inverted; if the convergence condition is not met, the updated inversion parameter data is used as the current parameter data to be inverted, and the partial derivative of the target inversion function is calculated based on the current parameter data to be inverted, and iterative execution is continued until the convergence condition is met; if the convergence condition is met, the updated inversion parameter data is used as the target inversion result of the parameter to be inverted. Furthermore, based on the target inversion results of the parameters to be inverted, the target inversion results of each target parameter in the target parameter set can be calculated, and the target parameter set includes anisotropic parameters. It can be seen that the embodiment of the present invention can realize joint inversion through PP wave seismic data and PS wave seismic data, which can effectively improve the stability and accuracy of anisotropic parameter prediction.

[0070] Based on the description of the above-mentioned embodiments of the joint inversion method, the embodiment of the present invention further proposes a joint inversion device, which may be a computer program (including program code) running in an electronic device; Fig.11 As shown, the joint inversion device may include an acquisition unit 1101 and a processing unit 1102. The joint inversion device may perform Figure 1 or Figure 2 The joint inversion method shown, that is, the joint inversion device can run the above units: An acquisition unit 1101 is used to acquire initial parameter data of a target area, and acquire PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data includes upper layer medium parameter data and lower layer medium parameter data; The processing unit 1102 is used to construct a forward operator based on the initial parameter data; and to construct a target inversion function based on the forward operator, the PP wave seismic data and the PS wave seismic data; The processing unit 1102 is further configured to determine a target inversion result of the parameter to be inverted based on the target inversion function; The processing unit 1102 is further configured to calculate a target inversion result of each target parameter in a target parameter set based on the target inversion result of the parameter to be inverted, wherein the target parameter set includes anisotropic parameters.

[0071] In one implementation, the forward operator includes a first forward operator and a second forward operator; when the processing unit 1102 constructs the forward operator based on the initial parameter data, it can be specifically used to: Based on the initial parameter data, a stiffness matrix is ​​determined, the stiffness matrix includes a plurality of stiffness elements, an initial element value of a stiffness element is determined based on an initial value of each target parameter in at least one target parameter in the target parameter set, the initial parameter data includes initial values ​​of each target parameter in the target parameter set, and an initial value of a target parameter includes an upper initial value of the corresponding target parameter in the upper medium parameter data and / or a lower initial value in the lower medium parameter data; Based on the stiffness matrix, a plurality of PP wave parameter perturbation coefficients are calculated; and based on the stiffness matrix, a plurality of PS wave parameter perturbation coefficients are calculated; wherein the plurality of PP wave parameter perturbation coefficients include the PP wave parameter perturbation coefficients corresponding to each stiffness element in the plurality of stiffness elements, and the plurality of PS wave parameter perturbation coefficients include the PS wave parameter perturbation coefficients corresponding to each stiffness element; The first forward operator is constructed based on the multiple PP wave parameter perturbation coefficients; and the second forward operator is constructed based on the multiple PS wave parameter perturbation coefficients, so as to realize the construction of the forward operator.

[0072] In another implementation manner, when the processing unit 1102 calculates a plurality of PP wave parameter disturbance coefficients based on the stiffness matrix, it may be specifically used to: Based on the stiffness matrix, a plurality of intermediate coefficients are determined; and based on the plurality of intermediate coefficients, the longitudinal wave polarization direction coefficient and the longitudinal wave vertical slowness, the PP wave diagonal parameters are calculated; Calculating a plurality of PP wave parameter disturbance coefficients based on the PP wave diagonal parameters, the stiffness matrix, and the plurality of intermediate coefficients; When the processing unit 1102 calculates a plurality of PS wave parameter disturbance coefficients based on the stiffness matrix, it may be specifically used to: Calculating PS wave diagonal parameters based on the multiple intermediate coefficients, the shear wave polarization direction coefficient, and the shear wave vertical slowness; Based on the PS wave diagonal parameters, the stiffness matrix and the multiple intermediate coefficients, multiple PS wave parameter disturbance coefficients are calculated.

[0073] In another embodiment, the forward operator includes a first forward operator and a second forward operator; when the processing unit 1102 constructs a target inversion function based on the forward operator, the PP wave seismic data, and the PS wave seismic data, it can be specifically used to: Determine the stiffness matrix and anisotropic parameter relationship constraints; constructing a PP wave seismic data error term based on the first forward operator and the PP wave seismic data; and constructing a PS wave seismic data error term based on the second forward operator and the PS wave seismic data; The target inversion function is constructed by using the PP wave seismic data error term, the PS wave seismic data error term, and the stiffness matrix and anisotropic parameter relationship constraint term.

[0074] In another implementation manner, when the processing unit 1102 determines the target inversion result of the parameter to be inverted based on the target inversion function, it can be specifically used to: Determining current parameter data to be inverted of the parameter to be inverted, and calculating the partial derivative of the target inversion function based on the current parameter data to be inverted; Based on the partial derivative of the target inversion function, the parameter to be inverted is updated to obtain updated inversion parameter data of the parameter to be inverted; Based on the updated inversion parameter data and the current parameter data to be inverted, determining whether a convergence condition is met; If the convergence condition is not met, the updated inversion parameter data is used as the current parameter data to be inverted, and the partial derivative of the target inversion function is calculated based on the current parameter data to be inverted is continued to be iteratively executed until the convergence condition is met; If the convergence condition is met, the updated inversion parameter data is used as the target inversion result of the parameter to be inverted.

[0075] In another implementation manner, when the processing unit 1102 calculates the target inversion result of each target parameter in the target parameter set based on the target inversion result of the parameter to be inverted, it can be specifically used to: Based on the nonlinear relationship between the stiffness matrix and the target parameter set, respectively determining the calculation method corresponding to each target parameter in the target parameter set; The target inversion results of the target parameters are calculated according to the calculation methods corresponding to the target parameters and the target inversion results of the parameters to be inverted.

[0076] In another implementation, the processing unit 1102 may also be used to: Determining a target stiffness matrix based on the target inversion results of the respective target parameters; Based on the target stiffness matrix, determining a target longitudinal wave polarization direction coefficient and a target shear wave polarization direction coefficient; Determining a target PP wave diagonal parameter based on the target longitudinal wave polarization direction coefficient, and determining a target PS wave diagonal parameter based on the target shear wave polarization direction coefficient; The target PP wave diagonal parameter, the target stiffness matrix and the target longitudinal wave polarization direction coefficient are used to calculate the target PP wave reflection coefficient; and the target PS wave diagonal parameter, the target stiffness matrix and the target shear wave polarization direction coefficient are used to calculate the target PS wave reflection coefficient.

[0077] According to one embodiment of the present invention, Fig.11 Each unit in the joint inversion device shown can be separately or completely combined into one or several other units to form a unit, or one (or some) of the units can be further divided into multiple functionally smaller units to form a unit, which can achieve the same operation without affecting the realization of the technical effects of the embodiments of the present invention. The above-mentioned units are divided based on logical functions. In practical applications, the function of one unit can also be implemented by multiple units, or the functions of multiple units can be implemented by one unit. In other embodiments of the present invention, any joint inversion device can also include other units. In practical applications, these functions can also be implemented with the assistance of other units, and can be implemented by the collaboration of multiple units.

[0078] According to another embodiment of the present invention, the program can be executed by running a program on a general electronic device such as a computer including a central processing unit (CPU), a random access memory medium (RAM), a read-only memory medium (ROM), and other processing elements and storage elements. Figure 1 or Figure 2 A computer program (including program code) for each step involved in the corresponding method shown in Fig.11 The joint inversion device shown in the embodiment of the present invention is used to implement the joint inversion method of the embodiment of the present invention. The computer program can be recorded on a computer storage medium, for example, and loaded into the above electronic device through the computer storage medium and run therein.

[0079] Based on the description of the above method embodiment and device embodiment, the exemplary embodiment of the present invention further provides an electronic device, including: at least one processor; and a memory connected to the at least one processor in communication. The memory stores a computer program that can be executed by the at least one processor, and the computer program is used to enable the electronic device to perform the method according to the embodiment of the present invention when executed by the at least one processor.

[0080] Exemplary embodiments of the present invention also provide a non-transitory computer-readable storage medium storing a computer program, wherein the computer program, when executed by a processor of a computer, is used to cause the computer to perform a method according to an embodiment of the present invention.

[0081] An exemplary embodiment of the present invention further provides a computer program product, comprising a computer program, wherein when the computer program is executed by a processor of a computer, the computer is used to enable the computer to perform a method according to an embodiment of the present invention.

[0082] refer to Fig.12 , a block diagram of an electronic device 1200 that can be used as a server or client of the present invention will now be described, 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, workbenches, 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 processing, cellular phones, smart phones, wearable devices, and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely examples and are not intended to limit the implementation of the present invention described and / or required herein.

[0083] like Fig.12 As shown, the electronic device 1200 includes a computing unit 1201, which can perform various appropriate actions and processes according to a computer program stored in a read-only memory (ROM) 1202 or a computer program loaded from a storage unit 1208 to a random access memory (RAM) 1203. In the RAM 1203, various programs and data required for the operation of the electronic device 1200 can also be stored. The computing unit 1201, the ROM 1202, and the RAM 1203 are connected to each other via a bus 1204. An input / output (I / O) interface 1205 is also connected to the bus 1204.

[0084] Multiple components in the electronic device 1200 are connected to the I / O interface 1205, including: an input unit 1206, an output unit 1207, a storage unit 1208, and a communication unit 1209. The input unit 1206 can be any type of device that can input information to the electronic device 1200, and the input unit 1206 can receive input digital or character information, and generate key signal input related to user settings and / or function control of the electronic device. The output unit 1207 can be any type of device that can present information, and can include but is not limited to a display, a speaker, a video / audio output terminal, a vibrator, and / or a printer. The storage unit 1208 can include but is not limited to a disk, an optical disk. The communication unit 1209 allows the electronic device 1200 to exchange information / data with other devices through a computer network such as the Internet and / or various telecommunication networks, and can include but is not limited to a modem, a network card, an infrared communication device, a wireless communication transceiver, and / or a chipset, such as a Bluetooth™ device, a WiFi device, a WiMax device, a cellular communication device, and / or the like.

[0085] The computing unit 1201 may be a variety of general and / or special processing components with processing and computing capabilities. Some examples of the computing unit 1201 include, but are not limited to, a central processing unit (CPU), a graphics processing unit (GPU), various dedicated artificial intelligence (AI) computing chips, various computing units running machine learning model algorithms, digital signal processors (DSPs), and any appropriate processors, controllers, microcontrollers, etc. The computing unit 1201 performs the various methods and processes described above. For example, in some embodiments, the joint inversion method may be implemented as a computer software program, which is tangibly contained in a machine-readable medium, such as a storage unit 1208. In some embodiments, part or all of the computer program may be loaded and / or installed on the electronic device 1200 via the ROM 1202 and / or the communication unit 1209. In some embodiments, the computing unit 1201 may be configured to perform the joint inversion method in any other appropriate manner (e.g., by means of firmware).

[0086] The program code for implementing the method of the present invention can be written in any combination of one or more programming languages. These program codes can be provided to a processor or controller of a general-purpose computer, a special-purpose computer or other programmable data processing device, so that the program code, when executed by the processor or controller, enables the functions / operations specified in the flow chart and / or block diagram to be implemented. The program code can be executed entirely on the machine, partially on the machine, partially on the machine as a stand-alone software package and partially on a remote machine, or entirely on a remote machine or server.

[0087] In the context of the present invention, a machine-readable medium may be a tangible medium that may contain or store a program for use by or in conjunction with an instruction execution system, device, or equipment. A machine-readable medium may be a machine-readable signal medium or a machine-readable storage medium. A machine-readable medium may include, but is not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, device, or equipment, or any suitable combination of the foregoing. A more specific example of a machine-readable storage medium may include an electrical connection based on one or more lines, a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.

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

[0089] To provide interaction with a user, the systems and techniques described herein can be implemented on a computer having: a display device (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor) for displaying information to the user; and a keyboard and pointing device (e.g., a mouse or trackball) through which the user can provide input to the computer. Other types of devices can also be used to provide interaction with the user; for example, the 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 acoustic input, voice input, or tactile input).

[0090] Furthermore, it should be understood that what is disclosed above is only a preferred embodiment of the present invention, and certainly cannot be used to limit the scope of rights of the present invention. Therefore, equivalent changes made according to the claims of the present invention are still within the scope of the present invention.

Claims

1. A joint inversion method, characterized in that: include: Acquiring initial parameter data of a target area, and acquiring PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data includes upper layer medium parameter data and lower layer medium parameter data; Based on the initial parameter data, a forward operator is constructed; and based on the forward operator, the PP wave seismic data and the PS wave seismic data, a target inversion function is constructed; Based on the target inversion function, determining a target inversion result of the parameter to be inverted; Based on the target inversion results of the parameters to be inverted, target inversion results of each target parameter in a target parameter set are calculated, and the target parameter set includes anisotropic parameters.

2. The method according to claim 1, characterized in that The forward operator includes a first forward operator and a second forward operator; the constructing of the forward operator based on the initial parameter data includes: Based on the initial parameter data, a stiffness matrix is ​​determined, the stiffness matrix includes a plurality of stiffness elements, an initial element value of a stiffness element is determined based on an initial value of each target parameter in at least one target parameter in the target parameter set, the initial parameter data includes initial values ​​of each target parameter in the target parameter set, and an initial value of a target parameter includes an upper initial value of the corresponding target parameter in the upper medium parameter data and / or a lower initial value in the lower medium parameter data; Based on the stiffness matrix, a plurality of PP wave parameter perturbation coefficients are calculated; and based on the stiffness matrix, a plurality of PS wave parameter perturbation coefficients are calculated; wherein the plurality of PP wave parameter perturbation coefficients include the PP wave parameter perturbation coefficients corresponding to each stiffness element in the plurality of stiffness elements, and the plurality of PS wave parameter perturbation coefficients include the PS wave parameter perturbation coefficients corresponding to each stiffness element; The first forward operator is constructed based on the multiple PP wave parameter perturbation coefficients; and the second forward operator is constructed based on the multiple PS wave parameter perturbation coefficients, so as to realize the construction of the forward operator.

3. The method according to claim 2, characterized in that The step of calculating a plurality of PP wave parameter disturbance coefficients based on the stiffness matrix includes: Based on the stiffness matrix, a plurality of intermediate coefficients are determined; and based on the plurality of intermediate coefficients, the longitudinal wave polarization direction coefficient and the longitudinal wave vertical slowness, the PP wave diagonal parameters are calculated; Calculating a plurality of PP wave parameter disturbance coefficients based on the PP wave diagonal parameters, the stiffness matrix, and the plurality of intermediate coefficients; The step of calculating a plurality of PS wave parameter disturbance coefficients based on the stiffness matrix comprises: Calculating PS wave diagonal parameters based on the multiple intermediate coefficients, the shear wave polarization direction coefficient, and the shear wave vertical slowness; Based on the PS wave diagonal parameters, the stiffness matrix and the multiple intermediate coefficients, multiple PS wave parameter disturbance coefficients are calculated.

4. The method according to any one of claims 1 to 3, characterized in that: The forward operator includes a first forward operator and a second forward operator; and constructing a target inversion function based on the forward operator, the PP wave seismic data, and the PS wave seismic data includes: Determine the stiffness matrix and anisotropic parameter relationship constraints; constructing a PP wave seismic data error term based on the first forward operator and the PP wave seismic data; and constructing a PS wave seismic data error term based on the second forward operator and the PS wave seismic data; The target inversion function is constructed by using the PP wave seismic data error term, the PS wave seismic data error term, and the stiffness matrix and anisotropic parameter relationship constraint term.

5. The method according to any one of claims 1 to 3, characterized in that: The step of determining a target inversion result of the parameter to be inverted based on the target inversion function comprises: Determining current parameter data to be inverted of the parameter to be inverted, and calculating the partial derivative of the target inversion function based on the current parameter data to be inverted; Based on the partial derivative of the target inversion function, the parameter to be inverted is updated to obtain updated inversion parameter data of the parameter to be inverted; Based on the updated inversion parameter data and the current parameter data to be inverted, determining whether a convergence condition is met; If the convergence condition is not met, the updated inversion parameter data is used as the current parameter data to be inverted, and the partial derivative of the target inversion function is calculated based on the current parameter data to be inverted is continued to be iteratively executed until the convergence condition is met; If the convergence condition is met, the updated inversion parameter data is used as the target inversion result of the parameter to be inverted.

6. The method according to any one of claims 1 to 3, characterized in that: The target inversion result based on the parameter to be inverted, calculating the target inversion result of each target parameter in the target parameter set, comprises: Based on the nonlinear relationship between the stiffness matrix and the target parameter set, respectively determining the calculation method corresponding to each target parameter in the target parameter set; The target inversion results of the target parameters are calculated according to the calculation methods corresponding to the target parameters and the target inversion results of the parameters to be inverted.

7. The method according to any one of claims 1 to 3, characterized in that: The method further comprises: Determining a target stiffness matrix based on the target inversion results of the respective target parameters; Based on the target stiffness matrix, determining a target longitudinal wave polarization direction coefficient and a target shear wave polarization direction coefficient; Determining a target PP wave diagonal parameter based on the target longitudinal wave polarization direction coefficient, and determining a target PS wave diagonal parameter based on the target shear wave polarization direction coefficient; The target PP wave diagonal parameters, the target stiffness matrix and the target longitudinal wave polarization direction coefficient are used to calculate the target PP wave reflection coefficient; and the target PS wave diagonal parameters, the target stiffness matrix and the target shear wave polarization direction coefficient are used to calculate the target PS wave reflection coefficient.

8. A joint inversion device, characterized in that: The device comprises: An acquisition unit, used to acquire initial parameter data of a target area, and acquire PP wave seismic data and PS wave seismic data of the target area, wherein the initial parameter data includes upper layer medium parameter data and lower layer medium parameter data; A processing unit, configured to construct a forward operator based on the initial parameter data; and to construct a target inversion function based on the forward operator, the PP wave seismic data, and the PS wave seismic data; The processing unit is further used to determine a target inversion result of the parameter to be inverted based on the target inversion function; The processing unit is further used to calculate the target inversion result of each target parameter in a target parameter set based on the target inversion result of the parameter to be inverted, and the target parameter set includes anisotropic parameters.

9. An electronic device, characterized in that: include: processor; as well as Memory for storing programs, The program includes instructions, which, when executed by the processor, cause the processor to perform the method according to any one of claims 1 to 7.

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

Images